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: #ifdef PETSC_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 defined(PETSC_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: #endif
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 *lnorms, *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(PetscCalloc2(numFields, &lnorms, 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) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
191: }
192: }
193: PetscCall(VecRestoreArrayRead(res, &r));
194: PetscCallMPI(MPIU_Allreduce(lnorms, 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(PetscFree2(lnorms, 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: /*@
326: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
328: Input Parameters:
329: + dm - The mesh
330: . X - Local solution
331: - ctx - The application context
333: Output Parameter:
334: . F - Local output vector
336: Level: developer
338: Note:
339: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
341: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
342: @*/
343: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, PetscCtx ctx)
344: {
345: DM plex;
346: IS allcellIS;
347: PetscInt Nds, s;
349: PetscFunctionBegin;
350: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
351: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
352: PetscCall(DMGetNumDS(dm, &Nds));
353: for (s = 0; s < Nds; ++s) {
354: PetscDS ds;
355: IS cellIS;
356: PetscFormKey key;
358: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
359: key.value = 0;
360: key.field = 0;
361: key.part = 0;
362: if (!key.label) {
363: PetscCall(PetscObjectReference((PetscObject)allcellIS));
364: cellIS = allcellIS;
365: } else {
366: IS pointIS;
368: key.value = 1;
369: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
370: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
371: PetscCall(ISDestroy(&pointIS));
372: }
373: PetscCall(DMPlexComputeResidualByKey(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, ctx));
374: PetscCall(ISDestroy(&cellIS));
375: }
376: PetscCall(ISDestroy(&allcellIS));
377: PetscCall(DMDestroy(&plex));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*@
382: DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
384: Input Parameters:
385: + dm - The mesh
386: . X - Local solution
387: - ctx - The application context
389: Output Parameter:
390: . F - Local output vector
392: Level: developer
394: Note:
395: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
397: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
398: @*/
399: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, PetscCtx ctx)
400: {
401: DM plex;
402: IS allcellIS;
403: PetscInt Nds, s;
405: PetscFunctionBegin;
406: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
407: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
408: PetscCall(DMGetNumDS(dm, &Nds));
409: for (s = 0; s < Nds; ++s) {
410: PetscDS ds;
411: DMLabel label;
412: IS cellIS;
414: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
415: {
416: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
417: PetscWeakForm wf;
418: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
419: PetscFormKey *reskeys;
421: /* Get unique residual keys */
422: for (m = 0; m < Nm; ++m) {
423: PetscInt Nkm;
424: PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
425: Nk += Nkm;
426: }
427: PetscCall(PetscMalloc1(Nk, &reskeys));
428: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
429: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
430: PetscCall(PetscFormKeySort(Nk, reskeys));
431: for (k = 0, kp = 1; kp < Nk; ++kp) {
432: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
433: ++k;
434: if (kp != k) reskeys[k] = reskeys[kp];
435: }
436: }
437: Nk = k;
439: PetscCall(PetscDSGetWeakForm(ds, &wf));
440: for (k = 0; k < Nk; ++k) {
441: DMLabel label = reskeys[k].label;
442: PetscInt val = reskeys[k].value;
444: if (!label) {
445: PetscCall(PetscObjectReference((PetscObject)allcellIS));
446: cellIS = allcellIS;
447: } else {
448: IS pointIS;
450: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
451: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
452: PetscCall(ISDestroy(&pointIS));
453: }
454: PetscCall(DMPlexComputeResidualByKey(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, ctx));
455: PetscCall(ISDestroy(&cellIS));
456: }
457: PetscCall(PetscFree(reskeys));
458: }
459: }
460: PetscCall(ISDestroy(&allcellIS));
461: PetscCall(DMDestroy(&plex));
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
468: Input Parameters:
469: + dm - The mesh
470: - ctx - The application context
472: Output Parameter:
473: . X - Local solution
475: Level: developer
477: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
478: @*/
479: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, PetscCtx ctx)
480: {
481: DM plex;
483: PetscFunctionBegin;
484: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
485: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
486: PetscCall(DMDestroy(&plex));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
493: Input Parameters:
494: + dm - The `DM`
495: . X - Local solution vector
496: . Y - Local input vector
497: - ctx - The application context
499: Output Parameter:
500: . F - local output vector
502: Level: developer
504: Note:
505: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
507: This only works with `DMPLEX`
509: Developer Note:
510: This should be called `DMPlexSNESComputeJacobianAction()`
512: .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
513: @*/
514: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, PetscCtx ctx)
515: {
516: DM plex;
517: IS allcellIS;
518: PetscInt Nds, s;
520: PetscFunctionBegin;
521: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
522: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
523: PetscCall(DMGetNumDS(dm, &Nds));
524: for (s = 0; s < Nds; ++s) {
525: PetscDS ds;
526: DMLabel label;
527: IS cellIS;
529: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
530: {
531: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
532: PetscWeakForm wf;
533: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
534: PetscFormKey *jackeys;
536: /* Get unique Jacobian keys */
537: for (m = 0; m < Nm; ++m) {
538: PetscInt Nkm;
539: PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
540: Nk += Nkm;
541: }
542: PetscCall(PetscMalloc1(Nk, &jackeys));
543: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
544: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
545: PetscCall(PetscFormKeySort(Nk, jackeys));
546: for (k = 0, kp = 1; kp < Nk; ++kp) {
547: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
548: ++k;
549: if (kp != k) jackeys[k] = jackeys[kp];
550: }
551: }
552: Nk = k;
554: PetscCall(PetscDSGetWeakForm(ds, &wf));
555: for (k = 0; k < Nk; ++k) {
556: DMLabel label = jackeys[k].label;
557: PetscInt val = jackeys[k].value;
559: if (!label) {
560: PetscCall(PetscObjectReference((PetscObject)allcellIS));
561: cellIS = allcellIS;
562: } else {
563: IS pointIS;
565: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
566: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
567: PetscCall(ISDestroy(&pointIS));
568: }
569: PetscCall(DMPlexComputeJacobianActionByKey(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, ctx));
570: PetscCall(ISDestroy(&cellIS));
571: }
572: PetscCall(PetscFree(jackeys));
573: }
574: }
575: PetscCall(ISDestroy(&allcellIS));
576: PetscCall(DMDestroy(&plex));
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
583: Input Parameters:
584: + dm - The `DM`
585: . X - Local input vector
586: - ctx - The application context
588: Output Parameters:
589: + Jac - Jacobian matrix
590: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
592: Level: developer
594: Note:
595: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
596: like a GPU, or vectorize on a multicore machine.
598: .seealso: [](ch_snes), `DMPLEX`, `Mat`
599: @*/
600: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, PetscCtx ctx)
601: {
602: DM plex;
603: IS allcellIS;
604: PetscBool hasJac, hasPrec;
605: PetscInt Nds;
607: PetscFunctionBegin;
608: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
609: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
610: PetscCall(DMGetNumDS(dm, &Nds));
611: for (PetscInt s = 0; s < Nds; ++s) {
612: PetscDS ds;
613: IS cellIS;
614: PetscFormKey key;
616: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
617: key.value = 0;
618: key.field = 0;
619: key.part = 0;
620: if (!key.label) {
621: PetscCall(PetscObjectReference((PetscObject)allcellIS));
622: cellIS = allcellIS;
623: } else {
624: IS pointIS;
626: key.value = 1;
627: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
628: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
629: PetscCall(ISDestroy(&pointIS));
630: }
631: if (!s) {
632: PetscCall(PetscDSHasJacobian(ds, &hasJac));
633: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
634: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
635: PetscCall(MatZeroEntries(JacP));
636: }
637: PetscCall(DMPlexComputeJacobianByKey(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, ctx));
638: PetscCall(ISDestroy(&cellIS));
639: }
640: PetscCall(ISDestroy(&allcellIS));
641: PetscCall(DMDestroy(&plex));
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: struct _DMSNESJacobianMFCtx {
646: DM dm;
647: Vec X;
648: PetscCtx ctx;
649: };
651: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
652: {
653: struct _DMSNESJacobianMFCtx *ctx;
655: PetscFunctionBegin;
656: PetscCall(MatShellGetContext(A, &ctx));
657: PetscCall(MatShellSetContext(A, NULL));
658: PetscCall(DMDestroy(&ctx->dm));
659: PetscCall(VecDestroy(&ctx->X));
660: PetscCall(PetscFree(ctx));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
665: {
666: struct _DMSNESJacobianMFCtx *ctx;
668: PetscFunctionBegin;
669: PetscCall(MatShellGetContext(A, &ctx));
670: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*@
675: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
677: Collective
679: Input Parameters:
680: + dm - The `DM`
681: . X - The evaluation point for the Jacobian
682: - ctx - An application context, or `NULL`
684: Output Parameter:
685: . J - The `Mat`
687: Level: advanced
689: Notes:
690: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
692: This only works for `DMPLEX`
694: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
695: @*/
696: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, PetscCtx ctx, Mat *J)
697: {
698: struct _DMSNESJacobianMFCtx *ictx;
699: PetscInt n, N;
701: PetscFunctionBegin;
702: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
703: PetscCall(MatSetType(*J, MATSHELL));
704: PetscCall(VecGetLocalSize(X, &n));
705: PetscCall(VecGetSize(X, &N));
706: PetscCall(MatSetSizes(*J, n, n, N, N));
707: PetscCall(PetscObjectReference((PetscObject)dm));
708: PetscCall(PetscObjectReference((PetscObject)X));
709: PetscCall(PetscMalloc1(1, &ictx));
710: ictx->dm = dm;
711: ictx->X = X;
712: ictx->ctx = ctx;
713: PetscCall(MatShellSetContext(*J, ictx));
714: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (PetscErrorCodeFn *)DMSNESJacobianMF_Destroy_Private));
715: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (PetscErrorCodeFn *)DMSNESJacobianMF_Mult_Private));
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, PetscCtx ctx)
720: {
721: SNES snes;
722: Mat pJ;
723: DM ovldm, origdm;
724: DMSNES sdm;
725: PetscErrorCode (*bfun)(DM, Vec, void *);
726: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
727: void *bctx, *jctx;
729: PetscFunctionBegin;
730: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
731: PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
732: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
733: PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
734: PetscCall(MatGetDM(pJ, &ovldm));
735: PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
736: PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
737: PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
738: PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
739: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
740: if (!snes) {
741: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
742: PetscCall(SNESSetDM(snes, ovldm));
743: PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
744: PetscCall(PetscObjectDereference((PetscObject)snes));
745: }
746: PetscCall(DMGetDMSNES(ovldm, &sdm));
747: PetscCall(VecLockReadPush(X));
748: {
749: PetscCtx ctx;
750: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
751: PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
752: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
753: }
754: PetscCall(VecLockReadPop(X));
755: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
756: {
757: Mat locpJ;
759: PetscCall(MatISGetLocalMat(pJ, &locpJ));
760: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
761: }
762: PetscFunctionReturn(PETSC_SUCCESS);
763: }
765: /*@
766: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
768: Input Parameters:
769: + dm - The `DM` object
770: . use_obj - Use the objective function callback
771: - ctx - The application context that will be passed to pointwise evaluation routines
773: Level: developer
775: .seealso: [](ch_snes), `DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
776: @*/
777: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, PetscCtx ctx)
778: {
779: PetscBool useCeed;
781: PetscFunctionBegin;
782: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
783: PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx));
784: if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx));
785: if (useCeed) {
786: #ifdef PETSC_HAVE_LIBCEED
787: PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx));
788: #else
789: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
790: #endif
791: } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx));
792: PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx));
793: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*@
798: DMSNESCheckDiscretization - Check the discretization error of the exact solution
800: Input Parameters:
801: + snes - the `SNES` object
802: . dm - the `DM`
803: . t - the time
804: . u - a `DM` vector
805: - tol - A tolerance for the check, or -1 to print the results instead
807: Output Parameter:
808: . error - An array which holds the discretization error in each field, or `NULL`
810: Level: developer
812: Note:
813: The user must call `PetscDSSetExactSolution()` beforehand
815: Developer Note:
816: How is this related to `PetscConvEst`?
818: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
819: @*/
820: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
821: {
822: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
823: void **ectxs;
824: PetscReal *err;
825: MPI_Comm comm;
826: PetscInt Nf;
828: PetscFunctionBegin;
832: if (error) PetscAssertPointer(error, 6);
834: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
835: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
837: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
838: PetscCall(DMGetNumFields(dm, &Nf));
839: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
840: {
841: PetscInt Nds;
843: PetscCall(DMGetNumDS(dm, &Nds));
844: for (PetscInt s = 0; s < Nds; ++s) {
845: PetscDS ds;
846: DMLabel label;
847: IS fieldIS;
848: const PetscInt *fields;
849: PetscInt dsNf;
851: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
852: PetscCall(PetscDSGetNumFields(ds, &dsNf));
853: PetscCall(ISGetIndices(fieldIS, &fields));
854: for (PetscInt f = 0; f < dsNf; ++f) {
855: const PetscInt field = fields[f];
856: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
857: }
858: PetscCall(ISRestoreIndices(fieldIS, &fields));
859: }
860: }
861: if (Nf > 1) {
862: PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
863: if (tol >= 0.0) {
864: 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);
865: } else if (error) {
866: for (PetscInt f = 0; f < Nf; ++f) error[f] = err[f];
867: } else {
868: PetscCall(PetscPrintf(comm, "L_2 Error: ["));
869: for (PetscInt f = 0; f < Nf; ++f) {
870: if (f) PetscCall(PetscPrintf(comm, ", "));
871: PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
872: }
873: PetscCall(PetscPrintf(comm, "]\n"));
874: }
875: } else {
876: PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
877: if (tol >= 0.0) {
878: PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
879: } else if (error) {
880: error[0] = err[0];
881: } else {
882: PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
883: }
884: }
885: PetscCall(PetscFree3(exacts, ectxs, err));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: /*@
890: DMSNESCheckResidual - Check the residual of the exact solution
892: Input Parameters:
893: + snes - the `SNES` object
894: . dm - the `DM`
895: . u - a `DM` vector
896: - tol - A tolerance for the check, or -1 to print the results instead
898: Output Parameter:
899: . residual - The residual norm of the exact solution, or `NULL`
901: Level: developer
903: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
904: @*/
905: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
906: {
907: MPI_Comm comm;
908: Vec r;
909: PetscReal res;
911: PetscFunctionBegin;
915: if (residual) PetscAssertPointer(residual, 5);
916: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
917: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
918: PetscCall(VecDuplicate(u, &r));
919: PetscCall(SNESComputeFunction(snes, u, r));
920: PetscCall(VecNorm(r, NORM_2, &res));
921: if (tol >= 0.0) {
922: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
923: } else if (residual) {
924: *residual = res;
925: } else {
926: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
927: PetscCall(VecFilter(r, 1.0e-10));
928: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
929: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
930: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)snes));
931: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
932: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL));
933: }
934: PetscCall(VecDestroy(&r));
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: /*@
939: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
941: Input Parameters:
942: + snes - the `SNES` object
943: . dm - the `DM`
944: . u - a `DM` vector
945: - tol - A tolerance for the check, or -1 to print the results instead
947: Output Parameters:
948: + isLinear - Flag indicaing that the function looks linear, or `NULL`
949: - convRate - The rate of convergence of the linear model, or `NULL`
951: Level: developer
953: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
954: @*/
955: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
956: {
957: MPI_Comm comm;
958: PetscDS ds;
959: Mat J, M;
960: MatNullSpace nullspace;
961: PetscReal slope, intercept;
962: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
964: PetscFunctionBegin;
968: if (isLinear) PetscAssertPointer(isLinear, 5);
969: if (convRate) PetscAssertPointer(convRate, 6);
970: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
971: if (!dm) PetscCall(SNESGetDM(snes, &dm));
972: if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
973: else PetscCall(SNESGetSolution(snes, &u));
974: /* Create and view matrices */
975: PetscCall(DMCreateMatrix(dm, &J));
976: PetscCall(DMGetDS(dm, &ds));
977: PetscCall(PetscDSHasJacobian(ds, &hasJac));
978: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
979: if (hasJac && hasPrec) {
980: PetscCall(DMCreateMatrix(dm, &M));
981: PetscCall(SNESComputeJacobian(snes, u, J, M));
982: PetscCall(PetscObjectSetName((PetscObject)M, "Matrix used to construct preconditioner"));
983: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
984: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
985: PetscCall(MatDestroy(&M));
986: } else {
987: PetscCall(SNESComputeJacobian(snes, u, J, J));
988: }
989: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
990: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
991: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
992: /* Check nullspace */
993: PetscCall(MatGetNullSpace(J, &nullspace));
994: if (nullspace) {
995: PetscBool isNull;
996: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
997: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
998: }
999: /* Taylor test */
1000: {
1001: PetscRandom rand;
1002: Vec du, uhat, r, rhat, df;
1003: PetscReal h;
1004: PetscReal *es, *hs, *errors;
1005: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1006: PetscInt Nv, v;
1008: /* Choose a perturbation direction */
1009: PetscCall(PetscRandomCreate(comm, &rand));
1010: PetscCall(VecDuplicate(u, &du));
1011: PetscCall(VecSetRandom(du, rand));
1012: PetscCall(PetscRandomDestroy(&rand));
1013: PetscCall(VecDuplicate(u, &df));
1014: PetscCall(MatMult(J, du, df));
1015: /* Evaluate residual at u, F(u), save in vector r */
1016: PetscCall(VecDuplicate(u, &r));
1017: PetscCall(SNESComputeFunction(snes, u, r));
1018: /* Look at the convergence of our Taylor approximation as we approach u */
1019: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1020: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1021: PetscCall(VecDuplicate(u, &uhat));
1022: PetscCall(VecDuplicate(u, &rhat));
1023: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1024: PetscCall(VecWAXPY(uhat, h, du, u));
1025: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1026: PetscCall(SNESComputeFunction(snes, uhat, rhat));
1027: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1028: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
1030: es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1031: hs[Nv] = PetscLog10Real(h);
1032: }
1033: PetscCall(VecDestroy(&uhat));
1034: PetscCall(VecDestroy(&rhat));
1035: PetscCall(VecDestroy(&df));
1036: PetscCall(VecDestroy(&r));
1037: PetscCall(VecDestroy(&du));
1038: for (v = 0; v < Nv; ++v) {
1039: if ((tol >= 0) && (errors[v] > tol)) break;
1040: else if (errors[v] > PETSC_SMALL) break;
1041: }
1042: if (v == Nv) isLin = PETSC_TRUE;
1043: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
1044: PetscCall(PetscFree3(es, hs, errors));
1045: /* Slope should be about 2 */
1046: if (tol >= 0) {
1047: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
1048: } else if (isLinear || convRate) {
1049: if (isLinear) *isLinear = isLin;
1050: if (convRate) *convRate = slope;
1051: } else {
1052: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
1053: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
1054: }
1055: }
1056: PetscCall(MatDestroy(&J));
1057: PetscFunctionReturn(PETSC_SUCCESS);
1058: }
1060: static PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1061: {
1062: PetscFunctionBegin;
1063: PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1064: PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1065: PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@
1070: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1072: Input Parameters:
1073: + snes - the `SNES` object
1074: - u - representative `SNES` vector
1076: Level: developer
1078: Note:
1079: The user must call `PetscDSSetExactSolution()` before this call
1081: .seealso: [](ch_snes), `SNES`, `DM`
1082: @*/
1083: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1084: {
1085: DM dm;
1086: Vec sol;
1087: PetscBool check;
1089: PetscFunctionBegin;
1090: PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1091: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1092: PetscCall(SNESGetDM(snes, &dm));
1093: PetscCall(VecDuplicate(u, &sol));
1094: PetscCall(SNESSetSolution(snes, sol));
1095: PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1096: PetscCall(VecDestroy(&sol));
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /*@
1101: DMPlexSetSNESVariableBounds - Compute upper and lower bounds for the solution using pointsie functions from the `PetscDS`
1103: Collective
1105: Input Parameters:
1106: + dm - The `DM` object
1107: - snes - the `SNES` object
1109: Level: intermediate
1111: Notes:
1112: This calls `SNESVISetVariableBounds()` after generating the bounds vectors, so it only applied to `SNESVI` solves.
1114: We project the actual bounds into the current finite element space so that they become more accurate with refinement.
1116: .seealso: `SNESVISetVariableBounds()`, `SNESVI`, [](ch_snes), `DM`
1117: @*/
1118: PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)
1119: {
1120: PetscDS ds;
1121: Vec lb, ub;
1122: PetscSimplePointFn **lfuncs, **ufuncs;
1123: void **lctxs, **uctxs;
1124: PetscBool hasBound, hasLower = PETSC_FALSE, hasUpper = PETSC_FALSE;
1125: PetscInt Nf;
1127: PetscFunctionBegin;
1128: PetscCall(DMHasBound(dm, &hasBound));
1129: if (!hasBound) PetscFunctionReturn(PETSC_SUCCESS);
1130: // TODO Generalize for multiple DSes
1131: PetscCall(DMGetDS(dm, &ds));
1132: PetscCall(PetscDSGetNumFields(ds, &Nf));
1133: PetscCall(PetscMalloc4(Nf, &lfuncs, Nf, &lctxs, Nf, &ufuncs, Nf, &uctxs));
1134: for (PetscInt f = 0; f < Nf; ++f) {
1135: PetscCall(PetscDSGetLowerBound(ds, f, &lfuncs[f], &lctxs[f]));
1136: PetscCall(PetscDSGetUpperBound(ds, f, &ufuncs[f], &uctxs[f]));
1137: if (lfuncs[f]) hasLower = PETSC_TRUE;
1138: if (ufuncs[f]) hasUpper = PETSC_TRUE;
1139: }
1140: PetscCall(DMCreateGlobalVector(dm, &lb));
1141: PetscCall(DMCreateGlobalVector(dm, &ub));
1142: PetscCall(PetscObjectSetName((PetscObject)lb, "Lower Bound"));
1143: PetscCall(PetscObjectSetName((PetscObject)ub, "Upper Bound"));
1144: PetscCall(VecSet(lb, PETSC_NINFINITY));
1145: PetscCall(VecSet(ub, PETSC_INFINITY));
1146: if (hasLower) PetscCall(DMProjectFunction(dm, 0., lfuncs, lctxs, INSERT_VALUES, lb));
1147: if (hasUpper) PetscCall(DMProjectFunction(dm, 0., ufuncs, uctxs, INSERT_VALUES, ub));
1148: PetscCall(DMPlexInsertBounds(dm, PETSC_TRUE, 0., lb));
1149: PetscCall(DMPlexInsertBounds(dm, PETSC_FALSE, 0., ub));
1150: PetscCall(VecViewFromOptions(lb, NULL, "-dm_plex_snes_lb_view"));
1151: PetscCall(VecViewFromOptions(ub, NULL, "-dm_plex_snes_ub_view"));
1152: PetscCall(SNESVISetVariableBounds(snes, lb, ub));
1153: PetscCall(VecDestroy(&lb));
1154: PetscCall(VecDestroy(&ub));
1155: PetscCall(PetscFree4(lfuncs, lctxs, ufuncs, uctxs));
1156: PetscFunctionReturn(PETSC_SUCCESS);
1157: }