Actual source code: dmplexsnes.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petscdraw.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/petscfeimpl.h>
8: #if PetscDefined(HAVE_LIBCEED)
9: #include <petscdmceed.h>
10: #include <petscdmplexceed.h>
11: #endif
13: static void pressure_Private(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 p[])
14: {
15: p[0] = u[uOff[1]];
16: }
18: /*
19: SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
20: This is normally used only to evaluate convergence rates for the pressure accurately.
22: Collective
24: Input Parameters:
25: + snes - The `SNES`
26: . pfield - The field number for pressure
27: . nullspace - The pressure nullspace
28: . u - The solution vector
29: - ctx - An optional application context
31: Output Parameter:
32: . u - The solution with a continuum pressure integral of zero
34: Level: developer
36: Note:
37: If int(u) = a and int(n) = b, then int(u - a/b n) = a - a/b b = 0. We assume that the nullspace is a single vector given explicitly.
39: .seealso: [](ch_snes), `SNESConvergedCorrectPressure()`
40: */
41: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, PetscCtx ctx)
42: {
43: DM dm;
44: PetscDS ds;
45: const Vec *nullvecs;
46: PetscScalar pintd, *intc, *intn;
47: MPI_Comm comm;
48: PetscInt Nf, Nv;
50: PetscFunctionBegin;
51: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
52: PetscCall(SNESGetDM(snes, &dm));
53: PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
54: PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
55: PetscCall(DMGetDS(dm, &ds));
56: PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private));
57: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs));
58: PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv);
59: PetscCall(VecDot(nullvecs[0], u, &pintd));
60: PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd));
61: PetscCall(PetscDSGetNumFields(ds, &Nf));
62: PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn));
63: PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx));
64: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
65: PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]));
66: if (PetscDefined(USE_DEBUG)) {
67: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
68: PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield]));
69: }
70: PetscCall(PetscFree2(intc, intn));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*@C
75: SNESConvergedCorrectPressure - The regular `SNES` convergence test that, up on convergence, adds a vector in the nullspace
76: to make the continuum integral of the pressure field equal to zero.
78: Logically Collective
80: Input Parameters:
81: + snes - the `SNES` context
82: . it - the iteration (0 indicates before any Newton steps)
83: . xnorm - 2-norm of current iterate
84: . gnorm - 2-norm of current step
85: . f - 2-norm of function at current iterate
86: - ctx - Optional application context
88: Output Parameter:
89: . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FUNCTION_NANORINF`
91: Options Database Key:
92: . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
94: Level: advanced
96: Notes:
97: In order to use this convergence test, you must set up several PETSc structures. First fields must be added to the `DM`, and a `PetscDS`
98: must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
99: The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
100: Last we must be able to integrate the pressure over the domain, so the `DM` attached to the SNES `must` be a `DMPLEX` at this time.
102: Developer Note:
103: This is a total misuse of the `SNES` convergence test handling system. It should be removed. Perhaps a `SNESSetPostSolve()` could
104: be constructed to handle this process.
106: .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`
107: @*/
108: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, PetscCtx ctx)
109: {
110: PetscBool monitorIntegral = PETSC_FALSE;
112: PetscFunctionBegin;
113: PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx));
114: if (monitorIntegral) {
115: Mat J;
116: Vec u;
117: MatNullSpace nullspace;
118: const Vec *nullvecs;
119: PetscScalar pintd;
121: PetscCall(SNESGetSolution(snes, &u));
122: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
123: PetscCall(MatGetNullSpace(J, &nullspace));
124: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs));
125: PetscCall(VecDot(nullvecs[0], u, &pintd));
126: PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd)));
127: }
128: if (*reason > 0) {
129: Mat J;
130: Vec u;
131: MatNullSpace nullspace;
132: PetscInt pfield = 1;
134: PetscCall(SNESGetSolution(snes, &u));
135: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
136: PetscCall(MatGetNullSpace(J, &nullspace));
137: PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx));
138: }
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
143: {
144: PetscBool isPlex;
146: PetscFunctionBegin;
147: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
148: if (isPlex) {
149: *plex = dm;
150: PetscCall(PetscObjectReference((PetscObject)dm));
151: } else {
152: PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
153: if (!*plex) {
154: PetscCall(DMConvert(dm, DMPLEX, plex));
155: PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
156: } else {
157: PetscCall(PetscObjectReference((PetscObject)*plex));
158: }
159: if (copy) {
160: PetscCall(DMCopyDMSNES(dm, *plex));
161: PetscCall(DMCopyAuxiliaryVec(dm, *plex));
162: }
163: }
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode SNESMonitorFields_ASCII(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewer viewer, PetscViewerFormat format)
168: {
169: Vec res;
170: DM dm;
171: PetscSection s;
172: const PetscScalar *r;
173: PetscReal *norms;
174: PetscInt numFields, f, pStart, pEnd, p;
176: PetscFunctionBegin;
177: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
178: PetscCall(SNESGetDM(snes, &dm));
179: PetscCall(DMGetLocalSection(dm, &s));
180: PetscCall(PetscSectionGetNumFields(s, &numFields));
181: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
182: PetscCall(PetscCalloc1(numFields, &norms));
183: PetscCall(VecGetArrayRead(res, &r));
184: for (p = pStart; p < pEnd; ++p) {
185: for (f = 0; f < numFields; ++f) {
186: PetscInt fdof, foff, d;
188: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
189: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
190: for (d = 0; d < fdof; ++d) norms[f] += PetscRealPart(PetscSqr(r[foff + d]));
191: }
192: }
193: PetscCall(VecRestoreArrayRead(res, &r));
194: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
195: PetscCall(PetscViewerPushFormat(viewer, format));
196: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
197: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
198: for (f = 0; f < numFields; ++f) {
199: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
200: PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
201: }
202: PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
203: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
204: PetscCall(PetscViewerPopFormat(viewer));
205: PetscCall(PetscFree(norms));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode SNESMonitorFields_Draw(SNES snes, PetscInt its, PetscViewer viewer, PetscViewerFormat format)
210: {
211: DM *subdm, dm;
212: Vec *subv, res;
213: IS *subidx;
214: PetscViewer *subview;
215: PetscSection s;
216: PetscInt Nf;
217: const char *prefix;
219: PetscFunctionBegin;
220: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)snes, &prefix));
221: PetscCall(SNESGetDM(snes, &dm));
222: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
223: PetscCall(DMGetLocalSection(dm, &s));
224: PetscCall(PetscSectionGetNumFields(s, &Nf));
225: PetscCall(PetscMalloc4(Nf, &subdm, Nf, &subv, Nf, &subidx, Nf, &subview));
227: for (PetscInt f = 0; f < Nf; ++f) {
228: PetscDraw draw;
229: const char *name;
231: PetscCall(DMCreateSubDM(dm, 1, &f, &subidx[f], &subdm[f]));
232: PetscCall(DMGetGlobalVector(subdm[f], &subv[f]));
233: PetscCall(PetscSectionGetFieldName(s, f, &name));
234: PetscCall(PetscObjectSetName((PetscObject)subv[f], name));
235: PetscCall(VecISCopy(res, subidx[f], SCATTER_REVERSE, subv[f]));
237: PetscCall(PetscViewerDrawOpen(PetscObjectComm((PetscObject)snes), NULL, name, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, PETSC_DETERMINE, &subview[f]));
238: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subview[f], prefix));
239: PetscCall(PetscViewerSetFromOptions(subview[f]));
240: PetscCall(PetscViewerDrawGetDraw(subview[f], 0, &draw));
241: PetscCall(PetscDrawSetPause(draw, -2.));
242: PetscCall(VecView(subv[f], subview[f]));
243: }
245: for (PetscInt f = 0; f < Nf; ++f) {
246: PetscCall(DMRestoreGlobalVector(subdm[f], &subv[f]));
247: PetscCall(ISDestroy(&subidx[f]));
248: PetscCall(DMDestroy(&subdm[f]));
249: PetscCall(PetscViewerDestroy(&subview[f]));
250: }
251: PetscCall(PetscFree4(subdm, subv, subidx, subview));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@C
256: SNESMonitorFields - Monitors the residual norm or draws the residual, for each field separately
258: Collective
260: Input Parameters:
261: + snes - the `SNES` context, must have an attached `DM`
262: . its - iteration number
263: . fgnorm - 2-norm of residual
264: - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII` or `PETSCVIEWERDRAW`
266: Level: intermediate
268: Note:
269: This routine prints the residual norm at each iteration.
271: .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
272: @*/
273: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
274: {
275: PetscViewer viewer = vf->viewer;
276: PetscBool isascii, isdraw;
278: PetscFunctionBegin;
280: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
281: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
282: if (isascii) PetscCall(SNESMonitorFields_ASCII(snes, its, fgnorm, viewer, vf->format));
283: else if (isdraw) PetscCall(SNESMonitorFields_Draw(snes, its, viewer, vf->format));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /********************* SNES callbacks **************************/
289: /*@
290: DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user
292: Input Parameters:
293: + dm - The mesh
294: . X - Local solution
295: - ctx - The application context
297: Output Parameter:
298: . obj - Local objective value
300: Level: developer
302: .seealso: `DM`, `DMPlexSNESComputeResidualFEM()`
303: @*/
304: PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, PetscCtx ctx)
305: {
306: PetscInt Nf, cellHeight, cStart, cEnd;
307: PetscScalar *cintegral;
309: PetscFunctionBegin;
310: PetscCall(DMGetNumFields(dm, &Nf));
311: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
312: PetscCall(DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd));
313: PetscCall(PetscCalloc1((cEnd - cStart) * Nf, &cintegral));
314: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
315: PetscCall(DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, ctx));
316: /* Sum up values */
317: *obj = 0;
318: for (PetscInt c = cStart; c < cEnd; ++c)
319: for (PetscInt f = 0; f < Nf; ++f) *obj += PetscRealPart(cintegral[(c - cStart) * Nf + f]);
320: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
321: PetscCall(PetscFree(cintegral));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: static PetscErrorCode CreateSurfaceCellIS_Private(DM dm, IS *cohesiveCells)
326: {
327: PetscInt cMax, cEnd;
328: PetscMPIInt size;
330: PetscFunctionBegin;
331: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
332: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, NULL, &cMax));
333: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
334: if (size > 1) {
335: PetscSF sf;
336: const PetscInt *leaves;
337: PetscInt *points;
338: PetscInt Nl, l, Ncoh = 0;
340: PetscCall(DMGetPointSF(dm, &sf));
341: PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
342: for (PetscInt c = cMax; c < cEnd; ++c) {
343: if (leaves) PetscCall(PetscFindInt(c, Nl, leaves, &l));
344: else l = (c >= 0 && c < Nl) ? c : -1;
345: if (l < 0) ++Ncoh;
346: }
347: PetscCall(PetscMalloc1(Ncoh, &points));
348: Ncoh = 0;
349: for (PetscInt c = cMax; c < cEnd; ++c) {
350: if (leaves) PetscCall(PetscFindInt(c, Nl, leaves, &l));
351: else l = (c >= 0 && c < Nl) ? c : -1;
352: if (l < 0) points[Ncoh++] = c;
353: }
354: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, Ncoh, points, PETSC_OWN_POINTER, cohesiveCells));
355: } else {
356: PetscCall(ISCreateStride(PETSC_COMM_SELF, cEnd - cMax, cMax, 1, cohesiveCells));
357: }
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*@
362: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `locF` from the local input `locX` using pointwise functions specified by the user
364: Collective
366: Input Parameters:
367: + dm - The mesh
368: . locX - Local solution
369: - ctx - The application context
371: Output Parameter:
372: . locF - Local output vector
374: Level: developer
376: Note:
377: The residual is summed into `locF`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `locF` is intentional.
379: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
380: @*/
381: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec locX, Vec locF, PetscCtx ctx)
382: {
383: DM plex;
384: IS allcellIS;
385: PetscInt Nds, s;
387: PetscFunctionBegin;
388: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
389: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
390: PetscCall(DMGetNumDS(dm, &Nds));
391: for (s = 0; s < Nds; ++s) {
392: PetscDS ds;
393: IS cellIS;
394: PetscWeakForm wf;
395: PetscFormKey key;
396: PetscFormKey keys[3];
397: PetscFormKey *bdf0keys, *bdf1keys;
398: PetscInt bdf0Nk, bdf1Nk;
400: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
401: key.value = 0;
402: key.field = 0;
403: key.part = 0;
404: if (!key.label) {
405: PetscCall(PetscObjectReference((PetscObject)allcellIS));
406: cellIS = allcellIS;
407: } else {
408: IS pointIS;
410: key.value = 1;
411: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
412: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
413: PetscCall(ISDestroy(&pointIS));
414: }
415: PetscCall(DMPlexComputeResidualByKey(plex, key, cellIS, PETSC_MIN_REAL, locX, NULL, 0.0, locF, ctx));
416: PetscCall(ISDestroy(&cellIS));
417: // Hybrid evaluation
418: PetscCall(PetscDSGetWeakForm(ds, &wf));
419: PetscCall(PetscWeakFormGetKeys(wf, PETSC_WF_BDF0, &bdf0Nk, &bdf0keys));
420: PetscCall(PetscWeakFormGetKeys(wf, PETSC_WF_BDF1, &bdf1Nk, &bdf1keys));
421: if (bdf0Nk + bdf1Nk > 1) {
422: IS cohesiveCells;
423: DMLabel label0 = NULL, label1 = NULL;
424: PetscInt value0A = 0, value0B = 0, value1 = 0;
426: // In the future, we need a way to construct the keys for both sides and the surface, and also the cellIS from the label
427: for (PetscInt i = 0; i < bdf0Nk; ++i) {
428: if (bdf0keys[i].field == 0) {
429: if (!label0) {
430: label0 = bdf0keys[i].label;
431: value0A = bdf0keys[i].value;
432: value0B = value0A;
433: } else if (bdf0keys[i].value != value0A) {
434: value0B = bdf0keys[i].value;
435: }
436: }
437: if (bdf0keys[i].field == 1) {
438: label1 = bdf0keys[i].label;
439: value1 = bdf0keys[i].value;
440: }
441: }
442: keys[0].label = label0;
443: keys[0].value = value0A;
444: keys[0].field = 0;
445: keys[0].part = 0;
446: keys[1].label = label0;
447: keys[1].value = value0B;
448: keys[1].field = 0;
449: keys[1].part = 1;
450: keys[2].label = label1;
451: keys[2].value = value1;
452: keys[2].field = 1;
453: keys[2].part = 2;
454: PetscCall(CreateSurfaceCellIS_Private(plex, &cohesiveCells));
455: PetscCall(DMPlexComputeResidualHybridByKey(plex, keys, cohesiveCells, PETSC_MIN_REAL, locX, NULL, 0.0, locF, ctx));
456: PetscCall(ISDestroy(&cohesiveCells));
457: }
458: PetscCall(PetscFree(bdf0keys));
459: PetscCall(PetscFree(bdf1keys));
460: }
461: PetscCall(ISDestroy(&allcellIS));
462: PetscCall(DMDestroy(&plex));
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /*@
467: DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
469: Input Parameters:
470: + dm - The mesh
471: . X - Local solution
472: - ctx - The application context
474: Output Parameter:
475: . F - Local output vector
477: Level: developer
479: Note:
480: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
482: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
483: @*/
484: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, PetscCtx ctx)
485: {
486: DM plex;
487: IS allcellIS;
488: PetscInt Nds, s;
490: PetscFunctionBegin;
491: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
492: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
493: PetscCall(DMGetNumDS(dm, &Nds));
494: for (s = 0; s < Nds; ++s) {
495: PetscDS ds;
496: DMLabel label;
497: IS cellIS;
499: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
500: {
501: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
502: PetscWeakForm wf;
503: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
504: PetscFormKey *reskeys;
506: /* Get unique residual keys */
507: for (m = 0; m < Nm; ++m) {
508: PetscInt Nkm;
509: PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
510: Nk += Nkm;
511: }
512: PetscCall(PetscMalloc1(Nk, &reskeys));
513: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
514: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
515: PetscCall(PetscFormKeySort(Nk, reskeys));
516: for (k = 0, kp = 1; kp < Nk; ++kp) {
517: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
518: ++k;
519: if (kp != k) reskeys[k] = reskeys[kp];
520: }
521: }
522: Nk = k;
524: PetscCall(PetscDSGetWeakForm(ds, &wf));
525: for (k = 0; k < Nk; ++k) {
526: DMLabel label = reskeys[k].label;
527: PetscInt val = reskeys[k].value;
529: if (!label) {
530: PetscCall(PetscObjectReference((PetscObject)allcellIS));
531: cellIS = allcellIS;
532: } else {
533: IS pointIS;
535: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
536: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
537: PetscCall(ISDestroy(&pointIS));
538: }
539: PetscCall(DMPlexComputeResidualByKey(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, ctx));
540: PetscCall(ISDestroy(&cellIS));
541: }
542: PetscCall(PetscFree(reskeys));
543: }
544: }
545: PetscCall(ISDestroy(&allcellIS));
546: PetscCall(DMDestroy(&plex));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: /*@
551: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
553: Input Parameters:
554: + dm - The mesh
555: - ctx - The application context
557: Output Parameter:
558: . X - Local solution
560: Level: developer
562: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
563: @*/
564: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, PetscCtx ctx)
565: {
566: DM plex;
568: PetscFunctionBegin;
569: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
570: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
571: PetscCall(DMDestroy(&plex));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: /*@
576: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
578: Input Parameters:
579: + dm - The `DM`
580: . X - Local solution vector
581: . Y - Local input vector
582: - ctx - The application context
584: Output Parameter:
585: . F - local output vector
587: Level: developer
589: Note:
590: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
592: This only works with `DMPLEX`
594: Developer Note:
595: This should be called `DMPlexSNESComputeJacobianAction()`
597: .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
598: @*/
599: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, PetscCtx ctx)
600: {
601: DM plex;
602: IS allcellIS;
603: PetscInt Nds, s;
605: PetscFunctionBegin;
606: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
607: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
608: PetscCall(DMGetNumDS(dm, &Nds));
609: for (s = 0; s < Nds; ++s) {
610: PetscDS ds;
611: DMLabel label;
612: IS cellIS;
614: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
615: {
616: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
617: PetscWeakForm wf;
618: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
619: PetscFormKey *jackeys;
621: /* Get unique Jacobian keys */
622: for (m = 0; m < Nm; ++m) {
623: PetscInt Nkm;
624: PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
625: Nk += Nkm;
626: }
627: PetscCall(PetscMalloc1(Nk, &jackeys));
628: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
629: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
630: PetscCall(PetscFormKeySort(Nk, jackeys));
631: for (k = 0, kp = 1; kp < Nk; ++kp) {
632: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
633: ++k;
634: if (kp != k) jackeys[k] = jackeys[kp];
635: }
636: }
637: Nk = k;
639: PetscCall(PetscDSGetWeakForm(ds, &wf));
640: for (k = 0; k < Nk; ++k) {
641: DMLabel label = jackeys[k].label;
642: PetscInt val = jackeys[k].value;
644: if (!label) {
645: PetscCall(PetscObjectReference((PetscObject)allcellIS));
646: cellIS = allcellIS;
647: } else {
648: IS pointIS;
650: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
651: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
652: PetscCall(ISDestroy(&pointIS));
653: }
654: PetscCall(DMPlexComputeJacobianActionByKey(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, ctx));
655: PetscCall(ISDestroy(&cellIS));
656: }
657: PetscCall(PetscFree(jackeys));
658: }
659: }
660: PetscCall(ISDestroy(&allcellIS));
661: PetscCall(DMDestroy(&plex));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: /*@
666: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
668: Input Parameters:
669: + dm - The `DM`
670: . X - Local input vector
671: - ctx - The application context
673: Output Parameters:
674: + Jac - Jacobian matrix
675: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
677: Level: developer
679: Note:
680: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
681: like a GPU, or vectorize on a multicore machine.
683: .seealso: [](ch_snes), `DMPLEX`, `Mat`
684: @*/
685: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, PetscCtx ctx)
686: {
687: DM plex;
688: IS allcellIS;
689: PetscBool hasJac, hasPrec;
690: PetscInt Nds;
692: PetscFunctionBegin;
693: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
694: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
695: PetscCall(DMGetNumDS(dm, &Nds));
696: for (PetscInt s = 0; s < Nds; ++s) {
697: PetscDS ds;
698: IS cellIS;
699: PetscFormKey key;
701: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
702: key.value = 0;
703: key.field = 0;
704: key.part = 0;
705: if (!key.label) {
706: PetscCall(PetscObjectReference((PetscObject)allcellIS));
707: cellIS = allcellIS;
708: } else {
709: IS pointIS;
711: key.value = 1;
712: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
713: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
714: PetscCall(ISDestroy(&pointIS));
715: }
716: if (!s) {
717: PetscCall(PetscDSHasJacobian(ds, &hasJac));
718: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
719: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
720: PetscCall(MatZeroEntries(JacP));
721: }
722: PetscCall(DMPlexComputeJacobianByKey(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, ctx));
723: PetscCall(ISDestroy(&cellIS));
724: }
725: PetscCall(ISDestroy(&allcellIS));
726: PetscCall(DMDestroy(&plex));
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: struct _DMSNESJacobianMFCtx {
731: DM dm;
732: Vec X;
733: PetscCtx ctx;
734: };
736: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
737: {
738: struct _DMSNESJacobianMFCtx *ctx;
740: PetscFunctionBegin;
741: PetscCall(MatShellGetContext(A, &ctx));
742: PetscCall(MatShellSetContext(A, NULL));
743: PetscCall(DMDestroy(&ctx->dm));
744: PetscCall(VecDestroy(&ctx->X));
745: PetscCall(PetscFree(ctx));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
750: {
751: struct _DMSNESJacobianMFCtx *ctx;
753: PetscFunctionBegin;
754: PetscCall(MatShellGetContext(A, &ctx));
755: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: /*@
760: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
762: Collective
764: Input Parameters:
765: + dm - The `DM`
766: . X - The evaluation point for the Jacobian
767: - ctx - An application context, or `NULL`
769: Output Parameter:
770: . J - The `Mat`
772: Level: advanced
774: Notes:
775: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
777: This only works for `DMPLEX`
779: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
780: @*/
781: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, PetscCtx ctx, Mat *J)
782: {
783: struct _DMSNESJacobianMFCtx *ictx;
784: PetscInt n, N;
786: PetscFunctionBegin;
787: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
788: PetscCall(MatSetType(*J, MATSHELL));
789: PetscCall(VecGetLocalSize(X, &n));
790: PetscCall(VecGetSize(X, &N));
791: PetscCall(MatSetSizes(*J, n, n, N, N));
792: PetscCall(PetscObjectReference((PetscObject)dm));
793: PetscCall(PetscObjectReference((PetscObject)X));
794: PetscCall(PetscMalloc1(1, &ictx));
795: ictx->dm = dm;
796: ictx->X = X;
797: ictx->ctx = ctx;
798: PetscCall(MatShellSetContext(*J, ictx));
799: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (PetscErrorCodeFn *)DMSNESJacobianMF_Destroy_Private));
800: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (PetscErrorCodeFn *)DMSNESJacobianMF_Mult_Private));
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, PetscCtx ctx)
805: {
806: SNES snes;
807: Mat pJ;
808: DM ovldm, origdm;
809: DMSNES sdm;
810: PetscErrorCode (*bfun)(DM, Vec, void *);
811: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
812: void *bctx, *jctx;
814: PetscFunctionBegin;
815: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
816: PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
817: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
818: PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
819: PetscCall(MatGetDM(pJ, &ovldm));
820: PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
821: PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
822: PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
823: PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
824: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
825: if (!snes) {
826: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
827: PetscCall(SNESSetDM(snes, ovldm));
828: PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
829: PetscCall(PetscObjectDereference((PetscObject)snes));
830: }
831: PetscCall(DMGetDMSNES(ovldm, &sdm));
832: PetscCall(VecLockReadPush(X));
833: {
834: PetscCtx ctx;
835: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
836: PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
837: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
838: }
839: PetscCall(VecLockReadPop(X));
840: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
841: {
842: Mat locpJ;
844: PetscCall(MatISGetLocalMat(pJ, &locpJ));
845: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
846: }
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@
851: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
853: Input Parameters:
854: + dm - The `DM` object
855: . use_obj - Use the objective function callback
856: - ctx - The application context that will be passed to pointwise evaluation routines
858: Level: developer
860: .seealso: [](ch_snes), `DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
861: @*/
862: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, PetscCtx ctx)
863: {
864: PetscBool useCeed;
866: PetscFunctionBegin;
867: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
868: PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx));
869: if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx));
870: if (useCeed) {
871: #if PetscDefined(HAVE_LIBCEED)
872: PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx));
873: #else
874: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
875: #endif
876: } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx));
877: PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx));
878: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: /*@
883: DMSNESCheckDiscretization - Check the discretization error of the exact solution
885: Input Parameters:
886: + snes - the `SNES` object
887: . dm - the `DM`
888: . t - the time
889: . u - a `DM` vector
890: - tol - A tolerance for the check, or -1 to print the results instead
892: Output Parameter:
893: . error - An array which holds the discretization error in each field, or `NULL`
895: Level: developer
897: Note:
898: The user must call `PetscDSSetExactSolution()` beforehand
900: Developer Note:
901: How is this related to `PetscConvEst`?
903: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
904: @*/
905: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
906: {
907: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
908: void **ectxs;
909: PetscReal *err;
910: MPI_Comm comm;
911: PetscInt Nf;
913: PetscFunctionBegin;
917: if (error) PetscAssertPointer(error, 6);
919: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
920: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
922: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
923: PetscCall(DMGetNumFields(dm, &Nf));
924: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
925: {
926: PetscInt Nds;
928: PetscCall(DMGetNumDS(dm, &Nds));
929: for (PetscInt s = 0; s < Nds; ++s) {
930: PetscDS ds;
931: DMLabel label;
932: IS fieldIS;
933: const PetscInt *fields;
934: PetscInt dsNf;
936: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
937: PetscCall(PetscDSGetNumFields(ds, &dsNf));
938: PetscCall(ISGetIndices(fieldIS, &fields));
939: for (PetscInt f = 0; f < dsNf; ++f) {
940: const PetscInt field = fields[f];
941: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
942: }
943: PetscCall(ISRestoreIndices(fieldIS, &fields));
944: }
945: }
946: if (Nf > 1) {
947: PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
948: if (tol >= 0.0) {
949: for (PetscInt f = 0; f < Nf; ++f) PetscCheck(err[f] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %" PetscInt_FMT " exceeds tolerance %g", (double)err[f], f, (double)tol);
950: } else if (error) {
951: for (PetscInt f = 0; f < Nf; ++f) error[f] = err[f];
952: } else {
953: PetscCall(PetscPrintf(comm, "L_2 Error: ["));
954: for (PetscInt f = 0; f < Nf; ++f) {
955: if (f) PetscCall(PetscPrintf(comm, ", "));
956: PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
957: }
958: PetscCall(PetscPrintf(comm, "]\n"));
959: }
960: } else {
961: PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
962: if (tol >= 0.0) {
963: PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
964: } else if (error) {
965: error[0] = err[0];
966: } else {
967: PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
968: }
969: }
970: PetscCall(PetscFree3(exacts, ectxs, err));
971: PetscFunctionReturn(PETSC_SUCCESS);
972: }
974: /*@
975: DMSNESCheckResidual - Check the residual of the exact solution
977: Input Parameters:
978: + snes - the `SNES` object
979: . dm - the `DM`
980: . u - a `DM` vector
981: - tol - A tolerance for the check, or -1 to print the results instead
983: Output Parameter:
984: . residual - The residual norm of the exact solution, or `NULL`
986: Level: developer
988: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
989: @*/
990: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
991: {
992: MPI_Comm comm;
993: Vec r;
994: PetscReal res;
996: PetscFunctionBegin;
1000: if (residual) PetscAssertPointer(residual, 5);
1001: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1002: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1003: PetscCall(VecDuplicate(u, &r));
1004: PetscCall(SNESComputeFunction(snes, u, r));
1005: PetscCall(VecNorm(r, NORM_2, &res));
1006: if (tol >= 0.0) {
1007: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
1008: } else if (residual) {
1009: *residual = res;
1010: } else {
1011: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
1012: PetscCall(VecFilter(r, 1.0e-10));
1013: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
1014: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
1015: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)snes));
1016: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
1017: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL));
1018: }
1019: PetscCall(VecDestroy(&r));
1020: PetscFunctionReturn(PETSC_SUCCESS);
1021: }
1023: /*@
1024: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
1026: Input Parameters:
1027: + snes - the `SNES` object
1028: . dm - the `DM`
1029: . u - a `DM` vector
1030: - tol - A tolerance for the check, or -1 to print the results instead
1032: Output Parameters:
1033: + isLinear - Flag indicaing that the function looks linear, or `NULL`
1034: - convRate - The rate of convergence of the linear model, or `NULL`
1036: Level: developer
1038: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1039: @*/
1040: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1041: {
1042: MPI_Comm comm;
1043: PetscDS ds;
1044: Mat J, M;
1045: MatNullSpace nullspace;
1046: PetscReal slope, intercept;
1047: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
1049: PetscFunctionBegin;
1053: if (isLinear) PetscAssertPointer(isLinear, 5);
1054: if (convRate) PetscAssertPointer(convRate, 6);
1055: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1056: if (!dm) PetscCall(SNESGetDM(snes, &dm));
1057: if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1058: else PetscCall(SNESGetSolution(snes, &u));
1059: /* Create and view matrices */
1060: PetscCall(DMCreateMatrix(dm, &J));
1061: PetscCall(DMGetDS(dm, &ds));
1062: PetscCall(PetscDSHasJacobian(ds, &hasJac));
1063: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1064: if (hasJac && hasPrec) {
1065: PetscCall(DMCreateMatrix(dm, &M));
1066: PetscCall(SNESComputeJacobian(snes, u, J, M));
1067: PetscCall(PetscObjectSetName((PetscObject)M, "Matrix used to construct preconditioner"));
1068: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
1069: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
1070: PetscCall(MatDestroy(&M));
1071: } else {
1072: PetscCall(SNESComputeJacobian(snes, u, J, J));
1073: }
1074: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
1075: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
1076: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
1077: /* Check nullspace */
1078: PetscCall(MatGetNullSpace(J, &nullspace));
1079: if (nullspace) {
1080: PetscBool isNull;
1081: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
1082: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1083: }
1084: /* Taylor test */
1085: {
1086: PetscRandom rand;
1087: Vec du, uhat, r, rhat, df;
1088: PetscReal h;
1089: PetscReal *es, *hs, *errors;
1090: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1091: PetscInt Nv, v;
1093: /* Choose a perturbation direction */
1094: PetscCall(PetscRandomCreate(comm, &rand));
1095: PetscCall(VecDuplicate(u, &du));
1096: PetscCall(VecSetRandom(du, rand));
1097: PetscCall(PetscRandomDestroy(&rand));
1098: PetscCall(VecDuplicate(u, &df));
1099: PetscCall(MatMult(J, du, df));
1100: /* Evaluate residual at u, F(u), save in vector r */
1101: PetscCall(VecDuplicate(u, &r));
1102: PetscCall(SNESComputeFunction(snes, u, r));
1103: /* Look at the convergence of our Taylor approximation as we approach u */
1104: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1105: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1106: PetscCall(VecDuplicate(u, &uhat));
1107: PetscCall(VecDuplicate(u, &rhat));
1108: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1109: PetscCall(VecWAXPY(uhat, h, du, u));
1110: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1111: PetscCall(SNESComputeFunction(snes, uhat, rhat));
1112: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1113: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
1115: es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1116: hs[Nv] = PetscLog10Real(h);
1117: }
1118: PetscCall(VecDestroy(&uhat));
1119: PetscCall(VecDestroy(&rhat));
1120: PetscCall(VecDestroy(&df));
1121: PetscCall(VecDestroy(&r));
1122: PetscCall(VecDestroy(&du));
1123: for (v = 0; v < Nv; ++v) {
1124: if ((tol >= 0) && (errors[v] > tol)) break;
1125: else if (errors[v] > PETSC_SMALL) break;
1126: }
1127: if (v == Nv) isLin = PETSC_TRUE;
1128: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
1129: PetscCall(PetscFree3(es, hs, errors));
1130: /* Slope should be about 2 */
1131: if (tol >= 0) {
1132: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
1133: } else if (isLinear || convRate) {
1134: if (isLinear) *isLinear = isLin;
1135: if (convRate) *convRate = slope;
1136: } else {
1137: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
1138: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
1139: }
1140: }
1141: PetscCall(MatDestroy(&J));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: static PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1146: {
1147: PetscFunctionBegin;
1148: PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1149: PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1150: PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: /*@
1155: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1157: Input Parameters:
1158: + snes - the `SNES` object
1159: - u - representative `SNES` vector
1161: Level: developer
1163: Note:
1164: The user must call `PetscDSSetExactSolution()` before this call
1166: .seealso: [](ch_snes), `SNES`, `DM`
1167: @*/
1168: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1169: {
1170: DM dm;
1171: Vec sol;
1172: PetscBool check;
1174: PetscFunctionBegin;
1175: PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1176: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1177: PetscCall(SNESGetDM(snes, &dm));
1178: PetscCall(VecDuplicate(u, &sol));
1179: PetscCall(SNESSetSolution(snes, sol));
1180: PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1181: PetscCall(VecDestroy(&sol));
1182: PetscFunctionReturn(PETSC_SUCCESS);
1183: }
1185: /*@
1186: DMPlexSetSNESVariableBounds - Compute upper and lower bounds for the solution using pointsie functions from the `PetscDS`
1188: Collective
1190: Input Parameters:
1191: + dm - The `DM` object
1192: - snes - the `SNES` object
1194: Level: intermediate
1196: Notes:
1197: This calls `SNESVISetVariableBounds()` after generating the bounds vectors, so it only applied to `SNESVI` solves.
1199: We project the actual bounds into the current finite element space so that they become more accurate with refinement.
1201: .seealso: `SNESVISetVariableBounds()`, `SNESVI`, [](ch_snes), `DM`
1202: @*/
1203: PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)
1204: {
1205: PetscDS ds;
1206: Vec lb, ub;
1207: PetscSimplePointFn **lfuncs, **ufuncs;
1208: void **lctxs, **uctxs;
1209: PetscBool hasBound, hasLower = PETSC_FALSE, hasUpper = PETSC_FALSE;
1210: PetscInt Nf;
1212: PetscFunctionBegin;
1213: PetscCall(DMHasBound(dm, &hasBound));
1214: if (!hasBound) PetscFunctionReturn(PETSC_SUCCESS);
1215: // TODO Generalize for multiple DSes
1216: PetscCall(DMGetDS(dm, &ds));
1217: PetscCall(PetscDSGetNumFields(ds, &Nf));
1218: PetscCall(PetscMalloc4(Nf, &lfuncs, Nf, &lctxs, Nf, &ufuncs, Nf, &uctxs));
1219: for (PetscInt f = 0; f < Nf; ++f) {
1220: PetscCall(PetscDSGetLowerBound(ds, f, &lfuncs[f], &lctxs[f]));
1221: PetscCall(PetscDSGetUpperBound(ds, f, &ufuncs[f], &uctxs[f]));
1222: if (lfuncs[f]) hasLower = PETSC_TRUE;
1223: if (ufuncs[f]) hasUpper = PETSC_TRUE;
1224: }
1225: PetscCall(DMCreateGlobalVector(dm, &lb));
1226: PetscCall(DMCreateGlobalVector(dm, &ub));
1227: PetscCall(PetscObjectSetName((PetscObject)lb, "Lower Bound"));
1228: PetscCall(PetscObjectSetName((PetscObject)ub, "Upper Bound"));
1229: if (hasLower) {
1230: Vec locb;
1232: PetscCall(DMGetLocalVector(dm, &locb));
1233: PetscCall(VecSet(locb, PETSC_NINFINITY));
1234: PetscCall(DMProjectFunctionLocal(dm, 0., lfuncs, lctxs, INSERT_VALUES, locb));
1235: PetscCall(DMPlexInsertBounds(dm, PETSC_TRUE, 0., locb));
1236: PetscCall(DMLocalToGlobalBegin(dm, locb, INSERT_VALUES, lb));
1237: PetscCall(DMLocalToGlobalEnd(dm, locb, INSERT_VALUES, lb));
1238: PetscCall(DMRestoreLocalVector(dm, &locb));
1239: } else {
1240: PetscCall(VecSet(lb, PETSC_NINFINITY));
1241: }
1242: if (hasUpper) {
1243: Vec locb;
1245: PetscCall(DMGetLocalVector(dm, &locb));
1246: PetscCall(VecSet(locb, PETSC_INFINITY));
1247: PetscCall(DMProjectFunctionLocal(dm, 0., ufuncs, uctxs, INSERT_VALUES, locb));
1248: PetscCall(DMPlexInsertBounds(dm, PETSC_FALSE, 0., locb));
1249: PetscCall(DMLocalToGlobalBegin(dm, locb, INSERT_VALUES, ub));
1250: PetscCall(DMLocalToGlobalEnd(dm, locb, INSERT_VALUES, ub));
1251: PetscCall(DMRestoreLocalVector(dm, &locb));
1252: } else {
1253: PetscCall(VecSet(ub, PETSC_INFINITY));
1254: }
1255: PetscCall(VecViewFromOptions(lb, NULL, "-dm_plex_snes_lb_view"));
1256: PetscCall(VecViewFromOptions(ub, NULL, "-dm_plex_snes_ub_view"));
1257: PetscCall(SNESVISetVariableBounds(snes, lb, ub));
1258: PetscCall(VecDestroy(&lb));
1259: PetscCall(VecDestroy(&ub));
1260: PetscCall(PetscFree4(lfuncs, lctxs, ufuncs, uctxs));
1261: PetscFunctionReturn(PETSC_SUCCESS);
1262: }