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