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: 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[])
  8: {
  9:   p[0] = u[uOff[1]];
 10: }

 12: /*
 13:   SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
 14:   This is normally used only to evaluate convergence rates for the pressure accurately.

 16:   Collective

 18:   Input Parameters:
 19: + snes      - The SNES
 20: . pfield    - The field number for pressure
 21: . nullspace - The pressure nullspace
 22: . u         - The solution vector
 23: - ctx       - An optional user context

 25:   Output Parameter:
 26: . u         - The solution with a continuum pressure integral of zero

 28:   Level: developer

 30:   Notes:
 31:   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.

 33: .seealso: `SNESConvergedCorrectPressure()`
 34: */
 35: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
 36: {
 37:   DM          dm;
 38:   PetscDS     ds;
 39:   const Vec  *nullvecs;
 40:   PetscScalar pintd, *intc, *intn;
 41:   MPI_Comm    comm;
 42:   PetscInt    Nf, Nv;

 44:   PetscFunctionBegin;
 45:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
 46:   PetscCall(SNESGetDM(snes, &dm));
 47:   PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
 48:   PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
 49:   PetscCall(DMGetDS(dm, &ds));
 50:   PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private));
 51:   PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs));
 52:   PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv);
 53:   PetscCall(VecDot(nullvecs[0], u, &pintd));
 54:   PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd));
 55:   PetscCall(PetscDSGetNumFields(ds, &Nf));
 56:   PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn));
 57:   PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx));
 58:   PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
 59:   PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]));
 60: #if defined(PETSC_USE_DEBUG)
 61:   PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
 62:   PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield]));
 63: #endif
 64:   PetscCall(PetscFree2(intc, intn));
 65:   PetscFunctionReturn(PETSC_SUCCESS);
 66: }

 68: /*@C
 69:    SNESConvergedCorrectPressure - Convergence test that adds a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
 70:    This is normally used only to evaluate convergence rates for the pressure accurately. The convergence test itself just mimics `SNESConvergedDefault()`.

 72:    Logically Collective

 74:    Input Parameters:
 75: +  snes - the `SNES` context
 76: .  it - the iteration (0 indicates before any Newton steps)
 77: .  xnorm - 2-norm of current iterate
 78: .  snorm - 2-norm of current step
 79: .  fnorm - 2-norm of function at current iterate
 80: -  ctx   - Optional user context

 82:    Output Parameter:
 83: .  reason  - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`

 85:    Options Database Key:
 86: .  -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`

 88:    Level: advanced

 90:    Notes:
 91:    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`
 92:    must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
 93:    The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
 94:    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.

 96: .seealso: `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`, `DMSetNullSpaceConstructor()`
 97: @*/
 98: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
 99: {
100:   PetscBool monitorIntegral = PETSC_FALSE;

102:   PetscFunctionBegin;
103:   PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx));
104:   if (monitorIntegral) {
105:     Mat          J;
106:     Vec          u;
107:     MatNullSpace nullspace;
108:     const Vec   *nullvecs;
109:     PetscScalar  pintd;

111:     PetscCall(SNESGetSolution(snes, &u));
112:     PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
113:     PetscCall(MatGetNullSpace(J, &nullspace));
114:     PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs));
115:     PetscCall(VecDot(nullvecs[0], u, &pintd));
116:     PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd)));
117:   }
118:   if (*reason > 0) {
119:     Mat          J;
120:     Vec          u;
121:     MatNullSpace nullspace;
122:     PetscInt     pfield = 1;

124:     PetscCall(SNESGetSolution(snes, &u));
125:     PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
126:     PetscCall(MatGetNullSpace(J, &nullspace));
127:     PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx));
128:   }
129:   PetscFunctionReturn(PETSC_SUCCESS);
130: }

132: /************************** Interpolation *******************************/

134: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
135: {
136:   PetscBool isPlex;

138:   PetscFunctionBegin;
139:   PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
140:   if (isPlex) {
141:     *plex = dm;
142:     PetscCall(PetscObjectReference((PetscObject)dm));
143:   } else {
144:     PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
145:     if (!*plex) {
146:       PetscCall(DMConvert(dm, DMPLEX, plex));
147:       PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
148:       if (copy) {
149:         PetscCall(DMCopyDMSNES(dm, *plex));
150:         PetscCall(DMCopyAuxiliaryVec(dm, *plex));
151:       }
152:     } else {
153:       PetscCall(PetscObjectReference((PetscObject)*plex));
154:     }
155:   }
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: /*@C
160:   DMInterpolationCreate - Creates a `DMInterpolationInfo` context

162:   Collective

164:   Input Parameter:
165: . comm - the communicator

167:   Output Parameter:
168: . ctx - the context

170:   Level: beginner

172: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
173: @*/
174: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
175: {
176:   PetscFunctionBegin;
178:   PetscCall(PetscNew(ctx));

180:   (*ctx)->comm   = comm;
181:   (*ctx)->dim    = -1;
182:   (*ctx)->nInput = 0;
183:   (*ctx)->points = NULL;
184:   (*ctx)->cells  = NULL;
185:   (*ctx)->n      = -1;
186:   (*ctx)->coords = NULL;
187:   PetscFunctionReturn(PETSC_SUCCESS);
188: }

190: /*@C
191:   DMInterpolationSetDim - Sets the spatial dimension for the interpolation context

193:   Not Collective

195:   Input Parameters:
196: + ctx - the context
197: - dim - the spatial dimension

199:   Level: intermediate

201: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
202: @*/
203: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
204: {
205:   PetscFunctionBegin;
206:   PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
207:   ctx->dim = dim;
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: /*@C
212:   DMInterpolationGetDim - Gets the spatial dimension for the interpolation context

214:   Not Collective

216:   Input Parameter:
217: . ctx - the context

219:   Output Parameter:
220: . dim - the spatial dimension

222:   Level: intermediate

224: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
225: @*/
226: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
227: {
228:   PetscFunctionBegin;
230:   *dim = ctx->dim;
231:   PetscFunctionReturn(PETSC_SUCCESS);
232: }

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

237:   Not Collective

239:   Input Parameters:
240: + ctx - the context
241: - dof - the number of fields

243:   Level: intermediate

245: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
246: @*/
247: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
248: {
249:   PetscFunctionBegin;
250:   PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
251:   ctx->dof = dof;
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

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

258:   Not Collective

260:   Input Parameter:
261: . ctx - the context

263:   Output Parameter:
264: . dof - the number of fields

266:   Level: intermediate

268: .seealso: DMInterpolationInfo, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
269: @*/
270: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
271: {
272:   PetscFunctionBegin;
274:   *dof = ctx->dof;
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

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

281:   Not Collective

283:   Input Parameters:
284: + ctx    - the context
285: . n      - the number of points
286: - points - the coordinates for each point, an array of size n * dim

288:   Level: intermediate

290:   Note:
291:   The coordinate information is copied.

293: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
294: @*/
295: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
296: {
297:   PetscFunctionBegin;
298:   PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
299:   PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
300:   ctx->nInput = n;

302:   PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
303:   PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: /*@C
308:   DMInterpolationSetUp - Compute spatial indices for point location during interpolation

310:   Collective

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

318:   Level: intermediate

320: .seealso: DMInterpolationInfo, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
321: @*/
322: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
323: {
324:   MPI_Comm           comm = ctx->comm;
325:   PetscScalar       *a;
326:   PetscInt           p, q, i;
327:   PetscMPIInt        rank, size;
328:   Vec                pointVec;
329:   PetscSF            cellSF;
330:   PetscLayout        layout;
331:   PetscReal         *globalPoints;
332:   PetscScalar       *globalPointsScalar;
333:   const PetscInt    *ranges;
334:   PetscMPIInt       *counts, *displs;
335:   const PetscSFNode *foundCells;
336:   const PetscInt    *foundPoints;
337:   PetscMPIInt       *foundProcs, *globalProcs;
338:   PetscInt           n, N, numFound;

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

408:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
409:       ctx->cells[q] = foundCells[q].index;
410:       ++q;
411:     }
412:     if (globalProcs[p] == size && rank == 0) {
413:       PetscInt d;

415:       for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
416:       ctx->cells[q] = -1;
417:       ++q;
418:     }
419:   }
420:   PetscCall(VecRestoreArray(ctx->coords, &a));
421: #if 0
422:   PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
423: #else
424:   PetscCall(PetscFree2(foundProcs, globalProcs));
425:   PetscCall(PetscSFDestroy(&cellSF));
426:   PetscCall(VecDestroy(&pointVec));
427: #endif
428:   if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
429:   if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
430:   PetscCall(PetscLayoutDestroy(&layout));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

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

437:   Collective

439:   Input Parameter:
440: . ctx - the context

442:   Output Parameter:
443: . coordinates  - the coordinates of interpolation points

445:   Level: intermediate

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

451: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
452: @*/
453: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
454: {
455:   PetscFunctionBegin;
457:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
458:   *coordinates = ctx->coords;
459:   PetscFunctionReturn(PETSC_SUCCESS);
460: }

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

465:   Collective

467:   Input Parameter:
468: . ctx - the context

470:   Output Parameter:
471: . v  - a vector capable of holding the interpolated field values

473:   Level: intermediate

475:   Note:
476:   This vector should be returned using `DMInterpolationRestoreVector()`.

478: .seealso: `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
479: @*/
480: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
481: {
482:   PetscFunctionBegin;
484:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
485:   PetscCall(VecCreate(ctx->comm, v));
486:   PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
487:   PetscCall(VecSetBlockSize(*v, ctx->dof));
488:   PetscCall(VecSetType(*v, VECSTANDARD));
489:   PetscFunctionReturn(PETSC_SUCCESS);
490: }

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

495:   Collective

497:   Input Parameters:
498: + ctx - the context
499: - v  - a vector capable of holding the interpolated field values

501:   Level: intermediate

503: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
504: @*/
505: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
506: {
507:   PetscFunctionBegin;
509:   PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
510:   PetscCall(VecDestroy(v));
511:   PetscFunctionReturn(PETSC_SUCCESS);
512: }

514: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
515: {
516:   PetscReal          v0, J, invJ, detJ;
517:   const PetscInt     dof = ctx->dof;
518:   const PetscScalar *coords;
519:   PetscScalar       *a;
520:   PetscInt           p;

522:   PetscFunctionBegin;
523:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
524:   PetscCall(VecGetArray(v, &a));
525:   for (p = 0; p < ctx->n; ++p) {
526:     PetscInt     c = ctx->cells[p];
527:     PetscScalar *x = NULL;
528:     PetscReal    xir[1];
529:     PetscInt     xSize, comp;

531:     PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
532:     PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
533:     xir[0] = invJ * PetscRealPart(coords[p] - v0);
534:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
535:     if (2 * dof == xSize) {
536:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
537:     } else if (dof == xSize) {
538:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
539:     } 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);
540:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
541:   }
542:   PetscCall(VecRestoreArray(v, &a));
543:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
548: {
549:   PetscReal         *v0, *J, *invJ, detJ;
550:   const PetscScalar *coords;
551:   PetscScalar       *a;
552:   PetscInt           p;

554:   PetscFunctionBegin;
555:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
556:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
557:   PetscCall(VecGetArray(v, &a));
558:   for (p = 0; p < ctx->n; ++p) {
559:     PetscInt     c = ctx->cells[p];
560:     PetscScalar *x = NULL;
561:     PetscReal    xi[4];
562:     PetscInt     d, f, comp;

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

569:     for (d = 0; d < ctx->dim; ++d) {
570:       xi[d] = 0.0;
571:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
572:       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];
573:     }
574:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
575:   }
576:   PetscCall(VecRestoreArray(v, &a));
577:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
578:   PetscCall(PetscFree3(v0, J, invJ));
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
583: {
584:   PetscReal         *v0, *J, *invJ, detJ;
585:   const PetscScalar *coords;
586:   PetscScalar       *a;
587:   PetscInt           p;

589:   PetscFunctionBegin;
590:   PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
591:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
592:   PetscCall(VecGetArray(v, &a));
593:   for (p = 0; p < ctx->n; ++p) {
594:     PetscInt       c        = ctx->cells[p];
595:     const PetscInt order[3] = {2, 1, 3};
596:     PetscScalar   *x        = NULL;
597:     PetscReal      xi[4];
598:     PetscInt       d, f, comp;

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

605:     for (d = 0; d < ctx->dim; ++d) {
606:       xi[d] = 0.0;
607:       for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
608:       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];
609:     }
610:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
611:   }
612:   PetscCall(VecRestoreArray(v, &a));
613:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
614:   PetscCall(PetscFree3(v0, J, invJ));
615:   PetscFunctionReturn(PETSC_SUCCESS);
616: }

618: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
619: {
620:   const PetscScalar *vertices = (const PetscScalar *)ctx;
621:   const PetscScalar  x0       = vertices[0];
622:   const PetscScalar  y0       = vertices[1];
623:   const PetscScalar  x1       = vertices[2];
624:   const PetscScalar  y1       = vertices[3];
625:   const PetscScalar  x2       = vertices[4];
626:   const PetscScalar  y2       = vertices[5];
627:   const PetscScalar  x3       = vertices[6];
628:   const PetscScalar  y3       = vertices[7];
629:   const PetscScalar  f_1      = x1 - x0;
630:   const PetscScalar  g_1      = y1 - y0;
631:   const PetscScalar  f_3      = x3 - x0;
632:   const PetscScalar  g_3      = y3 - y0;
633:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
634:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
635:   const PetscScalar *ref;
636:   PetscScalar       *real;

638:   PetscFunctionBegin;
639:   PetscCall(VecGetArrayRead(Xref, &ref));
640:   PetscCall(VecGetArray(Xreal, &real));
641:   {
642:     const PetscScalar p0 = ref[0];
643:     const PetscScalar p1 = ref[1];

645:     real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
646:     real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
647:   }
648:   PetscCall(PetscLogFlops(28));
649:   PetscCall(VecRestoreArrayRead(Xref, &ref));
650:   PetscCall(VecRestoreArray(Xreal, &real));
651:   PetscFunctionReturn(PETSC_SUCCESS);
652: }

654: #include <petsc/private/dmimpl.h>
655: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
656: {
657:   const PetscScalar *vertices = (const PetscScalar *)ctx;
658:   const PetscScalar  x0       = vertices[0];
659:   const PetscScalar  y0       = vertices[1];
660:   const PetscScalar  x1       = vertices[2];
661:   const PetscScalar  y1       = vertices[3];
662:   const PetscScalar  x2       = vertices[4];
663:   const PetscScalar  y2       = vertices[5];
664:   const PetscScalar  x3       = vertices[6];
665:   const PetscScalar  y3       = vertices[7];
666:   const PetscScalar  f_01     = x2 - x1 - x3 + x0;
667:   const PetscScalar  g_01     = y2 - y1 - y3 + y0;
668:   const PetscScalar *ref;

670:   PetscFunctionBegin;
671:   PetscCall(VecGetArrayRead(Xref, &ref));
672:   {
673:     const PetscScalar x       = ref[0];
674:     const PetscScalar y       = ref[1];
675:     const PetscInt    rows[2] = {0, 1};
676:     PetscScalar       values[4];

678:     values[0] = (x1 - x0 + f_01 * y) * 0.5;
679:     values[1] = (x3 - x0 + f_01 * x) * 0.5;
680:     values[2] = (y1 - y0 + g_01 * y) * 0.5;
681:     values[3] = (y3 - y0 + g_01 * x) * 0.5;
682:     PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
683:   }
684:   PetscCall(PetscLogFlops(30));
685:   PetscCall(VecRestoreArrayRead(Xref, &ref));
686:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
687:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
688:   PetscFunctionReturn(PETSC_SUCCESS);
689: }

691: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
692: {
693:   DM                 dmCoord;
694:   PetscFE            fem = NULL;
695:   SNES               snes;
696:   KSP                ksp;
697:   PC                 pc;
698:   Vec                coordsLocal, r, ref, real;
699:   Mat                J;
700:   PetscTabulation    T = NULL;
701:   const PetscScalar *coords;
702:   PetscScalar       *a;
703:   PetscReal          xir[2] = {0., 0.};
704:   PetscInt           Nf, p;
705:   const PetscInt     dof = ctx->dof;

707:   PetscFunctionBegin;
708:   PetscCall(DMGetNumFields(dm, &Nf));
709:   if (Nf) {
710:     PetscObject  obj;
711:     PetscClassId id;

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

740:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
741:   PetscCall(VecGetArray(v, &a));
742:   for (p = 0; p < ctx->n; ++p) {
743:     PetscScalar *x = NULL, *vertices = NULL;
744:     PetscScalar *xi;
745:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

747:     /* Can make this do all points at once */
748:     PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
749:     PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
750:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
751:     PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
752:     PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
753:     PetscCall(VecGetArray(real, &xi));
754:     xi[0] = coords[p * ctx->dim + 0];
755:     xi[1] = coords[p * ctx->dim + 1];
756:     PetscCall(VecRestoreArray(real, &xi));
757:     PetscCall(SNESSolve(snes, real, ref));
758:     PetscCall(VecGetArray(ref, &xi));
759:     xir[0] = PetscRealPart(xi[0]);
760:     xir[1] = PetscRealPart(xi[1]);
761:     if (4 * dof == xSize) {
762:       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];
763:     } else if (dof == xSize) {
764:       for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
765:     } else {
766:       PetscInt d;

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

785:   PetscCall(SNESDestroy(&snes));
786:   PetscCall(VecDestroy(&r));
787:   PetscCall(VecDestroy(&ref));
788:   PetscCall(VecDestroy(&real));
789:   PetscCall(MatDestroy(&J));
790:   PetscFunctionReturn(PETSC_SUCCESS);
791: }

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

844:   PetscFunctionBegin;
845:   PetscCall(VecGetArrayRead(Xref, &ref));
846:   PetscCall(VecGetArray(Xreal, &real));
847:   {
848:     const PetscScalar p0 = ref[0];
849:     const PetscScalar p1 = ref[1];
850:     const PetscScalar p2 = ref[2];

852:     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;
853:     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;
854:     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;
855:   }
856:   PetscCall(PetscLogFlops(114));
857:   PetscCall(VecRestoreArrayRead(Xref, &ref));
858:   PetscCall(VecRestoreArray(Xreal, &real));
859:   PetscFunctionReturn(PETSC_SUCCESS);
860: }

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

903:   PetscFunctionBegin;
904:   PetscCall(VecGetArrayRead(Xref, &ref));
905:   {
906:     const PetscScalar x       = ref[0];
907:     const PetscScalar y       = ref[1];
908:     const PetscScalar z       = ref[2];
909:     const PetscInt    rows[3] = {0, 1, 2};
910:     PetscScalar       values[9];

912:     values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
913:     values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
914:     values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
915:     values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
916:     values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
917:     values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
918:     values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
919:     values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
920:     values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;

922:     PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
923:   }
924:   PetscCall(PetscLogFlops(152));
925:   PetscCall(VecRestoreArrayRead(Xref, &ref));
926:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
927:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
928:   PetscFunctionReturn(PETSC_SUCCESS);
929: }

931: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
932: {
933:   DM                 dmCoord;
934:   SNES               snes;
935:   KSP                ksp;
936:   PC                 pc;
937:   Vec                coordsLocal, r, ref, real;
938:   Mat                J;
939:   const PetscScalar *coords;
940:   PetscScalar       *a;
941:   PetscInt           p;

943:   PetscFunctionBegin;
944:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
945:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
946:   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
947:   PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
948:   PetscCall(VecCreate(PETSC_COMM_SELF, &r));
949:   PetscCall(VecSetSizes(r, 3, 3));
950:   PetscCall(VecSetType(r, dm->vectype));
951:   PetscCall(VecDuplicate(r, &ref));
952:   PetscCall(VecDuplicate(r, &real));
953:   PetscCall(MatCreate(PETSC_COMM_SELF, &J));
954:   PetscCall(MatSetSizes(J, 3, 3, 3, 3));
955:   PetscCall(MatSetType(J, MATSEQDENSE));
956:   PetscCall(MatSetUp(J));
957:   PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
958:   PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
959:   PetscCall(SNESGetKSP(snes, &ksp));
960:   PetscCall(KSPGetPC(ksp, &pc));
961:   PetscCall(PCSetType(pc, PCLU));
962:   PetscCall(SNESSetFromOptions(snes));

964:   PetscCall(VecGetArrayRead(ctx->coords, &coords));
965:   PetscCall(VecGetArray(v, &a));
966:   for (p = 0; p < ctx->n; ++p) {
967:     PetscScalar *x = NULL, *vertices = NULL;
968:     PetscScalar *xi;
969:     PetscReal    xir[3];
970:     PetscInt     c = ctx->cells[p], comp, coordSize, xSize;

972:     /* Can make this do all points at once */
973:     PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
974:     PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
975:     PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
976:     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);
977:     PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
978:     PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
979:     PetscCall(VecGetArray(real, &xi));
980:     xi[0] = coords[p * ctx->dim + 0];
981:     xi[1] = coords[p * ctx->dim + 1];
982:     xi[2] = coords[p * ctx->dim + 2];
983:     PetscCall(VecRestoreArray(real, &xi));
984:     PetscCall(SNESSolve(snes, real, ref));
985:     PetscCall(VecGetArray(ref, &xi));
986:     xir[0] = PetscRealPart(xi[0]);
987:     xir[1] = PetscRealPart(xi[1]);
988:     xir[2] = PetscRealPart(xi[2]);
989:     if (8 * ctx->dof == xSize) {
990:       for (comp = 0; comp < ctx->dof; ++comp) {
991:         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]) +
992:                                  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];
993:       }
994:     } else {
995:       for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
996:     }
997:     PetscCall(VecRestoreArray(ref, &xi));
998:     PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
999:     PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
1000:   }
1001:   PetscCall(VecRestoreArray(v, &a));
1002:   PetscCall(VecRestoreArrayRead(ctx->coords, &coords));

1004:   PetscCall(SNESDestroy(&snes));
1005:   PetscCall(VecDestroy(&r));
1006:   PetscCall(VecDestroy(&ref));
1007:   PetscCall(VecDestroy(&real));
1008:   PetscCall(MatDestroy(&J));
1009:   PetscFunctionReturn(PETSC_SUCCESS);
1010: }

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

1015:   Input Parameters:
1016: + ctx - The `DMInterpolationInfo` context
1017: . dm  - The `DM`
1018: - x   - The local vector containing the field to be interpolated

1020:   Output Parameter:
1021: . v   - The vector containing the interpolated values

1023:   Level: beginner

1025:   Note:
1026:   A suitable `v` can be obtained using `DMInterpolationGetVector()`.

1028: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1029: @*/
1030: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1031: {
1032:   PetscDS   ds;
1033:   PetscInt  n, p, Nf, field;
1034:   PetscBool useDS = PETSC_FALSE;

1036:   PetscFunctionBegin;
1040:   PetscCall(VecGetLocalSize(v, &n));
1041:   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);
1042:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
1043:   PetscCall(DMGetDS(dm, &ds));
1044:   if (ds) {
1045:     useDS = PETSC_TRUE;
1046:     PetscCall(PetscDSGetNumFields(ds, &Nf));
1047:     for (field = 0; field < Nf; ++field) {
1048:       PetscObject  obj;
1049:       PetscClassId id;

1051:       PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1052:       PetscCall(PetscObjectGetClassId(obj, &id));
1053:       if (id != PETSCFE_CLASSID) {
1054:         useDS = PETSC_FALSE;
1055:         break;
1056:       }
1057:     }
1058:   }
1059:   if (useDS) {
1060:     const PetscScalar *coords;
1061:     PetscScalar       *interpolant;
1062:     PetscInt           cdim, d;

1064:     PetscCall(DMGetCoordinateDim(dm, &cdim));
1065:     PetscCall(VecGetArrayRead(ctx->coords, &coords));
1066:     PetscCall(VecGetArrayWrite(v, &interpolant));
1067:     for (p = 0; p < ctx->n; ++p) {
1068:       PetscReal    pcoords[3], xi[3];
1069:       PetscScalar *xa   = NULL;
1070:       PetscInt     coff = 0, foff = 0, clSize;

1072:       if (ctx->cells[p] < 0) continue;
1073:       for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
1074:       PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
1075:       PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1076:       for (field = 0; field < Nf; ++field) {
1077:         PetscTabulation T;
1078:         PetscFE         fe;

1080:         PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1081:         PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
1082:         {
1083:           const PetscReal *basis = T->T[0];
1084:           const PetscInt   Nb    = T->Nb;
1085:           const PetscInt   Nc    = T->Nc;
1086:           PetscInt         f, fc;

1088:           for (fc = 0; fc < Nc; ++fc) {
1089:             interpolant[p * ctx->dof + coff + fc] = 0.0;
1090:             for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
1091:           }
1092:           coff += Nc;
1093:           foff += Nb;
1094:         }
1095:         PetscCall(PetscTabulationDestroy(&T));
1096:       }
1097:       PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1098:       PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
1099:       PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
1100:     }
1101:     PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
1102:     PetscCall(VecRestoreArrayWrite(v, &interpolant));
1103:   } else {
1104:     DMPolytopeType ct;

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

1131: /*@C
1132:   DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context

1134:   Collective

1136:   Input Parameter:
1137: . ctx - the context

1139:   Level: beginner

1141: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1142: @*/
1143: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1144: {
1145:   PetscFunctionBegin;
1147:   PetscCall(VecDestroy(&(*ctx)->coords));
1148:   PetscCall(PetscFree((*ctx)->points));
1149:   PetscCall(PetscFree((*ctx)->cells));
1150:   PetscCall(PetscFree(*ctx));
1151:   *ctx = NULL;
1152:   PetscFunctionReturn(PETSC_SUCCESS);
1153: }

1155: /*@C
1156:   SNESMonitorFields - Monitors the residual for each field separately

1158:   Collective

1160:   Input Parameters:
1161: + snes   - the `SNES` context
1162: . its    - iteration number
1163: . fgnorm - 2-norm of residual
1164: - vf  - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`

1166:   Level: intermediate

1168:   Note:
1169:   This routine prints the residual norm at each iteration.

1171: .seealso: `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
1172: @*/
1173: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1174: {
1175:   PetscViewer        viewer = vf->viewer;
1176:   Vec                res;
1177:   DM                 dm;
1178:   PetscSection       s;
1179:   const PetscScalar *r;
1180:   PetscReal         *lnorms, *norms;
1181:   PetscInt           numFields, f, pStart, pEnd, p;

1183:   PetscFunctionBegin;
1185:   PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
1186:   PetscCall(SNESGetDM(snes, &dm));
1187:   PetscCall(DMGetLocalSection(dm, &s));
1188:   PetscCall(PetscSectionGetNumFields(s, &numFields));
1189:   PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1190:   PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
1191:   PetscCall(VecGetArrayRead(res, &r));
1192:   for (p = pStart; p < pEnd; ++p) {
1193:     for (f = 0; f < numFields; ++f) {
1194:       PetscInt fdof, foff, d;

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

1217: /********************* Residual Computation **************************/

1219: PetscErrorCode DMPlexGetAllCells_Internal(DM plex, IS *cellIS)
1220: {
1221:   PetscInt depth;

1223:   PetscFunctionBegin;
1224:   PetscCall(DMPlexGetDepth(plex, &depth));
1225:   PetscCall(DMGetStratumIS(plex, "dim", depth, cellIS));
1226:   if (!*cellIS) PetscCall(DMGetStratumIS(plex, "depth", depth, cellIS));
1227:   PetscFunctionReturn(PETSC_SUCCESS);
1228: }

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

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

1238:   Output Parameter:
1239: . F  - Local output vector

1241:   Level: developer

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

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

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

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

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

1286: PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user)
1287: {
1288:   DM       plex;
1289:   IS       allcellIS;
1290:   PetscInt Nds, s;

1292:   PetscFunctionBegin;
1293:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1294:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1295:   PetscCall(DMGetNumDS(dm, &Nds));
1296:   for (s = 0; s < Nds; ++s) {
1297:     PetscDS ds;
1298:     DMLabel label;
1299:     IS      cellIS;

1301:     PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1302:     {
1303:       PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1304:       PetscWeakForm     wf;
1305:       PetscInt          Nm = 2, m, Nk = 0, k, kp, off = 0;
1306:       PetscFormKey     *reskeys;

1308:       /* Get unique residual keys */
1309:       for (m = 0; m < Nm; ++m) {
1310:         PetscInt Nkm;
1311:         PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
1312:         Nk += Nkm;
1313:       }
1314:       PetscCall(PetscMalloc1(Nk, &reskeys));
1315:       for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
1316:       PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1317:       PetscCall(PetscFormKeySort(Nk, reskeys));
1318:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1319:         if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1320:           ++k;
1321:           if (kp != k) reskeys[k] = reskeys[kp];
1322:         }
1323:       }
1324:       Nk = k;

1326:       PetscCall(PetscDSGetWeakForm(ds, &wf));
1327:       for (k = 0; k < Nk; ++k) {
1328:         DMLabel  label = reskeys[k].label;
1329:         PetscInt val   = reskeys[k].value;

1331:         if (!label) {
1332:           PetscCall(PetscObjectReference((PetscObject)allcellIS));
1333:           cellIS = allcellIS;
1334:         } else {
1335:           IS pointIS;

1337:           PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1338:           PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1339:           PetscCall(ISDestroy(&pointIS));
1340:         }
1341:         PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
1342:         PetscCall(ISDestroy(&cellIS));
1343:       }
1344:       PetscCall(PetscFree(reskeys));
1345:     }
1346:   }
1347:   PetscCall(ISDestroy(&allcellIS));
1348:   PetscCall(DMDestroy(&plex));
1349:   PetscFunctionReturn(PETSC_SUCCESS);
1350: }

1352: /*@
1353:   DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X

1355:   Input Parameters:
1356: + dm - The mesh
1357: - user - The user context

1359:   Output Parameter:
1360: . X  - Local solution

1362:   Level: developer

1364: .seealso: `DMPLEX`, `DMPlexComputeJacobianAction()`
1365: @*/
1366: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1367: {
1368:   DM plex;

1370:   PetscFunctionBegin;
1371:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1372:   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
1373:   PetscCall(DMDestroy(&plex));
1374:   PetscFunctionReturn(PETSC_SUCCESS);
1375: }

1377: /*@
1378:   DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y

1380:   Input Parameters:
1381: + dm   - The `DM`
1382: . X    - Local solution vector
1383: . Y    - Local input vector
1384: - user - The user context

1386:   Output Parameter:
1387: . F    - local output vector

1389:   Level: developer

1391:   Notes:
1392:   Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.

1394: .seealso: `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
1395: @*/
1396: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1397: {
1398:   DM       plex;
1399:   IS       allcellIS;
1400:   PetscInt Nds, s;

1402:   PetscFunctionBegin;
1403:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1404:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1405:   PetscCall(DMGetNumDS(dm, &Nds));
1406:   for (s = 0; s < Nds; ++s) {
1407:     PetscDS ds;
1408:     DMLabel label;
1409:     IS      cellIS;

1411:     PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1412:     {
1413:       PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1414:       PetscWeakForm     wf;
1415:       PetscInt          Nm = 4, m, Nk = 0, k, kp, off = 0;
1416:       PetscFormKey     *jackeys;

1418:       /* Get unique Jacobian keys */
1419:       for (m = 0; m < Nm; ++m) {
1420:         PetscInt Nkm;
1421:         PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
1422:         Nk += Nkm;
1423:       }
1424:       PetscCall(PetscMalloc1(Nk, &jackeys));
1425:       for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
1426:       PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1427:       PetscCall(PetscFormKeySort(Nk, jackeys));
1428:       for (k = 0, kp = 1; kp < Nk; ++kp) {
1429:         if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1430:           ++k;
1431:           if (kp != k) jackeys[k] = jackeys[kp];
1432:         }
1433:       }
1434:       Nk = k;

1436:       PetscCall(PetscDSGetWeakForm(ds, &wf));
1437:       for (k = 0; k < Nk; ++k) {
1438:         DMLabel  label = jackeys[k].label;
1439:         PetscInt val   = jackeys[k].value;

1441:         if (!label) {
1442:           PetscCall(PetscObjectReference((PetscObject)allcellIS));
1443:           cellIS = allcellIS;
1444:         } else {
1445:           IS pointIS;

1447:           PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1448:           PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1449:           PetscCall(ISDestroy(&pointIS));
1450:         }
1451:         PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
1452:         PetscCall(ISDestroy(&cellIS));
1453:       }
1454:       PetscCall(PetscFree(jackeys));
1455:     }
1456:   }
1457:   PetscCall(ISDestroy(&allcellIS));
1458:   PetscCall(DMDestroy(&plex));
1459:   PetscFunctionReturn(PETSC_SUCCESS);
1460: }

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

1465:   Input Parameters:
1466: + dm - The `DM`
1467: . X  - Local input vector
1468: - user - The user context

1470:   Output Parameters:
1471: + Jac  - Jacobian matrix
1472: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`

1474:   Level: developer

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

1480: .seealso: `DMPLEX`, `Mat`
1481: @*/
1482: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
1483: {
1484:   DM        plex;
1485:   IS        allcellIS;
1486:   PetscBool hasJac, hasPrec;
1487:   PetscInt  Nds, s;

1489:   PetscFunctionBegin;
1490:   PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1491:   PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1492:   PetscCall(DMGetNumDS(dm, &Nds));
1493:   for (s = 0; s < Nds; ++s) {
1494:     PetscDS      ds;
1495:     IS           cellIS;
1496:     PetscFormKey key;

1498:     PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
1499:     key.value = 0;
1500:     key.field = 0;
1501:     key.part  = 0;
1502:     if (!key.label) {
1503:       PetscCall(PetscObjectReference((PetscObject)allcellIS));
1504:       cellIS = allcellIS;
1505:     } else {
1506:       IS pointIS;

1508:       key.value = 1;
1509:       PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
1510:       PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1511:       PetscCall(ISDestroy(&pointIS));
1512:     }
1513:     if (!s) {
1514:       PetscCall(PetscDSHasJacobian(ds, &hasJac));
1515:       PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1516:       if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
1517:       PetscCall(MatZeroEntries(JacP));
1518:     }
1519:     PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
1520:     PetscCall(ISDestroy(&cellIS));
1521:   }
1522:   PetscCall(ISDestroy(&allcellIS));
1523:   PetscCall(DMDestroy(&plex));
1524:   PetscFunctionReturn(PETSC_SUCCESS);
1525: }

1527: struct _DMSNESJacobianMFCtx {
1528:   DM    dm;
1529:   Vec   X;
1530:   void *ctx;
1531: };

1533: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1534: {
1535:   struct _DMSNESJacobianMFCtx *ctx;

1537:   PetscFunctionBegin;
1538:   PetscCall(MatShellGetContext(A, &ctx));
1539:   PetscCall(MatShellSetContext(A, NULL));
1540:   PetscCall(DMDestroy(&ctx->dm));
1541:   PetscCall(VecDestroy(&ctx->X));
1542:   PetscCall(PetscFree(ctx));
1543:   PetscFunctionReturn(PETSC_SUCCESS);
1544: }

1546: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1547: {
1548:   struct _DMSNESJacobianMFCtx *ctx;

1550:   PetscFunctionBegin;
1551:   PetscCall(MatShellGetContext(A, &ctx));
1552:   PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
1553:   PetscFunctionReturn(PETSC_SUCCESS);
1554: }

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

1559:   Collective

1561:   Input Parameters:
1562: + dm   - The `DM`
1563: . X    - The evaluation point for the Jacobian
1564: - user - A user context, or `NULL`

1566:   Output Parameter:
1567: . J    - The `Mat`

1569:   Level: advanced

1571:   Note:
1572:   Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.

1574: .seealso: `DM`, `DMSNESComputeJacobianAction()`
1575: @*/
1576: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1577: {
1578:   struct _DMSNESJacobianMFCtx *ctx;
1579:   PetscInt                     n, N;

1581:   PetscFunctionBegin;
1582:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
1583:   PetscCall(MatSetType(*J, MATSHELL));
1584:   PetscCall(VecGetLocalSize(X, &n));
1585:   PetscCall(VecGetSize(X, &N));
1586:   PetscCall(MatSetSizes(*J, n, n, N, N));
1587:   PetscCall(PetscObjectReference((PetscObject)dm));
1588:   PetscCall(PetscObjectReference((PetscObject)X));
1589:   PetscCall(PetscMalloc1(1, &ctx));
1590:   ctx->dm  = dm;
1591:   ctx->X   = X;
1592:   ctx->ctx = user;
1593:   PetscCall(MatShellSetContext(*J, ctx));
1594:   PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
1595:   PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
1596:   PetscFunctionReturn(PETSC_SUCCESS);
1597: }

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

1602:    Input Parameters:
1603: +     X - `SNES` linearization point
1604: .     ovl - index set of overlapping subdomains

1606:    Output Parameter:
1607: .     J - unassembled (Neumann) local matrix

1609:    Level: intermediate

1611: .seealso: `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()`
1612: */
1613: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1614: {
1615:   SNES   snes;
1616:   Mat    pJ;
1617:   DM     ovldm, origdm;
1618:   DMSNES sdm;
1619:   PetscErrorCode (*bfun)(DM, Vec, void *);
1620:   PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
1621:   void *bctx, *jctx;

1623:   PetscFunctionBegin;
1624:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
1625:   PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
1626:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
1627:   PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
1628:   PetscCall(MatGetDM(pJ, &ovldm));
1629:   PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
1630:   PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
1631:   PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
1632:   PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
1633:   PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
1634:   if (!snes) {
1635:     PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
1636:     PetscCall(SNESSetDM(snes, ovldm));
1637:     PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
1638:     PetscCall(PetscObjectDereference((PetscObject)snes));
1639:   }
1640:   PetscCall(DMGetDMSNES(ovldm, &sdm));
1641:   PetscCall(VecLockReadPush(X));
1642:   {
1643:     void *ctx;
1644:     PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
1645:     PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
1646:     PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
1647:   }
1648:   PetscCall(VecLockReadPop(X));
1649:   /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1650:   {
1651:     Mat locpJ;

1653:     PetscCall(MatISGetLocalMat(pJ, &locpJ));
1654:     PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
1655:   }
1656:   PetscFunctionReturn(PETSC_SUCCESS);
1657: }

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

1662:   Input Parameters:
1663: + dm - The `DM` object
1664: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`)
1665: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`)
1666: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`)

1668:   Level: developer

1670: .seealso: `DMPLEX`, `SNES`
1671: @*/
1672: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1673: {
1674:   PetscFunctionBegin;
1675:   PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx));
1676:   PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx));
1677:   PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx));
1678:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
1679:   PetscFunctionReturn(PETSC_SUCCESS);
1680: }

1682: /*@C
1683:   DMSNESCheckDiscretization - Check the discretization error of the exact solution

1685:   Input Parameters:
1686: + snes - the `SNES` object
1687: . dm   - the `DM`
1688: . t    - the time
1689: . u    - a `DM` vector
1690: - tol  - A tolerance for the check, or -1 to print the results instead

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

1695:   Level: developer

1697:   Note:
1698:   The user must call `PetscDSSetExactSolution()` beforehand

1700: .seealso: `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`, `PetscDSSetExactSolution()`
1701: @*/
1702: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1703: {
1704:   PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1705:   void     **ectxs;
1706:   PetscReal *err;
1707:   MPI_Comm   comm;
1708:   PetscInt   Nf, f;

1710:   PetscFunctionBegin;

1716:   PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1717:   PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));

1719:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1720:   PetscCall(DMGetNumFields(dm, &Nf));
1721:   PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
1722:   {
1723:     PetscInt Nds, s;

1725:     PetscCall(DMGetNumDS(dm, &Nds));
1726:     for (s = 0; s < Nds; ++s) {
1727:       PetscDS         ds;
1728:       DMLabel         label;
1729:       IS              fieldIS;
1730:       const PetscInt *fields;
1731:       PetscInt        dsNf, f;

1733:       PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
1734:       PetscCall(PetscDSGetNumFields(ds, &dsNf));
1735:       PetscCall(ISGetIndices(fieldIS, &fields));
1736:       for (f = 0; f < dsNf; ++f) {
1737:         const PetscInt field = fields[f];
1738:         PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
1739:       }
1740:       PetscCall(ISRestoreIndices(fieldIS, &fields));
1741:     }
1742:   }
1743:   if (Nf > 1) {
1744:     PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
1745:     if (tol >= 0.0) {
1746:       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);
1747:     } else if (error) {
1748:       for (f = 0; f < Nf; ++f) error[f] = err[f];
1749:     } else {
1750:       PetscCall(PetscPrintf(comm, "L_2 Error: ["));
1751:       for (f = 0; f < Nf; ++f) {
1752:         if (f) PetscCall(PetscPrintf(comm, ", "));
1753:         PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
1754:       }
1755:       PetscCall(PetscPrintf(comm, "]\n"));
1756:     }
1757:   } else {
1758:     PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
1759:     if (tol >= 0.0) {
1760:       PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
1761:     } else if (error) {
1762:       error[0] = err[0];
1763:     } else {
1764:       PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
1765:     }
1766:   }
1767:   PetscCall(PetscFree3(exacts, ectxs, err));
1768:   PetscFunctionReturn(PETSC_SUCCESS);
1769: }

1771: /*@C
1772:   DMSNESCheckResidual - Check the residual of the exact solution

1774:   Input Parameters:
1775: + snes - the `SNES` object
1776: . dm   - the `DM`
1777: . u    - a `DM` vector
1778: - tol  - A tolerance for the check, or -1 to print the results instead

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

1783:   Level: developer

1785: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
1786: @*/
1787: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1788: {
1789:   MPI_Comm  comm;
1790:   Vec       r;
1791:   PetscReal res;

1793:   PetscFunctionBegin;
1798:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1799:   PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1800:   PetscCall(VecDuplicate(u, &r));
1801:   PetscCall(SNESComputeFunction(snes, u, r));
1802:   PetscCall(VecNorm(r, NORM_2, &res));
1803:   if (tol >= 0.0) {
1804:     PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
1805:   } else if (residual) {
1806:     *residual = res;
1807:   } else {
1808:     PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
1809:     PetscCall(VecChop(r, 1.0e-10));
1810:     PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
1811:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
1812:     PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
1813:   }
1814:   PetscCall(VecDestroy(&r));
1815:   PetscFunctionReturn(PETSC_SUCCESS);
1816: }

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

1821:   Input Parameters:
1822: + snes - the `SNES` object
1823: . dm   - the `DM`
1824: . u    - a `DM` vector
1825: - tol  - A tolerance for the check, or -1 to print the results instead

1827:   Output Parameters:
1828: + isLinear - Flag indicaing that the function looks linear, or `NULL`
1829: - convRate - The rate of convergence of the linear model, or `NULL`

1831:   Level: developer

1833: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1834: @*/
1835: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1836: {
1837:   MPI_Comm     comm;
1838:   PetscDS      ds;
1839:   Mat          J, M;
1840:   MatNullSpace nullspace;
1841:   PetscReal    slope, intercept;
1842:   PetscBool    hasJac, hasPrec, isLin = PETSC_FALSE;

1844:   PetscFunctionBegin;
1850:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1851:   PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1852:   /* Create and view matrices */
1853:   PetscCall(DMCreateMatrix(dm, &J));
1854:   PetscCall(DMGetDS(dm, &ds));
1855:   PetscCall(PetscDSHasJacobian(ds, &hasJac));
1856:   PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1857:   if (hasJac && hasPrec) {
1858:     PetscCall(DMCreateMatrix(dm, &M));
1859:     PetscCall(SNESComputeJacobian(snes, u, J, M));
1860:     PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
1861:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
1862:     PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
1863:     PetscCall(MatDestroy(&M));
1864:   } else {
1865:     PetscCall(SNESComputeJacobian(snes, u, J, J));
1866:   }
1867:   PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
1868:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
1869:   PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
1870:   /* Check nullspace */
1871:   PetscCall(MatGetNullSpace(J, &nullspace));
1872:   if (nullspace) {
1873:     PetscBool isNull;
1874:     PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
1875:     PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1876:   }
1877:   /* Taylor test */
1878:   {
1879:     PetscRandom rand;
1880:     Vec         du, uhat, r, rhat, df;
1881:     PetscReal   h;
1882:     PetscReal  *es, *hs, *errors;
1883:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1884:     PetscInt    Nv, v;

1886:     /* Choose a perturbation direction */
1887:     PetscCall(PetscRandomCreate(comm, &rand));
1888:     PetscCall(VecDuplicate(u, &du));
1889:     PetscCall(VecSetRandom(du, rand));
1890:     PetscCall(PetscRandomDestroy(&rand));
1891:     PetscCall(VecDuplicate(u, &df));
1892:     PetscCall(MatMult(J, du, df));
1893:     /* Evaluate residual at u, F(u), save in vector r */
1894:     PetscCall(VecDuplicate(u, &r));
1895:     PetscCall(SNESComputeFunction(snes, u, r));
1896:     /* Look at the convergence of our Taylor approximation as we approach u */
1897:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
1898:       ;
1899:     PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1900:     PetscCall(VecDuplicate(u, &uhat));
1901:     PetscCall(VecDuplicate(u, &rhat));
1902:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1903:       PetscCall(VecWAXPY(uhat, h, du, u));
1904:       /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1905:       PetscCall(SNESComputeFunction(snes, uhat, rhat));
1906:       PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1907:       PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));

1909:       es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1910:       hs[Nv] = PetscLog10Real(h);
1911:     }
1912:     PetscCall(VecDestroy(&uhat));
1913:     PetscCall(VecDestroy(&rhat));
1914:     PetscCall(VecDestroy(&df));
1915:     PetscCall(VecDestroy(&r));
1916:     PetscCall(VecDestroy(&du));
1917:     for (v = 0; v < Nv; ++v) {
1918:       if ((tol >= 0) && (errors[v] > tol)) break;
1919:       else if (errors[v] > PETSC_SMALL) break;
1920:     }
1921:     if (v == Nv) isLin = PETSC_TRUE;
1922:     PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
1923:     PetscCall(PetscFree3(es, hs, errors));
1924:     /* Slope should be about 2 */
1925:     if (tol >= 0) {
1926:       PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
1927:     } else if (isLinear || convRate) {
1928:       if (isLinear) *isLinear = isLin;
1929:       if (convRate) *convRate = slope;
1930:     } else {
1931:       if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
1932:       else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
1933:     }
1934:   }
1935:   PetscCall(MatDestroy(&J));
1936:   PetscFunctionReturn(PETSC_SUCCESS);
1937: }

1939: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1940: {
1941:   PetscFunctionBegin;
1942:   PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1943:   PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1944:   PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1945:   PetscFunctionReturn(PETSC_SUCCESS);
1946: }

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

1951:   Input Parameters:
1952: + snes - the `SNES` object
1953: - u    - representative `SNES` vector

1955:   Level: developer

1957:   Note:
1958:   The user must call `PetscDSSetExactSolution()` beforehand

1960: .seealso: `SNES`, `DM`
1961: @*/
1962: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1963: {
1964:   DM        dm;
1965:   Vec       sol;
1966:   PetscBool check;

1968:   PetscFunctionBegin;
1969:   PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1970:   if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1971:   PetscCall(SNESGetDM(snes, &dm));
1972:   PetscCall(VecDuplicate(u, &sol));
1973:   PetscCall(SNESSetSolution(snes, sol));
1974:   PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1975:   PetscCall(VecDestroy(&sol));
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }