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:   PetscCall(MPIU_Allreduce(lnorms, norms, 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, &section));
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: /*@C
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: /*@C
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: /*@C
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: /*@C
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: }