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: /************************** Interpolation *******************************/

143: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
144: {
145:   PetscBool isPlex;

147:   PetscFunctionBegin;
148:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
149:   if (isPlex) {
150:     *plex = dm;
151:     PetscCall(PetscObjectReference((PetscObject)dm));
152:   } else {
153:     PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
154:     if (!*plex) {
155:       PetscCall(DMConvert(dm, DMPLEX, plex));
156:       PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
157:       if (copy) {
158:         PetscCall(DMCopyDMSNES(dm, *plex));
159:         PetscCall(DMCopyAuxiliaryVec(dm, *plex));
160:       }
161:     } else {
162:       PetscCall(PetscObjectReference((PetscObject)*plex));
163:     }
164:   }
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: /*@C
169:   DMInterpolationCreate - Creates a `DMInterpolationInfo` context

171:   Collective

173:   Input Parameter:
174: . comm - the communicator

176:   Output Parameter:
177: . ctx - the context

179:   Level: beginner

181: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
182: @*/
183: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
184: {
185:   PetscFunctionBegin;
186:   PetscAssertPointer(ctx, 2);
187:   PetscCall(PetscNew(ctx));

189:   (*ctx)->comm   = comm;
190:   (*ctx)->dim    = -1;
191:   (*ctx)->nInput = 0;
192:   (*ctx)->points = NULL;
193:   (*ctx)->cells  = NULL;
194:   (*ctx)->n      = -1;
195:   (*ctx)->coords = NULL;
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: /*@C
200:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

202:   Not Collective

204:   Input Parameters:
205: + ctx - the context
206: - dim - the spatial dimension

208:   Level: intermediate

210: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
211: @*/
212: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
213: {
214:   PetscFunctionBegin;
215:   PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
216:   ctx->dim = dim;
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: /*@C
221:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

223:   Not Collective

225:   Input Parameter:
226: . ctx - the context

228:   Output Parameter:
229: . dim - the spatial dimension

231:   Level: intermediate

233: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
234: @*/
235: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
236: {
237:   PetscFunctionBegin;
238:   PetscAssertPointer(dim, 2);
239:   *dim = ctx->dim;
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*@C
244:   DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context

246:   Not Collective

248:   Input Parameters:
249: + ctx - the context
250: - dof - the number of fields

252:   Level: intermediate

254: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
255: @*/
256: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
257: {
258:   PetscFunctionBegin;
259:   PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
260:   ctx->dof = dof;
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: /*@C
265:   DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context

267:   Not Collective

269:   Input Parameter:
270: . ctx - the context

272:   Output Parameter:
273: . dof - the number of fields

275:   Level: intermediate

277: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
278: @*/
279: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
280: {
281:   PetscFunctionBegin;
282:   PetscAssertPointer(dof, 2);
283:   *dof = ctx->dof;
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: /*@C
288:   DMInterpolationAddPoints - Add points at which we will interpolate the fields

290:   Not Collective

292:   Input Parameters:
293: + ctx    - the context
294: . n      - the number of points
295: - points - the coordinates for each point, an array of size n * dim

297:   Level: intermediate

299:   Note:
300:   The coordinate information is copied.

302: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
303: @*/
304: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
305: {
306:   PetscFunctionBegin;
307:   PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
308:   PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
309:   ctx->nInput = n;

311:   PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
312:   PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
313:   PetscFunctionReturn(PETSC_SUCCESS);
314: }

316: /*@C
317:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

319:   Collective

321:   Input Parameters:
322: + ctx                 - the context
323: . dm                  - the `DM` for the function space used for interpolation
324: . redundantPoints     - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
325: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error

327:   Level: intermediate

329: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
330: @*/
331: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
332: {
333:   MPI_Comm           comm = ctx->comm;
334:   PetscScalar       *a;
335:   PetscInt           p, q, i;
336:   PetscMPIInt        rank, size;
337:   Vec                pointVec;
338:   PetscSF            cellSF;
339:   PetscLayout        layout;
340:   PetscReal         *globalPoints;
341:   PetscScalar       *globalPointsScalar;
342:   const PetscInt    *ranges;
343:   PetscMPIInt       *counts, *displs;
344:   const PetscSFNode *foundCells;
345:   const PetscInt    *foundPoints;
346:   PetscMPIInt       *foundProcs, *globalProcs;
347:   PetscInt           n, N, numFound;

349:   PetscFunctionBegin;
351:   PetscCallMPI(MPI_Comm_size(comm, &size));
352:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
353:   PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
354:   /* Locate points */
355:   n = ctx->nInput;
356:   if (!redundantPoints) {
357:     PetscCall(PetscLayoutCreate(comm, &layout));
358:     PetscCall(PetscLayoutSetBlockSize(layout, 1));
359:     PetscCall(PetscLayoutSetLocalSize(layout, n));
360:     PetscCall(PetscLayoutSetUp(layout));
361:     PetscCall(PetscLayoutGetSize(layout, &N));
362:     /* Communicate all points to all processes */
363:     PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs));
364:     PetscCall(PetscLayoutGetRanges(layout, &ranges));
365:     for (p = 0; p < size; ++p) {
366:       counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
367:       displs[p] = ranges[p] * ctx->dim;
368:     }
369:     PetscCallMPI(MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
370:   } else {
371:     N            = n;
372:     globalPoints = ctx->points;
373:     counts = displs = NULL;
374:     layout          = NULL;
375:   }
376: #if 0
377:   PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
378:   /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
379: #else
380:   #if defined(PETSC_USE_COMPLEX)
381:   PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
382:   for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
383:   #else
384:   globalPointsScalar = globalPoints;
385:   #endif
386:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
387:   PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
388:   for (p = 0; p < N; ++p) foundProcs[p] = size;
389:   cellSF = NULL;
390:   PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
391:   PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
392: #endif
393:   for (p = 0; p < numFound; ++p) {
394:     if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
395:   }
396:   /* Let the lowest rank process own each point */
397:   PetscCall(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
398:   ctx->n = 0;
399:   for (p = 0; p < N; ++p) {
400:     if (globalProcs[p] == size) {
401:       PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0.0),
402:                  (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
403:       if (rank == 0) ++ctx->n;
404:     } else if (globalProcs[p] == rank) ++ctx->n;
405:   }
406:   /* Create coordinates vector and array of owned cells */
407:   PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
408:   PetscCall(VecCreate(comm, &ctx->coords));
409:   PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
410:   PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
411:   PetscCall(VecSetType(ctx->coords, VECSTANDARD));
412:   PetscCall(VecGetArray(ctx->coords, &a));
413:   for (p = 0, q = 0, i = 0; p < N; ++p) {
414:     if (globalProcs[p] == rank) {
415:       PetscInt d;

417:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
418:       ctx->cells[q] = foundCells[q].index;
419:       ++q;
420:     }
421:     if (globalProcs[p] == size && rank == 0) {
422:       PetscInt d;

424:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
425:       ctx->cells[q] = -1;
426:       ++q;
427:     }
428:   }
429:   PetscCall(VecRestoreArray(ctx->coords, &a));
430: #if 0
431:   PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
432: #else
433:   PetscCall(PetscFree2(foundProcs, globalProcs));
434:   PetscCall(PetscSFDestroy(&cellSF));
435:   PetscCall(VecDestroy(&pointVec));
436: #endif
437:   if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
438:   if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
439:   PetscCall(PetscLayoutDestroy(&layout));
440:   PetscFunctionReturn(PETSC_SUCCESS);
441: }

443: /*@C
444:   DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point

446:   Collective

448:   Input Parameter:
449: . ctx - the context

451:   Output Parameter:
452: . coordinates - the coordinates of interpolation points

454:   Level: intermediate

456:   Note:
457:   The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
458:   This is a borrowed vector that the user should not destroy.

460: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
461: @*/
462: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
463: {
464:   PetscFunctionBegin;
465:   PetscAssertPointer(coordinates, 2);
466:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
467:   *coordinates = ctx->coords;
468:   PetscFunctionReturn(PETSC_SUCCESS);
469: }

471: /*@C
472:   DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values

474:   Collective

476:   Input Parameter:
477: . ctx - the context

479:   Output Parameter:
480: . v - a vector capable of holding the interpolated field values

482:   Level: intermediate

484:   Note:
485:   This vector should be returned using `DMInterpolationRestoreVector()`.

487: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
488: @*/
489: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
490: {
491:   PetscFunctionBegin;
492:   PetscAssertPointer(v, 2);
493:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
494:   PetscCall(VecCreate(ctx->comm, v));
495:   PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
496:   PetscCall(VecSetBlockSize(*v, ctx->dof));
497:   PetscCall(VecSetType(*v, VECSTANDARD));
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: /*@C
502:   DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values

504:   Collective

506:   Input Parameters:
507: + ctx - the context
508: - v   - a vector capable of holding the interpolated field values

510:   Level: intermediate

512: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
513: @*/
514: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
515: {
516:   PetscFunctionBegin;
517:   PetscAssertPointer(v, 2);
518:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
519:   PetscCall(VecDestroy(v));
520:   PetscFunctionReturn(PETSC_SUCCESS);
521: }

523: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
524: {
525:   PetscReal          v0, J, invJ, detJ;
526:   const PetscInt     dof = ctx->dof;
527:   const PetscScalar *coords;
528:   PetscScalar       *a;
529:   PetscInt           p;

531:   PetscFunctionBegin;
532:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
533:   PetscCall(VecGetArray(v, &a));
534:   for (p = 0; p < ctx->n; ++p) {
535:     PetscInt     c = ctx->cells[p];
536:     PetscScalar *x = NULL;
537:     PetscReal    xir[1];
538:     PetscInt     xSize, comp;

540:     PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
541:     PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
542:     xir[0] = invJ * PetscRealPart(coords[p] - v0);
543:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
544:     if (2 * dof == xSize) {
545:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
546:     } else if (dof == xSize) {
547:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
548:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof);
549:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
550:   }
551:   PetscCall(VecRestoreArray(v, &a));
552:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
557: {
558:   PetscReal         *v0, *J, *invJ, detJ;
559:   const PetscScalar *coords;
560:   PetscScalar       *a;
561:   PetscInt           p;

563:   PetscFunctionBegin;
564:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
565:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
566:   PetscCall(VecGetArray(v, &a));
567:   for (p = 0; p < ctx->n; ++p) {
568:     PetscInt     c = ctx->cells[p];
569:     PetscScalar *x = NULL;
570:     PetscReal    xi[4];
571:     PetscInt     d, f, comp;

573:     PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
574:     PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
575:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
576:     for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];

578:     for (d = 0; d < ctx->dim; ++d) {
579:       xi[d] = 0.0;
580:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
581:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
582:     }
583:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
584:   }
585:   PetscCall(VecRestoreArray(v, &a));
586:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
587:   PetscCall(PetscFree3(v0, J, invJ));
588:   PetscFunctionReturn(PETSC_SUCCESS);
589: }

591: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
592: {
593:   PetscReal         *v0, *J, *invJ, detJ;
594:   const PetscScalar *coords;
595:   PetscScalar       *a;
596:   PetscInt           p;

598:   PetscFunctionBegin;
599:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
600:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
601:   PetscCall(VecGetArray(v, &a));
602:   for (p = 0; p < ctx->n; ++p) {
603:     PetscInt       c        = ctx->cells[p];
604:     const PetscInt order[3] = {2, 1, 3};
605:     PetscScalar   *x        = NULL;
606:     PetscReal      xi[4];
607:     PetscInt       d, f, comp;

609:     PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
610:     PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
611:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
612:     for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];

614:     for (d = 0; d < ctx->dim; ++d) {
615:       xi[d] = 0.0;
616:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
617:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
618:     }
619:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
620:   }
621:   PetscCall(VecRestoreArray(v, &a));
622:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
623:   PetscCall(PetscFree3(v0, J, invJ));
624:   PetscFunctionReturn(PETSC_SUCCESS);
625: }

627: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
628: {
629:   const PetscScalar *vertices = (const PetscScalar *)ctx;
630:   const PetscScalar  x0       = vertices[0];
631:   const PetscScalar  y0       = vertices[1];
632:   const PetscScalar  x1       = vertices[2];
633:   const PetscScalar  y1       = vertices[3];
634:   const PetscScalar  x2       = vertices[4];
635:   const PetscScalar  y2       = vertices[5];
636:   const PetscScalar  x3       = vertices[6];
637:   const PetscScalar  y3       = vertices[7];
638:   const PetscScalar  f_1      = x1 - x0;
639:   const PetscScalar  g_1      = y1 - y0;
640:   const PetscScalar  f_3      = x3 - x0;
641:   const PetscScalar  g_3      = y3 - y0;
642:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
643:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
644:   const PetscScalar *ref;
645:   PetscScalar       *real;

647:   PetscFunctionBegin;
648:   PetscCall(VecGetArrayRead(Xref, &ref));
649:   PetscCall(VecGetArray(Xreal, &real));
650:   {
651:     const PetscScalar p0 = ref[0];
652:     const PetscScalar p1 = ref[1];

654:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
655:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
656:   }
657:   PetscCall(PetscLogFlops(28));
658:   PetscCall(VecRestoreArrayRead(Xref, &ref));
659:   PetscCall(VecRestoreArray(Xreal, &real));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: #include <petsc/private/dmimpl.h>
664: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
665: {
666:   const PetscScalar *vertices = (const PetscScalar *)ctx;
667:   const PetscScalar  x0       = vertices[0];
668:   const PetscScalar  y0       = vertices[1];
669:   const PetscScalar  x1       = vertices[2];
670:   const PetscScalar  y1       = vertices[3];
671:   const PetscScalar  x2       = vertices[4];
672:   const PetscScalar  y2       = vertices[5];
673:   const PetscScalar  x3       = vertices[6];
674:   const PetscScalar  y3       = vertices[7];
675:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
676:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
677:   const PetscScalar *ref;

679:   PetscFunctionBegin;
680:   PetscCall(VecGetArrayRead(Xref, &ref));
681:   {
682:     const PetscScalar x       = ref[0];
683:     const PetscScalar y       = ref[1];
684:     const PetscInt    rows[2] = {0, 1};
685:     PetscScalar       values[4];

687:     values[0] = (x1 - x0 + f_01 * y) * 0.5;
688:     values[1] = (x3 - x0 + f_01 * x) * 0.5;
689:     values[2] = (y1 - y0 + g_01 * y) * 0.5;
690:     values[3] = (y3 - y0 + g_01 * x) * 0.5;
691:     PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
692:   }
693:   PetscCall(PetscLogFlops(30));
694:   PetscCall(VecRestoreArrayRead(Xref, &ref));
695:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
696:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
697:   PetscFunctionReturn(PETSC_SUCCESS);
698: }

700: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
701: {
702:   DM                 dmCoord;
703:   PetscFE            fem = NULL;
704:   SNES               snes;
705:   KSP                ksp;
706:   PC                 pc;
707:   Vec                coordsLocal, r, ref, real;
708:   Mat                J;
709:   PetscTabulation    T = NULL;
710:   const PetscScalar *coords;
711:   PetscScalar       *a;
712:   PetscReal          xir[2] = {0., 0.};
713:   PetscInt           Nf, p;
714:   const PetscInt     dof = ctx->dof;

716:   PetscFunctionBegin;
717:   PetscCall(DMGetNumFields(dm, &Nf));
718:   if (Nf) {
719:     PetscObject  obj;
720:     PetscClassId id;

722:     PetscCall(DMGetField(dm, 0, NULL, &obj));
723:     PetscCall(PetscObjectGetClassId(obj, &id));
724:     if (id == PETSCFE_CLASSID) {
725:       fem = (PetscFE)obj;
726:       PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T));
727:     }
728:   }
729:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
730:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
731:   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
732:   PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_"));
733:   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
734:   PetscCall(VecSetSizes(r, 2, 2));
735:   PetscCall(VecSetType(r, dm->vectype));
736:   PetscCall(VecDuplicate(r, &ref));
737:   PetscCall(VecDuplicate(r, &real));
738:   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
739:   PetscCall(MatSetSizes(J, 2, 2, 2, 2));
740:   PetscCall(MatSetType(J, MATSEQDENSE));
741:   PetscCall(MatSetUp(J));
742:   PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL));
743:   PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL));
744:   PetscCall(SNESGetKSP(snes, &ksp));
745:   PetscCall(KSPGetPC(ksp, &pc));
746:   PetscCall(PCSetType(pc, PCLU));
747:   PetscCall(SNESSetFromOptions(snes));

749:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
750:   PetscCall(VecGetArray(v, &a));
751:   for (p = 0; p < ctx->n; ++p) {
752:     PetscScalar *x = NULL, *vertices = NULL;
753:     PetscScalar *xi;
754:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

756:     /* Can make this do all points at once */
757:     PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
758:     PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
759:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
760:     PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
761:     PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
762:     PetscCall(VecGetArray(real, &xi));
763:     xi[0] = coords[p * ctx->dim + 0];
764:     xi[1] = coords[p * ctx->dim + 1];
765:     PetscCall(VecRestoreArray(real, &xi));
766:     PetscCall(SNESSolve(snes, real, ref));
767:     PetscCall(VecGetArray(ref, &xi));
768:     xir[0] = PetscRealPart(xi[0]);
769:     xir[1] = PetscRealPart(xi[1]);
770:     if (4 * dof == xSize) {
771:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1];
772:     } else if (dof == xSize) {
773:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
774:     } else {
775:       PetscInt d;

777:       PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE");
778:       xir[0] = 2.0 * xir[0] - 1.0;
779:       xir[1] = 2.0 * xir[1] - 1.0;
780:       PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T));
781:       for (comp = 0; comp < dof; ++comp) {
782:         a[p * dof + comp] = 0.0;
783:         for (d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
784:       }
785:     }
786:     PetscCall(VecRestoreArray(ref, &xi));
787:     PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
788:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
789:   }
790:   PetscCall(PetscTabulationDestroy(&T));
791:   PetscCall(VecRestoreArray(v, &a));
792:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));

794:   PetscCall(SNESDestroy(&snes));
795:   PetscCall(VecDestroy(&r));
796:   PetscCall(VecDestroy(&ref));
797:   PetscCall(VecDestroy(&real));
798:   PetscCall(MatDestroy(&J));
799:   PetscFunctionReturn(PETSC_SUCCESS);
800: }

802: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
803: {
804:   const PetscScalar *vertices = (const PetscScalar *)ctx;
805:   const PetscScalar  x0       = vertices[0];
806:   const PetscScalar  y0       = vertices[1];
807:   const PetscScalar  z0       = vertices[2];
808:   const PetscScalar  x1       = vertices[9];
809:   const PetscScalar  y1       = vertices[10];
810:   const PetscScalar  z1       = vertices[11];
811:   const PetscScalar  x2       = vertices[6];
812:   const PetscScalar  y2       = vertices[7];
813:   const PetscScalar  z2       = vertices[8];
814:   const PetscScalar  x3       = vertices[3];
815:   const PetscScalar  y3       = vertices[4];
816:   const PetscScalar  z3       = vertices[5];
817:   const PetscScalar  x4       = vertices[12];
818:   const PetscScalar  y4       = vertices[13];
819:   const PetscScalar  z4       = vertices[14];
820:   const PetscScalar  x5       = vertices[15];
821:   const PetscScalar  y5       = vertices[16];
822:   const PetscScalar  z5       = vertices[17];
823:   const PetscScalar  x6       = vertices[18];
824:   const PetscScalar  y6       = vertices[19];
825:   const PetscScalar  z6       = vertices[20];
826:   const PetscScalar  x7       = vertices[21];
827:   const PetscScalar  y7       = vertices[22];
828:   const PetscScalar  z7       = vertices[23];
829:   const PetscScalar  f_1      = x1 - x0;
830:   const PetscScalar  g_1      = y1 - y0;
831:   const PetscScalar  h_1      = z1 - z0;
832:   const PetscScalar  f_3      = x3 - x0;
833:   const PetscScalar  g_3      = y3 - y0;
834:   const PetscScalar  h_3      = z3 - z0;
835:   const PetscScalar  f_4      = x4 - x0;
836:   const PetscScalar  g_4      = y4 - y0;
837:   const PetscScalar  h_4      = z4 - z0;
838:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
839:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
840:   const PetscScalar  h_01     = z2 - z1 - z3 + z0;
841:   const PetscScalar  f_12     = x7 - x3 - x4 + x0;
842:   const PetscScalar  g_12     = y7 - y3 - y4 + y0;
843:   const PetscScalar  h_12     = z7 - z3 - z4 + z0;
844:   const PetscScalar  f_02     = x5 - x1 - x4 + x0;
845:   const PetscScalar  g_02     = y5 - y1 - y4 + y0;
846:   const PetscScalar  h_02     = z5 - z1 - z4 + z0;
847:   const PetscScalar  f_012    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
848:   const PetscScalar  g_012    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
849:   const PetscScalar  h_012    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
850:   const PetscScalar *ref;
851:   PetscScalar       *real;

853:   PetscFunctionBegin;
854:   PetscCall(VecGetArrayRead(Xref, &ref));
855:   PetscCall(VecGetArray(Xreal, &real));
856:   {
857:     const PetscScalar p0 = ref[0];
858:     const PetscScalar p1 = ref[1];
859:     const PetscScalar p2 = ref[2];

861:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2;
862:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2;
863:     real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2;
864:   }
865:   PetscCall(PetscLogFlops(114));
866:   PetscCall(VecRestoreArrayRead(Xref, &ref));
867:   PetscCall(VecRestoreArray(Xreal, &real));
868:   PetscFunctionReturn(PETSC_SUCCESS);
869: }

871: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
872: {
873:   const PetscScalar *vertices = (const PetscScalar *)ctx;
874:   const PetscScalar  x0       = vertices[0];
875:   const PetscScalar  y0       = vertices[1];
876:   const PetscScalar  z0       = vertices[2];
877:   const PetscScalar  x1       = vertices[9];
878:   const PetscScalar  y1       = vertices[10];
879:   const PetscScalar  z1       = vertices[11];
880:   const PetscScalar  x2       = vertices[6];
881:   const PetscScalar  y2       = vertices[7];
882:   const PetscScalar  z2       = vertices[8];
883:   const PetscScalar  x3       = vertices[3];
884:   const PetscScalar  y3       = vertices[4];
885:   const PetscScalar  z3       = vertices[5];
886:   const PetscScalar  x4       = vertices[12];
887:   const PetscScalar  y4       = vertices[13];
888:   const PetscScalar  z4       = vertices[14];
889:   const PetscScalar  x5       = vertices[15];
890:   const PetscScalar  y5       = vertices[16];
891:   const PetscScalar  z5       = vertices[17];
892:   const PetscScalar  x6       = vertices[18];
893:   const PetscScalar  y6       = vertices[19];
894:   const PetscScalar  z6       = vertices[20];
895:   const PetscScalar  x7       = vertices[21];
896:   const PetscScalar  y7       = vertices[22];
897:   const PetscScalar  z7       = vertices[23];
898:   const PetscScalar  f_xy     = x2 - x1 - x3 + x0;
899:   const PetscScalar  g_xy     = y2 - y1 - y3 + y0;
900:   const PetscScalar  h_xy     = z2 - z1 - z3 + z0;
901:   const PetscScalar  f_yz     = x7 - x3 - x4 + x0;
902:   const PetscScalar  g_yz     = y7 - y3 - y4 + y0;
903:   const PetscScalar  h_yz     = z7 - z3 - z4 + z0;
904:   const PetscScalar  f_xz     = x5 - x1 - x4 + x0;
905:   const PetscScalar  g_xz     = y5 - y1 - y4 + y0;
906:   const PetscScalar  h_xz     = z5 - z1 - z4 + z0;
907:   const PetscScalar  f_xyz    = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
908:   const PetscScalar  g_xyz    = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
909:   const PetscScalar  h_xyz    = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
910:   const PetscScalar *ref;

912:   PetscFunctionBegin;
913:   PetscCall(VecGetArrayRead(Xref, &ref));
914:   {
915:     const PetscScalar x       = ref[0];
916:     const PetscScalar y       = ref[1];
917:     const PetscScalar z       = ref[2];
918:     const PetscInt    rows[3] = {0, 1, 2};
919:     PetscScalar       values[9];

921:     values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
922:     values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
923:     values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
924:     values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
925:     values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
926:     values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
927:     values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
928:     values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
929:     values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;

931:     PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
932:   }
933:   PetscCall(PetscLogFlops(152));
934:   PetscCall(VecRestoreArrayRead(Xref, &ref));
935:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
936:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
937:   PetscFunctionReturn(PETSC_SUCCESS);
938: }

940: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
941: {
942:   DM                 dmCoord;
943:   SNES               snes;
944:   KSP                ksp;
945:   PC                 pc;
946:   Vec                coordsLocal, r, ref, real;
947:   Mat                J;
948:   const PetscScalar *coords;
949:   PetscScalar       *a;
950:   PetscInt           p;

952:   PetscFunctionBegin;
953:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
954:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
955:   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
956:   PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
957:   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
958:   PetscCall(VecSetSizes(r, 3, 3));
959:   PetscCall(VecSetType(r, dm->vectype));
960:   PetscCall(VecDuplicate(r, &ref));
961:   PetscCall(VecDuplicate(r, &real));
962:   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
963:   PetscCall(MatSetSizes(J, 3, 3, 3, 3));
964:   PetscCall(MatSetType(J, MATSEQDENSE));
965:   PetscCall(MatSetUp(J));
966:   PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
967:   PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
968:   PetscCall(SNESGetKSP(snes, &ksp));
969:   PetscCall(KSPGetPC(ksp, &pc));
970:   PetscCall(PCSetType(pc, PCLU));
971:   PetscCall(SNESSetFromOptions(snes));

973:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
974:   PetscCall(VecGetArray(v, &a));
975:   for (p = 0; p < ctx->n; ++p) {
976:     PetscScalar *x = NULL, *vertices = NULL;
977:     PetscScalar *xi;
978:     PetscReal    xir[3];
979:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

981:     /* Can make this do all points at once */
982:     PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
983:     PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
984:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
985:     PetscCheck((8 * ctx->dof == xSize) || (ctx->dof == xSize), ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input closure size %" PetscInt_FMT " should be %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 8 * ctx->dof, ctx->dof);
986:     PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
987:     PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
988:     PetscCall(VecGetArray(real, &xi));
989:     xi[0] = coords[p * ctx->dim + 0];
990:     xi[1] = coords[p * ctx->dim + 1];
991:     xi[2] = coords[p * ctx->dim + 2];
992:     PetscCall(VecRestoreArray(real, &xi));
993:     PetscCall(SNESSolve(snes, real, ref));
994:     PetscCall(VecGetArray(ref, &xi));
995:     xir[0] = PetscRealPart(xi[0]);
996:     xir[1] = PetscRealPart(xi[1]);
997:     xir[2] = PetscRealPart(xi[2]);
998:     if (8 * ctx->dof == xSize) {
999:       for (comp = 0; comp < ctx->dof; ++comp) {
1000:         a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) +
1001:                                  x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2];
1002:       }
1003:     } else {
1004:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
1005:     }
1006:     PetscCall(VecRestoreArray(ref, &xi));
1007:     PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
1008:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
1009:   }
1010:   PetscCall(VecRestoreArray(v, &a));
1011:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));

1013:   PetscCall(SNESDestroy(&snes));
1014:   PetscCall(VecDestroy(&r));
1015:   PetscCall(VecDestroy(&ref));
1016:   PetscCall(VecDestroy(&real));
1017:   PetscCall(MatDestroy(&J));
1018:   PetscFunctionReturn(PETSC_SUCCESS);
1019: }

1021: /*@C
1022:   DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.

1024:   Input Parameters:
1025: + ctx - The `DMInterpolationInfo` context
1026: . dm  - The `DM`
1027: - x   - The local vector containing the field to be interpolated

1029:   Output Parameter:
1030: . v - The vector containing the interpolated values

1032:   Level: beginner

1034:   Note:
1035:   A suitable `v` can be obtained using `DMInterpolationGetVector()`.

1037: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1038: @*/
1039: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1040: {
1041:   PetscDS   ds;
1042:   PetscInt  n, p, Nf, field;
1043:   PetscBool useDS = PETSC_FALSE;

1045:   PetscFunctionBegin;
1049:   PetscCall(VecGetLocalSize(v, &n));
1050:   PetscCheck(n == ctx->n * ctx->dof, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %" PetscInt_FMT " should be %" PetscInt_FMT, n, ctx->n * ctx->dof);
1051:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
1052:   PetscCall(DMGetDS(dm, &ds));
1053:   if (ds) {
1054:     useDS = PETSC_TRUE;
1055:     PetscCall(PetscDSGetNumFields(ds, &Nf));
1056:     for (field = 0; field < Nf; ++field) {
1057:       PetscObject  obj;
1058:       PetscClassId id;

1060:       PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1061:       PetscCall(PetscObjectGetClassId(obj, &id));
1062:       if (id != PETSCFE_CLASSID) {
1063:         useDS = PETSC_FALSE;
1064:         break;
1065:       }
1066:     }
1067:   }
1068:   if (useDS) {
1069:     const PetscScalar *coords;
1070:     PetscScalar       *interpolant;
1071:     PetscInt           cdim, d;

1073:     PetscCall(DMGetCoordinateDim(dm, &cdim));
1074:     PetscCall(VecGetArrayRead(ctx->coords, &coords));
1075:     PetscCall(VecGetArrayWrite(v, &interpolant));
1076:     for (p = 0; p < ctx->n; ++p) {
1077:       PetscReal    pcoords[3], xi[3];
1078:       PetscScalar *xa   = NULL;
1079:       PetscInt     coff = 0, foff = 0, clSize;

1081:       if (ctx->cells[p] < 0) continue;
1082:       for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
1083:       PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
1084:       PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1085:       for (field = 0; field < Nf; ++field) {
1086:         PetscTabulation T;
1087:         PetscFE         fe;

1089:         PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1090:         PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
1091:         {
1092:           const PetscReal *basis = T->T[0];
1093:           const PetscInt   Nb    = T->Nb;
1094:           const PetscInt   Nc    = T->Nc;
1095:           PetscInt         f, fc;

1097:           for (fc = 0; fc < Nc; ++fc) {
1098:             interpolant[p * ctx->dof + coff + fc] = 0.0;
1099:             for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
1100:           }
1101:           coff += Nc;
1102:           foff += Nb;
1103:         }
1104:         PetscCall(PetscTabulationDestroy(&T));
1105:       }
1106:       PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1107:       PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
1108:       PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
1109:     }
1110:     PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
1111:     PetscCall(VecRestoreArrayWrite(v, &interpolant));
1112:   } else {
1113:     DMPolytopeType ct;

1115:     /* TODO Check each cell individually */
1116:     PetscCall(DMPlexGetCellType(dm, ctx->cells[0], &ct));
1117:     switch (ct) {
1118:     case DM_POLYTOPE_SEGMENT:
1119:       PetscCall(DMInterpolate_Segment_Private(ctx, dm, x, v));
1120:       break;
1121:     case DM_POLYTOPE_TRIANGLE:
1122:       PetscCall(DMInterpolate_Triangle_Private(ctx, dm, x, v));
1123:       break;
1124:     case DM_POLYTOPE_QUADRILATERAL:
1125:       PetscCall(DMInterpolate_Quad_Private(ctx, dm, x, v));
1126:       break;
1127:     case DM_POLYTOPE_TETRAHEDRON:
1128:       PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, x, v));
1129:       break;
1130:     case DM_POLYTOPE_HEXAHEDRON:
1131:       PetscCall(DMInterpolate_Hex_Private(ctx, dm, x, v));
1132:       break;
1133:     default:
1134:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
1135:     }
1136:   }
1137:   PetscFunctionReturn(PETSC_SUCCESS);
1138: }

1140: /*@C
1141:   DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context

1143:   Collective

1145:   Input Parameter:
1146: . ctx - the context

1148:   Level: beginner

1150: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1151: @*/
1152: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1153: {
1154:   PetscFunctionBegin;
1155:   PetscAssertPointer(ctx, 1);
1156:   PetscCall(VecDestroy(&(*ctx)->coords));
1157:   PetscCall(PetscFree((*ctx)->points));
1158:   PetscCall(PetscFree((*ctx)->cells));
1159:   PetscCall(PetscFree(*ctx));
1160:   *ctx = NULL;
1161:   PetscFunctionReturn(PETSC_SUCCESS);
1162: }

1164: /*@C
1165:   SNESMonitorFields - Monitors the residual for each field separately

1167:   Collective

1169:   Input Parameters:
1170: + snes   - the `SNES` context, must have an attached `DM`
1171: . its    - iteration number
1172: . fgnorm - 2-norm of residual
1173: - vf     - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`

1175:   Level: intermediate

1177:   Note:
1178:   This routine prints the residual norm at each iteration.

1180: .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
1181: @*/
1182: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1183: {
1184:   PetscViewer        viewer = vf->viewer;
1185:   Vec                res;
1186:   DM                 dm;
1187:   PetscSection       s;
1188:   const PetscScalar *r;
1189:   PetscReal         *lnorms, *norms;
1190:   PetscInt           numFields, f, pStart, pEnd, p;

1192:   PetscFunctionBegin;
1194:   PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
1195:   PetscCall(SNESGetDM(snes, &dm));
1196:   PetscCall(DMGetLocalSection(dm, &s));
1197:   PetscCall(PetscSectionGetNumFields(s, &numFields));
1198:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1199:   PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
1200:   PetscCall(VecGetArrayRead(res, &r));
1201:   for (p = pStart; p < pEnd; ++p) {
1202:     for (f = 0; f < numFields; ++f) {
1203:       PetscInt fdof, foff, d;

1205:       PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
1206:       PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
1207:       for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
1208:     }
1209:   }
1210:   PetscCall(VecRestoreArrayRead(res, &r));
1211:   PetscCall(MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
1212:   PetscCall(PetscViewerPushFormat(viewer, vf->format));
1213:   PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
1214:   PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
1215:   for (f = 0; f < numFields; ++f) {
1216:     if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
1217:     PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
1218:   }
1219:   PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
1220:   PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
1221:   PetscCall(PetscViewerPopFormat(viewer));
1222:   PetscCall(PetscFree2(lnorms, norms));
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: /********************* Residual Computation **************************/

1228: /*@
1229:   DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user

1231:   Input Parameters:
1232: + dm   - The mesh
1233: . X    - Local solution
1234: - user - The user context

1236:   Output Parameter:
1237: . F - Local output vector

1239:   Level: developer

1241:   Note:
1242:   The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.

1244: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1245: @*/
1246: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1247: {
1248:   DM       plex;
1249:   IS       allcellIS;
1250:   PetscInt Nds, s;

1252:   PetscFunctionBegin;
1253:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1254:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1255:   PetscCall(DMGetNumDS(dm, &Nds));
1256:   for (s = 0; s < Nds; ++s) {
1257:     PetscDS      ds;
1258:     IS           cellIS;
1259:     PetscFormKey key;

1261:     PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
1262:     key.value = 0;
1263:     key.field = 0;
1264:     key.part  = 0;
1265:     if (!key.label) {
1266:       PetscCall(PetscObjectReference((PetscObject)allcellIS));
1267:       cellIS = allcellIS;
1268:     } else {
1269:       IS pointIS;

1271:       key.value = 1;
1272:       PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
1273:       PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1274:       PetscCall(ISDestroy(&pointIS));
1275:     }
1276:     PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
1277:     PetscCall(ISDestroy(&cellIS));
1278:   }
1279:   PetscCall(ISDestroy(&allcellIS));
1280:   PetscCall(DMDestroy(&plex));
1281:   PetscFunctionReturn(PETSC_SUCCESS);
1282: }

1284: /*@
1285:   DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`

1287:   Input Parameters:
1288: + dm   - The mesh
1289: . X    - Local solution
1290: - user - The user context

1292:   Output Parameter:
1293: . F - Local output vector

1295:   Level: developer

1297:   Note:
1298:   The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.

1300: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1301: @*/
1302: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user)
1303: {
1304:   DM       plex;
1305:   IS       allcellIS;
1306:   PetscInt Nds, s;

1308:   PetscFunctionBegin;
1309:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1310:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1311:   PetscCall(DMGetNumDS(dm, &Nds));
1312:   for (s = 0; s < Nds; ++s) {
1313:     PetscDS ds;
1314:     DMLabel label;
1315:     IS      cellIS;

1317:     PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1318:     {
1319:       PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1320:       PetscWeakForm     wf;
1321:       PetscInt          Nm = 2, m, Nk = 0, k, kp, off = 0;
1322:       PetscFormKey     *reskeys;

1324:       /* Get unique residual keys */
1325:       for (m = 0; m < Nm; ++m) {
1326:         PetscInt Nkm;
1327:         PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
1328:         Nk += Nkm;
1329:       }
1330:       PetscCall(PetscMalloc1(Nk, &reskeys));
1331:       for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
1332:       PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1333:       PetscCall(PetscFormKeySort(Nk, reskeys));
1334:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1335:         if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1336:           ++k;
1337:           if (kp != k) reskeys[k] = reskeys[kp];
1338:         }
1339:       }
1340:       Nk = k;

1342:       PetscCall(PetscDSGetWeakForm(ds, &wf));
1343:       for (k = 0; k < Nk; ++k) {
1344:         DMLabel  label = reskeys[k].label;
1345:         PetscInt val   = reskeys[k].value;

1347:         if (!label) {
1348:           PetscCall(PetscObjectReference((PetscObject)allcellIS));
1349:           cellIS = allcellIS;
1350:         } else {
1351:           IS pointIS;

1353:           PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1354:           PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1355:           PetscCall(ISDestroy(&pointIS));
1356:         }
1357:         PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
1358:         PetscCall(ISDestroy(&cellIS));
1359:       }
1360:       PetscCall(PetscFree(reskeys));
1361:     }
1362:   }
1363:   PetscCall(ISDestroy(&allcellIS));
1364:   PetscCall(DMDestroy(&plex));
1365:   PetscFunctionReturn(PETSC_SUCCESS);
1366: }

1368: #ifdef PETSC_HAVE_LIBCEED
1369: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user)
1370: {
1371:   Ceed       ceed;
1372:   DMCeed     sd = dm->dmceed;
1373:   CeedVector clocX, clocF;

1375:   PetscFunctionBegin;
1376:   PetscCall(DMGetCeed(dm, &ceed));
1377:   PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
1378:   PetscCall(DMCeedComputeGeometry(dm, sd));

1380:   PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
1381:   PetscCall(VecGetCeedVector(locF, ceed, &clocF));
1382:   PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
1383:   PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
1384:   PetscCall(VecRestoreCeedVector(locF, &clocF));

1386:   {
1387:     DM_Plex *mesh = (DM_Plex *)dm->data;

1389:     if (mesh->printFEM) {
1390:       PetscSection section;
1391:       Vec          locFbc;
1392:       PetscInt     pStart, pEnd, p, maxDof;
1393:       PetscScalar *zeroes;

1395:       PetscCall(DMGetLocalSection(dm, &section));
1396:       PetscCall(VecDuplicate(locF, &locFbc));
1397:       PetscCall(VecCopy(locF, locFbc));
1398:       PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
1399:       PetscCall(PetscSectionGetMaxDof(section, &maxDof));
1400:       PetscCall(PetscCalloc1(maxDof, &zeroes));
1401:       for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
1402:       PetscCall(PetscFree(zeroes));
1403:       PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
1404:       PetscCall(VecDestroy(&locFbc));
1405:     }
1406:   }
1407:   PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1409: #endif

1411: /*@
1412:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`

1414:   Input Parameters:
1415: + dm   - The mesh
1416: - user - The user context

1418:   Output Parameter:
1419: . X - Local solution

1421:   Level: developer

1423: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1424: @*/
1425: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1426: {
1427:   DM plex;

1429:   PetscFunctionBegin;
1430:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1431:   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
1432:   PetscCall(DMDestroy(&plex));
1433:   PetscFunctionReturn(PETSC_SUCCESS);
1434: }

1436: /*@
1437:   DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`

1439:   Input Parameters:
1440: + dm   - The `DM`
1441: . X    - Local solution vector
1442: . Y    - Local input vector
1443: - user - The user context

1445:   Output Parameter:
1446: . F - local output vector

1448:   Level: developer

1450:   Note:
1451:   Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.

1453:   This only works with `DMPLEX`

1455:   Developer Note:
1456:   This should be called `DMPlexSNESComputeJacobianAction()`

1458: .seealso: [](ch_snes), `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
1459: @*/
1460: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1461: {
1462:   DM       plex;
1463:   IS       allcellIS;
1464:   PetscInt Nds, s;

1466:   PetscFunctionBegin;
1467:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1468:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1469:   PetscCall(DMGetNumDS(dm, &Nds));
1470:   for (s = 0; s < Nds; ++s) {
1471:     PetscDS ds;
1472:     DMLabel label;
1473:     IS      cellIS;

1475:     PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1476:     {
1477:       PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1478:       PetscWeakForm     wf;
1479:       PetscInt          Nm = 4, m, Nk = 0, k, kp, off = 0;
1480:       PetscFormKey     *jackeys;

1482:       /* Get unique Jacobian keys */
1483:       for (m = 0; m < Nm; ++m) {
1484:         PetscInt Nkm;
1485:         PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
1486:         Nk += Nkm;
1487:       }
1488:       PetscCall(PetscMalloc1(Nk, &jackeys));
1489:       for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
1490:       PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1491:       PetscCall(PetscFormKeySort(Nk, jackeys));
1492:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1493:         if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1494:           ++k;
1495:           if (kp != k) jackeys[k] = jackeys[kp];
1496:         }
1497:       }
1498:       Nk = k;

1500:       PetscCall(PetscDSGetWeakForm(ds, &wf));
1501:       for (k = 0; k < Nk; ++k) {
1502:         DMLabel  label = jackeys[k].label;
1503:         PetscInt val   = jackeys[k].value;

1505:         if (!label) {
1506:           PetscCall(PetscObjectReference((PetscObject)allcellIS));
1507:           cellIS = allcellIS;
1508:         } else {
1509:           IS pointIS;

1511:           PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1512:           PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1513:           PetscCall(ISDestroy(&pointIS));
1514:         }
1515:         PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
1516:         PetscCall(ISDestroy(&cellIS));
1517:       }
1518:       PetscCall(PetscFree(jackeys));
1519:     }
1520:   }
1521:   PetscCall(ISDestroy(&allcellIS));
1522:   PetscCall(DMDestroy(&plex));
1523:   PetscFunctionReturn(PETSC_SUCCESS);
1524: }

1526: /*@
1527:   DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.

1529:   Input Parameters:
1530: + dm   - The `DM`
1531: . X    - Local input vector
1532: - user - The user context

1534:   Output Parameters:
1535: + Jac  - Jacobian matrix
1536: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`

1538:   Level: developer

1540:   Note:
1541:   We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1542:   like a GPU, or vectorize on a multicore machine.

1544: .seealso: [](ch_snes), `DMPLEX`, `Mat`
1545: @*/
1546: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
1547: {
1548:   DM        plex;
1549:   IS        allcellIS;
1550:   PetscBool hasJac, hasPrec;
1551:   PetscInt  Nds, s;

1553:   PetscFunctionBegin;
1554:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1555:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1556:   PetscCall(DMGetNumDS(dm, &Nds));
1557:   for (s = 0; s < Nds; ++s) {
1558:     PetscDS      ds;
1559:     IS           cellIS;
1560:     PetscFormKey key;

1562:     PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
1563:     key.value = 0;
1564:     key.field = 0;
1565:     key.part  = 0;
1566:     if (!key.label) {
1567:       PetscCall(PetscObjectReference((PetscObject)allcellIS));
1568:       cellIS = allcellIS;
1569:     } else {
1570:       IS pointIS;

1572:       key.value = 1;
1573:       PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
1574:       PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1575:       PetscCall(ISDestroy(&pointIS));
1576:     }
1577:     if (!s) {
1578:       PetscCall(PetscDSHasJacobian(ds, &hasJac));
1579:       PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1580:       if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
1581:       PetscCall(MatZeroEntries(JacP));
1582:     }
1583:     PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
1584:     PetscCall(ISDestroy(&cellIS));
1585:   }
1586:   PetscCall(ISDestroy(&allcellIS));
1587:   PetscCall(DMDestroy(&plex));
1588:   PetscFunctionReturn(PETSC_SUCCESS);
1589: }

1591: struct _DMSNESJacobianMFCtx {
1592:   DM    dm;
1593:   Vec   X;
1594:   void *ctx;
1595: };

1597: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1598: {
1599:   struct _DMSNESJacobianMFCtx *ctx;

1601:   PetscFunctionBegin;
1602:   PetscCall(MatShellGetContext(A, &ctx));
1603:   PetscCall(MatShellSetContext(A, NULL));
1604:   PetscCall(DMDestroy(&ctx->dm));
1605:   PetscCall(VecDestroy(&ctx->X));
1606:   PetscCall(PetscFree(ctx));
1607:   PetscFunctionReturn(PETSC_SUCCESS);
1608: }

1610: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1611: {
1612:   struct _DMSNESJacobianMFCtx *ctx;

1614:   PetscFunctionBegin;
1615:   PetscCall(MatShellGetContext(A, &ctx));
1616:   PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
1617:   PetscFunctionReturn(PETSC_SUCCESS);
1618: }

1620: /*@
1621:   DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free

1623:   Collective

1625:   Input Parameters:
1626: + dm   - The `DM`
1627: . X    - The evaluation point for the Jacobian
1628: - user - A user context, or `NULL`

1630:   Output Parameter:
1631: . J - The `Mat`

1633:   Level: advanced

1635:   Notes:
1636:   Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.

1638:   This only works for `DMPLEX`

1640: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
1641: @*/
1642: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1643: {
1644:   struct _DMSNESJacobianMFCtx *ctx;
1645:   PetscInt                     n, N;

1647:   PetscFunctionBegin;
1648:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
1649:   PetscCall(MatSetType(*J, MATSHELL));
1650:   PetscCall(VecGetLocalSize(X, &n));
1651:   PetscCall(VecGetSize(X, &N));
1652:   PetscCall(MatSetSizes(*J, n, n, N, N));
1653:   PetscCall(PetscObjectReference((PetscObject)dm));
1654:   PetscCall(PetscObjectReference((PetscObject)X));
1655:   PetscCall(PetscMalloc1(1, &ctx));
1656:   ctx->dm  = dm;
1657:   ctx->X   = X;
1658:   ctx->ctx = user;
1659:   PetscCall(MatShellSetContext(*J, ctx));
1660:   PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
1661:   PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
1662:   PetscFunctionReturn(PETSC_SUCCESS);
1663: }

1665: /*
1666:      MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.

1668:    Input Parameters:
1669: +     X - `SNES` linearization point
1670: .     ovl - index set of overlapping subdomains

1672:    Output Parameter:
1673: .     J - unassembled (Neumann) local matrix

1675:    Level: intermediate

1677: .seealso: [](ch_snes), `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()`
1678: */
1679: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1680: {
1681:   SNES   snes;
1682:   Mat    pJ;
1683:   DM     ovldm, origdm;
1684:   DMSNES sdm;
1685:   PetscErrorCode (*bfun)(DM, Vec, void *);
1686:   PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
1687:   void *bctx, *jctx;

1689:   PetscFunctionBegin;
1690:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
1691:   PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
1692:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
1693:   PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
1694:   PetscCall(MatGetDM(pJ, &ovldm));
1695:   PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
1696:   PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
1697:   PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
1698:   PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
1699:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
1700:   if (!snes) {
1701:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
1702:     PetscCall(SNESSetDM(snes, ovldm));
1703:     PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
1704:     PetscCall(PetscObjectDereference((PetscObject)snes));
1705:   }
1706:   PetscCall(DMGetDMSNES(ovldm, &sdm));
1707:   PetscCall(VecLockReadPush(X));
1708:   {
1709:     void *ctx;
1710:     PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
1711:     PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
1712:     PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
1713:   }
1714:   PetscCall(VecLockReadPop(X));
1715:   /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1716:   {
1717:     Mat locpJ;

1719:     PetscCall(MatISGetLocalMat(pJ, &locpJ));
1720:     PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
1721:   }
1722:   PetscFunctionReturn(PETSC_SUCCESS);
1723: }

1725: /*@
1726:   DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian.

1728:   Input Parameters:
1729: + dm          - The `DM` object
1730: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`)
1731: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`)
1732: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`)

1734:   Level: developer

1736: .seealso: [](ch_snes), `DMPLEX`, `SNES`
1737: @*/
1738: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1739: {
1740:   PetscBool useCeed;

1742:   PetscFunctionBegin;
1743:   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1744:   PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx));
1745:   if (useCeed) {
1746: #ifdef PETSC_HAVE_LIBCEED
1747:     PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, residualctx));
1748: #else
1749:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
1750: #endif
1751:   } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx));
1752:   PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx));
1753:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
1754:   PetscFunctionReturn(PETSC_SUCCESS);
1755: }

1757: /*@C
1758:   DMSNESCheckDiscretization - Check the discretization error of the exact solution

1760:   Input Parameters:
1761: + snes - the `SNES` object
1762: . dm   - the `DM`
1763: . t    - the time
1764: . u    - a `DM` vector
1765: - tol  - A tolerance for the check, or -1 to print the results instead

1767:   Output Parameter:
1768: . error - An array which holds the discretization error in each field, or `NULL`

1770:   Level: developer

1772:   Note:
1773:   The user must call `PetscDSSetExactSolution()` beforehand

1775:   Developer Note:
1776:   How is this related to `PetscConvEst`?

1778: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
1779: @*/
1780: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1781: {
1782:   PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1783:   void     **ectxs;
1784:   PetscReal *err;
1785:   MPI_Comm   comm;
1786:   PetscInt   Nf, f;

1788:   PetscFunctionBegin;
1792:   if (error) PetscAssertPointer(error, 6);

1794:   PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1795:   PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));

1797:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1798:   PetscCall(DMGetNumFields(dm, &Nf));
1799:   PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
1800:   {
1801:     PetscInt Nds, s;

1803:     PetscCall(DMGetNumDS(dm, &Nds));
1804:     for (s = 0; s < Nds; ++s) {
1805:       PetscDS         ds;
1806:       DMLabel         label;
1807:       IS              fieldIS;
1808:       const PetscInt *fields;
1809:       PetscInt        dsNf, f;

1811:       PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
1812:       PetscCall(PetscDSGetNumFields(ds, &dsNf));
1813:       PetscCall(ISGetIndices(fieldIS, &fields));
1814:       for (f = 0; f < dsNf; ++f) {
1815:         const PetscInt field = fields[f];
1816:         PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
1817:       }
1818:       PetscCall(ISRestoreIndices(fieldIS, &fields));
1819:     }
1820:   }
1821:   if (Nf > 1) {
1822:     PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
1823:     if (tol >= 0.0) {
1824:       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);
1825:     } else if (error) {
1826:       for (f = 0; f < Nf; ++f) error[f] = err[f];
1827:     } else {
1828:       PetscCall(PetscPrintf(comm, "L_2 Error: ["));
1829:       for (f = 0; f < Nf; ++f) {
1830:         if (f) PetscCall(PetscPrintf(comm, ", "));
1831:         PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
1832:       }
1833:       PetscCall(PetscPrintf(comm, "]\n"));
1834:     }
1835:   } else {
1836:     PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
1837:     if (tol >= 0.0) {
1838:       PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
1839:     } else if (error) {
1840:       error[0] = err[0];
1841:     } else {
1842:       PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
1843:     }
1844:   }
1845:   PetscCall(PetscFree3(exacts, ectxs, err));
1846:   PetscFunctionReturn(PETSC_SUCCESS);
1847: }

1849: /*@C
1850:   DMSNESCheckResidual - Check the residual of the exact solution

1852:   Input Parameters:
1853: + snes - the `SNES` object
1854: . dm   - the `DM`
1855: . u    - a `DM` vector
1856: - tol  - A tolerance for the check, or -1 to print the results instead

1858:   Output Parameter:
1859: . residual - The residual norm of the exact solution, or `NULL`

1861:   Level: developer

1863: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
1864: @*/
1865: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1866: {
1867:   MPI_Comm  comm;
1868:   Vec       r;
1869:   PetscReal res;

1871:   PetscFunctionBegin;
1875:   if (residual) PetscAssertPointer(residual, 5);
1876:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1877:   PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1878:   PetscCall(VecDuplicate(u, &r));
1879:   PetscCall(SNESComputeFunction(snes, u, r));
1880:   PetscCall(VecNorm(r, NORM_2, &res));
1881:   if (tol >= 0.0) {
1882:     PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
1883:   } else if (residual) {
1884:     *residual = res;
1885:   } else {
1886:     PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
1887:     PetscCall(VecFilter(r, 1.0e-10));
1888:     PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
1889:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
1890:     PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
1891:   }
1892:   PetscCall(VecDestroy(&r));
1893:   PetscFunctionReturn(PETSC_SUCCESS);
1894: }

1896: /*@C
1897:   DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test

1899:   Input Parameters:
1900: + snes - the `SNES` object
1901: . dm   - the `DM`
1902: . u    - a `DM` vector
1903: - tol  - A tolerance for the check, or -1 to print the results instead

1905:   Output Parameters:
1906: + isLinear - Flag indicaing that the function looks linear, or `NULL`
1907: - convRate - The rate of convergence of the linear model, or `NULL`

1909:   Level: developer

1911: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1912: @*/
1913: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1914: {
1915:   MPI_Comm     comm;
1916:   PetscDS      ds;
1917:   Mat          J, M;
1918:   MatNullSpace nullspace;
1919:   PetscReal    slope, intercept;
1920:   PetscBool    hasJac, hasPrec, isLin = PETSC_FALSE;

1922:   PetscFunctionBegin;
1926:   if (isLinear) PetscAssertPointer(isLinear, 5);
1927:   if (convRate) PetscAssertPointer(convRate, 6);
1928:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1929:   PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1930:   /* Create and view matrices */
1931:   PetscCall(DMCreateMatrix(dm, &J));
1932:   PetscCall(DMGetDS(dm, &ds));
1933:   PetscCall(PetscDSHasJacobian(ds, &hasJac));
1934:   PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1935:   if (hasJac && hasPrec) {
1936:     PetscCall(DMCreateMatrix(dm, &M));
1937:     PetscCall(SNESComputeJacobian(snes, u, J, M));
1938:     PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
1939:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
1940:     PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
1941:     PetscCall(MatDestroy(&M));
1942:   } else {
1943:     PetscCall(SNESComputeJacobian(snes, u, J, J));
1944:   }
1945:   PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
1946:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
1947:   PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
1948:   /* Check nullspace */
1949:   PetscCall(MatGetNullSpace(J, &nullspace));
1950:   if (nullspace) {
1951:     PetscBool isNull;
1952:     PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
1953:     PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1954:   }
1955:   /* Taylor test */
1956:   {
1957:     PetscRandom rand;
1958:     Vec         du, uhat, r, rhat, df;
1959:     PetscReal   h;
1960:     PetscReal  *es, *hs, *errors;
1961:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1962:     PetscInt    Nv, v;

1964:     /* Choose a perturbation direction */
1965:     PetscCall(PetscRandomCreate(comm, &rand));
1966:     PetscCall(VecDuplicate(u, &du));
1967:     PetscCall(VecSetRandom(du, rand));
1968:     PetscCall(PetscRandomDestroy(&rand));
1969:     PetscCall(VecDuplicate(u, &df));
1970:     PetscCall(MatMult(J, du, df));
1971:     /* Evaluate residual at u, F(u), save in vector r */
1972:     PetscCall(VecDuplicate(u, &r));
1973:     PetscCall(SNESComputeFunction(snes, u, r));
1974:     /* Look at the convergence of our Taylor approximation as we approach u */
1975:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
1976:       ;
1977:     PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1978:     PetscCall(VecDuplicate(u, &uhat));
1979:     PetscCall(VecDuplicate(u, &rhat));
1980:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1981:       PetscCall(VecWAXPY(uhat, h, du, u));
1982:       /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1983:       PetscCall(SNESComputeFunction(snes, uhat, rhat));
1984:       PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1985:       PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));

1987:       es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1988:       hs[Nv] = PetscLog10Real(h);
1989:     }
1990:     PetscCall(VecDestroy(&uhat));
1991:     PetscCall(VecDestroy(&rhat));
1992:     PetscCall(VecDestroy(&df));
1993:     PetscCall(VecDestroy(&r));
1994:     PetscCall(VecDestroy(&du));
1995:     for (v = 0; v < Nv; ++v) {
1996:       if ((tol >= 0) && (errors[v] > tol)) break;
1997:       else if (errors[v] > PETSC_SMALL) break;
1998:     }
1999:     if (v == Nv) isLin = PETSC_TRUE;
2000:     PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
2001:     PetscCall(PetscFree3(es, hs, errors));
2002:     /* Slope should be about 2 */
2003:     if (tol >= 0) {
2004:       PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
2005:     } else if (isLinear || convRate) {
2006:       if (isLinear) *isLinear = isLin;
2007:       if (convRate) *convRate = slope;
2008:     } else {
2009:       if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
2010:       else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
2011:     }
2012:   }
2013:   PetscCall(MatDestroy(&J));
2014:   PetscFunctionReturn(PETSC_SUCCESS);
2015: }

2017: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
2018: {
2019:   PetscFunctionBegin;
2020:   PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
2021:   PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
2022:   PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
2023:   PetscFunctionReturn(PETSC_SUCCESS);
2024: }

2026: /*@C
2027:   DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information

2029:   Input Parameters:
2030: + snes - the `SNES` object
2031: - u    - representative `SNES` vector

2033:   Level: developer

2035:   Note:
2036:   The user must call `PetscDSSetExactSolution()` before this call

2038: .seealso: [](ch_snes), `SNES`, `DM`
2039: @*/
2040: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
2041: {
2042:   DM        dm;
2043:   Vec       sol;
2044:   PetscBool check;

2046:   PetscFunctionBegin;
2047:   PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
2048:   if (!check) PetscFunctionReturn(PETSC_SUCCESS);
2049:   PetscCall(SNESGetDM(snes, &dm));
2050:   PetscCall(VecDuplicate(u, &sol));
2051:   PetscCall(SNESSetSolution(snes, sol));
2052:   PetscCall(DMSNESCheck_Internal(snes, dm, sol));
2053:   PetscCall(VecDestroy(&sol));
2054:   PetscFunctionReturn(PETSC_SUCCESS);
2055: }