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