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. This is normally used only to evaluate convergence rates for the pressure accurately.
15: Collective
17: Input Parameters:
18: + snes - The SNES
19: . pfield - The field number for pressure
20: . nullspace - The pressure nullspace
21: . u - The solution vector
22: - ctx - An optional user context
24: Output Parameter:
25: . u - The solution with a continuum pressure integral of zero
27: Notes:
28: 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.
30: Level: developer
32: .seealso: `SNESConvergedCorrectPressure()`
33: */
34: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
35: {
36: DM dm;
37: PetscDS ds;
38: const Vec *nullvecs;
39: PetscScalar pintd, *intc, *intn;
40: MPI_Comm comm;
41: PetscInt Nf, Nv;
43: PetscObjectGetComm((PetscObject)snes, &comm);
44: SNESGetDM(snes, &dm);
47: DMGetDS(dm, &ds);
48: PetscDSSetObjective(ds, pfield, pressure_Private);
49: MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs);
51: VecDot(nullvecs[0], u, &pintd);
53: PetscDSGetNumFields(ds, &Nf);
54: PetscMalloc2(Nf, &intc, Nf, &intn);
55: DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx);
56: DMPlexComputeIntegralFEM(dm, u, intc, ctx);
57: VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]);
58: #if defined(PETSC_USE_DEBUG)
59: DMPlexComputeIntegralFEM(dm, u, intc, ctx);
61: #endif
62: PetscFree2(intc, intn);
63: return 0;
64: }
66: /*@C
67: SNESConvergedCorrectPressure - Convergence test that adds a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
68: This is normally used only to evaluate convergence rates for the pressure accurately. The convergence test itself just mimics `SNESConvergedDefault()`.
70: Logically Collective
72: Input Parameters:
73: + snes - the `SNES` context
74: . it - the iteration (0 indicates before any Newton steps)
75: . xnorm - 2-norm of current iterate
76: . snorm - 2-norm of current step
77: . fnorm - 2-norm of function at current iterate
78: - ctx - Optional user context
80: Output Parameter:
81: . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`
83: Notes:
84: 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`
85: must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
86: The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
87: 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.
89: Options Database Key:
90: . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
92: Level: advanced
94: .seealso: `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`, `DMSetNullSpaceConstructor()`
95: @*/
96: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
97: {
98: PetscBool monitorIntegral = PETSC_FALSE;
100: SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx);
101: if (monitorIntegral) {
102: Mat J;
103: Vec u;
104: MatNullSpace nullspace;
105: const Vec *nullvecs;
106: PetscScalar pintd;
108: SNESGetSolution(snes, &u);
109: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
110: MatGetNullSpace(J, &nullspace);
111: MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs);
112: VecDot(nullvecs[0], u, &pintd);
113: PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd));
114: }
115: if (*reason > 0) {
116: Mat J;
117: Vec u;
118: MatNullSpace nullspace;
119: PetscInt pfield = 1;
121: SNESGetSolution(snes, &u);
122: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
123: MatGetNullSpace(J, &nullspace);
124: SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx);
125: }
126: return 0;
127: }
129: /************************** Interpolation *******************************/
131: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
132: {
133: PetscBool isPlex;
135: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
136: if (isPlex) {
137: *plex = dm;
138: PetscObjectReference((PetscObject)dm);
139: } else {
140: PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex);
141: if (!*plex) {
142: DMConvert(dm, DMPLEX, plex);
143: PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex);
144: if (copy) {
145: DMCopyDMSNES(dm, *plex);
146: DMCopyAuxiliaryVec(dm, *plex);
147: }
148: } else {
149: PetscObjectReference((PetscObject)*plex);
150: }
151: }
152: return 0;
153: }
155: /*@C
156: DMInterpolationCreate - Creates a `DMInterpolationInfo` context
158: Collective
160: Input Parameter:
161: . comm - the communicator
163: Output Parameter:
164: . ctx - the context
166: Level: beginner
168: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationDestroy()`
169: @*/
170: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
171: {
173: PetscNew(ctx);
175: (*ctx)->comm = comm;
176: (*ctx)->dim = -1;
177: (*ctx)->nInput = 0;
178: (*ctx)->points = NULL;
179: (*ctx)->cells = NULL;
180: (*ctx)->n = -1;
181: (*ctx)->coords = NULL;
182: return 0;
183: }
185: /*@C
186: DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
188: Not collective
190: Input Parameters:
191: + ctx - the context
192: - dim - the spatial dimension
194: Level: intermediate
196: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
197: @*/
198: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
199: {
201: ctx->dim = dim;
202: return 0;
203: }
205: /*@C
206: DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
208: Not collective
210: Input Parameter:
211: . ctx - the context
213: Output Parameter:
214: . dim - the spatial dimension
216: Level: intermediate
218: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
219: @*/
220: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
221: {
223: *dim = ctx->dim;
224: return 0;
225: }
227: /*@C
228: DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
230: Not collective
232: Input Parameters:
233: + ctx - the context
234: - dof - the number of fields
236: Level: intermediate
238: .seealso: `DMInterpolationInfo`, `DMInterpolationGetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
239: @*/
240: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
241: {
243: ctx->dof = dof;
244: return 0;
245: }
247: /*@C
248: DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
250: Not collective
252: Input Parameter:
253: . ctx - the context
255: Output Parameter:
256: . dof - the number of fields
258: Level: intermediate
260: .seealso: DMInterpolationInfo, `DMInterpolationSetDof()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`
261: @*/
262: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
263: {
265: *dof = ctx->dof;
266: return 0;
267: }
269: /*@C
270: DMInterpolationAddPoints - Add points at which we will interpolate the fields
272: Not collective
274: Input Parameters:
275: + ctx - the context
276: . n - the number of points
277: - points - the coordinates for each point, an array of size n * dim
279: Note:
280: The coordinate information is copied.
282: Level: intermediate
284: .seealso: `DMInterpolationInfo`, `DMInterpolationSetDim()`, `DMInterpolationEvaluate()`, `DMInterpolationCreate()`
285: @*/
286: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
287: {
290: ctx->nInput = n;
292: PetscMalloc1(n * ctx->dim, &ctx->points);
293: PetscArraycpy(ctx->points, points, n * ctx->dim);
294: return 0;
295: }
297: /*@C
298: DMInterpolationSetUp - Compute spatial indices for point location during interpolation
300: Collective
302: Input Parameters:
303: + ctx - the context
304: . dm - the `DM` for the function space used for interpolation
305: . redundantPoints - If `PETSC_TRUE`, all processes are passing in the same array of points. Otherwise, points need to be communicated among processes.
306: - ignoreOutsideDomain - If `PETSC_TRUE`, ignore points outside the domain, otherwise return an error
308: Level: intermediate
310: .seealso: DMInterpolationInfo, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
311: @*/
312: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints, PetscBool ignoreOutsideDomain)
313: {
314: MPI_Comm comm = ctx->comm;
315: PetscScalar *a;
316: PetscInt p, q, i;
317: PetscMPIInt rank, size;
318: Vec pointVec;
319: PetscSF cellSF;
320: PetscLayout layout;
321: PetscReal *globalPoints;
322: PetscScalar *globalPointsScalar;
323: const PetscInt *ranges;
324: PetscMPIInt *counts, *displs;
325: const PetscSFNode *foundCells;
326: const PetscInt *foundPoints;
327: PetscMPIInt *foundProcs, *globalProcs;
328: PetscInt n, N, numFound;
331: MPI_Comm_size(comm, &size);
332: MPI_Comm_rank(comm, &rank);
334: /* Locate points */
335: n = ctx->nInput;
336: if (!redundantPoints) {
337: PetscLayoutCreate(comm, &layout);
338: PetscLayoutSetBlockSize(layout, 1);
339: PetscLayoutSetLocalSize(layout, n);
340: PetscLayoutSetUp(layout);
341: PetscLayoutGetSize(layout, &N);
342: /* Communicate all points to all processes */
343: PetscMalloc3(N * ctx->dim, &globalPoints, size, &counts, size, &displs);
344: PetscLayoutGetRanges(layout, &ranges);
345: for (p = 0; p < size; ++p) {
346: counts[p] = (ranges[p + 1] - ranges[p]) * ctx->dim;
347: displs[p] = ranges[p] * ctx->dim;
348: }
349: MPI_Allgatherv(ctx->points, n * ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
350: } else {
351: N = n;
352: globalPoints = ctx->points;
353: counts = displs = NULL;
354: layout = NULL;
355: }
356: #if 0
357: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
358: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
359: #else
360: #if defined(PETSC_USE_COMPLEX)
361: PetscMalloc1(N * ctx->dim, &globalPointsScalar);
362: for (i = 0; i < N * ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
363: #else
364: globalPointsScalar = globalPoints;
365: #endif
366: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N * ctx->dim, globalPointsScalar, &pointVec);
367: PetscMalloc2(N, &foundProcs, N, &globalProcs);
368: for (p = 0; p < N; ++p) foundProcs[p] = size;
369: cellSF = NULL;
370: DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
371: PetscSFGetGraph(cellSF, NULL, &numFound, &foundPoints, &foundCells);
372: #endif
373: for (p = 0; p < numFound; ++p) {
374: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
375: }
376: /* Let the lowest rank process own each point */
377: MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
378: ctx->n = 0;
379: for (p = 0; p < N; ++p) {
380: if (globalProcs[p] == size) {
382: (double)(ctx->dim > 2 ? globalPoints[p * ctx->dim + 2] : 0.0));
383: if (rank == 0) ++ctx->n;
384: } else if (globalProcs[p] == rank) ++ctx->n;
385: }
386: /* Create coordinates vector and array of owned cells */
387: PetscMalloc1(ctx->n, &ctx->cells);
388: VecCreate(comm, &ctx->coords);
389: VecSetSizes(ctx->coords, ctx->n * ctx->dim, PETSC_DECIDE);
390: VecSetBlockSize(ctx->coords, ctx->dim);
391: VecSetType(ctx->coords, VECSTANDARD);
392: VecGetArray(ctx->coords, &a);
393: for (p = 0, q = 0, i = 0; p < N; ++p) {
394: if (globalProcs[p] == rank) {
395: PetscInt d;
397: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p * ctx->dim + d];
398: ctx->cells[q] = foundCells[q].index;
399: ++q;
400: }
401: if (globalProcs[p] == size && rank == 0) {
402: PetscInt d;
404: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = 0.;
405: ctx->cells[q] = -1;
406: ++q;
407: }
408: }
409: VecRestoreArray(ctx->coords, &a);
410: #if 0
411: PetscFree3(foundCells,foundProcs,globalProcs);
412: #else
413: PetscFree2(foundProcs, globalProcs);
414: PetscSFDestroy(&cellSF);
415: VecDestroy(&pointVec);
416: #endif
417: if ((void *)globalPointsScalar != (void *)globalPoints) PetscFree(globalPointsScalar);
418: if (!redundantPoints) PetscFree3(globalPoints, counts, displs);
419: PetscLayoutDestroy(&layout);
420: return 0;
421: }
423: /*@C
424: DMInterpolationGetCoordinates - Gets a `Vec` with the coordinates of each interpolation point
426: Collective
428: Input Parameter:
429: . ctx - the context
431: Output Parameter:
432: . coordinates - the coordinates of interpolation points
434: Note:
435: The local vector entries correspond to interpolation points lying on this process, according to the associated `DM`.
436: This is a borrowed vector that the user should not destroy.
438: Level: intermediate
440: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
441: @*/
442: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
443: {
446: *coordinates = ctx->coords;
447: return 0;
448: }
450: /*@C
451: DMInterpolationGetVector - Gets a `Vec` which can hold all the interpolated field values
453: Collective
455: Input Parameter:
456: . ctx - the context
458: Output Parameter:
459: . v - a vector capable of holding the interpolated field values
461: Note:
462: This vector should be returned using `DMInterpolationRestoreVector()`.
464: Level: intermediate
466: .seealso: `DMInterpolationInfo`, `DMInterpolationRestoreVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
467: @*/
468: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
469: {
472: VecCreate(ctx->comm, v);
473: VecSetSizes(*v, ctx->n * ctx->dof, PETSC_DECIDE);
474: VecSetBlockSize(*v, ctx->dof);
475: VecSetType(*v, VECSTANDARD);
476: return 0;
477: }
479: /*@C
480: DMInterpolationRestoreVector - Returns a `Vec` which can hold all the interpolated field values
482: Collective
484: Input Parameters:
485: + ctx - the context
486: - v - a vector capable of holding the interpolated field values
488: Level: intermediate
490: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
491: @*/
492: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
493: {
496: VecDestroy(v);
497: return 0;
498: }
500: static inline PetscErrorCode DMInterpolate_Segment_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
501: {
502: PetscReal v0, J, invJ, detJ;
503: const PetscInt dof = ctx->dof;
504: const PetscScalar *coords;
505: PetscScalar *a;
506: PetscInt p;
508: VecGetArrayRead(ctx->coords, &coords);
509: VecGetArray(v, &a);
510: for (p = 0; p < ctx->n; ++p) {
511: PetscInt c = ctx->cells[p];
512: PetscScalar *x = NULL;
513: PetscReal xir[1];
514: PetscInt xSize, comp;
516: DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ);
518: xir[0] = invJ * PetscRealPart(coords[p] - v0);
519: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
520: if (2 * dof == xSize) {
521: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp] * (1 - xir[0]) + x[1 * dof + comp] * xir[0];
522: } else if (dof == xSize) {
523: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
524: } 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);
525: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
526: }
527: VecRestoreArray(v, &a);
528: VecRestoreArrayRead(ctx->coords, &coords);
529: return 0;
530: }
532: static inline PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
533: {
534: PetscReal *v0, *J, *invJ, detJ;
535: const PetscScalar *coords;
536: PetscScalar *a;
537: PetscInt p;
539: PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ);
540: VecGetArrayRead(ctx->coords, &coords);
541: VecGetArray(v, &a);
542: for (p = 0; p < ctx->n; ++p) {
543: PetscInt c = ctx->cells[p];
544: PetscScalar *x = NULL;
545: PetscReal xi[4];
546: PetscInt d, f, comp;
548: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
550: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
551: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
553: for (d = 0; d < ctx->dim; ++d) {
554: xi[d] = 0.0;
555: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
556: 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];
557: }
558: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
559: }
560: VecRestoreArray(v, &a);
561: VecRestoreArrayRead(ctx->coords, &coords);
562: PetscFree3(v0, J, invJ);
563: return 0;
564: }
566: static inline PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
567: {
568: PetscReal *v0, *J, *invJ, detJ;
569: const PetscScalar *coords;
570: PetscScalar *a;
571: PetscInt p;
573: PetscMalloc3(ctx->dim, &v0, ctx->dim * ctx->dim, &J, ctx->dim * ctx->dim, &invJ);
574: VecGetArrayRead(ctx->coords, &coords);
575: VecGetArray(v, &a);
576: for (p = 0; p < ctx->n; ++p) {
577: PetscInt c = ctx->cells[p];
578: const PetscInt order[3] = {2, 1, 3};
579: PetscScalar *x = NULL;
580: PetscReal xi[4];
581: PetscInt d, f, comp;
583: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
585: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
586: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
588: for (d = 0; d < ctx->dim; ++d) {
589: xi[d] = 0.0;
590: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d * ctx->dim + f] * 0.5 * PetscRealPart(coords[p * ctx->dim + f] - v0[f]);
591: 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];
592: }
593: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
594: }
595: VecRestoreArray(v, &a);
596: VecRestoreArrayRead(ctx->coords, &coords);
597: PetscFree3(v0, J, invJ);
598: return 0;
599: }
601: static inline PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
602: {
603: const PetscScalar *vertices = (const PetscScalar *)ctx;
604: const PetscScalar x0 = vertices[0];
605: const PetscScalar y0 = vertices[1];
606: const PetscScalar x1 = vertices[2];
607: const PetscScalar y1 = vertices[3];
608: const PetscScalar x2 = vertices[4];
609: const PetscScalar y2 = vertices[5];
610: const PetscScalar x3 = vertices[6];
611: const PetscScalar y3 = vertices[7];
612: const PetscScalar f_1 = x1 - x0;
613: const PetscScalar g_1 = y1 - y0;
614: const PetscScalar f_3 = x3 - x0;
615: const PetscScalar g_3 = y3 - y0;
616: const PetscScalar f_01 = x2 - x1 - x3 + x0;
617: const PetscScalar g_01 = y2 - y1 - y3 + y0;
618: const PetscScalar *ref;
619: PetscScalar *real;
621: VecGetArrayRead(Xref, &ref);
622: VecGetArray(Xreal, &real);
623: {
624: const PetscScalar p0 = ref[0];
625: const PetscScalar p1 = ref[1];
627: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
628: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
629: }
630: PetscLogFlops(28);
631: VecRestoreArrayRead(Xref, &ref);
632: VecRestoreArray(Xreal, &real);
633: return 0;
634: }
636: #include <petsc/private/dmimpl.h>
637: static inline PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
638: {
639: const PetscScalar *vertices = (const PetscScalar *)ctx;
640: const PetscScalar x0 = vertices[0];
641: const PetscScalar y0 = vertices[1];
642: const PetscScalar x1 = vertices[2];
643: const PetscScalar y1 = vertices[3];
644: const PetscScalar x2 = vertices[4];
645: const PetscScalar y2 = vertices[5];
646: const PetscScalar x3 = vertices[6];
647: const PetscScalar y3 = vertices[7];
648: const PetscScalar f_01 = x2 - x1 - x3 + x0;
649: const PetscScalar g_01 = y2 - y1 - y3 + y0;
650: const PetscScalar *ref;
652: VecGetArrayRead(Xref, &ref);
653: {
654: const PetscScalar x = ref[0];
655: const PetscScalar y = ref[1];
656: const PetscInt rows[2] = {0, 1};
657: PetscScalar values[4];
659: values[0] = (x1 - x0 + f_01 * y) * 0.5;
660: values[1] = (x3 - x0 + f_01 * x) * 0.5;
661: values[2] = (y1 - y0 + g_01 * y) * 0.5;
662: values[3] = (y3 - y0 + g_01 * x) * 0.5;
663: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
664: }
665: PetscLogFlops(30);
666: VecRestoreArrayRead(Xref, &ref);
667: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
668: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
669: return 0;
670: }
672: static inline PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
673: {
674: DM dmCoord;
675: PetscFE fem = NULL;
676: SNES snes;
677: KSP ksp;
678: PC pc;
679: Vec coordsLocal, r, ref, real;
680: Mat J;
681: PetscTabulation T = NULL;
682: const PetscScalar *coords;
683: PetscScalar *a;
684: PetscReal xir[2] = {0., 0.};
685: PetscInt Nf, p;
686: const PetscInt dof = ctx->dof;
688: DMGetNumFields(dm, &Nf);
689: if (Nf) {
690: PetscObject obj;
691: PetscClassId id;
693: DMGetField(dm, 0, NULL, &obj);
694: PetscObjectGetClassId(obj, &id);
695: if (id == PETSCFE_CLASSID) {
696: fem = (PetscFE)obj;
697: PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);
698: }
699: }
700: DMGetCoordinatesLocal(dm, &coordsLocal);
701: DMGetCoordinateDM(dm, &dmCoord);
702: SNESCreate(PETSC_COMM_SELF, &snes);
703: SNESSetOptionsPrefix(snes, "quad_interp_");
704: VecCreate(PETSC_COMM_SELF, &r);
705: VecSetSizes(r, 2, 2);
706: VecSetType(r, dm->vectype);
707: VecDuplicate(r, &ref);
708: VecDuplicate(r, &real);
709: MatCreate(PETSC_COMM_SELF, &J);
710: MatSetSizes(J, 2, 2, 2, 2);
711: MatSetType(J, MATSEQDENSE);
712: MatSetUp(J);
713: SNESSetFunction(snes, r, QuadMap_Private, NULL);
714: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
715: SNESGetKSP(snes, &ksp);
716: KSPGetPC(ksp, &pc);
717: PCSetType(pc, PCLU);
718: SNESSetFromOptions(snes);
720: VecGetArrayRead(ctx->coords, &coords);
721: VecGetArray(v, &a);
722: for (p = 0; p < ctx->n; ++p) {
723: PetscScalar *x = NULL, *vertices = NULL;
724: PetscScalar *xi;
725: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
727: /* Can make this do all points at once */
728: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
730: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
731: SNESSetFunction(snes, NULL, NULL, vertices);
732: SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
733: VecGetArray(real, &xi);
734: xi[0] = coords[p * ctx->dim + 0];
735: xi[1] = coords[p * ctx->dim + 1];
736: VecRestoreArray(real, &xi);
737: SNESSolve(snes, real, ref);
738: VecGetArray(ref, &xi);
739: xir[0] = PetscRealPart(xi[0]);
740: xir[1] = PetscRealPart(xi[1]);
741: if (4 * dof == xSize) {
742: 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];
743: } else if (dof == xSize) {
744: for (comp = 0; comp < dof; ++comp) a[p * dof + comp] = x[0 * dof + comp];
745: } else {
746: PetscInt d;
749: xir[0] = 2.0 * xir[0] - 1.0;
750: xir[1] = 2.0 * xir[1] - 1.0;
751: PetscFEComputeTabulation(fem, 1, xir, 0, T);
752: for (comp = 0; comp < dof; ++comp) {
753: a[p * dof + comp] = 0.0;
754: for (d = 0; d < xSize / dof; ++d) a[p * dof + comp] += x[d * dof + comp] * T->T[0][d * dof + comp];
755: }
756: }
757: VecRestoreArray(ref, &xi);
758: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
759: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
760: }
761: PetscTabulationDestroy(&T);
762: VecRestoreArray(v, &a);
763: VecRestoreArrayRead(ctx->coords, &coords);
765: SNESDestroy(&snes);
766: VecDestroy(&r);
767: VecDestroy(&ref);
768: VecDestroy(&real);
769: MatDestroy(&J);
770: return 0;
771: }
773: static inline PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
774: {
775: const PetscScalar *vertices = (const PetscScalar *)ctx;
776: const PetscScalar x0 = vertices[0];
777: const PetscScalar y0 = vertices[1];
778: const PetscScalar z0 = vertices[2];
779: const PetscScalar x1 = vertices[9];
780: const PetscScalar y1 = vertices[10];
781: const PetscScalar z1 = vertices[11];
782: const PetscScalar x2 = vertices[6];
783: const PetscScalar y2 = vertices[7];
784: const PetscScalar z2 = vertices[8];
785: const PetscScalar x3 = vertices[3];
786: const PetscScalar y3 = vertices[4];
787: const PetscScalar z3 = vertices[5];
788: const PetscScalar x4 = vertices[12];
789: const PetscScalar y4 = vertices[13];
790: const PetscScalar z4 = vertices[14];
791: const PetscScalar x5 = vertices[15];
792: const PetscScalar y5 = vertices[16];
793: const PetscScalar z5 = vertices[17];
794: const PetscScalar x6 = vertices[18];
795: const PetscScalar y6 = vertices[19];
796: const PetscScalar z6 = vertices[20];
797: const PetscScalar x7 = vertices[21];
798: const PetscScalar y7 = vertices[22];
799: const PetscScalar z7 = vertices[23];
800: const PetscScalar f_1 = x1 - x0;
801: const PetscScalar g_1 = y1 - y0;
802: const PetscScalar h_1 = z1 - z0;
803: const PetscScalar f_3 = x3 - x0;
804: const PetscScalar g_3 = y3 - y0;
805: const PetscScalar h_3 = z3 - z0;
806: const PetscScalar f_4 = x4 - x0;
807: const PetscScalar g_4 = y4 - y0;
808: const PetscScalar h_4 = z4 - z0;
809: const PetscScalar f_01 = x2 - x1 - x3 + x0;
810: const PetscScalar g_01 = y2 - y1 - y3 + y0;
811: const PetscScalar h_01 = z2 - z1 - z3 + z0;
812: const PetscScalar f_12 = x7 - x3 - x4 + x0;
813: const PetscScalar g_12 = y7 - y3 - y4 + y0;
814: const PetscScalar h_12 = z7 - z3 - z4 + z0;
815: const PetscScalar f_02 = x5 - x1 - x4 + x0;
816: const PetscScalar g_02 = y5 - y1 - y4 + y0;
817: const PetscScalar h_02 = z5 - z1 - z4 + z0;
818: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
819: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
820: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
821: const PetscScalar *ref;
822: PetscScalar *real;
824: VecGetArrayRead(Xref, &ref);
825: VecGetArray(Xreal, &real);
826: {
827: const PetscScalar p0 = ref[0];
828: const PetscScalar p1 = ref[1];
829: const PetscScalar p2 = ref[2];
831: 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;
832: 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;
833: 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;
834: }
835: PetscLogFlops(114);
836: VecRestoreArrayRead(Xref, &ref);
837: VecRestoreArray(Xreal, &real);
838: return 0;
839: }
841: static inline PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
842: {
843: const PetscScalar *vertices = (const PetscScalar *)ctx;
844: const PetscScalar x0 = vertices[0];
845: const PetscScalar y0 = vertices[1];
846: const PetscScalar z0 = vertices[2];
847: const PetscScalar x1 = vertices[9];
848: const PetscScalar y1 = vertices[10];
849: const PetscScalar z1 = vertices[11];
850: const PetscScalar x2 = vertices[6];
851: const PetscScalar y2 = vertices[7];
852: const PetscScalar z2 = vertices[8];
853: const PetscScalar x3 = vertices[3];
854: const PetscScalar y3 = vertices[4];
855: const PetscScalar z3 = vertices[5];
856: const PetscScalar x4 = vertices[12];
857: const PetscScalar y4 = vertices[13];
858: const PetscScalar z4 = vertices[14];
859: const PetscScalar x5 = vertices[15];
860: const PetscScalar y5 = vertices[16];
861: const PetscScalar z5 = vertices[17];
862: const PetscScalar x6 = vertices[18];
863: const PetscScalar y6 = vertices[19];
864: const PetscScalar z6 = vertices[20];
865: const PetscScalar x7 = vertices[21];
866: const PetscScalar y7 = vertices[22];
867: const PetscScalar z7 = vertices[23];
868: const PetscScalar f_xy = x2 - x1 - x3 + x0;
869: const PetscScalar g_xy = y2 - y1 - y3 + y0;
870: const PetscScalar h_xy = z2 - z1 - z3 + z0;
871: const PetscScalar f_yz = x7 - x3 - x4 + x0;
872: const PetscScalar g_yz = y7 - y3 - y4 + y0;
873: const PetscScalar h_yz = z7 - z3 - z4 + z0;
874: const PetscScalar f_xz = x5 - x1 - x4 + x0;
875: const PetscScalar g_xz = y5 - y1 - y4 + y0;
876: const PetscScalar h_xz = z5 - z1 - z4 + z0;
877: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
878: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
879: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
880: const PetscScalar *ref;
882: VecGetArrayRead(Xref, &ref);
883: {
884: const PetscScalar x = ref[0];
885: const PetscScalar y = ref[1];
886: const PetscScalar z = ref[2];
887: const PetscInt rows[3] = {0, 1, 2};
888: PetscScalar values[9];
890: values[0] = (x1 - x0 + f_xy * y + f_xz * z + f_xyz * y * z) / 2.0;
891: values[1] = (x3 - x0 + f_xy * x + f_yz * z + f_xyz * x * z) / 2.0;
892: values[2] = (x4 - x0 + f_yz * y + f_xz * x + f_xyz * x * y) / 2.0;
893: values[3] = (y1 - y0 + g_xy * y + g_xz * z + g_xyz * y * z) / 2.0;
894: values[4] = (y3 - y0 + g_xy * x + g_yz * z + g_xyz * x * z) / 2.0;
895: values[5] = (y4 - y0 + g_yz * y + g_xz * x + g_xyz * x * y) / 2.0;
896: values[6] = (z1 - z0 + h_xy * y + h_xz * z + h_xyz * y * z) / 2.0;
897: values[7] = (z3 - z0 + h_xy * x + h_yz * z + h_xyz * x * z) / 2.0;
898: values[8] = (z4 - z0 + h_yz * y + h_xz * x + h_xyz * x * y) / 2.0;
900: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
901: }
902: PetscLogFlops(152);
903: VecRestoreArrayRead(Xref, &ref);
904: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
905: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
906: return 0;
907: }
909: static inline PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
910: {
911: DM dmCoord;
912: SNES snes;
913: KSP ksp;
914: PC pc;
915: Vec coordsLocal, r, ref, real;
916: Mat J;
917: const PetscScalar *coords;
918: PetscScalar *a;
919: PetscInt p;
921: DMGetCoordinatesLocal(dm, &coordsLocal);
922: DMGetCoordinateDM(dm, &dmCoord);
923: SNESCreate(PETSC_COMM_SELF, &snes);
924: SNESSetOptionsPrefix(snes, "hex_interp_");
925: VecCreate(PETSC_COMM_SELF, &r);
926: VecSetSizes(r, 3, 3);
927: VecSetType(r, dm->vectype);
928: VecDuplicate(r, &ref);
929: VecDuplicate(r, &real);
930: MatCreate(PETSC_COMM_SELF, &J);
931: MatSetSizes(J, 3, 3, 3, 3);
932: MatSetType(J, MATSEQDENSE);
933: MatSetUp(J);
934: SNESSetFunction(snes, r, HexMap_Private, NULL);
935: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
936: SNESGetKSP(snes, &ksp);
937: KSPGetPC(ksp, &pc);
938: PCSetType(pc, PCLU);
939: SNESSetFromOptions(snes);
941: VecGetArrayRead(ctx->coords, &coords);
942: VecGetArray(v, &a);
943: for (p = 0; p < ctx->n; ++p) {
944: PetscScalar *x = NULL, *vertices = NULL;
945: PetscScalar *xi;
946: PetscReal xir[3];
947: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
949: /* Can make this do all points at once */
950: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
952: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
954: SNESSetFunction(snes, NULL, NULL, vertices);
955: SNESSetJacobian(snes, NULL, NULL, NULL, vertices);
956: VecGetArray(real, &xi);
957: xi[0] = coords[p * ctx->dim + 0];
958: xi[1] = coords[p * ctx->dim + 1];
959: xi[2] = coords[p * ctx->dim + 2];
960: VecRestoreArray(real, &xi);
961: SNESSolve(snes, real, ref);
962: VecGetArray(ref, &xi);
963: xir[0] = PetscRealPart(xi[0]);
964: xir[1] = PetscRealPart(xi[1]);
965: xir[2] = PetscRealPart(xi[2]);
966: if (8 * ctx->dof == xSize) {
967: for (comp = 0; comp < ctx->dof; ++comp) {
968: 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]) +
969: 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];
970: }
971: } else {
972: for (comp = 0; comp < ctx->dof; ++comp) a[p * ctx->dof + comp] = x[0 * ctx->dof + comp];
973: }
974: VecRestoreArray(ref, &xi);
975: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
976: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
977: }
978: VecRestoreArray(v, &a);
979: VecRestoreArrayRead(ctx->coords, &coords);
981: SNESDestroy(&snes);
982: VecDestroy(&r);
983: VecDestroy(&ref);
984: VecDestroy(&real);
985: MatDestroy(&J);
986: return 0;
987: }
989: /*@C
990: DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.
992: Input Parameters:
993: + ctx - The `DMInterpolationInfo` context
994: . dm - The `DM`
995: - x - The local vector containing the field to be interpolated
997: Output Parameters:
998: . v - The vector containing the interpolated values
1000: Note:
1001: A suitable v can be obtained using `DMInterpolationGetVector()`.
1003: Level: beginner
1005: .seealso: `DMInterpolationInfo`, `DMInterpolationGetVector()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1006: @*/
1007: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1008: {
1009: PetscDS ds;
1010: PetscInt n, p, Nf, field;
1011: PetscBool useDS = PETSC_FALSE;
1016: VecGetLocalSize(v, &n);
1018: if (!n) return 0;
1019: DMGetDS(dm, &ds);
1020: if (ds) {
1021: useDS = PETSC_TRUE;
1022: PetscDSGetNumFields(ds, &Nf);
1023: for (field = 0; field < Nf; ++field) {
1024: PetscObject obj;
1025: PetscClassId id;
1027: PetscDSGetDiscretization(ds, field, &obj);
1028: PetscObjectGetClassId(obj, &id);
1029: if (id != PETSCFE_CLASSID) {
1030: useDS = PETSC_FALSE;
1031: break;
1032: }
1033: }
1034: }
1035: if (useDS) {
1036: const PetscScalar *coords;
1037: PetscScalar *interpolant;
1038: PetscInt cdim, d;
1040: DMGetCoordinateDim(dm, &cdim);
1041: VecGetArrayRead(ctx->coords, &coords);
1042: VecGetArrayWrite(v, &interpolant);
1043: for (p = 0; p < ctx->n; ++p) {
1044: PetscReal pcoords[3], xi[3];
1045: PetscScalar *xa = NULL;
1046: PetscInt coff = 0, foff = 0, clSize;
1048: if (ctx->cells[p] < 0) continue;
1049: for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p * cdim + d]);
1050: DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi);
1051: DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1052: for (field = 0; field < Nf; ++field) {
1053: PetscTabulation T;
1054: PetscFE fe;
1056: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
1057: PetscFECreateTabulation(fe, 1, 1, xi, 0, &T);
1058: {
1059: const PetscReal *basis = T->T[0];
1060: const PetscInt Nb = T->Nb;
1061: const PetscInt Nc = T->Nc;
1062: PetscInt f, fc;
1064: for (fc = 0; fc < Nc; ++fc) {
1065: interpolant[p * ctx->dof + coff + fc] = 0.0;
1066: for (f = 0; f < Nb; ++f) interpolant[p * ctx->dof + coff + fc] += xa[foff + f] * basis[(0 * Nb + f) * Nc + fc];
1067: }
1068: coff += Nc;
1069: foff += Nb;
1070: }
1071: PetscTabulationDestroy(&T);
1072: }
1073: DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], &clSize, &xa);
1076: }
1077: VecRestoreArrayRead(ctx->coords, &coords);
1078: VecRestoreArrayWrite(v, &interpolant);
1079: } else {
1080: DMPolytopeType ct;
1082: /* TODO Check each cell individually */
1083: DMPlexGetCellType(dm, ctx->cells[0], &ct);
1084: switch (ct) {
1085: case DM_POLYTOPE_SEGMENT:
1086: DMInterpolate_Segment_Private(ctx, dm, x, v);
1087: break;
1088: case DM_POLYTOPE_TRIANGLE:
1089: DMInterpolate_Triangle_Private(ctx, dm, x, v);
1090: break;
1091: case DM_POLYTOPE_QUADRILATERAL:
1092: DMInterpolate_Quad_Private(ctx, dm, x, v);
1093: break;
1094: case DM_POLYTOPE_TETRAHEDRON:
1095: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
1096: break;
1097: case DM_POLYTOPE_HEXAHEDRON:
1098: DMInterpolate_Hex_Private(ctx, dm, x, v);
1099: break;
1100: default:
1101: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for cell type %s", DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
1102: }
1103: }
1104: return 0;
1105: }
1107: /*@C
1108: DMInterpolationDestroy - Destroys a `DMInterpolationInfo` context
1110: Collective
1112: Input Parameter:
1113: . ctx - the context
1115: Level: beginner
1117: .seealso: `DMInterpolationInfo`, `DMInterpolationEvaluate()`, `DMInterpolationAddPoints()`, `DMInterpolationCreate()`
1118: @*/
1119: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1120: {
1122: VecDestroy(&(*ctx)->coords);
1123: PetscFree((*ctx)->points);
1124: PetscFree((*ctx)->cells);
1125: PetscFree(*ctx);
1126: *ctx = NULL;
1127: return 0;
1128: }
1130: /*@C
1131: SNESMonitorFields - Monitors the residual for each field separately
1133: Collective
1135: Input Parameters:
1136: + snes - the `SNES` context
1137: . its - iteration number
1138: . fgnorm - 2-norm of residual
1139: - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`
1141: Note:
1142: This routine prints the residual norm at each iteration.
1144: Level: intermediate
1146: .seealso: `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
1147: @*/
1148: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1149: {
1150: PetscViewer viewer = vf->viewer;
1151: Vec res;
1152: DM dm;
1153: PetscSection s;
1154: const PetscScalar *r;
1155: PetscReal *lnorms, *norms;
1156: PetscInt numFields, f, pStart, pEnd, p;
1159: SNESGetFunction(snes, &res, NULL, NULL);
1160: SNESGetDM(snes, &dm);
1161: DMGetLocalSection(dm, &s);
1162: PetscSectionGetNumFields(s, &numFields);
1163: PetscSectionGetChart(s, &pStart, &pEnd);
1164: PetscCalloc2(numFields, &lnorms, numFields, &norms);
1165: VecGetArrayRead(res, &r);
1166: for (p = pStart; p < pEnd; ++p) {
1167: for (f = 0; f < numFields; ++f) {
1168: PetscInt fdof, foff, d;
1170: PetscSectionGetFieldDof(s, p, f, &fdof);
1171: PetscSectionGetFieldOffset(s, p, f, &foff);
1172: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
1173: }
1174: }
1175: VecRestoreArrayRead(res, &r);
1176: MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));
1177: PetscViewerPushFormat(viewer, vf->format);
1178: PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
1179: PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm);
1180: for (f = 0; f < numFields; ++f) {
1181: if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
1182: PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f]));
1183: }
1184: PetscViewerASCIIPrintf(viewer, "]\n");
1185: PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
1186: PetscViewerPopFormat(viewer);
1187: PetscFree2(lnorms, norms);
1188: return 0;
1189: }
1191: /********************* Residual Computation **************************/
1193: PetscErrorCode DMPlexGetAllCells_Internal(DM plex, IS *cellIS)
1194: {
1195: PetscInt depth;
1197: DMPlexGetDepth(plex, &depth);
1198: DMGetStratumIS(plex, "dim", depth, cellIS);
1199: if (!*cellIS) DMGetStratumIS(plex, "depth", depth, cellIS);
1200: return 0;
1201: }
1203: /*@
1204: DMPlexSNESComputeResidualFEM - Sums the local residual into vector F from the local input X using pointwise functions specified by the user
1206: Input Parameters:
1207: + dm - The mesh
1208: . X - Local solution
1209: - user - The user context
1211: Output Parameter:
1212: . F - Local output vector
1214: Note:
1215: The residual is summed into F; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in F is intentional.
1217: Level: developer
1219: .seealso: `DM`, `DMPlexComputeJacobianAction()`
1220: @*/
1221: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1222: {
1223: DM plex;
1224: IS allcellIS;
1225: PetscInt Nds, s;
1227: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1228: DMPlexGetAllCells_Internal(plex, &allcellIS);
1229: DMGetNumDS(dm, &Nds);
1230: for (s = 0; s < Nds; ++s) {
1231: PetscDS ds;
1232: IS cellIS;
1233: PetscFormKey key;
1235: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1236: key.value = 0;
1237: key.field = 0;
1238: key.part = 0;
1239: if (!key.label) {
1240: PetscObjectReference((PetscObject)allcellIS);
1241: cellIS = allcellIS;
1242: } else {
1243: IS pointIS;
1245: key.value = 1;
1246: DMLabelGetStratumIS(key.label, key.value, &pointIS);
1247: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1248: ISDestroy(&pointIS);
1249: }
1250: DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1251: ISDestroy(&cellIS);
1252: }
1253: ISDestroy(&allcellIS);
1254: DMDestroy(&plex);
1255: return 0;
1256: }
1258: PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user)
1259: {
1260: DM plex;
1261: IS allcellIS;
1262: PetscInt Nds, s;
1264: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1265: DMPlexGetAllCells_Internal(plex, &allcellIS);
1266: DMGetNumDS(dm, &Nds);
1267: for (s = 0; s < Nds; ++s) {
1268: PetscDS ds;
1269: DMLabel label;
1270: IS cellIS;
1272: DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1273: {
1274: PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
1275: PetscWeakForm wf;
1276: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
1277: PetscFormKey *reskeys;
1279: /* Get unique residual keys */
1280: for (m = 0; m < Nm; ++m) {
1281: PetscInt Nkm;
1282: PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm);
1283: Nk += Nkm;
1284: }
1285: PetscMalloc1(Nk, &reskeys);
1286: for (m = 0; m < Nm; ++m) PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys);
1288: PetscFormKeySort(Nk, reskeys);
1289: for (k = 0, kp = 1; kp < Nk; ++kp) {
1290: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1291: ++k;
1292: if (kp != k) reskeys[k] = reskeys[kp];
1293: }
1294: }
1295: Nk = k;
1297: PetscDSGetWeakForm(ds, &wf);
1298: for (k = 0; k < Nk; ++k) {
1299: DMLabel label = reskeys[k].label;
1300: PetscInt val = reskeys[k].value;
1302: if (!label) {
1303: PetscObjectReference((PetscObject)allcellIS);
1304: cellIS = allcellIS;
1305: } else {
1306: IS pointIS;
1308: DMLabelGetStratumIS(label, val, &pointIS);
1309: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1310: ISDestroy(&pointIS);
1311: }
1312: DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1313: ISDestroy(&cellIS);
1314: }
1315: PetscFree(reskeys);
1316: }
1317: }
1318: ISDestroy(&allcellIS);
1319: DMDestroy(&plex);
1320: return 0;
1321: }
1323: /*@
1324: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
1326: Input Parameters:
1327: + dm - The mesh
1328: - user - The user context
1330: Output Parameter:
1331: . X - Local solution
1333: Level: developer
1335: .seealso: `DMPLEX`, `DMPlexComputeJacobianAction()`
1336: @*/
1337: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1338: {
1339: DM plex;
1341: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1342: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1343: DMDestroy(&plex);
1344: return 0;
1345: }
1347: /*@
1348: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(X) on Y
1350: Input Parameters:
1351: + dm - The `DM`
1352: . X - Local solution vector
1353: . Y - Local input vector
1354: - user - The user context
1356: Output Parameter:
1357: . F - lcoal output vector
1359: Level: developer
1361: Notes:
1362: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
1364: .seealso: `DM`, ``DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
1365: @*/
1366: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
1367: {
1368: DM plex;
1369: IS allcellIS;
1370: PetscInt Nds, s;
1372: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1373: DMPlexGetAllCells_Internal(plex, &allcellIS);
1374: DMGetNumDS(dm, &Nds);
1375: for (s = 0; s < Nds; ++s) {
1376: PetscDS ds;
1377: DMLabel label;
1378: IS cellIS;
1380: DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1381: {
1382: PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
1383: PetscWeakForm wf;
1384: PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
1385: PetscFormKey *jackeys;
1387: /* Get unique Jacobian keys */
1388: for (m = 0; m < Nm; ++m) {
1389: PetscInt Nkm;
1390: PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm);
1391: Nk += Nkm;
1392: }
1393: PetscMalloc1(Nk, &jackeys);
1394: for (m = 0; m < Nm; ++m) PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys);
1396: PetscFormKeySort(Nk, jackeys);
1397: for (k = 0, kp = 1; kp < Nk; ++kp) {
1398: if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
1399: ++k;
1400: if (kp != k) jackeys[k] = jackeys[kp];
1401: }
1402: }
1403: Nk = k;
1405: PetscDSGetWeakForm(ds, &wf);
1406: for (k = 0; k < Nk; ++k) {
1407: DMLabel label = jackeys[k].label;
1408: PetscInt val = jackeys[k].value;
1410: if (!label) {
1411: PetscObjectReference((PetscObject)allcellIS);
1412: cellIS = allcellIS;
1413: } else {
1414: IS pointIS;
1416: DMLabelGetStratumIS(label, val, &pointIS);
1417: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1418: ISDestroy(&pointIS);
1419: }
1420: DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user);
1421: ISDestroy(&cellIS);
1422: }
1423: PetscFree(jackeys);
1424: }
1425: }
1426: ISDestroy(&allcellIS);
1427: DMDestroy(&plex);
1428: return 0;
1429: }
1431: /*@
1432: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1434: Input Parameters:
1435: + dm - The mesh
1436: . X - Local input vector
1437: - user - The user context
1439: Output Parameter:
1440: . Jac - Jacobian matrix
1442: Note:
1443: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1444: like a GPU, or vectorize on a multicore machine.
1446: Level: developer
1448: .seealso: `DMPLEX`, `Mat`
1449: @*/
1450: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
1451: {
1452: DM plex;
1453: IS allcellIS;
1454: PetscBool hasJac, hasPrec;
1455: PetscInt Nds, s;
1457: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1458: DMPlexGetAllCells_Internal(plex, &allcellIS);
1459: DMGetNumDS(dm, &Nds);
1460: for (s = 0; s < Nds; ++s) {
1461: PetscDS ds;
1462: IS cellIS;
1463: PetscFormKey key;
1465: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1466: key.value = 0;
1467: key.field = 0;
1468: key.part = 0;
1469: if (!key.label) {
1470: PetscObjectReference((PetscObject)allcellIS);
1471: cellIS = allcellIS;
1472: } else {
1473: IS pointIS;
1475: key.value = 1;
1476: DMLabelGetStratumIS(key.label, key.value, &pointIS);
1477: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1478: ISDestroy(&pointIS);
1479: }
1480: if (!s) {
1481: PetscDSHasJacobian(ds, &hasJac);
1482: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1483: if (hasJac && hasPrec) MatZeroEntries(Jac);
1484: MatZeroEntries(JacP);
1485: }
1486: DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);
1487: ISDestroy(&cellIS);
1488: }
1489: ISDestroy(&allcellIS);
1490: DMDestroy(&plex);
1491: return 0;
1492: }
1494: struct _DMSNESJacobianMFCtx {
1495: DM dm;
1496: Vec X;
1497: void *ctx;
1498: };
1500: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
1501: {
1502: struct _DMSNESJacobianMFCtx *ctx;
1504: MatShellGetContext(A, &ctx);
1505: MatShellSetContext(A, NULL);
1506: DMDestroy(&ctx->dm);
1507: VecDestroy(&ctx->X);
1508: PetscFree(ctx);
1509: return 0;
1510: }
1512: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
1513: {
1514: struct _DMSNESJacobianMFCtx *ctx;
1516: MatShellGetContext(A, &ctx);
1517: DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx);
1518: return 0;
1519: }
1521: /*@
1522: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
1524: Collective
1526: Input Parameters:
1527: + dm - The `DM`
1528: . X - The evaluation point for the Jacobian
1529: - user - A user context, or NULL
1531: Output Parameter:
1532: . J - The `Mat`
1534: Level: advanced
1536: Note:
1537: Vec X is kept in `Mat` J, so updating X then updates the evaluation point.
1539: .seealso: `DM`, `DMSNESComputeJacobianAction()`
1540: @*/
1541: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
1542: {
1543: struct _DMSNESJacobianMFCtx *ctx;
1544: PetscInt n, N;
1546: MatCreate(PetscObjectComm((PetscObject)dm), J);
1547: MatSetType(*J, MATSHELL);
1548: VecGetLocalSize(X, &n);
1549: VecGetSize(X, &N);
1550: MatSetSizes(*J, n, n, N, N);
1551: PetscObjectReference((PetscObject)dm);
1552: PetscObjectReference((PetscObject)X);
1553: PetscMalloc1(1, &ctx);
1554: ctx->dm = dm;
1555: ctx->X = X;
1556: ctx->ctx = user;
1557: MatShellSetContext(*J, ctx);
1558: MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private);
1559: MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private);
1560: return 0;
1561: }
1563: /*
1564: MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.
1566: Input Parameters:
1567: + X - SNES linearization point
1568: . ovl - index set of overlapping subdomains
1570: Output Parameter:
1571: . J - unassembled (Neumann) local matrix
1573: Level: intermediate
1575: .seealso: `DMCreateNeumannOverlap()`, `MATIS`, `PCHPDDMSetAuxiliaryMat()`
1576: */
1577: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1578: {
1579: SNES snes;
1580: Mat pJ;
1581: DM ovldm, origdm;
1582: DMSNES sdm;
1583: PetscErrorCode (*bfun)(DM, Vec, void *);
1584: PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
1585: void *bctx, *jctx;
1587: PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ);
1589: PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm);
1591: MatGetDM(pJ, &ovldm);
1592: DMSNESGetBoundaryLocal(origdm, &bfun, &bctx);
1593: DMSNESSetBoundaryLocal(ovldm, bfun, bctx);
1594: DMSNESGetJacobianLocal(origdm, &jfun, &jctx);
1595: DMSNESSetJacobianLocal(ovldm, jfun, jctx);
1596: PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes);
1597: if (!snes) {
1598: SNESCreate(PetscObjectComm((PetscObject)ovl), &snes);
1599: SNESSetDM(snes, ovldm);
1600: PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes);
1601: PetscObjectDereference((PetscObject)snes);
1602: }
1603: DMGetDMSNES(ovldm, &sdm);
1604: VecLockReadPush(X);
1605: {
1606: void *ctx;
1607: PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
1608: DMSNESGetJacobian(ovldm, &J, &ctx);
1609: PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
1610: }
1611: VecLockReadPop(X);
1612: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1613: {
1614: Mat locpJ;
1616: MatISGetLocalMat(pJ, &locpJ);
1617: MatCopy(locpJ, J, SAME_NONZERO_PATTERN);
1618: }
1619: return 0;
1620: }
1622: /*@
1623: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, residual, and Jacobian.
1625: Input Parameters:
1626: + dm - The `DM` object
1627: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see `PetscDSAddBoundary()`)
1628: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see `PetscDSSetResidual()`)
1629: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see `PetscDSSetJacobian()`)
1631: Level: developer
1633: .seealso: `DMPLEX`, `SNES`
1634: @*/
1635: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1636: {
1637: DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, boundaryctx);
1638: DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, residualctx);
1639: DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, jacobianctx);
1640: PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex);
1641: return 0;
1642: }
1644: /*@C
1645: DMSNESCheckDiscretization - Check the discretization error of the exact solution
1647: Input Parameters:
1648: + snes - the `SNES` object
1649: . dm - the `DM`
1650: . t - the time
1651: . u - a `DM` vector
1652: - tol - A tolerance for the check, or -1 to print the results instead
1654: Output Parameters:
1655: . error - An array which holds the discretization error in each field, or NULL
1657: Note:
1658: The user must call `PetscDSSetExactSolution()` beforehand
1660: Level: developer
1662: .seealso: `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`, `PetscDSSetExactSolution()`
1663: @*/
1664: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1665: {
1666: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1667: void **ectxs;
1668: PetscReal *err;
1669: MPI_Comm comm;
1670: PetscInt Nf, f;
1677: DMComputeExactSolution(dm, t, u, NULL);
1678: VecViewFromOptions(u, NULL, "-vec_view");
1680: PetscObjectGetComm((PetscObject)snes, &comm);
1681: DMGetNumFields(dm, &Nf);
1682: PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);
1683: {
1684: PetscInt Nds, s;
1686: DMGetNumDS(dm, &Nds);
1687: for (s = 0; s < Nds; ++s) {
1688: PetscDS ds;
1689: DMLabel label;
1690: IS fieldIS;
1691: const PetscInt *fields;
1692: PetscInt dsNf, f;
1694: DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1695: PetscDSGetNumFields(ds, &dsNf);
1696: ISGetIndices(fieldIS, &fields);
1697: for (f = 0; f < dsNf; ++f) {
1698: const PetscInt field = fields[f];
1699: PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);
1700: }
1701: ISRestoreIndices(fieldIS, &fields);
1702: }
1703: }
1704: if (Nf > 1) {
1705: DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);
1706: if (tol >= 0.0) {
1708: } else if (error) {
1709: for (f = 0; f < Nf; ++f) error[f] = err[f];
1710: } else {
1711: PetscPrintf(comm, "L_2 Error: [");
1712: for (f = 0; f < Nf; ++f) {
1713: if (f) PetscPrintf(comm, ", ");
1714: PetscPrintf(comm, "%g", (double)err[f]);
1715: }
1716: PetscPrintf(comm, "]\n");
1717: }
1718: } else {
1719: DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);
1720: if (tol >= 0.0) {
1722: } else if (error) {
1723: error[0] = err[0];
1724: } else {
1725: PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]);
1726: }
1727: }
1728: PetscFree3(exacts, ectxs, err);
1729: return 0;
1730: }
1732: /*@C
1733: DMSNESCheckResidual - Check the residual of the exact solution
1735: Input Parameters:
1736: + snes - the `SNES` object
1737: . dm - the `DM`
1738: . u - a `DM` vector
1739: - tol - A tolerance for the check, or -1 to print the results instead
1741: Output Parameter:
1742: . residual - The residual norm of the exact solution, or NULL
1744: Level: developer
1746: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
1747: @*/
1748: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1749: {
1750: MPI_Comm comm;
1751: Vec r;
1752: PetscReal res;
1758: PetscObjectGetComm((PetscObject)snes, &comm);
1759: DMComputeExactSolution(dm, 0.0, u, NULL);
1760: VecDuplicate(u, &r);
1761: SNESComputeFunction(snes, u, r);
1762: VecNorm(r, NORM_2, &res);
1763: if (tol >= 0.0) {
1765: } else if (residual) {
1766: *residual = res;
1767: } else {
1768: PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
1769: VecChop(r, 1.0e-10);
1770: PetscObjectSetName((PetscObject)r, "Initial Residual");
1771: PetscObjectSetOptionsPrefix((PetscObject)r, "res_");
1772: VecViewFromOptions(r, NULL, "-vec_view");
1773: }
1774: VecDestroy(&r);
1775: return 0;
1776: }
1778: /*@C
1779: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
1781: Input Parameters:
1782: + snes - the `SNES` object
1783: . dm - the `DM`
1784: . u - a `DM` vector
1785: - tol - A tolerance for the check, or -1 to print the results instead
1787: Output Parameters:
1788: + isLinear - Flag indicaing that the function looks linear, or NULL
1789: - convRate - The rate of convergence of the linear model, or NULL
1791: Level: developer
1793: .seealso: `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
1794: @*/
1795: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1796: {
1797: MPI_Comm comm;
1798: PetscDS ds;
1799: Mat J, M;
1800: MatNullSpace nullspace;
1801: PetscReal slope, intercept;
1802: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
1809: PetscObjectGetComm((PetscObject)snes, &comm);
1810: DMComputeExactSolution(dm, 0.0, u, NULL);
1811: /* Create and view matrices */
1812: DMCreateMatrix(dm, &J);
1813: DMGetDS(dm, &ds);
1814: PetscDSHasJacobian(ds, &hasJac);
1815: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1816: if (hasJac && hasPrec) {
1817: DMCreateMatrix(dm, &M);
1818: SNESComputeJacobian(snes, u, J, M);
1819: PetscObjectSetName((PetscObject)M, "Preconditioning Matrix");
1820: PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_");
1821: MatViewFromOptions(M, NULL, "-mat_view");
1822: MatDestroy(&M);
1823: } else {
1824: SNESComputeJacobian(snes, u, J, J);
1825: }
1826: PetscObjectSetName((PetscObject)J, "Jacobian");
1827: PetscObjectSetOptionsPrefix((PetscObject)J, "jac_");
1828: MatViewFromOptions(J, NULL, "-mat_view");
1829: /* Check nullspace */
1830: MatGetNullSpace(J, &nullspace);
1831: if (nullspace) {
1832: PetscBool isNull;
1833: MatNullSpaceTest(nullspace, J, &isNull);
1835: }
1836: /* Taylor test */
1837: {
1838: PetscRandom rand;
1839: Vec du, uhat, r, rhat, df;
1840: PetscReal h;
1841: PetscReal *es, *hs, *errors;
1842: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1843: PetscInt Nv, v;
1845: /* Choose a perturbation direction */
1846: PetscRandomCreate(comm, &rand);
1847: VecDuplicate(u, &du);
1848: VecSetRandom(du, rand);
1849: PetscRandomDestroy(&rand);
1850: VecDuplicate(u, &df);
1851: MatMult(J, du, df);
1852: /* Evaluate residual at u, F(u), save in vector r */
1853: VecDuplicate(u, &r);
1854: SNESComputeFunction(snes, u, r);
1855: /* Look at the convergence of our Taylor approximation as we approach u */
1856: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
1857: ;
1858: PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
1859: VecDuplicate(u, &uhat);
1860: VecDuplicate(u, &rhat);
1861: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1862: VecWAXPY(uhat, h, du, u);
1863: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1864: SNESComputeFunction(snes, uhat, rhat);
1865: VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
1866: VecNorm(rhat, NORM_2, &errors[Nv]);
1868: es[Nv] = PetscLog10Real(errors[Nv]);
1869: hs[Nv] = PetscLog10Real(h);
1870: }
1871: VecDestroy(&uhat);
1872: VecDestroy(&rhat);
1873: VecDestroy(&df);
1874: VecDestroy(&r);
1875: VecDestroy(&du);
1876: for (v = 0; v < Nv; ++v) {
1877: if ((tol >= 0) && (errors[v] > tol)) break;
1878: else if (errors[v] > PETSC_SMALL) break;
1879: }
1880: if (v == Nv) isLin = PETSC_TRUE;
1881: PetscLinearRegression(Nv, hs, es, &slope, &intercept);
1882: PetscFree3(es, hs, errors);
1883: /* Slope should be about 2 */
1884: if (tol >= 0) {
1886: } else if (isLinear || convRate) {
1887: if (isLinear) *isLinear = isLin;
1888: if (convRate) *convRate = slope;
1889: } else {
1890: if (!isLin) PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope);
1891: else PetscPrintf(comm, "Function appears to be linear\n");
1892: }
1893: }
1894: MatDestroy(&J);
1895: return 0;
1896: }
1898: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1899: {
1900: DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);
1901: DMSNESCheckResidual(snes, dm, u, -1.0, NULL);
1902: DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);
1903: return 0;
1904: }
1906: /*@C
1907: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1909: Input Parameters:
1910: + snes - the `SNES` object
1911: - u - representative `SNES` vector
1913: Note:
1914: The user must call `PetscDSSetExactSolution()` beforehand
1916: Level: developer
1917: @*/
1918: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1919: {
1920: DM dm;
1921: Vec sol;
1922: PetscBool check;
1924: PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check);
1925: if (!check) return 0;
1926: SNESGetDM(snes, &dm);
1927: VecDuplicate(u, &sol);
1928: SNESSetSolution(snes, sol);
1929: DMSNESCheck_Internal(snes, dm, sol);
1930: VecDestroy(&sol);
1931: return 0;
1932: }