Actual source code: ex13.c
1: static char help[] = "Poisson Problem in 2d and 3d with finite elements.\n\
2: We solve the Poisson problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\
4: This example supports automatic convergence estimation\n\
5: and eventually adaptivity.\n\n\n";
7: #include <petscdmplex.h>
8: #include <petscsnes.h>
9: #include <petscds.h>
10: #include <petscconvest.h>
12: typedef struct {
13: /* Domain and mesh definition */
14: PetscBool spectral; /* Look at the spectrum along planes in the solution */
15: PetscBool shear; /* Shear the domain */
16: PetscBool adjoint; /* Solve the adjoint problem */
17: PetscBool homogeneous; /* Use homogeneous boundary conditions */
18: PetscBool viewError; /* Output the solution error */
19: } AppCtx;
21: static PetscErrorCode zero(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
22: {
23: *u = 0.0;
24: return 0;
25: }
27: static PetscErrorCode trig_inhomogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
28: {
29: PetscInt d;
30: *u = 0.0;
31: for (d = 0; d < dim; ++d) *u += PetscSinReal(2.0 * PETSC_PI * x[d]);
32: return 0;
33: }
35: static PetscErrorCode trig_homogeneous_u(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
36: {
37: PetscInt d;
38: *u = 1.0;
39: for (d = 0; d < dim; ++d) *u *= PetscSinReal(2.0 * PETSC_PI * x[d]);
40: return 0;
41: }
43: /* Compute integral of (residual of solution)*(adjoint solution - projection of adjoint solution) */
44: static void obj_error_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[])
45: {
46: obj[0] = a[aOff[0]] * (u[0] - a[aOff[1]]);
47: }
49: static void f0_trig_inhomogeneous_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
50: {
51: PetscInt d;
52: for (d = 0; d < dim; ++d) f0[0] += -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
53: }
55: static void f0_trig_homogeneous_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
56: {
57: PetscInt d;
58: for (d = 0; d < dim; ++d) {
59: PetscScalar v = 1.;
60: for (PetscInt e = 0; e < dim; e++) {
61: if (e == d) {
62: v *= -4.0 * PetscSqr(PETSC_PI) * PetscSinReal(2.0 * PETSC_PI * x[d]);
63: } else {
64: v *= PetscSinReal(2.0 * PETSC_PI * x[d]);
65: }
66: }
67: f0[0] += v;
68: }
69: }
71: static void f0_unity_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
72: {
73: f0[0] = 1.0;
74: }
76: static void f0_identityaux_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
77: {
78: f0[0] = a[0];
79: }
81: static void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
82: {
83: PetscInt d;
84: for (d = 0; d < dim; ++d) f1[d] = u_x[d];
85: }
87: static void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
88: {
89: PetscInt d;
90: for (d = 0; d < dim; ++d) g3[d * dim + d] = 1.0;
91: }
93: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
94: {
96: options->shear = PETSC_FALSE;
97: options->spectral = PETSC_FALSE;
98: options->adjoint = PETSC_FALSE;
99: options->homogeneous = PETSC_FALSE;
100: options->viewError = PETSC_FALSE;
102: PetscOptionsBegin(comm, "", "Poisson Problem Options", "DMPLEX");
103: PetscOptionsBool("-shear", "Shear the domain", "ex13.c", options->shear, &options->shear, NULL);
104: PetscOptionsBool("-spectral", "Look at the spectrum along planes of the solution", "ex13.c", options->spectral, &options->spectral, NULL);
105: PetscOptionsBool("-adjoint", "Solve the adjoint problem", "ex13.c", options->adjoint, &options->adjoint, NULL);
106: PetscOptionsBool("-homogeneous", "Use homogeneous boundary conditions", "ex13.c", options->homogeneous, &options->homogeneous, NULL);
107: PetscOptionsBool("-error_view", "Output the solution error", "ex13.c", options->viewError, &options->viewError, NULL);
108: PetscOptionsEnd();
109: return 0;
110: }
112: static PetscErrorCode CreateSpectralPlanes(DM dm, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user)
113: {
114: PetscSection coordSection;
115: Vec coordinates;
116: const PetscScalar *coords;
117: PetscInt dim, p, vStart, vEnd, v;
120: DMGetCoordinateDim(dm, &dim);
121: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
122: DMGetCoordinatesLocal(dm, &coordinates);
123: DMGetCoordinateSection(dm, &coordSection);
124: VecGetArrayRead(coordinates, &coords);
125: for (p = 0; p < numPlanes; ++p) {
126: DMLabel label;
127: char name[PETSC_MAX_PATH_LEN];
129: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p);
130: DMCreateLabel(dm, name);
131: DMGetLabel(dm, name, &label);
132: DMLabelAddStratum(label, 1);
133: for (v = vStart; v < vEnd; ++v) {
134: PetscInt off;
136: PetscSectionGetOffset(coordSection, v, &off);
137: if (PetscAbsReal(planeCoord[p] - PetscRealPart(coords[off + planeDir[p]])) < PETSC_SMALL) DMLabelSetValue(label, v, 1);
138: }
139: }
140: VecRestoreArrayRead(coordinates, &coords);
141: return 0;
142: }
144: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
145: {
147: DMCreate(comm, dm);
148: DMSetType(*dm, DMPLEX);
149: DMSetFromOptions(*dm);
150: if (user->shear) DMPlexShearGeometry(*dm, DM_X, NULL);
151: DMSetApplicationContext(*dm, user);
152: DMViewFromOptions(*dm, NULL, "-dm_view");
153: if (user->spectral) {
154: PetscInt planeDir[2] = {0, 1};
155: PetscReal planeCoord[2] = {0., 1.};
157: CreateSpectralPlanes(*dm, 2, planeDir, planeCoord, user);
158: }
159: return 0;
160: }
162: static PetscErrorCode SetupPrimalProblem(DM dm, AppCtx *user)
163: {
164: PetscDS ds;
165: DMLabel label;
166: const PetscInt id = 1;
167: PetscPointFunc f0 = user->homogeneous ? f0_trig_homogeneous_u : f0_trig_inhomogeneous_u;
168: PetscErrorCode (*ex)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = user->homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u;
171: DMGetDS(dm, &ds);
172: PetscDSSetResidual(ds, 0, f0, f1_u);
173: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
174: PetscDSSetExactSolution(ds, 0, ex, user);
175: DMGetLabel(dm, "marker", &label);
176: if (label) DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))ex, NULL, user, NULL);
177: return 0;
178: }
180: static PetscErrorCode SetupAdjointProblem(DM dm, AppCtx *user)
181: {
182: PetscDS ds;
183: DMLabel label;
184: const PetscInt id = 1;
187: DMGetDS(dm, &ds);
188: PetscDSSetResidual(ds, 0, f0_unity_u, f1_u);
189: PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_uu);
190: PetscDSSetObjective(ds, 0, obj_error_u);
191: DMGetLabel(dm, "marker", &label);
192: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void))zero, NULL, user, NULL);
193: return 0;
194: }
196: static PetscErrorCode SetupErrorProblem(DM dm, AppCtx *user)
197: {
198: PetscDS prob;
201: DMGetDS(dm, &prob);
202: return 0;
203: }
205: static PetscErrorCode SetupDiscretization(DM dm, const char name[], PetscErrorCode (*setup)(DM, AppCtx *), AppCtx *user)
206: {
207: DM cdm = dm;
208: PetscFE fe;
209: DMPolytopeType ct;
210: PetscBool simplex;
211: PetscInt dim, cStart;
212: char prefix[PETSC_MAX_PATH_LEN];
215: DMGetDimension(dm, &dim);
216: DMPlexGetHeightStratum(dm, 0, &cStart, NULL);
217: DMPlexGetCellType(dm, cStart, &ct);
218: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;
219: /* Create finite element */
220: PetscSNPrintf(prefix, PETSC_MAX_PATH_LEN, "%s_", name);
221: PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, name ? prefix : NULL, -1, &fe);
222: PetscObjectSetName((PetscObject)fe, name);
223: /* Set discretization and boundary conditions for each mesh */
224: DMSetField(dm, 0, NULL, (PetscObject)fe);
225: DMCreateDS(dm);
226: (*setup)(dm, user);
227: while (cdm) {
228: DMCopyDisc(dm, cdm);
229: /* TODO: Check whether the boundary of coarse meshes is marked */
230: DMGetCoarseDM(cdm, &cdm);
231: }
232: PetscFEDestroy(&fe);
233: return 0;
234: }
236: static PetscErrorCode ComputeSpectral(DM dm, Vec u, PetscInt numPlanes, const PetscInt planeDir[], const PetscReal planeCoord[], AppCtx *user)
237: {
238: MPI_Comm comm;
239: PetscSection coordSection, section;
240: Vec coordinates, uloc;
241: const PetscScalar *coords, *array;
242: PetscInt p;
243: PetscMPIInt size, rank;
246: PetscObjectGetComm((PetscObject)dm, &comm);
247: MPI_Comm_size(comm, &size);
248: MPI_Comm_rank(comm, &rank);
249: DMGetLocalVector(dm, &uloc);
250: DMGlobalToLocalBegin(dm, u, INSERT_VALUES, uloc);
251: DMGlobalToLocalEnd(dm, u, INSERT_VALUES, uloc);
252: DMPlexInsertBoundaryValues(dm, PETSC_TRUE, uloc, 0.0, NULL, NULL, NULL);
253: VecViewFromOptions(uloc, NULL, "-sol_view");
254: DMGetLocalSection(dm, §ion);
255: VecGetArrayRead(uloc, &array);
256: DMGetCoordinatesLocal(dm, &coordinates);
257: DMGetCoordinateSection(dm, &coordSection);
258: VecGetArrayRead(coordinates, &coords);
259: for (p = 0; p < numPlanes; ++p) {
260: DMLabel label;
261: char name[PETSC_MAX_PATH_LEN];
262: Mat F;
263: Vec x, y;
264: IS stratum;
265: PetscReal *ray, *gray;
266: PetscScalar *rvals, *svals, *gsvals;
267: PetscInt *perm, *nperm;
268: PetscInt n, N, i, j, off, offu;
269: const PetscInt *points;
271: PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "spectral_plane_%" PetscInt_FMT, p);
272: DMGetLabel(dm, name, &label);
273: DMLabelGetStratumIS(label, 1, &stratum);
274: ISGetLocalSize(stratum, &n);
275: ISGetIndices(stratum, &points);
276: PetscMalloc2(n, &ray, n, &svals);
277: for (i = 0; i < n; ++i) {
278: PetscSectionGetOffset(coordSection, points[i], &off);
279: PetscSectionGetOffset(section, points[i], &offu);
280: ray[i] = PetscRealPart(coords[off + ((planeDir[p] + 1) % 2)]);
281: svals[i] = array[offu];
282: }
283: /* Gather the ray data to proc 0 */
284: if (size > 1) {
285: PetscMPIInt *cnt, *displs, p;
287: PetscCalloc2(size, &cnt, size, &displs);
288: MPI_Gather(&n, 1, MPIU_INT, cnt, 1, MPIU_INT, 0, comm);
289: for (p = 1; p < size; ++p) displs[p] = displs[p - 1] + cnt[p - 1];
290: N = displs[size - 1] + cnt[size - 1];
291: PetscMalloc2(N, &gray, N, &gsvals);
292: MPI_Gatherv(ray, n, MPIU_REAL, gray, cnt, displs, MPIU_REAL, 0, comm);
293: MPI_Gatherv(svals, n, MPIU_SCALAR, gsvals, cnt, displs, MPIU_SCALAR, 0, comm);
294: PetscFree2(cnt, displs);
295: } else {
296: N = n;
297: gray = ray;
298: gsvals = svals;
299: }
300: if (rank == 0) {
301: /* Sort point along ray */
302: PetscMalloc2(N, &perm, N, &nperm);
303: for (i = 0; i < N; ++i) perm[i] = i;
304: PetscSortRealWithPermutation(N, gray, perm);
305: /* Count duplicates and squish mapping */
306: nperm[0] = perm[0];
307: for (i = 1, j = 1; i < N; ++i) {
308: if (PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) > PETSC_SMALL) nperm[j++] = perm[i];
309: }
310: /* Create FFT structs */
311: MatCreateFFT(PETSC_COMM_SELF, 1, &j, MATFFTW, &F);
312: MatCreateVecs(F, &x, &y);
313: PetscObjectSetName((PetscObject)y, name);
314: VecGetArray(x, &rvals);
315: for (i = 0, j = 0; i < N; ++i) {
316: if (i > 0 && PetscAbsReal(gray[perm[i]] - gray[perm[i - 1]]) < PETSC_SMALL) continue;
317: rvals[j] = gsvals[nperm[j]];
318: ++j;
319: }
320: PetscFree2(perm, nperm);
321: if (size > 1) PetscFree2(gray, gsvals);
322: VecRestoreArray(x, &rvals);
323: /* Do FFT along the ray */
324: MatMult(F, x, y);
325: /* Chop FFT */
326: VecChop(y, PETSC_SMALL);
327: VecViewFromOptions(x, NULL, "-real_view");
328: VecViewFromOptions(y, NULL, "-fft_view");
329: VecDestroy(&x);
330: VecDestroy(&y);
331: MatDestroy(&F);
332: }
333: ISRestoreIndices(stratum, &points);
334: ISDestroy(&stratum);
335: PetscFree2(ray, svals);
336: }
337: VecRestoreArrayRead(coordinates, &coords);
338: VecRestoreArrayRead(uloc, &array);
339: DMRestoreLocalVector(dm, &uloc);
340: return 0;
341: }
343: int main(int argc, char **argv)
344: {
345: DM dm; /* Problem specification */
346: SNES snes; /* Nonlinear solver */
347: Vec u; /* Solutions */
348: AppCtx user; /* User-defined work context */
351: PetscInitialize(&argc, &argv, NULL, help);
352: ProcessOptions(PETSC_COMM_WORLD, &user);
353: /* Primal system */
354: SNESCreate(PETSC_COMM_WORLD, &snes);
355: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
356: SNESSetDM(snes, dm);
357: SetupDiscretization(dm, "potential", SetupPrimalProblem, &user);
358: DMCreateGlobalVector(dm, &u);
359: VecSet(u, 0.0);
360: PetscObjectSetName((PetscObject)u, "potential");
361: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
362: SNESSetFromOptions(snes);
363: SNESSolve(snes, NULL, u);
364: SNESGetSolution(snes, &u);
365: VecViewFromOptions(u, NULL, "-potential_view");
366: if (user.viewError) {
367: PetscErrorCode (*sol)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
368: void *ctx;
369: PetscDS ds;
370: PetscReal error;
371: PetscInt N;
373: DMGetDS(dm, &ds);
374: PetscDSGetExactSolution(ds, 0, &sol, &ctx);
375: VecGetSize(u, &N);
376: DMComputeL2Diff(dm, 0.0, &sol, &ctx, u, &error);
377: PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g\n", N, (double)error);
378: }
379: if (user.spectral) {
380: PetscInt planeDir[2] = {0, 1};
381: PetscReal planeCoord[2] = {0., 1.};
383: ComputeSpectral(dm, u, 2, planeDir, planeCoord, &user);
384: }
385: /* Adjoint system */
386: if (user.adjoint) {
387: DM dmAdj;
388: SNES snesAdj;
389: Vec uAdj;
391: SNESCreate(PETSC_COMM_WORLD, &snesAdj);
392: PetscObjectSetOptionsPrefix((PetscObject)snesAdj, "adjoint_");
393: DMClone(dm, &dmAdj);
394: SNESSetDM(snesAdj, dmAdj);
395: SetupDiscretization(dmAdj, "adjoint", SetupAdjointProblem, &user);
396: DMCreateGlobalVector(dmAdj, &uAdj);
397: VecSet(uAdj, 0.0);
398: PetscObjectSetName((PetscObject)uAdj, "adjoint");
399: DMPlexSetSNESLocalFEM(dmAdj, &user, &user, &user);
400: SNESSetFromOptions(snesAdj);
401: SNESSolve(snesAdj, NULL, uAdj);
402: SNESGetSolution(snesAdj, &uAdj);
403: VecViewFromOptions(uAdj, NULL, "-adjoint_view");
404: /* Error representation */
405: {
406: DM dmErr, dmErrAux, dms[2];
407: Vec errorEst, errorL2, uErr, uErrLoc, uAdjLoc, uAdjProj;
408: IS *subis;
409: PetscReal errorEstTot, errorL2Norm, errorL2Tot;
410: PetscInt N, i;
411: PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *) = {user.homogeneous ? trig_homogeneous_u : trig_inhomogeneous_u};
412: void (*identity[1])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]) = {f0_identityaux_u};
413: void *ctxs[1] = {0};
415: ctxs[0] = &user;
416: DMClone(dm, &dmErr);
417: SetupDiscretization(dmErr, "error", SetupErrorProblem, &user);
418: DMGetGlobalVector(dmErr, &errorEst);
419: DMGetGlobalVector(dmErr, &errorL2);
420: /* Compute auxiliary data (solution and projection of adjoint solution) */
421: DMGetLocalVector(dmAdj, &uAdjLoc);
422: DMGlobalToLocalBegin(dmAdj, uAdj, INSERT_VALUES, uAdjLoc);
423: DMGlobalToLocalEnd(dmAdj, uAdj, INSERT_VALUES, uAdjLoc);
424: DMGetGlobalVector(dm, &uAdjProj);
425: DMSetAuxiliaryVec(dm, NULL, 0, 0, uAdjLoc);
426: DMProjectField(dm, 0.0, u, identity, INSERT_VALUES, uAdjProj);
427: DMSetAuxiliaryVec(dm, NULL, 0, 0, NULL);
428: DMRestoreLocalVector(dmAdj, &uAdjLoc);
429: /* Attach auxiliary data */
430: dms[0] = dm;
431: dms[1] = dm;
432: DMCreateSuperDM(dms, 2, &subis, &dmErrAux);
433: if (0) {
434: PetscSection sec;
436: DMGetLocalSection(dms[0], &sec);
437: PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
438: DMGetLocalSection(dms[1], &sec);
439: PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
440: DMGetLocalSection(dmErrAux, &sec);
441: PetscSectionView(sec, PETSC_VIEWER_STDOUT_WORLD);
442: }
443: DMViewFromOptions(dmErrAux, NULL, "-dm_err_view");
444: ISViewFromOptions(subis[0], NULL, "-super_is_view");
445: ISViewFromOptions(subis[1], NULL, "-super_is_view");
446: DMGetGlobalVector(dmErrAux, &uErr);
447: VecViewFromOptions(u, NULL, "-map_vec_view");
448: VecViewFromOptions(uAdjProj, NULL, "-map_vec_view");
449: VecViewFromOptions(uErr, NULL, "-map_vec_view");
450: VecISCopy(uErr, subis[0], SCATTER_FORWARD, u);
451: VecISCopy(uErr, subis[1], SCATTER_FORWARD, uAdjProj);
452: DMRestoreGlobalVector(dm, &uAdjProj);
453: for (i = 0; i < 2; ++i) ISDestroy(&subis[i]);
454: PetscFree(subis);
455: DMGetLocalVector(dmErrAux, &uErrLoc);
456: DMGlobalToLocalBegin(dm, uErr, INSERT_VALUES, uErrLoc);
457: DMGlobalToLocalEnd(dm, uErr, INSERT_VALUES, uErrLoc);
458: DMRestoreGlobalVector(dmErrAux, &uErr);
459: DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, uErrLoc);
460: /* Compute cellwise error estimate */
461: VecSet(errorEst, 0.0);
462: DMPlexComputeCellwiseIntegralFEM(dmAdj, uAdj, errorEst, &user);
463: DMSetAuxiliaryVec(dmAdj, NULL, 0, 0, NULL);
464: DMRestoreLocalVector(dmErrAux, &uErrLoc);
465: DMDestroy(&dmErrAux);
466: /* Plot cellwise error vector */
467: VecViewFromOptions(errorEst, NULL, "-error_view");
468: /* Compute ratio of estimate (sum over cells) with actual L_2 error */
469: DMComputeL2Diff(dm, 0.0, funcs, ctxs, u, &errorL2Norm);
470: DMPlexComputeL2DiffVec(dm, 0.0, funcs, ctxs, u, errorL2);
471: VecViewFromOptions(errorL2, NULL, "-l2_error_view");
472: VecNorm(errorL2, NORM_INFINITY, &errorL2Tot);
473: VecNorm(errorEst, NORM_INFINITY, &errorEstTot);
474: VecGetSize(errorEst, &N);
475: VecPointwiseDivide(errorEst, errorEst, errorL2);
476: PetscObjectSetName((PetscObject)errorEst, "Error ratio");
477: VecViewFromOptions(errorEst, NULL, "-error_ratio_view");
478: PetscPrintf(PETSC_COMM_WORLD, "N: %" PetscInt_FMT " L2 error: %g Error Ratio: %g/%g = %g\n", N, (double)errorL2Norm, (double)errorEstTot, (double)PetscSqrtReal(errorL2Tot), (double)(errorEstTot / PetscSqrtReal(errorL2Tot)));
479: DMRestoreGlobalVector(dmErr, &errorEst);
480: DMRestoreGlobalVector(dmErr, &errorL2);
481: DMDestroy(&dmErr);
482: }
483: DMDestroy(&dmAdj);
484: VecDestroy(&uAdj);
485: SNESDestroy(&snesAdj);
486: }
487: /* Cleanup */
488: VecDestroy(&u);
489: SNESDestroy(&snes);
490: DMDestroy(&dm);
491: PetscFinalize();
492: return 0;
493: }
495: /*TEST
497: test:
498: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
499: suffix: 2d_p1_conv
500: requires: triangle
501: args: -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
502: test:
503: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
504: suffix: 2d_p2_conv
505: requires: triangle
506: args: -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
507: test:
508: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
509: suffix: 2d_p3_conv
510: requires: triangle
511: args: -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
512: test:
513: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
514: suffix: 2d_q1_conv
515: args: -dm_plex_simplex 0 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
516: test:
517: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
518: suffix: 2d_q2_conv
519: args: -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
520: test:
521: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
522: suffix: 2d_q3_conv
523: args: -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
524: test:
525: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.9
526: suffix: 2d_q1_shear_conv
527: args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 2
528: test:
529: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.9
530: suffix: 2d_q2_shear_conv
531: args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 2
532: test:
533: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 3.9
534: suffix: 2d_q3_shear_conv
535: args: -dm_plex_simplex 0 -shear -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 2
536: test:
537: # Using -convest_num_refine 3 we get L_2 convergence rate: 1.7
538: suffix: 3d_p1_conv
539: requires: ctetgen
540: args: -dm_plex_dim 3 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
541: test:
542: # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 2.8
543: suffix: 3d_p2_conv
544: requires: ctetgen
545: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
546: test:
547: # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 4.0
548: suffix: 3d_p3_conv
549: requires: ctetgen
550: args: -dm_plex_dim 3 -dm_plex_box_faces 2,2,2 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
551: test:
552: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 1.8
553: suffix: 3d_q1_conv
554: args: -dm_plex_dim 3 -dm_plex_simplex 0 -dm_refine 1 -potential_petscspace_degree 1 -snes_convergence_estimate -convest_num_refine 1
555: test:
556: # Using -dm_refine 2 -convest_num_refine 3 we get L_2 convergence rate: 2.8
557: suffix: 3d_q2_conv
558: args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 2 -snes_convergence_estimate -convest_num_refine 1
559: test:
560: # Using -dm_refine 1 -convest_num_refine 3 we get L_2 convergence rate: 3.8
561: suffix: 3d_q3_conv
562: args: -dm_plex_dim 3 -dm_plex_simplex 0 -potential_petscspace_degree 3 -snes_convergence_estimate -convest_num_refine 1
563: test:
564: suffix: 2d_p1_fas_full
565: requires: triangle
566: args: -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
567: -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full -snes_fas_full_total \
568: -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
569: -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
570: -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
571: -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10
572: test:
573: suffix: 2d_p1_fas_full_homogeneous
574: requires: triangle
575: args: -homogeneous -potential_petscspace_degree 1 -dm_refine_hierarchy 5 \
576: -snes_max_it 1 -snes_type fas -snes_fas_levels 5 -snes_fas_type full \
577: -fas_coarse_snes_monitor -fas_coarse_snes_max_it 1 -fas_coarse_ksp_atol 1.e-13 \
578: -fas_levels_snes_monitor -fas_levels_snes_max_it 1 -fas_levels_snes_type newtonls \
579: -fas_levels_pc_type none -fas_levels_ksp_max_it 2 -fas_levels_ksp_converged_maxits -fas_levels_ksp_type chebyshev \
580: -fas_levels_esteig_ksp_type cg -fas_levels_ksp_chebyshev_esteig 0,0.25,0,1.1 -fas_levels_esteig_ksp_max_it 10
582: test:
583: suffix: 2d_p1_scalable
584: requires: triangle
585: args: -potential_petscspace_degree 1 -dm_refine 3 \
586: -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned \
587: -pc_type gamg -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \
588: -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 \
589: -pc_gamg_coarse_eq_limit 1000 \
590: -pc_gamg_threshold 0.05 \
591: -pc_gamg_threshold_scale .0 \
592: -mg_levels_ksp_type chebyshev \
593: -mg_levels_ksp_max_it 1 \
594: -mg_levels_pc_type jacobi \
595: -matptap_via scalable
596: test:
597: suffix: 2d_p1_gmg_vcycle
598: requires: triangle
599: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
600: -ksp_rtol 5e-10 -pc_type mg \
601: -mg_levels_ksp_max_it 1 \
602: -mg_levels_esteig_ksp_type cg \
603: -mg_levels_esteig_ksp_max_it 10 \
604: -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
605: -mg_levels_pc_type jacobi
606: test:
607: suffix: 2d_p1_gmg_fcycle
608: requires: triangle
609: args: -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
610: -ksp_rtol 5e-10 -pc_type mg -pc_mg_type full \
611: -mg_levels_ksp_max_it 2 \
612: -mg_levels_esteig_ksp_type cg \
613: -mg_levels_esteig_ksp_max_it 10 \
614: -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
615: -mg_levels_pc_type jacobi
616: test:
617: suffix: 2d_p1_gmg_vcycle_adapt
618: requires: triangle
619: args: -petscpartitioner_type simple -potential_petscspace_degree 1 -dm_plex_box_faces 2,2 -dm_refine_hierarchy 3 \
620: -ksp_rtol 5e-10 -pc_type mg -pc_mg_galerkin -pc_mg_adapt_interp_coarse_space harmonic -pc_mg_adapt_interp_n 8 \
621: -mg_levels_ksp_max_it 1 \
622: -mg_levels_esteig_ksp_type cg \
623: -mg_levels_esteig_ksp_max_it 10 \
624: -mg_levels_ksp_chebyshev_esteig 0,0.1,0,1.1 \
625: -mg_levels_pc_type jacobi
626: test:
627: suffix: 2d_p1_spectral_0
628: requires: triangle fftw !complex
629: args: -dm_plex_box_faces 1,1 -potential_petscspace_degree 1 -dm_refine 6 -spectral -fft_view
630: test:
631: suffix: 2d_p1_spectral_1
632: requires: triangle fftw !complex
633: nsize: 2
634: args: -dm_plex_box_faces 4,4 -potential_petscspace_degree 1 -spectral -fft_view
635: test:
636: suffix: 2d_p1_adj_0
637: requires: triangle
638: args: -potential_petscspace_degree 1 -dm_refine 1 -adjoint -adjoint_petscspace_degree 1 -error_petscspace_degree 0
639: test:
640: nsize: 2
641: requires: kokkos_kernels
642: suffix: kokkos
643: args: -dm_plex_dim 3 -dm_plex_box_faces 2,3,6 -petscpartitioner_type simple -dm_plex_simplex 0 -potential_petscspace_degree 1 \
644: -dm_refine 0 -ksp_type cg -ksp_rtol 1.e-11 -ksp_norm_type unpreconditioned -pc_type gamg -pc_gamg_coarse_eq_limit 1000 -pc_gamg_threshold 0.0 \
645: -pc_gamg_threshold_scale .5 -mg_levels_ksp_type chebyshev -mg_levels_ksp_max_it 2 -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10 \
646: -ksp_monitor -snes_monitor -dm_view -dm_mat_type aijkokkos -dm_vec_type kokkos
648: TEST*/