Actual source code: dmplexsnes.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petsc/private/petscimpl.h>
5: #include <petsc/private/petscfeimpl.h>
7: #ifdef PETSC_HAVE_LIBCEED
8: #include <petscdmceed.h>
9: #include <petscdmplexceed.h>
10: #endif
12: 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[])
13: {
14: p[0] = u[uOff[1]];
15: }
17: /*
18: SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
19: This is normally used only to evaluate convergence rates for the pressure accurately.
21: Collective
23: Input Parameters:
24: + snes - The `SNES`
25: . pfield - The field number for pressure
26: . nullspace - The pressure nullspace
27: . u - The solution vector
28: - ctx - An optional user context
30: Output Parameter:
31: . u - The solution with a continuum pressure integral of zero
33: Level: developer
35: Note:
36: 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.
38: .seealso: [](ch_snes), `SNESConvergedCorrectPressure()`
39: */
40: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
41: {
42: DM dm;
43: PetscDS ds;
44: const Vec *nullvecs;
45: PetscScalar pintd, *intc, *intn;
46: MPI_Comm comm;
47: PetscInt Nf, Nv;
49: PetscFunctionBegin;
50: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
51: PetscCall(SNESGetDM(snes, &dm));
52: PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
53: PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
54: PetscCall(DMGetDS(dm, &ds));
55: PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private));
56: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs));
57: PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv);
58: PetscCall(VecDot(nullvecs[0], u, &pintd));
59: PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd));
60: PetscCall(PetscDSGetNumFields(ds, &Nf));
61: PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn));
62: PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx));
63: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
64: PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]));
65: #if defined(PETSC_USE_DEBUG)
66: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
67: PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield]));
68: #endif
69: PetscCall(PetscFree2(intc, intn));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@C
74: SNESConvergedCorrectPressure - The regular `SNES` convergence test that, up on convergence, adds a vector in the nullspace
75: to make the continuum integral of the pressure field equal to zero.
77: Logically Collective
79: Input Parameters:
80: + snes - the `SNES` context
81: . it - the iteration (0 indicates before any Newton steps)
82: . xnorm - 2-norm of current iterate
83: . gnorm - 2-norm of current step
84: . f - 2-norm of function at current iterate
85: - ctx - Optional user context
87: Output Parameter:
88: . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`
90: Options Database Key:
91: . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
93: Level: advanced
95: Notes:
96: 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`
97: must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
98: The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
99: 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.
101: Developer Note:
102: This is a total misuse of the `SNES` convergence test handling system. It should be removed. Perhaps a `SNESSetPostSolve()` could
103: be constructed to handle this process.
105: .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`
106: @*/
107: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
108: {
109: PetscBool monitorIntegral = PETSC_FALSE;
111: PetscFunctionBegin;
112: PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx));
113: if (monitorIntegral) {
114: Mat J;
115: Vec u;
116: MatNullSpace nullspace;
117: const Vec *nullvecs;
118: PetscScalar pintd;
120: PetscCall(SNESGetSolution(snes, &u));
121: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
122: PetscCall(MatGetNullSpace(J, &nullspace));
123: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs));
124: PetscCall(VecDot(nullvecs[0], u, &pintd));
125: PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd)));
126: }
127: if (*reason > 0) {
128: Mat J;
129: Vec u;
130: MatNullSpace nullspace;
131: PetscInt pfield = 1;
133: PetscCall(SNESGetSolution(snes, &u));
134: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
135: PetscCall(MatGetNullSpace(J, &nullspace));
136: PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx));
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
142: {
143: PetscBool isPlex;
145: PetscFunctionBegin;
146: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
147: if (isPlex) {
148: *plex = dm;
149: PetscCall(PetscObjectReference((PetscObject)dm));
150: } else {
151: PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
152: if (!*plex) {
153: PetscCall(DMConvert(dm, DMPLEX, plex));
154: PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
155: } else {
156: PetscCall(PetscObjectReference((PetscObject)*plex));
157: }
158: if (copy) {
159: PetscCall(DMCopyDMSNES(dm, *plex));
160: PetscCall(DMCopyAuxiliaryVec(dm, *plex));
161: }
162: }
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /*@C
167: SNESMonitorFields - Monitors the residual for each field separately
169: Collective
171: Input Parameters:
172: + snes - the `SNES` context, must have an attached `DM`
173: . its - iteration number
174: . fgnorm - 2-norm of residual
175: - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`
177: Level: intermediate
179: Note:
180: This routine prints the residual norm at each iteration.
182: .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
183: @*/
184: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
185: {
186: PetscViewer viewer = vf->viewer;
187: Vec res;
188: DM dm;
189: PetscSection s;
190: const PetscScalar *r;
191: PetscReal *lnorms, *norms;
192: PetscInt numFields, f, pStart, pEnd, p;
194: PetscFunctionBegin;
196: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
197: PetscCall(SNESGetDM(snes, &dm));
198: PetscCall(DMGetLocalSection(dm, &s));
199: PetscCall(PetscSectionGetNumFields(s, &numFields));
200: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
201: PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
202: PetscCall(VecGetArrayRead(res, &r));
203: for (p = pStart; p < pEnd; ++p) {
204: for (f = 0; f < numFields; ++f) {
205: PetscInt fdof, foff, d;
207: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
208: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
209: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
210: }
211: }
212: PetscCall(VecRestoreArrayRead(res, &r));
213: PetscCallMPI(MPIU_Allreduce(lnorms, norms, (PetscMPIInt)numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
214: PetscCall(PetscViewerPushFormat(viewer, vf->format));
215: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
216: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
217: for (f = 0; f < numFields; ++f) {
218: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
219: PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
220: }
221: PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
222: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
223: PetscCall(PetscViewerPopFormat(viewer));
224: PetscCall(PetscFree2(lnorms, norms));
225: PetscFunctionReturn(PETSC_SUCCESS);
226: }
228: /********************* SNES callbacks **************************/
230: /*@
231: DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user
233: Input Parameters:
234: + dm - The mesh
235: . X - Local solution
236: - user - The user context
238: Output Parameter:
239: . obj - Local objective value
241: Level: developer
243: .seealso: `DM`, `DMPlexSNESComputeResidualFEM()`
244: @*/
245: PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, void *user)
246: {
247: PetscInt Nf, cellHeight, cStart, cEnd;
248: PetscScalar *cintegral;
250: PetscFunctionBegin;
251: PetscCall(DMGetNumFields(dm, &Nf));
252: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
253: PetscCall(DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd));
254: PetscCall(PetscCalloc1((cEnd - cStart) * Nf, &cintegral));
255: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
256: PetscCall(DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user));
257: /* Sum up values */
258: *obj = 0;
259: for (PetscInt c = cStart; c < cEnd; ++c)
260: for (PetscInt f = 0; f < Nf; ++f) *obj += PetscRealPart(cintegral[(c - cStart) * Nf + f]);
261: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
262: PetscCall(PetscFree(cintegral));
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: /*@
267: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
269: Input Parameters:
270: + dm - The mesh
271: . X - Local solution
272: - user - The user context
274: Output Parameter:
275: . F - Local output vector
277: Level: developer
279: Note:
280: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
282: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
283: @*/
284: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
285: {
286: DM plex;
287: IS allcellIS;
288: PetscInt Nds, s;
290: PetscFunctionBegin;
291: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
292: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
293: PetscCall(DMGetNumDS(dm, &Nds));
294: for (s = 0; s < Nds; ++s) {
295: PetscDS ds;
296: IS cellIS;
297: PetscFormKey key;
299: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
300: key.value = 0;
301: key.field = 0;
302: key.part = 0;
303: if (!key.label) {
304: PetscCall(PetscObjectReference((PetscObject)allcellIS));
305: cellIS = allcellIS;
306: } else {
307: IS pointIS;
309: key.value = 1;
310: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
311: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
312: PetscCall(ISDestroy(&pointIS));
313: }
314: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
315: PetscCall(ISDestroy(&cellIS));
316: }
317: PetscCall(ISDestroy(&allcellIS));
318: PetscCall(DMDestroy(&plex));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*@
323: DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
325: Input Parameters:
326: + dm - The mesh
327: . X - Local solution
328: - user - The user context
330: Output Parameter:
331: . F - Local output vector
333: Level: developer
335: Note:
336: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
338: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
339: @*/
340: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user)
341: {
342: DM plex;
343: IS allcellIS;
344: PetscInt Nds, s;
346: PetscFunctionBegin;
347: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
348: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
349: PetscCall(DMGetNumDS(dm, &Nds));
350: for (s = 0; s < Nds; ++s) {
351: PetscDS ds;
352: DMLabel label;
353: IS cellIS;
355: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
356: {
357: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
358: PetscWeakForm wf;
359: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
360: PetscFormKey *reskeys;
362: /* Get unique residual keys */
363: for (m = 0; m < Nm; ++m) {
364: PetscInt Nkm;
365: PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
366: Nk += Nkm;
367: }
368: PetscCall(PetscMalloc1(Nk, &reskeys));
369: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
370: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
371: PetscCall(PetscFormKeySort(Nk, reskeys));
372: for (k = 0, kp = 1; kp < Nk; ++kp) {
373: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
374: ++k;
375: if (kp != k) reskeys[k] = reskeys[kp];
376: }
377: }
378: Nk = k;
380: PetscCall(PetscDSGetWeakForm(ds, &wf));
381: for (k = 0; k < Nk; ++k) {
382: DMLabel label = reskeys[k].label;
383: PetscInt val = reskeys[k].value;
385: if (!label) {
386: PetscCall(PetscObjectReference((PetscObject)allcellIS));
387: cellIS = allcellIS;
388: } else {
389: IS pointIS;
391: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
392: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
393: PetscCall(ISDestroy(&pointIS));
394: }
395: PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
396: PetscCall(ISDestroy(&cellIS));
397: }
398: PetscCall(PetscFree(reskeys));
399: }
400: }
401: PetscCall(ISDestroy(&allcellIS));
402: PetscCall(DMDestroy(&plex));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: #ifdef PETSC_HAVE_LIBCEED
407: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user)
408: {
409: Ceed ceed;
410: DMCeed sd = dm->dmceed;
411: CeedVector clocX, clocF;
413: PetscFunctionBegin;
414: PetscCall(DMGetCeed(dm, &ceed));
415: PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
416: PetscCall(DMCeedComputeGeometry(dm, sd));
418: PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
419: PetscCall(VecGetCeedVector(locF, ceed, &clocF));
420: PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
421: PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
422: PetscCall(VecRestoreCeedVector(locF, &clocF));
424: {
425: DM_Plex *mesh = (DM_Plex *)dm->data;
427: if (mesh->printFEM) {
428: PetscSection section;
429: Vec locFbc;
430: PetscInt pStart, pEnd, p, maxDof;
431: PetscScalar *zeroes;
433: PetscCall(DMGetLocalSection(dm, §ion));
434: PetscCall(VecDuplicate(locF, &locFbc));
435: PetscCall(VecCopy(locF, locFbc));
436: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
437: PetscCall(PetscSectionGetMaxDof(section, &maxDof));
438: PetscCall(PetscCalloc1(maxDof, &zeroes));
439: for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
440: PetscCall(PetscFree(zeroes));
441: PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
442: PetscCall(VecDestroy(&locFbc));
443: }
444: }
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
447: #endif
449: /*@
450: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
452: Input Parameters:
453: + dm - The mesh
454: - user - The user context
456: Output Parameter:
457: . X - Local solution
459: Level: developer
461: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
462: @*/
463: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
464: {
465: DM plex;
467: PetscFunctionBegin;
468: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
469: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
470: PetscCall(DMDestroy(&plex));
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /*@
475: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
477: Input Parameters:
478: + dm - The `DM`
479: . X - Local solution vector
480: . Y - Local input vector
481: - user - The user context
483: Output Parameter:
484: . F - local output vector
486: Level: developer
488: Note:
489: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
491: This only works with `DMPLEX`
493: Developer Note:
494: This should be called `DMPlexSNESComputeJacobianAction()`
496: .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
497: @*/
498: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
499: {
500: DM plex;
501: IS allcellIS;
502: PetscInt Nds, s;
504: PetscFunctionBegin;
505: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
506: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
507: PetscCall(DMGetNumDS(dm, &Nds));
508: for (s = 0; s < Nds; ++s) {
509: PetscDS ds;
510: DMLabel label;
511: IS cellIS;
513: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
514: {
515: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
516: PetscWeakForm wf;
517: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
518: PetscFormKey *jackeys;
520: /* Get unique Jacobian keys */
521: for (m = 0; m < Nm; ++m) {
522: PetscInt Nkm;
523: PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
524: Nk += Nkm;
525: }
526: PetscCall(PetscMalloc1(Nk, &jackeys));
527: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
528: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
529: PetscCall(PetscFormKeySort(Nk, jackeys));
530: for (k = 0, kp = 1; kp < Nk; ++kp) {
531: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
532: ++k;
533: if (kp != k) jackeys[k] = jackeys[kp];
534: }
535: }
536: Nk = k;
538: PetscCall(PetscDSGetWeakForm(ds, &wf));
539: for (k = 0; k < Nk; ++k) {
540: DMLabel label = jackeys[k].label;
541: PetscInt val = jackeys[k].value;
543: if (!label) {
544: PetscCall(PetscObjectReference((PetscObject)allcellIS));
545: cellIS = allcellIS;
546: } else {
547: IS pointIS;
549: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
550: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
551: PetscCall(ISDestroy(&pointIS));
552: }
553: PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
554: PetscCall(ISDestroy(&cellIS));
555: }
556: PetscCall(PetscFree(jackeys));
557: }
558: }
559: PetscCall(ISDestroy(&allcellIS));
560: PetscCall(DMDestroy(&plex));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: /*@
565: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
567: Input Parameters:
568: + dm - The `DM`
569: . X - Local input vector
570: - user - The user context
572: Output Parameters:
573: + Jac - Jacobian matrix
574: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
576: Level: developer
578: Note:
579: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
580: like a GPU, or vectorize on a multicore machine.
582: .seealso: [](ch_snes), `DMPLEX`, `Mat`
583: @*/
584: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
585: {
586: DM plex;
587: IS allcellIS;
588: PetscBool hasJac, hasPrec;
589: PetscInt Nds, s;
591: PetscFunctionBegin;
592: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
593: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
594: PetscCall(DMGetNumDS(dm, &Nds));
595: for (s = 0; s < Nds; ++s) {
596: PetscDS ds;
597: IS cellIS;
598: PetscFormKey key;
600: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
601: key.value = 0;
602: key.field = 0;
603: key.part = 0;
604: if (!key.label) {
605: PetscCall(PetscObjectReference((PetscObject)allcellIS));
606: cellIS = allcellIS;
607: } else {
608: IS pointIS;
610: key.value = 1;
611: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
612: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
613: PetscCall(ISDestroy(&pointIS));
614: }
615: if (!s) {
616: PetscCall(PetscDSHasJacobian(ds, &hasJac));
617: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
618: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
619: PetscCall(MatZeroEntries(JacP));
620: }
621: PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
622: PetscCall(ISDestroy(&cellIS));
623: }
624: PetscCall(ISDestroy(&allcellIS));
625: PetscCall(DMDestroy(&plex));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: struct _DMSNESJacobianMFCtx {
630: DM dm;
631: Vec X;
632: void *ctx;
633: };
635: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
636: {
637: struct _DMSNESJacobianMFCtx *ctx;
639: PetscFunctionBegin;
640: PetscCall(MatShellGetContext(A, &ctx));
641: PetscCall(MatShellSetContext(A, NULL));
642: PetscCall(DMDestroy(&ctx->dm));
643: PetscCall(VecDestroy(&ctx->X));
644: PetscCall(PetscFree(ctx));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
649: {
650: struct _DMSNESJacobianMFCtx *ctx;
652: PetscFunctionBegin;
653: PetscCall(MatShellGetContext(A, &ctx));
654: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*@
659: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
661: Collective
663: Input Parameters:
664: + dm - The `DM`
665: . X - The evaluation point for the Jacobian
666: - user - A user context, or `NULL`
668: Output Parameter:
669: . J - The `Mat`
671: Level: advanced
673: Notes:
674: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
676: This only works for `DMPLEX`
678: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
679: @*/
680: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
681: {
682: struct _DMSNESJacobianMFCtx *ctx;
683: PetscInt n, N;
685: PetscFunctionBegin;
686: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
687: PetscCall(MatSetType(*J, MATSHELL));
688: PetscCall(VecGetLocalSize(X, &n));
689: PetscCall(VecGetSize(X, &N));
690: PetscCall(MatSetSizes(*J, n, n, N, N));
691: PetscCall(PetscObjectReference((PetscObject)dm));
692: PetscCall(PetscObjectReference((PetscObject)X));
693: PetscCall(PetscMalloc1(1, &ctx));
694: ctx->dm = dm;
695: ctx->X = X;
696: ctx->ctx = user;
697: PetscCall(MatShellSetContext(*J, ctx));
698: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
699: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
704: {
705: SNES snes;
706: Mat pJ;
707: DM ovldm, origdm;
708: DMSNES sdm;
709: PetscErrorCode (*bfun)(DM, Vec, void *);
710: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
711: void *bctx, *jctx;
713: PetscFunctionBegin;
714: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
715: PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
716: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
717: PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
718: PetscCall(MatGetDM(pJ, &ovldm));
719: PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
720: PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
721: PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
722: PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
723: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
724: if (!snes) {
725: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
726: PetscCall(SNESSetDM(snes, ovldm));
727: PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
728: PetscCall(PetscObjectDereference((PetscObject)snes));
729: }
730: PetscCall(DMGetDMSNES(ovldm, &sdm));
731: PetscCall(VecLockReadPush(X));
732: {
733: void *ctx;
734: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
735: PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
736: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
737: }
738: PetscCall(VecLockReadPop(X));
739: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
740: {
741: Mat locpJ;
743: PetscCall(MatISGetLocalMat(pJ, &locpJ));
744: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
745: }
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /*@
750: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
752: Input Parameters:
753: + dm - The `DM` object
754: . use_obj - Use the objective function callback
755: - ctx - The user context that will be passed to pointwise evaluation routines
757: Level: developer
759: .seealso: [](ch_snes),`DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
760: @*/
761: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, void *ctx)
762: {
763: PetscBool useCeed;
765: PetscFunctionBegin;
766: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
767: PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx));
768: if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx));
769: if (useCeed) {
770: #ifdef PETSC_HAVE_LIBCEED
771: PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx));
772: #else
773: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
774: #endif
775: } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx));
776: PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx));
777: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: /*@
782: DMSNESCheckDiscretization - Check the discretization error of the exact solution
784: Input Parameters:
785: + snes - the `SNES` object
786: . dm - the `DM`
787: . t - the time
788: . u - a `DM` vector
789: - tol - A tolerance for the check, or -1 to print the results instead
791: Output Parameter:
792: . error - An array which holds the discretization error in each field, or `NULL`
794: Level: developer
796: Note:
797: The user must call `PetscDSSetExactSolution()` beforehand
799: Developer Note:
800: How is this related to `PetscConvEst`?
802: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
803: @*/
804: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
805: {
806: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
807: void **ectxs;
808: PetscReal *err;
809: MPI_Comm comm;
810: PetscInt Nf, f;
812: PetscFunctionBegin;
816: if (error) PetscAssertPointer(error, 6);
818: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
819: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
821: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
822: PetscCall(DMGetNumFields(dm, &Nf));
823: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
824: {
825: PetscInt Nds, s;
827: PetscCall(DMGetNumDS(dm, &Nds));
828: for (s = 0; s < Nds; ++s) {
829: PetscDS ds;
830: DMLabel label;
831: IS fieldIS;
832: const PetscInt *fields;
833: PetscInt dsNf, f;
835: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
836: PetscCall(PetscDSGetNumFields(ds, &dsNf));
837: PetscCall(ISGetIndices(fieldIS, &fields));
838: for (f = 0; f < dsNf; ++f) {
839: const PetscInt field = fields[f];
840: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
841: }
842: PetscCall(ISRestoreIndices(fieldIS, &fields));
843: }
844: }
845: if (Nf > 1) {
846: PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
847: if (tol >= 0.0) {
848: for (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);
849: } else if (error) {
850: for (f = 0; f < Nf; ++f) error[f] = err[f];
851: } else {
852: PetscCall(PetscPrintf(comm, "L_2 Error: ["));
853: for (f = 0; f < Nf; ++f) {
854: if (f) PetscCall(PetscPrintf(comm, ", "));
855: PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
856: }
857: PetscCall(PetscPrintf(comm, "]\n"));
858: }
859: } else {
860: PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
861: if (tol >= 0.0) {
862: PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
863: } else if (error) {
864: error[0] = err[0];
865: } else {
866: PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
867: }
868: }
869: PetscCall(PetscFree3(exacts, ectxs, err));
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: DMSNESCheckResidual - Check the residual of the exact solution
876: Input Parameters:
877: + snes - the `SNES` object
878: . dm - the `DM`
879: . u - a `DM` vector
880: - tol - A tolerance for the check, or -1 to print the results instead
882: Output Parameter:
883: . residual - The residual norm of the exact solution, or `NULL`
885: Level: developer
887: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
888: @*/
889: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
890: {
891: MPI_Comm comm;
892: Vec r;
893: PetscReal res;
895: PetscFunctionBegin;
899: if (residual) PetscAssertPointer(residual, 5);
900: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
901: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
902: PetscCall(VecDuplicate(u, &r));
903: PetscCall(SNESComputeFunction(snes, u, r));
904: PetscCall(VecNorm(r, NORM_2, &res));
905: if (tol >= 0.0) {
906: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
907: } else if (residual) {
908: *residual = res;
909: } else {
910: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
911: PetscCall(VecFilter(r, 1.0e-10));
912: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
913: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
914: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
915: }
916: PetscCall(VecDestroy(&r));
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: /*@
921: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
923: Input Parameters:
924: + snes - the `SNES` object
925: . dm - the `DM`
926: . u - a `DM` vector
927: - tol - A tolerance for the check, or -1 to print the results instead
929: Output Parameters:
930: + isLinear - Flag indicaing that the function looks linear, or `NULL`
931: - convRate - The rate of convergence of the linear model, or `NULL`
933: Level: developer
935: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
936: @*/
937: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
938: {
939: MPI_Comm comm;
940: PetscDS ds;
941: Mat J, M;
942: MatNullSpace nullspace;
943: PetscReal slope, intercept;
944: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
946: PetscFunctionBegin;
950: if (isLinear) PetscAssertPointer(isLinear, 5);
951: if (convRate) PetscAssertPointer(convRate, 6);
952: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
953: if (!dm) PetscCall(SNESGetDM(snes, &dm));
954: if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
955: else PetscCall(SNESGetSolution(snes, &u));
956: /* Create and view matrices */
957: PetscCall(DMCreateMatrix(dm, &J));
958: PetscCall(DMGetDS(dm, &ds));
959: PetscCall(PetscDSHasJacobian(ds, &hasJac));
960: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
961: if (hasJac && hasPrec) {
962: PetscCall(DMCreateMatrix(dm, &M));
963: PetscCall(SNESComputeJacobian(snes, u, J, M));
964: PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
965: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
966: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
967: PetscCall(MatDestroy(&M));
968: } else {
969: PetscCall(SNESComputeJacobian(snes, u, J, J));
970: }
971: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
972: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
973: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
974: /* Check nullspace */
975: PetscCall(MatGetNullSpace(J, &nullspace));
976: if (nullspace) {
977: PetscBool isNull;
978: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
979: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
980: }
981: /* Taylor test */
982: {
983: PetscRandom rand;
984: Vec du, uhat, r, rhat, df;
985: PetscReal h;
986: PetscReal *es, *hs, *errors;
987: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
988: PetscInt Nv, v;
990: /* Choose a perturbation direction */
991: PetscCall(PetscRandomCreate(comm, &rand));
992: PetscCall(VecDuplicate(u, &du));
993: PetscCall(VecSetRandom(du, rand));
994: PetscCall(PetscRandomDestroy(&rand));
995: PetscCall(VecDuplicate(u, &df));
996: PetscCall(MatMult(J, du, df));
997: /* Evaluate residual at u, F(u), save in vector r */
998: PetscCall(VecDuplicate(u, &r));
999: PetscCall(SNESComputeFunction(snes, u, r));
1000: /* Look at the convergence of our Taylor approximation as we approach u */
1001: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1002: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1003: PetscCall(VecDuplicate(u, &uhat));
1004: PetscCall(VecDuplicate(u, &rhat));
1005: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1006: PetscCall(VecWAXPY(uhat, h, du, u));
1007: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1008: PetscCall(SNESComputeFunction(snes, uhat, rhat));
1009: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1010: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
1012: es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1013: hs[Nv] = PetscLog10Real(h);
1014: }
1015: PetscCall(VecDestroy(&uhat));
1016: PetscCall(VecDestroy(&rhat));
1017: PetscCall(VecDestroy(&df));
1018: PetscCall(VecDestroy(&r));
1019: PetscCall(VecDestroy(&du));
1020: for (v = 0; v < Nv; ++v) {
1021: if ((tol >= 0) && (errors[v] > tol)) break;
1022: else if (errors[v] > PETSC_SMALL) break;
1023: }
1024: if (v == Nv) isLin = PETSC_TRUE;
1025: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
1026: PetscCall(PetscFree3(es, hs, errors));
1027: /* Slope should be about 2 */
1028: if (tol >= 0) {
1029: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
1030: } else if (isLinear || convRate) {
1031: if (isLinear) *isLinear = isLin;
1032: if (convRate) *convRate = slope;
1033: } else {
1034: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
1035: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
1036: }
1037: }
1038: PetscCall(MatDestroy(&J));
1039: PetscFunctionReturn(PETSC_SUCCESS);
1040: }
1042: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1043: {
1044: PetscFunctionBegin;
1045: PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1046: PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1047: PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1054: Input Parameters:
1055: + snes - the `SNES` object
1056: - u - representative `SNES` vector
1058: Level: developer
1060: Note:
1061: The user must call `PetscDSSetExactSolution()` before this call
1063: .seealso: [](ch_snes), `SNES`, `DM`
1064: @*/
1065: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1066: {
1067: DM dm;
1068: Vec sol;
1069: PetscBool check;
1071: PetscFunctionBegin;
1072: PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1073: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1074: PetscCall(SNESGetDM(snes, &dm));
1075: PetscCall(VecDuplicate(u, &sol));
1076: PetscCall(SNESSetSolution(snes, sol));
1077: PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1078: PetscCall(VecDestroy(&sol));
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }