Actual source code: dmplexsnes.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petsc/private/petscimpl.h>
5: #include <petsc/private/petscfeimpl.h>
7: #ifdef PETSC_HAVE_LIBCEED
8: #include <petscdmceed.h>
9: #include <petscdmplexceed.h>
10: #endif
12: static void pressure_Private(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar p[])
13: {
14: p[0] = u[uOff[1]];
15: }
17: /*
18: SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
19: This is normally used only to evaluate convergence rates for the pressure accurately.
21: Collective
23: Input Parameters:
24: + snes - The `SNES`
25: . pfield - The field number for pressure
26: . nullspace - The pressure nullspace
27: . u - The solution vector
28: - ctx - An optional user context
30: Output Parameter:
31: . u - The solution with a continuum pressure integral of zero
33: Level: developer
35: Note:
36: If int(u) = a and int(n) = b, then int(u - a/b n) = a - a/b b = 0. We assume that the nullspace is a single vector given explicitly.
38: .seealso: [](ch_snes), `SNESConvergedCorrectPressure()`
39: */
40: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
41: {
42: DM dm;
43: PetscDS ds;
44: const Vec *nullvecs;
45: PetscScalar pintd, *intc, *intn;
46: MPI_Comm comm;
47: PetscInt Nf, Nv;
49: PetscFunctionBegin;
50: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
51: PetscCall(SNESGetDM(snes, &dm));
52: PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
53: PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
54: PetscCall(DMGetDS(dm, &ds));
55: PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private));
56: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs));
57: PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv);
58: PetscCall(VecDot(nullvecs[0], u, &pintd));
59: PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd));
60: PetscCall(PetscDSGetNumFields(ds, &Nf));
61: PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn));
62: PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx));
63: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
64: PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]));
65: #if defined(PETSC_USE_DEBUG)
66: PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
67: PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield]));
68: #endif
69: PetscCall(PetscFree2(intc, intn));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@C
74: SNESConvergedCorrectPressure - The regular `SNES` convergence test that, up on convergence, adds a vector in the nullspace
75: to make the continuum integral of the pressure field equal to zero.
77: Logically Collective
79: Input Parameters:
80: + snes - the `SNES` context
81: . it - the iteration (0 indicates before any Newton steps)
82: . xnorm - 2-norm of current iterate
83: . gnorm - 2-norm of current step
84: . f - 2-norm of function at current iterate
85: - ctx - Optional user context
87: Output Parameter:
88: . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`
90: Options Database Key:
91: . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
93: Level: advanced
95: Notes:
96: In order to use this convergence test, you must set up several PETSc structures. First fields must be added to the `DM`, and a `PetscDS`
97: must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
98: The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
99: Last we must be able to integrate the pressure over the domain, so the `DM` attached to the SNES `must` be a `DMPLEX` at this time.
101: Developer Note:
102: This is a total misuse of the `SNES` convergence test handling system. It should be removed. Perhaps a `SNESSetPostSolve()` could
103: be constructed to handle this process.
105: .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`
106: @*/
107: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
108: {
109: PetscBool monitorIntegral = PETSC_FALSE;
111: PetscFunctionBegin;
112: PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx));
113: if (monitorIntegral) {
114: Mat J;
115: Vec u;
116: MatNullSpace nullspace;
117: const Vec *nullvecs;
118: PetscScalar pintd;
120: PetscCall(SNESGetSolution(snes, &u));
121: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
122: PetscCall(MatGetNullSpace(J, &nullspace));
123: PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs));
124: PetscCall(VecDot(nullvecs[0], u, &pintd));
125: PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd)));
126: }
127: if (*reason > 0) {
128: Mat J;
129: Vec u;
130: MatNullSpace nullspace;
131: PetscInt pfield = 1;
133: PetscCall(SNESGetSolution(snes, &u));
134: PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
135: PetscCall(MatGetNullSpace(J, &nullspace));
136: PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx));
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /************************** Interpolation *******************************/
143: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
144: {
145: PetscBool isPlex;
147: PetscFunctionBegin;
148: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
149: if (isPlex) {
150: *plex = dm;
151: PetscCall(PetscObjectReference((PetscObject)dm));
152: } else {
153: PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
154: if (!*plex) {
155: PetscCall(DMConvert(dm, DMPLEX, plex));
156: PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
157: if (copy) {
158: PetscCall(DMCopyDMSNES(dm, *plex));
159: PetscCall(DMCopyAuxiliaryVec(dm, *plex));
160: }
161: } else {
162: PetscCall(PetscObjectReference((PetscObject)*plex));
163: }
164: }
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /*@C
169: DMInterpolationCreate - Creates a `DMInterpolationInfo` context
171: Collective
173: Input Parameter:
174: . comm - the communicator
176: Output Parameter:
177: . ctx - the context
179: Level: beginner
181: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
182: @*/
183: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
184: {
185: PetscFunctionBegin;
186: PetscAssertPointer(ctx, 2);
187: PetscCall(PetscNew(ctx));
189: (*ctx)->comm = comm;
190: (*ctx)->dim = -1;
191: (*ctx)->nInput = 0;
192: (*ctx)->points = NULL;
193: (*ctx)->cells = NULL;
194: (*ctx)->n = -1;
195: (*ctx)->coords = NULL;
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*@C
200: DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
202: Not Collective
204: Input Parameters:
205: + ctx - the context
206: - dim - the spatial dimension
208: Level: intermediate
210: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
211: @*/
212: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
213: {
214: PetscFunctionBegin;
215: PetscCheck(!(dim < 1) && !(dim > 3), ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %" PetscInt_FMT, dim);
216: ctx->dim = dim;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@C
221: DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
223: Not Collective
225: Input Parameter:
226: . ctx - the context
228: Output Parameter:
229: . dim - the spatial dimension
231: Level: intermediate
233: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
234: @*/
235: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
236: {
237: PetscFunctionBegin;
238: PetscAssertPointer(dim, 2);
239: *dim = ctx->dim;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@C
244: DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
246: Not Collective
248: Input Parameters:
249: + ctx - the context
250: - dof - the number of fields
252: Level: intermediate
254: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
255: @*/
256: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
257: {
258: PetscFunctionBegin;
259: PetscCheck(dof >= 1, ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %" PetscInt_FMT, dof);
260: ctx->dof = dof;
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@C
265: DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
267: Not Collective
269: Input Parameter:
270: . ctx - the context
272: Output Parameter:
273: . dof - the number of fields
275: Level: intermediate
277: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
278: @*/
279: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
280: {
281: PetscFunctionBegin;
282: PetscAssertPointer(dof, 2);
283: *dof = ctx->dof;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@C
288: DMInterpolationAddPoints - Add points at which we will interpolate the fields
290: Not Collective
292: Input Parameters:
293: + ctx - the context
294: . n - the number of points
295: - points - the coordinates for each point, an array of size n * dim
297: Level: intermediate
299: Note:
300: The coordinate information is copied.
302: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
303: @*/
304: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
305: {
306: PetscFunctionBegin;
307: PetscCheck(ctx->dim >= 0, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
308: PetscCheck(!ctx->points, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
309: ctx->nInput = n;
311: PetscCall(PetscMalloc1(n * ctx->dim, &ctx->points));
312: PetscCall(PetscArraycpy(ctx->points, points, n * ctx->dim));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: /*@C
317: DMInterpolationSetUp - Compute spatial indices for point location during interpolation
319: Collective
321: Input Parameters:
322: + ctx - the context
323: . dm - the `DM` for the function space used for interpolation
324: . redundantPoints - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
325: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error
327: Level: intermediate
329: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
330: @*/
331: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
332: {
333: MPI_Comm comm = ctx->comm;
334: PetscScalar *a;
335: PetscInt p, q, i;
336: PetscMPIInt rank, size;
337: Vec pointVec;
338: PetscSF cellSF;
339: PetscLayout layout;
340: PetscReal *globalPoints;
341: PetscScalar *globalPointsScalar;
342: const PetscInt *ranges;
343: PetscMPIInt *counts, *displs;
344: const PetscSFNode *foundCells;
345: const PetscInt *foundPoints;
346: PetscMPIInt *foundProcs, *globalProcs;
347: PetscInt n, N, numFound;
349: PetscFunctionBegin;
351: PetscCallMPI(MPI_Comm_size(comm, &size));
352: PetscCallMPI(MPI_Comm_rank(comm, &rank));
353: PetscCheck(ctx->dim >= 0, comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
354: /* Locate points */
355: n = ctx->nInput;
356: if (!redundantPoints) {
357: PetscCall(PetscLayoutCreate(comm, &layout));
358: PetscCall(PetscLayoutSetBlockSize(layout, 1));
359: PetscCall(PetscLayoutSetLocalSize(layout, n));
360: PetscCall(PetscLayoutSetUp(layout));
361: PetscCall(PetscLayoutGetSize(layout, &N));
362: /* Communicate all points to all processes */
363: PetscCall(PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs));
364: PetscCall(PetscLayoutGetRanges(layout, &ranges));
365: for (p = 0; p < size; ++p) {
366: counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
367: displs[p] = ranges[p] * ctx->dim;
368: }
369: PetscCallMPI(MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm));
370: } else {
371: N = n;
372: globalPoints = ctx->points;
373: counts = displs = NULL;
374: layout = NULL;
375: }
376: #if 0
377: PetscCall(PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs));
378: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
379: #else
380: #if defined(PETSC_USE_COMPLEX)
381: PetscCall(PetscMalloc1(N * ctx->dim, &globalPointsScalar));
382: for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
383: #else
384: globalPointsScalar = globalPoints;
385: #endif
386: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec));
387: PetscCall(PetscMalloc2(N, &foundProcs, N, &globalProcs));
388: for (p = 0; p < N; ++p) foundProcs[p] = size;
389: cellSF = NULL;
390: PetscCall(DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF));
391: PetscCall(PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells));
392: #endif
393: for (p = 0; p < numFound; ++p) {
394: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
395: }
396: /* Let the lowest rank process own each point */
397: PetscCall(MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm));
398: ctx->n = 0;
399: for (p = 0; p < N; ++p) {
400: if (globalProcs[p] == size) {
401: PetscCheck(ignoreOutsideDomain, comm, PETSC_ERR_PLIB, "Point %" PetscInt_FMT ": %g %g %g not located in mesh", p, (double)globalPoints[p * ctx->dim + 0], (double)(ctx->dim > 1 ? globalPoints[p * ctx->dim + 1] : 0.0),
402: (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
403: if (rank == 0) ++ctx->n;
404: } else if (globalProcs[p] == rank) ++ctx->n;
405: }
406: /* Create coordinates vector and array of owned cells */
407: PetscCall(PetscMalloc1(ctx->n, &ctx->cells));
408: PetscCall(VecCreate(comm, &ctx->coords));
409: PetscCall(VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE));
410: PetscCall(VecSetBlockSize(ctx->coords, ctx->dim));
411: PetscCall(VecSetType(ctx->coords, VECSTANDARD));
412: PetscCall(VecGetArray(ctx->coords, &a));
413: for (p = 0, q = 0, i = 0; p < N; ++p) {
414: if (globalProcs[p] == rank) {
415: PetscInt d;
417: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
418: ctx->cells[q] = foundCells[q].index;
419: ++q;
420: }
421: if (globalProcs[p] == size && rank == 0) {
422: PetscInt d;
424: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
425: ctx->cells[q] = -1;
426: ++q;
427: }
428: }
429: PetscCall(VecRestoreArray(ctx->coords, &a));
430: #if 0
431: PetscCall(PetscFree3(foundCells,foundProcs,globalProcs));
432: #else
433: PetscCall(PetscFree2(foundProcs, globalProcs));
434: PetscCall(PetscSFDestroy(&cellSF));
435: PetscCall(VecDestroy(&pointVec));
436: #endif
437: if ((void *)globalPointsScalar != (void *)globalPoints) PetscCall(PetscFree(globalPointsScalar));
438: if (!redundantPoints) PetscCall(PetscFree3(globalPoints, counts, displs));
439: PetscCall(PetscLayoutDestroy(&layout));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@C
444: DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
446: Collective
448: Input Parameter:
449: . ctx - the context
451: Output Parameter:
452: . coordinates - the coordinates of interpolation points
454: Level: intermediate
456: Note:
457: The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
458: This is a borrowed vector that the user should not destroy.
460: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
461: @*/
462: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
463: {
464: PetscFunctionBegin;
465: PetscAssertPointer(coordinates, 2);
466: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
467: *coordinates = ctx->coords;
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@C
472: DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
474: Collective
476: Input Parameter:
477: . ctx - the context
479: Output Parameter:
480: . v - a vector capable of holding the interpolated field values
482: Level: intermediate
484: Note:
485: This vector should be returned using `DMInterpolationRestoreVector()`.
487: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
488: @*/
489: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
490: {
491: PetscFunctionBegin;
492: PetscAssertPointer(v, 2);
493: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
494: PetscCall(VecCreate(ctx->comm, v));
495: PetscCall(VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE));
496: PetscCall(VecSetBlockSize(*v, ctx->dof));
497: PetscCall(VecSetType(*v, VECSTANDARD));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@C
502: DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
504: Collective
506: Input Parameters:
507: + ctx - the context
508: - v - a vector capable of holding the interpolated field values
510: Level: intermediate
512: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
513: @*/
514: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
515: {
516: PetscFunctionBegin;
517: PetscAssertPointer(v, 2);
518: PetscCheck(ctx->coords, ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
519: PetscCall(VecDestroy(v));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
524: {
525: PetscReal v0, J, invJ, detJ;
526: const PetscInt dof = ctx->dof;
527: const PetscScalar *coords;
528: PetscScalar *a;
529: PetscInt p;
531: PetscFunctionBegin;
532: PetscCall(VecGetArrayRead(ctx->coords, &coords));
533: PetscCall(VecGetArray(v, &a));
534: for (p = 0; p < ctx->n; ++p) {
535: PetscInt c = ctx->cells[p];
536: PetscScalar *x = NULL;
537: PetscReal xir[1];
538: PetscInt xSize, comp;
540: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
541: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
542: xir[0] = invJ * PetscRealPart(coords[p] - v0);
543: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
544: if (2 * dof == xSize) {
545: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
546: } else if (dof == xSize) {
547: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
548: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Input closure size %" PetscInt_FMT " must be either %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 2 * dof, dof);
549: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
550: }
551: PetscCall(VecRestoreArray(v, &a));
552: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
557: {
558: PetscReal *v0, *J, *invJ, detJ;
559: const PetscScalar *coords;
560: PetscScalar *a;
561: PetscInt p;
563: PetscFunctionBegin;
564: PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
565: PetscCall(VecGetArrayRead(ctx->coords, &coords));
566: PetscCall(VecGetArray(v, &a));
567: for (p = 0; p < ctx->n; ++p) {
568: PetscInt c = ctx->cells[p];
569: PetscScalar *x = NULL;
570: PetscReal xi[4];
571: PetscInt d, f, comp;
573: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
574: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
575: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
576: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
578: for (d = 0; d < ctx->dim; ++d) {
579: xi[d] = 0.0;
580: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
581: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[(d + 1) * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
582: }
583: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
584: }
585: PetscCall(VecRestoreArray(v, &a));
586: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
587: PetscCall(PetscFree3(v0, J, invJ));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
592: {
593: PetscReal *v0, *J, *invJ, detJ;
594: const PetscScalar *coords;
595: PetscScalar *a;
596: PetscInt p;
598: PetscFunctionBegin;
599: PetscCall(PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ));
600: PetscCall(VecGetArrayRead(ctx->coords, &coords));
601: PetscCall(VecGetArray(v, &a));
602: for (p = 0; p < ctx->n; ++p) {
603: PetscInt c = ctx->cells[p];
604: const PetscInt order[3] = {2, 1, 3};
605: PetscScalar *x = NULL;
606: PetscReal xi[4];
607: PetscInt d, f, comp;
609: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
610: PetscCheck(detJ > 0.0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %" PetscInt_FMT, (double)detJ, c);
611: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x));
612: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
614: for (d = 0; d < ctx->dim; ++d) {
615: xi[d] = 0.0;
616: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
617: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] += PetscRealPart(x[order[d] * ctx->dof + comp] - x[0 * ctx->dof + comp]) * xi[d];
618: }
619: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x));
620: }
621: PetscCall(VecRestoreArray(v, &a));
622: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
623: PetscCall(PetscFree3(v0, J, invJ));
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
628: {
629: const PetscScalar *vertices = (const PetscScalar *)ctx;
630: const PetscScalar x0 = vertices[0];
631: const PetscScalar y0 = vertices[1];
632: const PetscScalar x1 = vertices[2];
633: const PetscScalar y1 = vertices[3];
634: const PetscScalar x2 = vertices[4];
635: const PetscScalar y2 = vertices[5];
636: const PetscScalar x3 = vertices[6];
637: const PetscScalar y3 = vertices[7];
638: const PetscScalar f_1 = x1 - x0;
639: const PetscScalar g_1 = y1 - y0;
640: const PetscScalar f_3 = x3 - x0;
641: const PetscScalar g_3 = y3 - y0;
642: const PetscScalar f_01 = x2 - x1 - x3 + x0;
643: const PetscScalar g_01 = y2 - y1 - y3 + y0;
644: const PetscScalar *ref;
645: PetscScalar *real;
647: PetscFunctionBegin;
648: PetscCall(VecGetArrayRead(Xref, &ref));
649: PetscCall(VecGetArray(Xreal, &real));
650: {
651: const PetscScalar p0 = ref[0];
652: const PetscScalar p1 = ref[1];
654: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
655: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
656: }
657: PetscCall(PetscLogFlops(28));
658: PetscCall(VecRestoreArrayRead(Xref, &ref));
659: PetscCall(VecRestoreArray(Xreal, &real));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: #include <petsc/private/dmimpl.h>
664: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
665: {
666: const PetscScalar *vertices = (const PetscScalar *)ctx;
667: const PetscScalar x0 = vertices[0];
668: const PetscScalar y0 = vertices[1];
669: const PetscScalar x1 = vertices[2];
670: const PetscScalar y1 = vertices[3];
671: const PetscScalar x2 = vertices[4];
672: const PetscScalar y2 = vertices[5];
673: const PetscScalar x3 = vertices[6];
674: const PetscScalar y3 = vertices[7];
675: const PetscScalar f_01 = x2 - x1 - x3 + x0;
676: const PetscScalar g_01 = y2 - y1 - y3 + y0;
677: const PetscScalar *ref;
679: PetscFunctionBegin;
680: PetscCall(VecGetArrayRead(Xref, &ref));
681: {
682: const PetscScalar x = ref[0];
683: const PetscScalar y = ref[1];
684: const PetscInt rows[2] = {0, 1};
685: PetscScalar values[4];
687: values[0] = (x1 - x0 + f_01 * y) * 0.5;
688: values[1] = (x3 - x0 + f_01 * x) * 0.5;
689: values[2] = (y1 - y0 + g_01 * y) * 0.5;
690: values[3] = (y3 - y0 + g_01 * x) * 0.5;
691: PetscCall(MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES));
692: }
693: PetscCall(PetscLogFlops(30));
694: PetscCall(VecRestoreArrayRead(Xref, &ref));
695: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
696: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
701: {
702: DM dmCoord;
703: PetscFE fem = NULL;
704: SNES snes;
705: KSP ksp;
706: PC pc;
707: Vec coordsLocal, r, ref, real;
708: Mat J;
709: PetscTabulation T = NULL;
710: const PetscScalar *coords;
711: PetscScalar *a;
712: PetscReal xir[2] = {0., 0.};
713: PetscInt Nf, p;
714: const PetscInt dof = ctx->dof;
716: PetscFunctionBegin;
717: PetscCall(DMGetNumFields(dm, &Nf));
718: if (Nf) {
719: PetscObject obj;
720: PetscClassId id;
722: PetscCall(DMGetField(dm, 0, NULL, &obj));
723: PetscCall(PetscObjectGetClassId(obj, &id));
724: if (id == PETSCFE_CLASSID) {
725: fem = (PetscFE)obj;
726: PetscCall(PetscFECreateTabulation(fem, 1, 1, xir, 0, &T));
727: }
728: }
729: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
730: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
731: PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
732: PetscCall(SNESSetOptionsPrefix(snes, "quad_interp_"));
733: PetscCall(VecCreate(PETSC_COMM_SELF, &r));
734: PetscCall(VecSetSizes(r, 2, 2));
735: PetscCall(VecSetType(r, dm->vectype));
736: PetscCall(VecDuplicate(r, &ref));
737: PetscCall(VecDuplicate(r, &real));
738: PetscCall(MatCreate(PETSC_COMM_SELF, &J));
739: PetscCall(MatSetSizes(J, 2, 2, 2, 2));
740: PetscCall(MatSetType(J, MATSEQDENSE));
741: PetscCall(MatSetUp(J));
742: PetscCall(SNESSetFunction(snes, r, QuadMap_Private, NULL));
743: PetscCall(SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL));
744: PetscCall(SNESGetKSP(snes, &ksp));
745: PetscCall(KSPGetPC(ksp, &pc));
746: PetscCall(PCSetType(pc, PCLU));
747: PetscCall(SNESSetFromOptions(snes));
749: PetscCall(VecGetArrayRead(ctx->coords, &coords));
750: PetscCall(VecGetArray(v, &a));
751: for (p = 0; p < ctx->n; ++p) {
752: PetscScalar *x = NULL, *vertices = NULL;
753: PetscScalar *xi;
754: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
756: /* Can make this do all points at once */
757: PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
758: PetscCheck(4 * 2 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %" PetscInt_FMT " should be %d", coordSize, 4 * 2);
759: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
760: PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
761: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
762: PetscCall(VecGetArray(real, &xi));
763: xi[0] = coords[p * ctx->dim + 0];
764: xi[1] = coords[p * ctx->dim + 1];
765: PetscCall(VecRestoreArray(real, &xi));
766: PetscCall(SNESSolve(snes, real, ref));
767: PetscCall(VecGetArray(ref, &xi));
768: xir[0] = PetscRealPart(xi[0]);
769: xir[1] = PetscRealPart(xi[1]);
770: if (4 * dof == xSize) {
771: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) * (1 - xir[1]) + x[1 * dof + comp] * xir[0] * (1 - xir[1]) + x[2 * dof + comp] * xir[0] * xir[1] + x[3 * dof + comp] * (1 - xir[0]) * xir[1];
772: } else if (dof == xSize) {
773: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
774: } else {
775: PetscInt d;
777: PetscCheck(fem, ctx->comm, PETSC_ERR_ARG_WRONG, "Cannot have a higher order interpolant if the discretization is not PetscFE");
778: xir[0] = 2.0 * xir[0] - 1.0;
779: xir[1] = 2.0 * xir[1] - 1.0;
780: PetscCall(PetscFEComputeTabulation(fem, 1, xir, 0, T));
781: for (comp = 0; comp < dof; ++comp) {
782: a[p * dof + comp] = 0.0;
783: for (d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
784: }
785: }
786: PetscCall(VecRestoreArray(ref, &xi));
787: PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
788: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
789: }
790: PetscCall(PetscTabulationDestroy(&T));
791: PetscCall(VecRestoreArray(v, &a));
792: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
794: PetscCall(SNESDestroy(&snes));
795: PetscCall(VecDestroy(&r));
796: PetscCall(VecDestroy(&ref));
797: PetscCall(VecDestroy(&real));
798: PetscCall(MatDestroy(&J));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
803: {
804: const PetscScalar *vertices = (const PetscScalar *)ctx;
805: const PetscScalar x0 = vertices[0];
806: const PetscScalar y0 = vertices[1];
807: const PetscScalar z0 = vertices[2];
808: const PetscScalar x1 = vertices[9];
809: const PetscScalar y1 = vertices[10];
810: const PetscScalar z1 = vertices[11];
811: const PetscScalar x2 = vertices[6];
812: const PetscScalar y2 = vertices[7];
813: const PetscScalar z2 = vertices[8];
814: const PetscScalar x3 = vertices[3];
815: const PetscScalar y3 = vertices[4];
816: const PetscScalar z3 = vertices[5];
817: const PetscScalar x4 = vertices[12];
818: const PetscScalar y4 = vertices[13];
819: const PetscScalar z4 = vertices[14];
820: const PetscScalar x5 = vertices[15];
821: const PetscScalar y5 = vertices[16];
822: const PetscScalar z5 = vertices[17];
823: const PetscScalar x6 = vertices[18];
824: const PetscScalar y6 = vertices[19];
825: const PetscScalar z6 = vertices[20];
826: const PetscScalar x7 = vertices[21];
827: const PetscScalar y7 = vertices[22];
828: const PetscScalar z7 = vertices[23];
829: const PetscScalar f_1 = x1 - x0;
830: const PetscScalar g_1 = y1 - y0;
831: const PetscScalar h_1 = z1 - z0;
832: const PetscScalar f_3 = x3 - x0;
833: const PetscScalar g_3 = y3 - y0;
834: const PetscScalar h_3 = z3 - z0;
835: const PetscScalar f_4 = x4 - x0;
836: const PetscScalar g_4 = y4 - y0;
837: const PetscScalar h_4 = z4 - z0;
838: const PetscScalar f_01 = x2 - x1 - x3 + x0;
839: const PetscScalar g_01 = y2 - y1 - y3 + y0;
840: const PetscScalar h_01 = z2 - z1 - z3 + z0;
841: const PetscScalar f_12 = x7 - x3 - x4 + x0;
842: const PetscScalar g_12 = y7 - y3 - y4 + y0;
843: const PetscScalar h_12 = z7 - z3 - z4 + z0;
844: const PetscScalar f_02 = x5 - x1 - x4 + x0;
845: const PetscScalar g_02 = y5 - y1 - y4 + y0;
846: const PetscScalar h_02 = z5 - z1 - z4 + z0;
847: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
848: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
849: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
850: const PetscScalar *ref;
851: PetscScalar *real;
853: PetscFunctionBegin;
854: PetscCall(VecGetArrayRead(Xref, &ref));
855: PetscCall(VecGetArray(Xreal, &real));
856: {
857: const PetscScalar p0 = ref[0];
858: const PetscScalar p1 = ref[1];
859: const PetscScalar p2 = ref[2];
861: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_4 * p2 + f_01 * p0 * p1 + f_12 * p1 * p2 + f_02 * p0 * p2 + f_012 * p0 * p1 * p2;
862: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_4 * p2 + g_01 * p0 * p1 + g_01 * p0 * p1 + g_12 * p1 * p2 + g_02 * p0 * p2 + g_012 * p0 * p1 * p2;
863: real[2] = z0 + h_1 * p0 + h_3 * p1 + h_4 * p2 + h_01 * p0 * p1 + h_01 * p0 * p1 + h_12 * p1 * p2 + h_02 * p0 * p2 + h_012 * p0 * p1 * p2;
864: }
865: PetscCall(PetscLogFlops(114));
866: PetscCall(VecRestoreArrayRead(Xref, &ref));
867: PetscCall(VecRestoreArray(Xreal, &real));
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
872: {
873: const PetscScalar *vertices = (const PetscScalar *)ctx;
874: const PetscScalar x0 = vertices[0];
875: const PetscScalar y0 = vertices[1];
876: const PetscScalar z0 = vertices[2];
877: const PetscScalar x1 = vertices[9];
878: const PetscScalar y1 = vertices[10];
879: const PetscScalar z1 = vertices[11];
880: const PetscScalar x2 = vertices[6];
881: const PetscScalar y2 = vertices[7];
882: const PetscScalar z2 = vertices[8];
883: const PetscScalar x3 = vertices[3];
884: const PetscScalar y3 = vertices[4];
885: const PetscScalar z3 = vertices[5];
886: const PetscScalar x4 = vertices[12];
887: const PetscScalar y4 = vertices[13];
888: const PetscScalar z4 = vertices[14];
889: const PetscScalar x5 = vertices[15];
890: const PetscScalar y5 = vertices[16];
891: const PetscScalar z5 = vertices[17];
892: const PetscScalar x6 = vertices[18];
893: const PetscScalar y6 = vertices[19];
894: const PetscScalar z6 = vertices[20];
895: const PetscScalar x7 = vertices[21];
896: const PetscScalar y7 = vertices[22];
897: const PetscScalar z7 = vertices[23];
898: const PetscScalar f_xy = x2 - x1 - x3 + x0;
899: const PetscScalar g_xy = y2 - y1 - y3 + y0;
900: const PetscScalar h_xy = z2 - z1 - z3 + z0;
901: const PetscScalar f_yz = x7 - x3 - x4 + x0;
902: const PetscScalar g_yz = y7 - y3 - y4 + y0;
903: const PetscScalar h_yz = z7 - z3 - z4 + z0;
904: const PetscScalar f_xz = x5 - x1 - x4 + x0;
905: const PetscScalar g_xz = y5 - y1 - y4 + y0;
906: const PetscScalar h_xz = z5 - z1 - z4 + z0;
907: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
908: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
909: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
910: const PetscScalar *ref;
912: PetscFunctionBegin;
913: PetscCall(VecGetArrayRead(Xref, &ref));
914: {
915: const PetscScalar x = ref[0];
916: const PetscScalar y = ref[1];
917: const PetscScalar z = ref[2];
918: const PetscInt rows[3] = {0, 1, 2};
919: PetscScalar values[9];
921: values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
922: values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
923: values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
924: values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
925: values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
926: values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
927: values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
928: values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
929: values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;
931: PetscCall(MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES));
932: }
933: PetscCall(PetscLogFlops(152));
934: PetscCall(VecRestoreArrayRead(Xref, &ref));
935: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
936: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
941: {
942: DM dmCoord;
943: SNES snes;
944: KSP ksp;
945: PC pc;
946: Vec coordsLocal, r, ref, real;
947: Mat J;
948: const PetscScalar *coords;
949: PetscScalar *a;
950: PetscInt p;
952: PetscFunctionBegin;
953: PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
954: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
955: PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
956: PetscCall(SNESSetOptionsPrefix(snes, "hex_interp_"));
957: PetscCall(VecCreate(PETSC_COMM_SELF, &r));
958: PetscCall(VecSetSizes(r, 3, 3));
959: PetscCall(VecSetType(r, dm->vectype));
960: PetscCall(VecDuplicate(r, &ref));
961: PetscCall(VecDuplicate(r, &real));
962: PetscCall(MatCreate(PETSC_COMM_SELF, &J));
963: PetscCall(MatSetSizes(J, 3, 3, 3, 3));
964: PetscCall(MatSetType(J, MATSEQDENSE));
965: PetscCall(MatSetUp(J));
966: PetscCall(SNESSetFunction(snes, r, HexMap_Private, NULL));
967: PetscCall(SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL));
968: PetscCall(SNESGetKSP(snes, &ksp));
969: PetscCall(KSPGetPC(ksp, &pc));
970: PetscCall(PCSetType(pc, PCLU));
971: PetscCall(SNESSetFromOptions(snes));
973: PetscCall(VecGetArrayRead(ctx->coords, &coords));
974: PetscCall(VecGetArray(v, &a));
975: for (p = 0; p < ctx->n; ++p) {
976: PetscScalar *x = NULL, *vertices = NULL;
977: PetscScalar *xi;
978: PetscReal xir[3];
979: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
981: /* Can make this do all points at once */
982: PetscCall(DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
983: PetscCheck(8 * 3 == coordSize, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid coordinate closure size %" PetscInt_FMT " should be %d", coordSize, 8 * 3);
984: PetscCall(DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x));
985: PetscCheck((8 * ctx->dof == xSize) || (ctx->dof == xSize), ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input closure size %" PetscInt_FMT " should be %" PetscInt_FMT " or %" PetscInt_FMT, xSize, 8 * ctx->dof, ctx->dof);
986: PetscCall(SNESSetFunction(snes, NULL, NULL, vertices));
987: PetscCall(SNESSetJacobian(snes, NULL, NULL, NULL, vertices));
988: PetscCall(VecGetArray(real, &xi));
989: xi[0] = coords[p * ctx->dim + 0];
990: xi[1] = coords[p * ctx->dim + 1];
991: xi[2] = coords[p * ctx->dim + 2];
992: PetscCall(VecRestoreArray(real, &xi));
993: PetscCall(SNESSolve(snes, real, ref));
994: PetscCall(VecGetArray(ref, &xi));
995: xir[0] = PetscRealPart(xi[0]);
996: xir[1] = PetscRealPart(xi[1]);
997: xir[2] = PetscRealPart(xi[2]);
998: if (8 * ctx->dof == xSize) {
999: for (comp = 0; comp < ctx->dof; ++comp) {
1000: a[p * ctx->dof + comp] = x[0 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * (1 - xir[2]) + x[3 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * (1 - xir[2]) + x[2 * ctx->dof + comp] * xir[0] * xir[1] * (1 - xir[2]) + x[1 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * (1 - xir[2]) +
1001: x[4 * ctx->dof + comp] * (1 - xir[0]) * (1 - xir[1]) * xir[2] + x[5 * ctx->dof + comp] * xir[0] * (1 - xir[1]) * xir[2] + x[6 * ctx->dof + comp] * xir[0] * xir[1] * xir[2] + x[7 * ctx->dof + comp] * (1 - xir[0]) * xir[1] * xir[2];
1002: }
1003: } else {
1004: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
1005: }
1006: PetscCall(VecRestoreArray(ref, &xi));
1007: PetscCall(DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices));
1008: PetscCall(DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x));
1009: }
1010: PetscCall(VecRestoreArray(v, &a));
1011: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
1013: PetscCall(SNESDestroy(&snes));
1014: PetscCall(VecDestroy(&r));
1015: PetscCall(VecDestroy(&ref));
1016: PetscCall(VecDestroy(&real));
1017: PetscCall(MatDestroy(&J));
1018: PetscFunctionReturn(PETSC_SUCCESS);
1019: }
1021: /*@C
1022: DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.
1024: Input Parameters:
1025: + ctx - The `DMInterpolationInfo` context
1026: . dm - The `DM`
1027: - x - The local vector containing the field to be interpolated
1029: Output Parameter:
1030: . v - The vector containing the interpolated values
1032: Level: beginner
1034: Note:
1035: A suitable `v` can be obtained using `DMInterpolationGetVector()`.
1037: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1038: @*/
1039: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1040: {
1041: PetscDS ds;
1042: PetscInt n, p, Nf, field;
1043: PetscBool useDS = PETSC_FALSE;
1045: PetscFunctionBegin;
1049: PetscCall(VecGetLocalSize(v, &n));
1050: PetscCheck(n == ctx->n * ctx->dof, ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %" PetscInt_FMT " should be %" PetscInt_FMT, n, ctx->n * ctx->dof);
1051: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
1052: PetscCall(DMGetDS(dm, &ds));
1053: if (ds) {
1054: useDS = PETSC_TRUE;
1055: PetscCall(PetscDSGetNumFields(ds, &Nf));
1056: for (field = 0; field < Nf; ++field) {
1057: PetscObject obj;
1058: PetscClassId id;
1060: PetscCall(PetscDSGetDiscretization(ds, field, &obj));
1061: PetscCall(PetscObjectGetClassId(obj, &id));
1062: if (id != PETSCFE_CLASSID) {
1063: useDS = PETSC_FALSE;
1064: break;
1065: }
1066: }
1067: }
1068: if (useDS) {
1069: const PetscScalar *coords;
1070: PetscScalar *interpolant;
1071: PetscInt cdim, d;
1073: PetscCall(DMGetCoordinateDim(dm, &cdim));
1074: PetscCall(VecGetArrayRead(ctx->coords, &coords));
1075: PetscCall(VecGetArrayWrite(v, &interpolant));
1076: for (p = 0; p < ctx->n; ++p) {
1077: PetscReal pcoords[3], xi[3];
1078: PetscScalar *xa = NULL;
1079: PetscInt coff = 0, foff = 0, clSize;
1081: if (ctx->cells[p] < 0) continue;
1082: for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
1083: PetscCall(DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi));
1084: PetscCall(DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1085: for (field = 0; field < Nf; ++field) {
1086: PetscTabulation T;
1087: PetscFE fe;
1089: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1090: PetscCall(PetscFECreateTabulation(fe, 1, 1, xi, 0, &T));
1091: {
1092: const PetscReal *basis = T->T[0];
1093: const PetscInt Nb = T->Nb;
1094: const PetscInt Nc = T->Nc;
1095: PetscInt f, fc;
1097: for (fc = 0; fc < Nc; ++fc) {
1098: interpolant[p * ctx->dof + coff + fc] = 0.0;
1099: for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
1100: }
1101: coff += Nc;
1102: foff += Nb;
1103: }
1104: PetscCall(PetscTabulationDestroy(&T));
1105: }
1106: PetscCall(DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa));
1107: PetscCheck(coff == ctx->dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %" PetscInt_FMT " != %" PetscInt_FMT " dof specified for interpolation", coff, ctx->dof);
1108: PetscCheck(foff == clSize, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total FE space size %" PetscInt_FMT " != %" PetscInt_FMT " closure size", foff, clSize);
1109: }
1110: PetscCall(VecRestoreArrayRead(ctx->coords, &coords));
1111: PetscCall(VecRestoreArrayWrite(v, &interpolant));
1112: } else {
1113: DMPolytopeType ct;
1115: /* TODO Check each cell individually */
1116: PetscCall(DMPlexGetCellType(dm, ctx->cells[0], &ct));
1117: switch (ct) {
1118: case DM_POLYTOPE_SEGMENT:
1119: PetscCall(DMInterpolate_Segment_Private(ctx, dm, x, v));
1120: break;
1121: case DM_POLYTOPE_TRIANGLE:
1122: PetscCall(DMInterpolate_Triangle_Private(ctx, dm, x, v));
1123: break;
1124: case DM_POLYTOPE_QUADRILATERAL:
1125: PetscCall(DMInterpolate_Quad_Private(ctx, dm, x, v));
1126: break;
1127: case DM_POLYTOPE_TETRAHEDRON:
1128: PetscCall(DMInterpolate_Tetrahedron_Private(ctx, dm, x, v));
1129: break;
1130: case DM_POLYTOPE_HEXAHEDRON:
1131: PetscCall(DMInterpolate_Hex_Private(ctx, dm, x, v));
1132: break;
1133: default:
1134: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
1135: }
1136: }
1137: PetscFunctionReturn(PETSC_SUCCESS);
1138: }
1140: /*@C
1141: DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
1143: Collective
1145: Input Parameter:
1146: . ctx - the context
1148: Level: beginner
1150: .seealso: [](ch_snes), `DM`, `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1151: @*/
1152: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1153: {
1154: PetscFunctionBegin;
1155: PetscAssertPointer(ctx, 1);
1156: PetscCall(VecDestroy(&(*ctx)->coords));
1157: PetscCall(PetscFree((*ctx)->points));
1158: PetscCall(PetscFree((*ctx)->cells));
1159: PetscCall(PetscFree(*ctx));
1160: *ctx = NULL;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: /*@C
1165: SNESMonitorFields - Monitors the residual for each field separately
1167: Collective
1169: Input Parameters:
1170: + snes - the `SNES` context, must have an attached `DM`
1171: . its - iteration number
1172: . fgnorm - 2-norm of residual
1173: - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`
1175: Level: intermediate
1177: Note:
1178: This routine prints the residual norm at each iteration.
1180: .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
1181: @*/
1182: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1183: {
1184: PetscViewer viewer = vf->viewer;
1185: Vec res;
1186: DM dm;
1187: PetscSection s;
1188: const PetscScalar *r;
1189: PetscReal *lnorms, *norms;
1190: PetscInt numFields, f, pStart, pEnd, p;
1192: PetscFunctionBegin;
1194: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
1195: PetscCall(SNESGetDM(snes, &dm));
1196: PetscCall(DMGetLocalSection(dm, &s));
1197: PetscCall(PetscSectionGetNumFields(s, &numFields));
1198: PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
1199: PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
1200: PetscCall(VecGetArrayRead(res, &r));
1201: for (p = pStart; p < pEnd; ++p) {
1202: for (f = 0; f < numFields; ++f) {
1203: PetscInt fdof, foff, d;
1205: PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
1206: PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
1207: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
1208: }
1209: }
1210: PetscCall(VecRestoreArrayRead(res, &r));
1211: PetscCall(MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
1212: PetscCall(PetscViewerPushFormat(viewer, vf->format));
1213: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
1214: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
1215: for (f = 0; f < numFields; ++f) {
1216: if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
1217: PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
1218: }
1219: PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
1220: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
1221: PetscCall(PetscViewerPopFormat(viewer));
1222: PetscCall(PetscFree2(lnorms, norms));
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: /********************* Residual Computation **************************/
1228: /*@
1229: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
1231: Input Parameters:
1232: + dm - The mesh
1233: . X - Local solution
1234: - user - The user context
1236: Output Parameter:
1237: . F - Local output vector
1239: Level: developer
1241: Note:
1242: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
1244: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1245: @*/
1246: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1247: {
1248: DM plex;
1249: IS allcellIS;
1250: PetscInt Nds, s;
1252: PetscFunctionBegin;
1253: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1254: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1255: PetscCall(DMGetNumDS(dm, &Nds));
1256: for (s = 0; s < Nds; ++s) {
1257: PetscDS ds;
1258: IS cellIS;
1259: PetscFormKey key;
1261: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
1262: key.value = 0;
1263: key.field = 0;
1264: key.part = 0;
1265: if (!key.label) {
1266: PetscCall(PetscObjectReference((PetscObject)allcellIS));
1267: cellIS = allcellIS;
1268: } else {
1269: IS pointIS;
1271: key.value = 1;
1272: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
1273: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1274: PetscCall(ISDestroy(&pointIS));
1275: }
1276: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
1277: PetscCall(ISDestroy(&cellIS));
1278: }
1279: PetscCall(ISDestroy(&allcellIS));
1280: PetscCall(DMDestroy(&plex));
1281: PetscFunctionReturn(PETSC_SUCCESS);
1282: }
1284: /*@
1285: DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
1287: Input Parameters:
1288: + dm - The mesh
1289: . X - Local solution
1290: - user - The user context
1292: Output Parameter:
1293: . F - Local output vector
1295: Level: developer
1297: Note:
1298: The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
1300: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1301: @*/
1302: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user)
1303: {
1304: DM plex;
1305: IS allcellIS;
1306: PetscInt Nds, s;
1308: PetscFunctionBegin;
1309: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1310: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1311: PetscCall(DMGetNumDS(dm, &Nds));
1312: for (s = 0; s < Nds; ++s) {
1313: PetscDS ds;
1314: DMLabel label;
1315: IS cellIS;
1317: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1318: {
1319: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1320: PetscWeakForm wf;
1321: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
1322: PetscFormKey *reskeys;
1324: /* Get unique residual keys */
1325: for (m = 0; m < Nm; ++m) {
1326: PetscInt Nkm;
1327: PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
1328: Nk += Nkm;
1329: }
1330: PetscCall(PetscMalloc1(Nk, &reskeys));
1331: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
1332: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1333: PetscCall(PetscFormKeySort(Nk, reskeys));
1334: for (k = 0, kp = 1; kp < Nk; ++kp) {
1335: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1336: ++k;
1337: if (kp != k) reskeys[k] = reskeys[kp];
1338: }
1339: }
1340: Nk = k;
1342: PetscCall(PetscDSGetWeakForm(ds, &wf));
1343: for (k = 0; k < Nk; ++k) {
1344: DMLabel label = reskeys[k].label;
1345: PetscInt val = reskeys[k].value;
1347: if (!label) {
1348: PetscCall(PetscObjectReference((PetscObject)allcellIS));
1349: cellIS = allcellIS;
1350: } else {
1351: IS pointIS;
1353: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1354: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1355: PetscCall(ISDestroy(&pointIS));
1356: }
1357: PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
1358: PetscCall(ISDestroy(&cellIS));
1359: }
1360: PetscCall(PetscFree(reskeys));
1361: }
1362: }
1363: PetscCall(ISDestroy(&allcellIS));
1364: PetscCall(DMDestroy(&plex));
1365: PetscFunctionReturn(PETSC_SUCCESS);
1366: }
1368: #ifdef PETSC_HAVE_LIBCEED
1369: PetscErrorCode DMPlexSNESComputeResidualCEED(DM dm, Vec locX, Vec locF, void *user)
1370: {
1371: Ceed ceed;
1372: DMCeed sd = dm->dmceed;
1373: CeedVector clocX, clocF;
1375: PetscFunctionBegin;
1376: PetscCall(DMGetCeed(dm, &ceed));
1377: PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
1378: PetscCall(DMCeedComputeGeometry(dm, sd));
1380: PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
1381: PetscCall(VecGetCeedVector(locF, ceed, &clocF));
1382: PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
1383: PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
1384: PetscCall(VecRestoreCeedVector(locF, &clocF));
1386: {
1387: DM_Plex *mesh = (DM_Plex *)dm->data;
1389: if (mesh->printFEM) {
1390: PetscSection section;
1391: Vec locFbc;
1392: PetscInt pStart, pEnd, p, maxDof;
1393: PetscScalar *zeroes;
1395: PetscCall(DMGetLocalSection(dm, §ion));
1396: PetscCall(VecDuplicate(locF, &locFbc));
1397: PetscCall(VecCopy(locF, locFbc));
1398: PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
1399: PetscCall(PetscSectionGetMaxDof(section, &maxDof));
1400: PetscCall(PetscCalloc1(maxDof, &zeroes));
1401: for (p = pStart; p < pEnd; ++p) PetscCall(VecSetValuesSection(locFbc, section, p, zeroes, INSERT_BC_VALUES));
1402: PetscCall(PetscFree(zeroes));
1403: PetscCall(DMPrintLocalVec(dm, "Residual", mesh->printTol, locFbc));
1404: PetscCall(VecDestroy(&locFbc));
1405: }
1406: }
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1409: #endif
1411: /*@
1412: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
1414: Input Parameters:
1415: + dm - The mesh
1416: - user - The user context
1418: Output Parameter:
1419: . X - Local solution
1421: Level: developer
1423: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
1424: @*/
1425: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1426: {
1427: DM plex;
1429: PetscFunctionBegin;
1430: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1431: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
1432: PetscCall(DMDestroy(&plex));
1433: PetscFunctionReturn(PETSC_SUCCESS);
1434: }
1436: /*@
1437: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
1439: Input Parameters:
1440: + dm - The `DM`
1441: . X - Local solution vector
1442: . Y - Local input vector
1443: - user - The user context
1445: Output Parameter:
1446: . F - local output vector
1448: Level: developer
1450: Note:
1451: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
1453: This only works with `DMPLEX`
1455: Developer Note:
1456: This should be called `DMPlexSNESComputeJacobianAction()`
1458: .seealso: [](ch_snes), `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
1459: @*/
1460: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1461: {
1462: DM plex;
1463: IS allcellIS;
1464: PetscInt Nds, s;
1466: PetscFunctionBegin;
1467: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1468: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1469: PetscCall(DMGetNumDS(dm, &Nds));
1470: for (s = 0; s < Nds; ++s) {
1471: PetscDS ds;
1472: DMLabel label;
1473: IS cellIS;
1475: PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
1476: {
1477: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1478: PetscWeakForm wf;
1479: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
1480: PetscFormKey *jackeys;
1482: /* Get unique Jacobian keys */
1483: for (m = 0; m < Nm; ++m) {
1484: PetscInt Nkm;
1485: PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
1486: Nk += Nkm;
1487: }
1488: PetscCall(PetscMalloc1(Nk, &jackeys));
1489: for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
1490: PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
1491: PetscCall(PetscFormKeySort(Nk, jackeys));
1492: for (k = 0, kp = 1; kp < Nk; ++kp) {
1493: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1494: ++k;
1495: if (kp != k) jackeys[k] = jackeys[kp];
1496: }
1497: }
1498: Nk = k;
1500: PetscCall(PetscDSGetWeakForm(ds, &wf));
1501: for (k = 0; k < Nk; ++k) {
1502: DMLabel label = jackeys[k].label;
1503: PetscInt val = jackeys[k].value;
1505: if (!label) {
1506: PetscCall(PetscObjectReference((PetscObject)allcellIS));
1507: cellIS = allcellIS;
1508: } else {
1509: IS pointIS;
1511: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
1512: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1513: PetscCall(ISDestroy(&pointIS));
1514: }
1515: PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
1516: PetscCall(ISDestroy(&cellIS));
1517: }
1518: PetscCall(PetscFree(jackeys));
1519: }
1520: }
1521: PetscCall(ISDestroy(&allcellIS));
1522: PetscCall(DMDestroy(&plex));
1523: PetscFunctionReturn(PETSC_SUCCESS);
1524: }
1526: /*@
1527: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
1529: Input Parameters:
1530: + dm - The `DM`
1531: . X - Local input vector
1532: - user - The user context
1534: Output Parameters:
1535: + Jac - Jacobian matrix
1536: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
1538: Level: developer
1540: Note:
1541: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1542: like a GPU, or vectorize on a multicore machine.
1544: .seealso: [](ch_snes), `DMPLEX`, `Mat`
1545: @*/
1546: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
1547: {
1548: DM plex;
1549: IS allcellIS;
1550: PetscBool hasJac, hasPrec;
1551: PetscInt Nds, s;
1553: PetscFunctionBegin;
1554: PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
1555: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
1556: PetscCall(DMGetNumDS(dm, &Nds));
1557: for (s = 0; s < Nds; ++s) {
1558: PetscDS ds;
1559: IS cellIS;
1560: PetscFormKey key;
1562: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
1563: key.value = 0;
1564: key.field = 0;
1565: key.part = 0;
1566: if (!key.label) {
1567: PetscCall(PetscObjectReference((PetscObject)allcellIS));
1568: cellIS = allcellIS;
1569: } else {
1570: IS pointIS;
1572: key.value = 1;
1573: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
1574: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
1575: PetscCall(ISDestroy(&pointIS));
1576: }
1577: if (!s) {
1578: PetscCall(PetscDSHasJacobian(ds, &hasJac));
1579: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1580: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
1581: PetscCall(MatZeroEntries(JacP));
1582: }
1583: PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
1584: PetscCall(ISDestroy(&cellIS));
1585: }
1586: PetscCall(ISDestroy(&allcellIS));
1587: PetscCall(DMDestroy(&plex));
1588: PetscFunctionReturn(PETSC_SUCCESS);
1589: }
1591: struct _DMSNESJacobianMFCtx {
1592: DM dm;
1593: Vec X;
1594: void *ctx;
1595: };
1597: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1598: {
1599: struct _DMSNESJacobianMFCtx *ctx;
1601: PetscFunctionBegin;
1602: PetscCall(MatShellGetContext(A, &ctx));
1603: PetscCall(MatShellSetContext(A, NULL));
1604: PetscCall(DMDestroy(&ctx->dm));
1605: PetscCall(VecDestroy(&ctx->X));
1606: PetscCall(PetscFree(ctx));
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1611: {
1612: struct _DMSNESJacobianMFCtx *ctx;
1614: PetscFunctionBegin;
1615: PetscCall(MatShellGetContext(A, &ctx));
1616: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
1617: PetscFunctionReturn(PETSC_SUCCESS);
1618: }
1620: /*@
1621: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
1623: Collective
1625: Input Parameters:
1626: + dm - The `DM`
1627: . X - The evaluation point for the Jacobian
1628: - user - A user context, or `NULL`
1630: Output Parameter:
1631: . J - The `Mat`
1633: Level: advanced
1635: Notes:
1636: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
1638: This only works for `DMPLEX`
1640: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
1641: @*/
1642: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1643: {
1644: struct _DMSNESJacobianMFCtx *ctx;
1645: PetscInt n, N;
1647: PetscFunctionBegin;
1648: PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
1649: PetscCall(MatSetType(*J, MATSHELL));
1650: PetscCall(VecGetLocalSize(X, &n));
1651: PetscCall(VecGetSize(X, &N));
1652: PetscCall(MatSetSizes(*J, n, n, N, N));
1653: PetscCall(PetscObjectReference((PetscObject)dm));
1654: PetscCall(PetscObjectReference((PetscObject)X));
1655: PetscCall(PetscMalloc1(1, &ctx));
1656: ctx->dm = dm;
1657: ctx->X = X;
1658: ctx->ctx = user;
1659: PetscCall(MatShellSetContext(*J, ctx));
1660: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
1661: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
1662: PetscFunctionReturn(PETSC_SUCCESS);
1663: }
1665: /*
1666: MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.
1668: Input Parameters:
1669: + X - `SNES` linearization point
1670: . ovl - index set of overlapping subdomains
1672: Output Parameter:
1673: . J - unassembled (Neumann) local matrix
1675: Level: intermediate
1677: .seealso: [](ch_snes), `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()`
1678: */
1679: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1680: {
1681: SNES snes;
1682: Mat pJ;
1683: DM ovldm, origdm;
1684: DMSNES sdm;
1685: PetscErrorCode (*bfun)(DM, Vec, void *);
1686: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
1687: void *bctx, *jctx;
1689: PetscFunctionBegin;
1690: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
1691: PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
1692: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
1693: PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
1694: PetscCall(MatGetDM(pJ, &ovldm));
1695: PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
1696: PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
1697: PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
1698: PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
1699: PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
1700: if (!snes) {
1701: PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
1702: PetscCall(SNESSetDM(snes, ovldm));
1703: PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
1704: PetscCall(PetscObjectDereference((PetscObject)snes));
1705: }
1706: PetscCall(DMGetDMSNES(ovldm, &sdm));
1707: PetscCall(VecLockReadPush(X));
1708: {
1709: void *ctx;
1710: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
1711: PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
1712: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
1713: }
1714: PetscCall(VecLockReadPop(X));
1715: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1716: {
1717: Mat locpJ;
1719: PetscCall(MatISGetLocalMat(pJ, &locpJ));
1720: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
1721: }
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }
1725: /*@
1726: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian.
1728: Input Parameters:
1729: + dm - The `DM` object
1730: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`)
1731: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`)
1732: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`)
1734: Level: developer
1736: .seealso: [](ch_snes), `DMPLEX`, `SNES`
1737: @*/
1738: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1739: {
1740: PetscBool useCeed;
1742: PetscFunctionBegin;
1743: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1744: PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx));
1745: if (useCeed) {
1746: #ifdef PETSC_HAVE_LIBCEED
1747: PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, residualctx));
1748: #else
1749: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
1750: #endif
1751: } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx));
1752: PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx));
1753: PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: /*@C
1758: DMSNESCheckDiscretization - Check the discretization error of the exact solution
1760: Input Parameters:
1761: + snes - the `SNES` object
1762: . dm - the `DM`
1763: . t - the time
1764: . u - a `DM` vector
1765: - tol - A tolerance for the check, or -1 to print the results instead
1767: Output Parameter:
1768: . error - An array which holds the discretization error in each field, or `NULL`
1770: Level: developer
1772: Note:
1773: The user must call `PetscDSSetExactSolution()` beforehand
1775: Developer Note:
1776: How is this related to `PetscConvEst`?
1778: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
1779: @*/
1780: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1781: {
1782: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1783: void **ectxs;
1784: PetscReal *err;
1785: MPI_Comm comm;
1786: PetscInt Nf, f;
1788: PetscFunctionBegin;
1792: if (error) PetscAssertPointer(error, 6);
1794: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
1795: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
1797: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1798: PetscCall(DMGetNumFields(dm, &Nf));
1799: PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
1800: {
1801: PetscInt Nds, s;
1803: PetscCall(DMGetNumDS(dm, &Nds));
1804: for (s = 0; s < Nds; ++s) {
1805: PetscDS ds;
1806: DMLabel label;
1807: IS fieldIS;
1808: const PetscInt *fields;
1809: PetscInt dsNf, f;
1811: PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
1812: PetscCall(PetscDSGetNumFields(ds, &dsNf));
1813: PetscCall(ISGetIndices(fieldIS, &fields));
1814: for (f = 0; f < dsNf; ++f) {
1815: const PetscInt field = fields[f];
1816: PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
1817: }
1818: PetscCall(ISRestoreIndices(fieldIS, &fields));
1819: }
1820: }
1821: if (Nf > 1) {
1822: PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
1823: if (tol >= 0.0) {
1824: for (f = 0; f < Nf; ++f) PetscCheck(err[f] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %" PetscInt_FMT " exceeds tolerance %g", (double)err[f], f, (double)tol);
1825: } else if (error) {
1826: for (f = 0; f < Nf; ++f) error[f] = err[f];
1827: } else {
1828: PetscCall(PetscPrintf(comm, "L_2 Error: ["));
1829: for (f = 0; f < Nf; ++f) {
1830: if (f) PetscCall(PetscPrintf(comm, ", "));
1831: PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
1832: }
1833: PetscCall(PetscPrintf(comm, "]\n"));
1834: }
1835: } else {
1836: PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
1837: if (tol >= 0.0) {
1838: PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
1839: } else if (error) {
1840: error[0] = err[0];
1841: } else {
1842: PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
1843: }
1844: }
1845: PetscCall(PetscFree3(exacts, ectxs, err));
1846: PetscFunctionReturn(PETSC_SUCCESS);
1847: }
1849: /*@C
1850: DMSNESCheckResidual - Check the residual of the exact solution
1852: Input Parameters:
1853: + snes - the `SNES` object
1854: . dm - the `DM`
1855: . u - a `DM` vector
1856: - tol - A tolerance for the check, or -1 to print the results instead
1858: Output Parameter:
1859: . residual - The residual norm of the exact solution, or `NULL`
1861: Level: developer
1863: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
1864: @*/
1865: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1866: {
1867: MPI_Comm comm;
1868: Vec r;
1869: PetscReal res;
1871: PetscFunctionBegin;
1875: if (residual) PetscAssertPointer(residual, 5);
1876: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1877: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1878: PetscCall(VecDuplicate(u, &r));
1879: PetscCall(SNESComputeFunction(snes, u, r));
1880: PetscCall(VecNorm(r, NORM_2, &res));
1881: if (tol >= 0.0) {
1882: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
1883: } else if (residual) {
1884: *residual = res;
1885: } else {
1886: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
1887: PetscCall(VecFilter(r, 1.0e-10));
1888: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
1889: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
1890: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
1891: }
1892: PetscCall(VecDestroy(&r));
1893: PetscFunctionReturn(PETSC_SUCCESS);
1894: }
1896: /*@C
1897: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
1899: Input Parameters:
1900: + snes - the `SNES` object
1901: . dm - the `DM`
1902: . u - a `DM` vector
1903: - tol - A tolerance for the check, or -1 to print the results instead
1905: Output Parameters:
1906: + isLinear - Flag indicaing that the function looks linear, or `NULL`
1907: - convRate - The rate of convergence of the linear model, or `NULL`
1909: Level: developer
1911: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1912: @*/
1913: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1914: {
1915: MPI_Comm comm;
1916: PetscDS ds;
1917: Mat J, M;
1918: MatNullSpace nullspace;
1919: PetscReal slope, intercept;
1920: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
1922: PetscFunctionBegin;
1926: if (isLinear) PetscAssertPointer(isLinear, 5);
1927: if (convRate) PetscAssertPointer(convRate, 6);
1928: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1929: PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
1930: /* Create and view matrices */
1931: PetscCall(DMCreateMatrix(dm, &J));
1932: PetscCall(DMGetDS(dm, &ds));
1933: PetscCall(PetscDSHasJacobian(ds, &hasJac));
1934: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
1935: if (hasJac && hasPrec) {
1936: PetscCall(DMCreateMatrix(dm, &M));
1937: PetscCall(SNESComputeJacobian(snes, u, J, M));
1938: PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
1939: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
1940: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
1941: PetscCall(MatDestroy(&M));
1942: } else {
1943: PetscCall(SNESComputeJacobian(snes, u, J, J));
1944: }
1945: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
1946: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
1947: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
1948: /* Check nullspace */
1949: PetscCall(MatGetNullSpace(J, &nullspace));
1950: if (nullspace) {
1951: PetscBool isNull;
1952: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
1953: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1954: }
1955: /* Taylor test */
1956: {
1957: PetscRandom rand;
1958: Vec du, uhat, r, rhat, df;
1959: PetscReal h;
1960: PetscReal *es, *hs, *errors;
1961: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1962: PetscInt Nv, v;
1964: /* Choose a perturbation direction */
1965: PetscCall(PetscRandomCreate(comm, &rand));
1966: PetscCall(VecDuplicate(u, &du));
1967: PetscCall(VecSetRandom(du, rand));
1968: PetscCall(PetscRandomDestroy(&rand));
1969: PetscCall(VecDuplicate(u, &df));
1970: PetscCall(MatMult(J, du, df));
1971: /* Evaluate residual at u, F(u), save in vector r */
1972: PetscCall(VecDuplicate(u, &r));
1973: PetscCall(SNESComputeFunction(snes, u, r));
1974: /* Look at the convergence of our Taylor approximation as we approach u */
1975: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
1976: ;
1977: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
1978: PetscCall(VecDuplicate(u, &uhat));
1979: PetscCall(VecDuplicate(u, &rhat));
1980: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1981: PetscCall(VecWAXPY(uhat, h, du, u));
1982: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1983: PetscCall(SNESComputeFunction(snes, uhat, rhat));
1984: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
1985: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
1987: es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
1988: hs[Nv] = PetscLog10Real(h);
1989: }
1990: PetscCall(VecDestroy(&uhat));
1991: PetscCall(VecDestroy(&rhat));
1992: PetscCall(VecDestroy(&df));
1993: PetscCall(VecDestroy(&r));
1994: PetscCall(VecDestroy(&du));
1995: for (v = 0; v < Nv; ++v) {
1996: if ((tol >= 0) && (errors[v] > tol)) break;
1997: else if (errors[v] > PETSC_SMALL) break;
1998: }
1999: if (v == Nv) isLin = PETSC_TRUE;
2000: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
2001: PetscCall(PetscFree3(es, hs, errors));
2002: /* Slope should be about 2 */
2003: if (tol >= 0) {
2004: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
2005: } else if (isLinear || convRate) {
2006: if (isLinear) *isLinear = isLin;
2007: if (convRate) *convRate = slope;
2008: } else {
2009: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
2010: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
2011: }
2012: }
2013: PetscCall(MatDestroy(&J));
2014: PetscFunctionReturn(PETSC_SUCCESS);
2015: }
2017: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
2018: {
2019: PetscFunctionBegin;
2020: PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
2021: PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
2022: PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
2023: PetscFunctionReturn(PETSC_SUCCESS);
2024: }
2026: /*@C
2027: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
2029: Input Parameters:
2030: + snes - the `SNES` object
2031: - u - representative `SNES` vector
2033: Level: developer
2035: Note:
2036: The user must call `PetscDSSetExactSolution()` before this call
2038: .seealso: [](ch_snes), `SNES`, `DM`
2039: @*/
2040: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
2041: {
2042: DM dm;
2043: Vec sol;
2044: PetscBool check;
2046: PetscFunctionBegin;
2047: PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
2048: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
2049: PetscCall(SNESGetDM(snes, &dm));
2050: PetscCall(VecDuplicate(u, &sol));
2051: PetscCall(SNESSetSolution(snes, sol));
2052: PetscCall(DMSNESCheck_Internal(snes, dm, sol));
2053: PetscCall(VecDestroy(&sol));
2054: PetscFunctionReturn(PETSC_SUCCESS);
2055: }