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;
193: PetscMPIInt numFieldsi;
195: PetscFunctionBegin;
197: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
198: PetscCall(SNESGetDM(snes, &dm));
199: PetscCall(DMGetLocalSection(dm, &s));
200: PetscCall(PetscSectionGetNumFields(s, &numFields));
201: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
202: PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
203: PetscCall(VecGetArrayRead(res, &r));
204: for (p = pStart; p < pEnd; ++p) {
205: for (f = 0; f < numFields; ++f) {
206: PetscInt fdof, foff, d;
208: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
209: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
210: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
211: }
212: }
213: PetscCall(VecRestoreArrayRead(res, &r));
214: PetscCall(PetscMPIIntCast(numFields, &numFieldsi));
215: PetscCallMPI(MPIU_Allreduce(lnorms, norms, numFieldsi, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
216: PetscCall(PetscViewerPushFormat(viewer, vf->format));
217: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
218: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
219: for (f = 0; f < numFields; ++f) {
220: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
221: PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
222: }
223: PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
224: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
225: PetscCall(PetscViewerPopFormat(viewer));
226: PetscCall(PetscFree2(lnorms, norms));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: /********************* SNES callbacks **************************/
232: /*@
233: DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user
235: Input Parameters:
236: + dm - The mesh
237: . X - Local solution
238: - user - The user context
240: Output Parameter:
241: . obj - Local objective value
243: Level: developer
245: .seealso: `DM`, `DMPlexSNESComputeResidualFEM()`
246: @*/
247: PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, void *user)
248: {
249: PetscInt Nf, cellHeight, cStart, cEnd;
250: PetscScalar *cintegral;
252: PetscFunctionBegin;
253: PetscCall(DMGetNumFields(dm, &Nf));
254: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
255: PetscCall(DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd));
256: PetscCall(PetscCalloc1((cEnd - cStart) * Nf, &cintegral));
257: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
258: PetscCall(DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, user));
259: /* Sum up values */
260: *obj = 0;
261: for (PetscInt c = cStart; c < cEnd; ++c)
262: for (PetscInt f = 0; f < Nf; ++f) *obj += PetscRealPart(cintegral[(c - cStart) * Nf + f]);
263: PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
264: PetscCall(PetscFree(cintegral));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: /*@
269: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
271: Input Parameters:
272: + dm - The mesh
273: . X - Local solution
274: - user - The user context
276: Output Parameter:
277: . F - Local output vector
279: Level: developer
281: Note:
282: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
284: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
285: @*/
286: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
287: {
288: DM plex;
289: IS allcellIS;
290: PetscInt Nds, s;
292: PetscFunctionBegin;
293: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
294: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
295: PetscCall(DMGetNumDS(dm, &Nds));
296: for (s = 0; s < Nds; ++s) {
297: PetscDS ds;
298: IS cellIS;
299: PetscFormKey key;
301: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
302: key.value = 0;
303: key.field = 0;
304: key.part = 0;
305: if (!key.label) {
306: PetscCall(PetscObjectReference((PetscObject)allcellIS));
307: cellIS = allcellIS;
308: } else {
309: IS pointIS;
311: key.value = 1;
312: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
313: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
314: PetscCall(ISDestroy(&pointIS));
315: }
316: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
317: PetscCall(ISDestroy(&cellIS));
318: }
319: PetscCall(ISDestroy(&allcellIS));
320: PetscCall(DMDestroy(&plex));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
327: Input Parameters:
328: + dm - The mesh
329: . X - Local solution
330: - user - The user context
332: Output Parameter:
333: . F - Local output vector
335: Level: developer
337: Note:
338: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
340: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
341: @*/
342: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user)
343: {
344: DM plex;
345: IS allcellIS;
346: PetscInt Nds, s;
348: PetscFunctionBegin;
349: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
350: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
351: PetscCall(DMGetNumDS(dm, &Nds));
352: for (s = 0; s < Nds; ++s) {
353: PetscDS ds;
354: DMLabel label;
355: IS cellIS;
357: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
358: {
359: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
360: PetscWeakForm wf;
361: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
362: PetscFormKey *reskeys;
364: /* Get unique residual keys */
365: for (m = 0; m < Nm; ++m) {
366: PetscInt Nkm;
367: PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
368: Nk += Nkm;
369: }
370: PetscCall(PetscMalloc1(Nk, &reskeys));
371: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
372: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
373: PetscCall(PetscFormKeySort(Nk, reskeys));
374: for (k = 0, kp = 1; kp < Nk; ++kp) {
375: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
376: ++k;
377: if (kp != k) reskeys[k] = reskeys[kp];
378: }
379: }
380: Nk = k;
382: PetscCall(PetscDSGetWeakForm(ds, &wf));
383: for (k = 0; k < Nk; ++k) {
384: DMLabel label = reskeys[k].label;
385: PetscInt val = reskeys[k].value;
387: if (!label) {
388: PetscCall(PetscObjectReference((PetscObject)allcellIS));
389: cellIS = allcellIS;
390: } else {
391: IS pointIS;
393: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
394: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
395: PetscCall(ISDestroy(&pointIS));
396: }
397: PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
398: PetscCall(ISDestroy(&cellIS));
399: }
400: PetscCall(PetscFree(reskeys));
401: }
402: }
403: PetscCall(ISDestroy(&allcellIS));
404: PetscCall(DMDestroy(&plex));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: #ifdef PETSC_HAVE_LIBCEED
409: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user)
410: {
411: Ceed ceed;
412: DMCeed sd = dm->dmceed;
413: CeedVector clocX, clocF;
415: PetscFunctionBegin;
416: PetscCall(DMGetCeed(dm, &ceed));
417: PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
418: PetscCall(DMCeedComputeGeometry(dm, sd));
420: PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
421: PetscCall(VecGetCeedVector(locF, ceed, &clocF));
422: PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
423: PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
424: PetscCall(VecRestoreCeedVector(locF, &clocF));
426: {
427: DM_Plex *mesh = (DM_Plex *)dm->data;
429: if (mesh->printFEM) {
430: PetscSection section;
431: Vec locFbc;
432: PetscInt pStart, pEnd, p, maxDof;
433: PetscScalar *zeroes;
435: PetscCall(DMGetLocalSection(dm, §ion));
436: PetscCall(VecDuplicate(locF, &locFbc));
437: PetscCall(VecCopy(locF, locFbc));
438: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
439: PetscCall(PetscSectionGetMaxDof(section, &maxDof));
440: PetscCall(PetscCalloc1(maxDof, &zeroes));
441: for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
442: PetscCall(PetscFree(zeroes));
443: PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
444: PetscCall(VecDestroy(&locFbc));
445: }
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
449: #endif
451: /*@
452: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
454: Input Parameters:
455: + dm - The mesh
456: - user - The user context
458: Output Parameter:
459: . X - Local solution
461: Level: developer
463: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
464: @*/
465: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
466: {
467: DM plex;
469: PetscFunctionBegin;
470: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
471: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
472: PetscCall(DMDestroy(&plex));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: /*@
477: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
479: Input Parameters:
480: + dm - The `DM`
481: . X - Local solution vector
482: . Y - Local input vector
483: - user - The user context
485: Output Parameter:
486: . F - local output vector
488: Level: developer
490: Note:
491: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
493: This only works with `DMPLEX`
495: Developer Note:
496: This should be called `DMPlexSNESComputeJacobianAction()`
498: .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
499: @*/
500: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
501: {
502: DM plex;
503: IS allcellIS;
504: PetscInt Nds, s;
506: PetscFunctionBegin;
507: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
508: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
509: PetscCall(DMGetNumDS(dm, &Nds));
510: for (s = 0; s < Nds; ++s) {
511: PetscDS ds;
512: DMLabel label;
513: IS cellIS;
515: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
516: {
517: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
518: PetscWeakForm wf;
519: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
520: PetscFormKey *jackeys;
522: /* Get unique Jacobian keys */
523: for (m = 0; m < Nm; ++m) {
524: PetscInt Nkm;
525: PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
526: Nk += Nkm;
527: }
528: PetscCall(PetscMalloc1(Nk, &jackeys));
529: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
530: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
531: PetscCall(PetscFormKeySort(Nk, jackeys));
532: for (k = 0, kp = 1; kp < Nk; ++kp) {
533: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
534: ++k;
535: if (kp != k) jackeys[k] = jackeys[kp];
536: }
537: }
538: Nk = k;
540: PetscCall(PetscDSGetWeakForm(ds, &wf));
541: for (k = 0; k < Nk; ++k) {
542: DMLabel label = jackeys[k].label;
543: PetscInt val = jackeys[k].value;
545: if (!label) {
546: PetscCall(PetscObjectReference((PetscObject)allcellIS));
547: cellIS = allcellIS;
548: } else {
549: IS pointIS;
551: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
552: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
553: PetscCall(ISDestroy(&pointIS));
554: }
555: PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
556: PetscCall(ISDestroy(&cellIS));
557: }
558: PetscCall(PetscFree(jackeys));
559: }
560: }
561: PetscCall(ISDestroy(&allcellIS));
562: PetscCall(DMDestroy(&plex));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*@
567: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
569: Input Parameters:
570: + dm - The `DM`
571: . X - Local input vector
572: - user - The user context
574: Output Parameters:
575: + Jac - Jacobian matrix
576: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
578: Level: developer
580: Note:
581: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
582: like a GPU, or vectorize on a multicore machine.
584: .seealso: [](ch_snes), `DMPLEX`, `Mat`
585: @*/
586: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
587: {
588: DM plex;
589: IS allcellIS;
590: PetscBool hasJac, hasPrec;
591: PetscInt Nds, s;
593: PetscFunctionBegin;
594: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
595: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
596: PetscCall(DMGetNumDS(dm, &Nds));
597: for (s = 0; s < Nds; ++s) {
598: PetscDS ds;
599: IS cellIS;
600: PetscFormKey key;
602: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
603: key.value = 0;
604: key.field = 0;
605: key.part = 0;
606: if (!key.label) {
607: PetscCall(PetscObjectReference((PetscObject)allcellIS));
608: cellIS = allcellIS;
609: } else {
610: IS pointIS;
612: key.value = 1;
613: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
614: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
615: PetscCall(ISDestroy(&pointIS));
616: }
617: if (!s) {
618: PetscCall(PetscDSHasJacobian(ds, &hasJac));
619: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
620: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
621: PetscCall(MatZeroEntries(JacP));
622: }
623: PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
624: PetscCall(ISDestroy(&cellIS));
625: }
626: PetscCall(ISDestroy(&allcellIS));
627: PetscCall(DMDestroy(&plex));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: struct _DMSNESJacobianMFCtx {
632: DM dm;
633: Vec X;
634: void *ctx;
635: };
637: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
638: {
639: struct _DMSNESJacobianMFCtx *ctx;
641: PetscFunctionBegin;
642: PetscCall(MatShellGetContext(A, &ctx));
643: PetscCall(MatShellSetContext(A, NULL));
644: PetscCall(DMDestroy(&ctx->dm));
645: PetscCall(VecDestroy(&ctx->X));
646: PetscCall(PetscFree(ctx));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
651: {
652: struct _DMSNESJacobianMFCtx *ctx;
654: PetscFunctionBegin;
655: PetscCall(MatShellGetContext(A, &ctx));
656: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: /*@
661: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
663: Collective
665: Input Parameters:
666: + dm - The `DM`
667: . X - The evaluation point for the Jacobian
668: - user - A user context, or `NULL`
670: Output Parameter:
671: . J - The `Mat`
673: Level: advanced
675: Notes:
676: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
678: This only works for `DMPLEX`
680: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
681: @*/
682: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
683: {
684: struct _DMSNESJacobianMFCtx *ctx;
685: PetscInt n, N;
687: PetscFunctionBegin;
688: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
689: PetscCall(MatSetType(*J, MATSHELL));
690: PetscCall(VecGetLocalSize(X, &n));
691: PetscCall(VecGetSize(X, &N));
692: PetscCall(MatSetSizes(*J, n, n, N, N));
693: PetscCall(PetscObjectReference((PetscObject)dm));
694: PetscCall(PetscObjectReference((PetscObject)X));
695: PetscCall(PetscMalloc1(1, &ctx));
696: ctx->dm = dm;
697: ctx->X = X;
698: ctx->ctx = user;
699: PetscCall(MatShellSetContext(*J, ctx));
700: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
701: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
706: {
707: SNES snes;
708: Mat pJ;
709: DM ovldm, origdm;
710: DMSNES sdm;
711: PetscErrorCode (*bfun)(DM, Vec, void *);
712: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
713: void *bctx, *jctx;
715: PetscFunctionBegin;
716: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
717: PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
718: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
719: PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
720: PetscCall(MatGetDM(pJ, &ovldm));
721: PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
722: PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
723: PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
724: PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
725: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
726: if (!snes) {
727: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
728: PetscCall(SNESSetDM(snes, ovldm));
729: PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
730: PetscCall(PetscObjectDereference((PetscObject)snes));
731: }
732: PetscCall(DMGetDMSNES(ovldm, &sdm));
733: PetscCall(VecLockReadPush(X));
734: {
735: void *ctx;
736: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
737: PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
738: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
739: }
740: PetscCall(VecLockReadPop(X));
741: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
742: {
743: Mat locpJ;
745: PetscCall(MatISGetLocalMat(pJ, &locpJ));
746: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
747: }
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: /*@
752: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
754: Input Parameters:
755: + dm - The `DM` object
756: . use_obj - Use the objective function callback
757: - ctx - The user context that will be passed to pointwise evaluation routines
759: Level: developer
761: .seealso: [](ch_snes),`DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
762: @*/
763: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, void *ctx)
764: {
765: PetscBool useCeed;
767: PetscFunctionBegin;
768: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
769: PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx));
770: if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx));
771: if (useCeed) {
772: #ifdef PETSC_HAVE_LIBCEED
773: PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx));
774: #else
775: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
776: #endif
777: } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx));
778: PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx));
779: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: /*@
784: DMSNESCheckDiscretization - Check the discretization error of the exact solution
786: Input Parameters:
787: + snes - the `SNES` object
788: . dm - the `DM`
789: . t - the time
790: . u - a `DM` vector
791: - tol - A tolerance for the check, or -1 to print the results instead
793: Output Parameter:
794: . error - An array which holds the discretization error in each field, or `NULL`
796: Level: developer
798: Note:
799: The user must call `PetscDSSetExactSolution()` beforehand
801: Developer Note:
802: How is this related to `PetscConvEst`?
804: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
805: @*/
806: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
807: {
808: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
809: void **ectxs;
810: PetscReal *err;
811: MPI_Comm comm;
812: PetscInt Nf, f;
814: PetscFunctionBegin;
818: if (error) PetscAssertPointer(error, 6);
820: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
821: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
823: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
824: PetscCall(DMGetNumFields(dm, &Nf));
825: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
826: {
827: PetscInt Nds, s;
829: PetscCall(DMGetNumDS(dm, &Nds));
830: for (s = 0; s < Nds; ++s) {
831: PetscDS ds;
832: DMLabel label;
833: IS fieldIS;
834: const PetscInt *fields;
835: PetscInt dsNf, f;
837: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
838: PetscCall(PetscDSGetNumFields(ds, &dsNf));
839: PetscCall(ISGetIndices(fieldIS, &fields));
840: for (f = 0; f < dsNf; ++f) {
841: const PetscInt field = fields[f];
842: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
843: }
844: PetscCall(ISRestoreIndices(fieldIS, &fields));
845: }
846: }
847: if (Nf > 1) {
848: PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
849: if (tol >= 0.0) {
850: 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);
851: } else if (error) {
852: for (f = 0; f < Nf; ++f) error[f] = err[f];
853: } else {
854: PetscCall(PetscPrintf(comm, "L_2 Error: ["));
855: for (f = 0; f < Nf; ++f) {
856: if (f) PetscCall(PetscPrintf(comm, ", "));
857: PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
858: }
859: PetscCall(PetscPrintf(comm, "]\n"));
860: }
861: } else {
862: PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
863: if (tol >= 0.0) {
864: PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
865: } else if (error) {
866: error[0] = err[0];
867: } else {
868: PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
869: }
870: }
871: PetscCall(PetscFree3(exacts, ectxs, err));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: /*@
876: DMSNESCheckResidual - Check the residual of the exact solution
878: Input Parameters:
879: + snes - the `SNES` object
880: . dm - the `DM`
881: . u - a `DM` vector
882: - tol - A tolerance for the check, or -1 to print the results instead
884: Output Parameter:
885: . residual - The residual norm of the exact solution, or `NULL`
887: Level: developer
889: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
890: @*/
891: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
892: {
893: MPI_Comm comm;
894: Vec r;
895: PetscReal res;
897: PetscFunctionBegin;
901: if (residual) PetscAssertPointer(residual, 5);
902: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
903: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
904: PetscCall(VecDuplicate(u, &r));
905: PetscCall(SNESComputeFunction(snes, u, r));
906: PetscCall(VecNorm(r, NORM_2, &res));
907: if (tol >= 0.0) {
908: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
909: } else if (residual) {
910: *residual = res;
911: } else {
912: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
913: PetscCall(VecFilter(r, 1.0e-10));
914: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
915: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
916: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
917: }
918: PetscCall(VecDestroy(&r));
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: /*@
923: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
925: Input Parameters:
926: + snes - the `SNES` object
927: . dm - the `DM`
928: . u - a `DM` vector
929: - tol - A tolerance for the check, or -1 to print the results instead
931: Output Parameters:
932: + isLinear - Flag indicaing that the function looks linear, or `NULL`
933: - convRate - The rate of convergence of the linear model, or `NULL`
935: Level: developer
937: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
938: @*/
939: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
940: {
941: MPI_Comm comm;
942: PetscDS ds;
943: Mat J, M;
944: MatNullSpace nullspace;
945: PetscReal slope, intercept;
946: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
948: PetscFunctionBegin;
952: if (isLinear) PetscAssertPointer(isLinear, 5);
953: if (convRate) PetscAssertPointer(convRate, 6);
954: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
955: if (!dm) PetscCall(SNESGetDM(snes, &dm));
956: if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
957: else PetscCall(SNESGetSolution(snes, &u));
958: /* Create and view matrices */
959: PetscCall(DMCreateMatrix(dm, &J));
960: PetscCall(DMGetDS(dm, &ds));
961: PetscCall(PetscDSHasJacobian(ds, &hasJac));
962: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
963: if (hasJac && hasPrec) {
964: PetscCall(DMCreateMatrix(dm, &M));
965: PetscCall(SNESComputeJacobian(snes, u, J, M));
966: PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
967: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
968: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
969: PetscCall(MatDestroy(&M));
970: } else {
971: PetscCall(SNESComputeJacobian(snes, u, J, J));
972: }
973: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
974: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
975: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
976: /* Check nullspace */
977: PetscCall(MatGetNullSpace(J, &nullspace));
978: if (nullspace) {
979: PetscBool isNull;
980: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
981: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
982: }
983: /* Taylor test */
984: {
985: PetscRandom rand;
986: Vec du, uhat, r, rhat, df;
987: PetscReal h;
988: PetscReal *es, *hs, *errors;
989: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
990: PetscInt Nv, v;
992: /* Choose a perturbation direction */
993: PetscCall(PetscRandomCreate(comm, &rand));
994: PetscCall(VecDuplicate(u, &du));
995: PetscCall(VecSetRandom(du, rand));
996: PetscCall(PetscRandomDestroy(&rand));
997: PetscCall(VecDuplicate(u, &df));
998: PetscCall(MatMult(J, du, df));
999: /* Evaluate residual at u, F(u), save in vector r */
1000: PetscCall(VecDuplicate(u, &r));
1001: PetscCall(SNESComputeFunction(snes, u, r));
1002: /* Look at the convergence of our Taylor approximation as we approach u */
1003: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1004: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1005: PetscCall(VecDuplicate(u, &uhat));
1006: PetscCall(VecDuplicate(u, &rhat));
1007: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1008: PetscCall(VecWAXPY(uhat, h, du, u));
1009: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1010: PetscCall(SNESComputeFunction(snes, uhat, rhat));
1011: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1012: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
1014: es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1015: hs[Nv] = PetscLog10Real(h);
1016: }
1017: PetscCall(VecDestroy(&uhat));
1018: PetscCall(VecDestroy(&rhat));
1019: PetscCall(VecDestroy(&df));
1020: PetscCall(VecDestroy(&r));
1021: PetscCall(VecDestroy(&du));
1022: for (v = 0; v < Nv; ++v) {
1023: if ((tol >= 0) && (errors[v] > tol)) break;
1024: else if (errors[v] > PETSC_SMALL) break;
1025: }
1026: if (v == Nv) isLin = PETSC_TRUE;
1027: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
1028: PetscCall(PetscFree3(es, hs, errors));
1029: /* Slope should be about 2 */
1030: if (tol >= 0) {
1031: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
1032: } else if (isLinear || convRate) {
1033: if (isLinear) *isLinear = isLin;
1034: if (convRate) *convRate = slope;
1035: } else {
1036: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
1037: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
1038: }
1039: }
1040: PetscCall(MatDestroy(&J));
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1045: {
1046: PetscFunctionBegin;
1047: PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1048: PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1049: PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1050: PetscFunctionReturn(PETSC_SUCCESS);
1051: }
1053: /*@
1054: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1056: Input Parameters:
1057: + snes - the `SNES` object
1058: - u - representative `SNES` vector
1060: Level: developer
1062: Note:
1063: The user must call `PetscDSSetExactSolution()` before this call
1065: .seealso: [](ch_snes), `SNES`, `DM`
1066: @*/
1067: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1068: {
1069: DM dm;
1070: Vec sol;
1071: PetscBool check;
1073: PetscFunctionBegin;
1074: PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1075: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1076: PetscCall(SNESGetDM(snes, &dm));
1077: PetscCall(VecDuplicate(u, &sol));
1078: PetscCall(SNESSetSolution(snes, sol));
1079: PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1080: PetscCall(VecDestroy(&sol));
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }