Actual source code: ex68.c
2: #include <petscdt.h>
3: #include <petscdraw.h>
4: #include <petscviewer.h>
5: #include <petscksp.h>
7: /*
8: Solves -Laplacian u = f, u(-1) = u(1) = 0 with a single spectral element for n = 4 to N GLL points
10: Plots the L_2 norm of the error (evaluated via the GLL nodes and weights) as a function of n.
12: */
13: PetscErrorCode ComputeSolution(PetscInt n, PetscReal *nodes, PetscReal *weights, Vec x)
14: {
15: PetscInt i, m;
16: PetscScalar *xx;
17: PetscReal xd;
19: VecGetSize(x, &m);
20: VecGetArray(x, &xx);
21: for (i = 0; i < m; i++) {
22: xd = nodes[i];
23: xx[i] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
24: }
25: VecRestoreArray(x, &xx);
26: return 0;
27: }
29: /*
30: Evaluates \integral_{-1}^{1} f*v_i where v_i is the ith basis polynomial via the GLL nodes and weights, since the v_i
31: basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
32: */
33: PetscErrorCode ComputeRhs(PetscInt n, PetscReal *nodes, PetscReal *weights, Vec b)
34: {
35: PetscInt i, m;
36: PetscScalar *bb;
37: PetscReal xd;
39: VecGetSize(b, &m);
40: VecGetArray(b, &bb);
41: for (i = 0; i < m; i++) {
42: xd = nodes[i];
43: bb[i] = -weights[i] * (-20. * PETSC_PI * xd * PetscSinReal(5. * PETSC_PI * xd) + (2. - (5. * PETSC_PI) * (5. * PETSC_PI) * (xd * xd - 1.)) * PetscCosReal(5. * PETSC_PI * xd));
44: }
45: VecRestoreArray(b, &bb);
46: return 0;
47: }
49: int main(int argc, char **args)
50: {
51: PetscReal *nodes;
52: PetscReal *weights;
53: PetscInt N = 80, n;
54: PetscReal **A;
55: Mat K;
56: KSP ksp;
57: PC pc;
58: Vec x, b;
59: PetscInt rows[2];
60: PetscReal norm, xc, yc;
61: PetscScalar *f;
62: PetscDraw draw;
63: PetscDrawLG lg;
64: PetscDrawAxis axis;
67: PetscInitialize(&argc, &args, NULL, NULL);
68: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
70: PetscDrawCreate(PETSC_COMM_SELF, NULL, "Log(Error norm) vs Number of GLL points", 0, 0, 500, 500, &draw);
71: PetscDrawSetFromOptions(draw);
72: PetscDrawLGCreate(draw, 1, &lg);
73: PetscDrawLGSetUseMarkers(lg, PETSC_TRUE);
74: PetscDrawLGGetAxis(lg, &axis);
75: PetscDrawAxisSetLabels(axis, NULL, "Number of GLL points", "Log(Error Norm)");
77: for (n = 4; n < N; n += 2) {
78: /*
79: compute GLL node and weight values
80: */
81: PetscMalloc2(n, &nodes, n, &weights);
82: PetscDTGaussLobattoLegendreQuadrature(n, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, nodes, weights);
83: /*
84: Creates the element stiffness matrix for the given gll
85: */
86: PetscGaussLobattoLegendreElementLaplacianCreate(n, nodes, weights, &A);
87: MatCreateSeqDense(PETSC_COMM_SELF, n, n, &A[0][0], &K);
88: rows[0] = 0;
89: rows[1] = n - 1;
90: KSPCreate(PETSC_COMM_SELF, &ksp);
91: MatCreateVecs(K, &x, &b);
92: ComputeRhs(n, nodes, weights, b);
93: /*
94: Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
95: */
96: MatZeroRowsColumns(K, 2, rows, 1.0, x, b);
97: KSPSetOperators(ksp, K, K);
98: KSPGetPC(ksp, &pc);
99: PCSetType(pc, PCLU);
100: KSPSetFromOptions(ksp);
101: KSPSolve(ksp, b, x);
103: /* compute the error to the continium problem */
104: ComputeSolution(n, nodes, weights, b);
105: VecAXPY(x, -1.0, b);
107: /* compute the L^2 norm of the error */
108: VecGetArray(x, &f);
109: PetscGaussLobattoLegendreIntegrate(n, nodes, weights, f, &norm);
110: VecRestoreArray(x, &f);
111: norm = PetscSqrtReal(norm);
112: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "L^2 norm of the error %" PetscInt_FMT " %g\n", n, (double)norm);
113: xc = (PetscReal)n;
114: yc = PetscLog10Real(norm);
115: PetscDrawLGAddPoint(lg, &xc, &yc);
116: PetscDrawLGDraw(lg);
118: VecDestroy(&b);
119: VecDestroy(&x);
120: KSPDestroy(&ksp);
121: MatDestroy(&K);
122: PetscGaussLobattoLegendreElementLaplacianDestroy(n, nodes, weights, &A);
123: PetscFree2(nodes, weights);
124: }
125: PetscDrawSetPause(draw, -2);
126: PetscDrawLGDestroy(&lg);
127: PetscDrawDestroy(&draw);
128: PetscFinalize();
129: return 0;
130: }
132: /*TEST
134: build:
135: requires: !complex
137: test:
139: TEST*/