Actual source code: ex69.c
2: #include <petscdt.h>
3: #include <petscdraw.h>
4: #include <petscviewer.h>
5: #include <petscksp.h>
6: #include <petscdmda.h>
8: /*
9: Solves -Laplacian u = f, u(-1) = u(1) = 0 with multiple spectral elements
11: Uses DMDA to manage the parallelization of the elements
13: This is not intended to be highly optimized in either memory usage or time, but strifes for simplicity.
15: */
17: typedef struct {
18: PetscInt n; /* number of nodes */
19: PetscReal *nodes; /* GLL nodes */
20: PetscReal *weights; /* GLL weights */
21: } PetscGLL;
23: PetscErrorCode ComputeSolution(DM da, PetscGLL *gll, Vec u)
24: {
25: PetscInt j, xs, xn;
26: PetscScalar *uu, *xx;
27: PetscReal xd;
28: Vec x;
30: DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
31: DMGetCoordinates(da, &x);
32: DMDAVecGetArray(da, x, &xx);
33: DMDAVecGetArray(da, u, &uu);
34: /* loop over local nodes */
35: for (j = xs; j < xs + xn; j++) {
36: xd = xx[j];
37: uu[j] = (xd * xd - 1.0) * PetscCosReal(5. * PETSC_PI * xd);
38: }
39: DMDAVecRestoreArray(da, x, &xx);
40: DMDAVecRestoreArray(da, u, &uu);
41: return 0;
42: }
44: /*
45: 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
46: basis function is zero at all nodes except the ith one the integral is simply the weight_i * f(node_i)
47: */
48: PetscErrorCode ComputeRhs(DM da, PetscGLL *gll, Vec b)
49: {
50: PetscInt i, j, xs, xn, n = gll->n;
51: PetscScalar *bb, *xx;
52: PetscReal xd;
53: Vec blocal, xlocal;
55: DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
56: xs = xs / (n - 1);
57: xn = xn / (n - 1);
58: DMGetLocalVector(da, &blocal);
59: VecZeroEntries(blocal);
60: DMDAVecGetArray(da, blocal, &bb);
61: DMGetCoordinatesLocal(da, &xlocal);
62: DMDAVecGetArray(da, xlocal, &xx);
63: /* loop over local spectral elements */
64: for (j = xs; j < xs + xn; j++) {
65: /* loop over GLL points in each element */
66: for (i = 0; i < n; i++) {
67: xd = xx[j * (n - 1) + i];
68: bb[j * (n - 1) + i] += -gll->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));
69: }
70: }
71: DMDAVecRestoreArray(da, xlocal, &xx);
72: DMDAVecRestoreArray(da, blocal, &bb);
73: VecZeroEntries(b);
74: DMLocalToGlobalBegin(da, blocal, ADD_VALUES, b);
75: DMLocalToGlobalEnd(da, blocal, ADD_VALUES, b);
76: DMRestoreLocalVector(da, &blocal);
77: return 0;
78: }
80: /*
81: Run with -build_twosided allreduce -pc_type bjacobi -sub_pc_type lu -q 16 -ksp_rtol 1.e-34 (or 1.e-14 for double precision)
83: -q <q> number of spectral elements to use
84: -N <N> maximum number of GLL points per element
86: */
87: int main(int argc, char **args)
88: {
89: PetscGLL gll;
90: PetscInt N = 80, n, q = 8, xs, xn, j, l;
91: PetscReal **A;
92: Mat K;
93: KSP ksp;
94: PC pc;
95: Vec x, b;
96: PetscInt *rows;
97: PetscReal norm, xc, yc, h;
98: PetscScalar *f;
99: PetscDraw draw;
100: PetscDrawLG lg;
101: PetscDrawAxis axis;
102: DM da;
103: PetscMPIInt rank, size;
106: PetscInitialize(&argc, &args, NULL, NULL);
107: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
108: MPI_Comm_size(PETSC_COMM_WORLD, &size);
109: PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
110: PetscOptionsGetInt(NULL, NULL, "-q", &q, NULL);
112: PetscDrawCreate(PETSC_COMM_WORLD, NULL, "Log(Error norm) vs Number of GLL points", 0, 0, 500, 500, &draw);
113: PetscDrawSetFromOptions(draw);
114: PetscDrawLGCreate(draw, 1, &lg);
115: PetscDrawLGSetUseMarkers(lg, PETSC_TRUE);
116: PetscDrawLGGetAxis(lg, &axis);
117: PetscDrawAxisSetLabels(axis, NULL, "Number of GLL points", "Log(Error Norm)");
119: for (n = 4; n < N; n += 2) {
120: /*
121: da contains the information about the parallel layout of the elements
122: */
123: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, q * (n - 1) + 1, 1, 1, NULL, &da);
124: DMSetFromOptions(da);
125: DMSetUp(da);
126: DMDAGetInfo(da, NULL, &q, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
127: q = (q - 1) / (n - 1); /* number of spectral elements */
129: /*
130: gll simply contains the GLL node and weight values
131: */
132: PetscMalloc2(n, &gll.nodes, n, &gll.weights);
133: PetscDTGaussLobattoLegendreQuadrature(n, PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA, gll.nodes, gll.weights);
134: gll.n = n;
135: DMDASetGLLCoordinates(da, gll.n, gll.nodes);
137: /*
138: Creates the element stiffness matrix for the given gll
139: */
140: PetscGaussLobattoLegendreElementLaplacianCreate(gll.n, gll.nodes, gll.weights, &A);
142: /*
143: Scale the element stiffness and weights by the size of the element
144: */
145: h = 2.0 / q;
146: for (j = 0; j < n; j++) {
147: gll.weights[j] *= .5 * h;
148: for (l = 0; l < n; l++) A[j][l] = 2. * A[j][l] / h;
149: }
151: /*
152: Create the global stiffness matrix and add the element stiffness for each local element
153: */
154: DMCreateMatrix(da, &K);
155: MatSetOption(K, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
156: DMDAGetCorners(da, &xs, NULL, NULL, &xn, NULL, NULL);
157: xs = xs / (n - 1);
158: xn = xn / (n - 1);
159: PetscMalloc1(n, &rows);
160: /*
161: loop over local elements
162: */
163: for (j = xs; j < xs + xn; j++) {
164: for (l = 0; l < n; l++) rows[l] = j * (n - 1) + l;
165: MatSetValues(K, n, rows, n, rows, &A[0][0], ADD_VALUES);
166: }
167: MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY);
170: MatCreateVecs(K, &x, &b);
171: ComputeRhs(da, &gll, b);
173: /*
174: Replace the first and last rows/columns of the matrix with the identity to obtain the zero Dirichlet boundary conditions
175: */
176: rows[0] = 0;
177: rows[1] = q * (n - 1);
178: MatZeroRowsColumns(K, 2, rows, 1.0, x, b);
179: PetscFree(rows);
181: KSPCreate(PETSC_COMM_WORLD, &ksp);
182: KSPSetOperators(ksp, K, K);
183: KSPGetPC(ksp, &pc);
184: PCSetType(pc, PCLU);
185: KSPSetFromOptions(ksp);
186: KSPSolve(ksp, b, x);
188: /* compute the error to the continium problem */
189: ComputeSolution(da, &gll, b);
190: VecAXPY(x, -1.0, b);
192: /* compute the L^2 norm of the error */
193: VecGetArray(x, &f);
194: PetscGaussLobattoLegendreIntegrate(gll.n, gll.nodes, gll.weights, f, &norm);
195: VecRestoreArray(x, &f);
196: norm = PetscSqrtReal(norm);
197: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "L^2 norm of the error %" PetscInt_FMT " %g\n", n, (double)norm);
199: xc = (PetscReal)n;
200: yc = PetscLog10Real(norm);
201: PetscDrawLGAddPoint(lg, &xc, &yc);
202: PetscDrawLGDraw(lg);
204: VecDestroy(&b);
205: VecDestroy(&x);
206: KSPDestroy(&ksp);
207: MatDestroy(&K);
208: PetscGaussLobattoLegendreElementLaplacianDestroy(gll.n, gll.nodes, gll.weights, &A);
209: PetscFree2(gll.nodes, gll.weights);
210: DMDestroy(&da);
211: }
212: PetscDrawLGDestroy(&lg);
213: PetscDrawDestroy(&draw);
214: PetscFinalize();
215: return 0;
216: }
218: /*TEST
220: build:
221: requires: !complex
223: test:
224: requires: !single
226: test:
227: suffix: 2
228: nsize: 2
229: requires: superlu_dist
231: TEST*/