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: }