Actual source code: dmplexsnes.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/petscfeimpl.h>
8: static void pressure_Private(PetscInt dim, PetscInt Nf, PetscInt NfAux,
9: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
10: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
11: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar p[])
12: {
13: p[0] = u[uOff[1]];
14: }
16: /*
17: 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.
19: Collective on SNES
21: Input Parameters:
22: + snes - The SNES
23: . pfield - The field number for pressure
24: . nullspace - The pressure nullspace
25: . u - The solution vector
26: - ctx - An optional user context
28: Output Parameter:
29: . u - The solution with a continuum pressure integral of zero
31: Notes:
32: 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.
34: Level: developer
36: .seealso: SNESConvergedCorrectPressure()
37: */
38: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *ctx)
39: {
40: DM dm;
41: PetscDS ds;
42: const Vec *nullvecs;
43: PetscScalar pintd, *intc, *intn;
44: MPI_Comm comm;
45: PetscInt Nf, Nv;
49: PetscObjectGetComm((PetscObject) snes, &comm);
50: SNESGetDM(snes, &dm);
51: if (!dm) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
52: if (!nullspace) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
53: DMGetDS(dm, &ds);
54: PetscDSSetObjective(ds, pfield, pressure_Private);
55: MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs);
56: if (Nv != 1) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %D", Nv);
57: VecDot(nullvecs[0], u, &pintd);
58: if (PetscAbsScalar(pintd) > PETSC_SMALL) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g\n", (double) PetscRealPart(pintd));
59: PetscDSGetNumFields(ds, &Nf);
60: PetscMalloc2(Nf, &intc, Nf, &intn);
61: DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx);
62: DMPlexComputeIntegralFEM(dm, u, intc, ctx);
63: VecAXPY(u, -intc[pfield]/intn[pfield], nullvecs[0]);
64: #if defined (PETSC_USE_DEBUG)
65: DMPlexComputeIntegralFEM(dm, u, intc, ctx);
66: if (PetscAbsScalar(intc[pfield]) > PETSC_SMALL) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g\n", (double) PetscRealPart(intc[pfield]));
67: #endif
68: PetscFree2(intc, intn);
69: return(0);
70: }
72: /*@C
73: SNESConvergedCorrectPressure - Convergence test that adds 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. The convergence test itself just mimics SNESConvergedDefault().
75: Logically Collective on SNES
77: Input Parameters:
78: + snes - the SNES context
79: . it - the iteration (0 indicates before any Newton steps)
80: . xnorm - 2-norm of current iterate
81: . snorm - 2-norm of current step
82: . fnorm - 2-norm of function at current iterate
83: - ctx - Optional user context
85: Output Parameter:
86: . reason - SNES_CONVERGED_ITERATING, SNES_CONVERGED_ITS, or SNES_DIVERGED_FNORM_NAN
88: Notes:
89: In order to use this monitor, you must setup several PETSc structures. First fields must be added to the DM, and a PetscDS must be created with discretizations of those fields. We currently assume that the pressure field has index 1. The pressure field must have a nullspace, likely created using the DMSetNullSpaceConstructor() interface. Last we must be able to integrate the pressure over the domain, so the DM attached to the SNES must be a Plex at this time.
91: Level: advanced
93: .seealso: SNESConvergedDefault(), SNESSetConvergenceTest(), DMSetNullSpaceConstructor()
94: @*/
95: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
96: {
97: PetscBool monitorIntegral = PETSC_FALSE;
101: SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx);
102: if (monitorIntegral) {
103: Mat J;
104: Vec u;
105: MatNullSpace nullspace;
106: const Vec *nullvecs;
107: PetscScalar pintd;
109: SNESGetSolution(snes, &u);
110: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
111: MatGetNullSpace(J, &nullspace);
112: MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs);
113: VecDot(nullvecs[0], u, &pintd);
114: PetscInfo1(snes, "SNES: Discrete integral of pressure: %g\n", (double) PetscRealPart(pintd));
115: }
116: if (*reason > 0) {
117: Mat J;
118: Vec u;
119: MatNullSpace nullspace;
120: PetscInt pfield = 1;
122: SNESGetSolution(snes, &u);
123: SNESGetJacobian(snes, &J, NULL, NULL, NULL);
124: MatGetNullSpace(J, &nullspace);
125: SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx);
126: }
127: return(0);
128: }
130: /************************** Interpolation *******************************/
132: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
133: {
134: PetscBool isPlex;
138: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
139: if (isPlex) {
140: *plex = dm;
141: PetscObjectReference((PetscObject) dm);
142: } else {
143: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
144: if (!*plex) {
145: DMConvert(dm,DMPLEX,plex);
146: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
147: if (copy) {
148: PetscInt i;
149: PetscObject obj;
150: const char *comps[3] = {"A","dmAux","dmCh"};
152: DMCopyDMSNES(dm, *plex);
153: for (i = 0; i < 3; i++) {
154: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
155: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
156: }
157: }
158: } else {
159: PetscObjectReference((PetscObject) *plex);
160: }
161: }
162: return(0);
163: }
165: /*@C
166: DMInterpolationCreate - Creates a DMInterpolationInfo context
168: Collective
170: Input Parameter:
171: . comm - the communicator
173: Output Parameter:
174: . ctx - the context
176: Level: beginner
178: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationDestroy()
179: @*/
180: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
181: {
186: PetscNew(ctx);
188: (*ctx)->comm = comm;
189: (*ctx)->dim = -1;
190: (*ctx)->nInput = 0;
191: (*ctx)->points = NULL;
192: (*ctx)->cells = NULL;
193: (*ctx)->n = -1;
194: (*ctx)->coords = NULL;
195: return(0);
196: }
198: /*@C
199: DMInterpolationSetDim - Sets the spatial dimension for the interpolation context
201: Not collective
203: Input Parameters:
204: + ctx - the context
205: - dim - the spatial dimension
207: Level: intermediate
209: .seealso: DMInterpolationGetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
210: @*/
211: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
212: {
214: if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %D", dim);
215: ctx->dim = dim;
216: return(0);
217: }
219: /*@C
220: DMInterpolationGetDim - Gets the spatial dimension for the interpolation context
222: Not collective
224: Input Parameter:
225: . ctx - the context
227: Output Parameter:
228: . dim - the spatial dimension
230: Level: intermediate
232: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
233: @*/
234: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
235: {
238: *dim = ctx->dim;
239: return(0);
240: }
242: /*@C
243: DMInterpolationSetDof - Sets the number of fields interpolated at a point for the interpolation context
245: Not collective
247: Input Parameters:
248: + ctx - the context
249: - dof - the number of fields
251: Level: intermediate
253: .seealso: DMInterpolationGetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
254: @*/
255: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
256: {
258: if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %D", dof);
259: ctx->dof = dof;
260: return(0);
261: }
263: /*@C
264: DMInterpolationGetDof - Gets the number of fields interpolated at a point for the interpolation context
266: Not collective
268: Input Parameter:
269: . ctx - the context
271: Output Parameter:
272: . dof - the number of fields
274: Level: intermediate
276: .seealso: DMInterpolationSetDof(), DMInterpolationEvaluate(), DMInterpolationAddPoints()
277: @*/
278: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
279: {
282: *dof = ctx->dof;
283: return(0);
284: }
286: /*@C
287: DMInterpolationAddPoints - Add points at which we will interpolate the fields
289: Not collective
291: Input Parameters:
292: + ctx - the context
293: . n - the number of points
294: - points - the coordinates for each point, an array of size n * dim
296: Note: The coordinate information is copied.
298: Level: intermediate
300: .seealso: DMInterpolationSetDim(), DMInterpolationEvaluate(), DMInterpolationCreate()
301: @*/
302: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
303: {
307: if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
308: if (ctx->points) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
309: ctx->nInput = n;
311: PetscMalloc1(n*ctx->dim, &ctx->points);
312: PetscArraycpy(ctx->points, points, n*ctx->dim);
313: return(0);
314: }
316: /*@C
317: DMInterpolationSetUp - Compute spatial indices for point location during interpolation
319: Collective on ctx
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: 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: PetscErrorCode ierr;
338: Vec pointVec;
339: PetscSF cellSF;
340: PetscLayout layout;
341: PetscReal *globalPoints;
342: PetscScalar *globalPointsScalar;
343: const PetscInt *ranges;
344: PetscMPIInt *counts, *displs;
345: const PetscSFNode *foundCells;
346: const PetscInt *foundPoints;
347: PetscMPIInt *foundProcs, *globalProcs;
348: PetscInt n, N, numFound;
352: MPI_Comm_size(comm, &size);
353: MPI_Comm_rank(comm, &rank);
354: if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
355: /* Locate points */
356: n = ctx->nInput;
357: if (!redundantPoints) {
358: PetscLayoutCreate(comm, &layout);
359: PetscLayoutSetBlockSize(layout, 1);
360: PetscLayoutSetLocalSize(layout, n);
361: PetscLayoutSetUp(layout);
362: PetscLayoutGetSize(layout, &N);
363: /* Communicate all points to all processes */
364: PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
365: PetscLayoutGetRanges(layout, &ranges);
366: for (p = 0; p < size; ++p) {
367: counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
368: displs[p] = ranges[p]*ctx->dim;
369: }
370: MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
371: } else {
372: N = n;
373: globalPoints = ctx->points;
374: counts = displs = NULL;
375: layout = NULL;
376: }
377: #if 0
378: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
379: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
380: #else
381: #if defined(PETSC_USE_COMPLEX)
382: PetscMalloc1(N*ctx->dim,&globalPointsScalar);
383: for (i=0; i<N*ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
384: #else
385: globalPointsScalar = globalPoints;
386: #endif
387: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
388: PetscMalloc2(N,&foundProcs,N,&globalProcs);
389: for (p = 0; p < N; ++p) {foundProcs[p] = size;}
390: cellSF = NULL;
391: DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
392: PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
393: #endif
394: for (p = 0; p < numFound; ++p) {
395: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
396: }
397: /* Let the lowest rank process own each point */
398: MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
399: ctx->n = 0;
400: for (p = 0; p < N; ++p) {
401: if (globalProcs[p] == size) {
402: if (!ignoreOutsideDomain) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %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), (double)(ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0));
403: else if (!rank) ++ctx->n;
404: } else if (globalProcs[p] == rank) ++ctx->n;
405: }
406: /* Create coordinates vector and array of owned cells */
407: PetscMalloc1(ctx->n, &ctx->cells);
408: VecCreate(comm, &ctx->coords);
409: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
410: VecSetBlockSize(ctx->coords, ctx->dim);
411: VecSetType(ctx->coords,VECSTANDARD);
412: 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) {
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: VecRestoreArray(ctx->coords, &a);
430: #if 0
431: PetscFree3(foundCells,foundProcs,globalProcs);
432: #else
433: PetscFree2(foundProcs,globalProcs);
434: PetscSFDestroy(&cellSF);
435: VecDestroy(&pointVec);
436: #endif
437: if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
438: if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
439: PetscLayoutDestroy(&layout);
440: return(0);
441: }
443: /*@C
444: DMInterpolationGetCoordinates - Gets a Vec with the coordinates of each interpolation point
446: Collective on ctx
448: Input Parameter:
449: . ctx - the context
451: Output Parameter:
452: . coordinates - the coordinates of interpolation points
454: Note: The local vector entries correspond to interpolation points lying on this process, according to the associated DM. This is a borrowed vector that the user should not destroy.
456: Level: intermediate
458: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
459: @*/
460: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
461: {
464: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
465: *coordinates = ctx->coords;
466: return(0);
467: }
469: /*@C
470: DMInterpolationGetVector - Gets a Vec which can hold all the interpolated field values
472: Collective on ctx
474: Input Parameter:
475: . ctx - the context
477: Output Parameter:
478: . v - a vector capable of holding the interpolated field values
480: Note: This vector should be returned using DMInterpolationRestoreVector().
482: Level: intermediate
484: .seealso: DMInterpolationRestoreVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
485: @*/
486: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
487: {
492: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
493: VecCreate(ctx->comm, v);
494: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
495: VecSetBlockSize(*v, ctx->dof);
496: VecSetType(*v,VECSTANDARD);
497: return(0);
498: }
500: /*@C
501: DMInterpolationRestoreVector - Returns a Vec which can hold all the interpolated field values
503: Collective on ctx
505: Input Parameters:
506: + ctx - the context
507: - v - a vector capable of holding the interpolated field values
509: Level: intermediate
511: .seealso: DMInterpolationGetVector(), DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
512: @*/
513: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
514: {
519: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
520: VecDestroy(v);
521: return(0);
522: }
524: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
525: {
526: PetscReal *v0, *J, *invJ, detJ;
527: const PetscScalar *coords;
528: PetscScalar *a;
529: PetscInt p;
533: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
534: VecGetArrayRead(ctx->coords, &coords);
535: VecGetArray(v, &a);
536: for (p = 0; p < ctx->n; ++p) {
537: PetscInt c = ctx->cells[p];
538: PetscScalar *x = NULL;
539: PetscReal xi[4];
540: PetscInt d, f, comp;
542: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
543: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
544: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
545: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
547: for (d = 0; d < ctx->dim; ++d) {
548: xi[d] = 0.0;
549: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
550: 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];
551: }
552: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
553: }
554: VecRestoreArray(v, &a);
555: VecRestoreArrayRead(ctx->coords, &coords);
556: PetscFree3(v0, J, invJ);
557: return(0);
558: }
560: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
561: {
562: PetscReal *v0, *J, *invJ, detJ;
563: const PetscScalar *coords;
564: PetscScalar *a;
565: PetscInt p;
569: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
570: VecGetArrayRead(ctx->coords, &coords);
571: VecGetArray(v, &a);
572: for (p = 0; p < ctx->n; ++p) {
573: PetscInt c = ctx->cells[p];
574: const PetscInt order[3] = {2, 1, 3};
575: PetscScalar *x = NULL;
576: PetscReal xi[4];
577: PetscInt d, f, comp;
579: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
580: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double)detJ, c);
581: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
582: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
584: for (d = 0; d < ctx->dim; ++d) {
585: xi[d] = 0.0;
586: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
587: 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];
588: }
589: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
590: }
591: VecRestoreArray(v, &a);
592: VecRestoreArrayRead(ctx->coords, &coords);
593: PetscFree3(v0, J, invJ);
594: return(0);
595: }
597: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
598: {
599: const PetscScalar *vertices = (const PetscScalar*) ctx;
600: const PetscScalar x0 = vertices[0];
601: const PetscScalar y0 = vertices[1];
602: const PetscScalar x1 = vertices[2];
603: const PetscScalar y1 = vertices[3];
604: const PetscScalar x2 = vertices[4];
605: const PetscScalar y2 = vertices[5];
606: const PetscScalar x3 = vertices[6];
607: const PetscScalar y3 = vertices[7];
608: const PetscScalar f_1 = x1 - x0;
609: const PetscScalar g_1 = y1 - y0;
610: const PetscScalar f_3 = x3 - x0;
611: const PetscScalar g_3 = y3 - y0;
612: const PetscScalar f_01 = x2 - x1 - x3 + x0;
613: const PetscScalar g_01 = y2 - y1 - y3 + y0;
614: const PetscScalar *ref;
615: PetscScalar *real;
616: PetscErrorCode ierr;
619: VecGetArrayRead(Xref, &ref);
620: VecGetArray(Xreal, &real);
621: {
622: const PetscScalar p0 = ref[0];
623: const PetscScalar p1 = ref[1];
625: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
626: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
627: }
628: PetscLogFlops(28);
629: VecRestoreArrayRead(Xref, &ref);
630: VecRestoreArray(Xreal, &real);
631: return(0);
632: }
634: #include <petsc/private/dmimpl.h>
635: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
636: {
637: const PetscScalar *vertices = (const PetscScalar*) ctx;
638: const PetscScalar x0 = vertices[0];
639: const PetscScalar y0 = vertices[1];
640: const PetscScalar x1 = vertices[2];
641: const PetscScalar y1 = vertices[3];
642: const PetscScalar x2 = vertices[4];
643: const PetscScalar y2 = vertices[5];
644: const PetscScalar x3 = vertices[6];
645: const PetscScalar y3 = vertices[7];
646: const PetscScalar f_01 = x2 - x1 - x3 + x0;
647: const PetscScalar g_01 = y2 - y1 - y3 + y0;
648: const PetscScalar *ref;
649: PetscErrorCode ierr;
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; values[1] = (x3 - x0 + f_01*x) * 0.5;
660: values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
661: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
662: }
663: PetscLogFlops(30);
664: VecRestoreArrayRead(Xref, &ref);
665: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
666: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
667: return(0);
668: }
670: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
671: {
672: DM dmCoord;
673: PetscFE fem = NULL;
674: SNES snes;
675: KSP ksp;
676: PC pc;
677: Vec coordsLocal, r, ref, real;
678: Mat J;
679: PetscTabulation T;
680: const PetscScalar *coords;
681: PetscScalar *a;
682: PetscReal xir[2];
683: PetscInt Nf, p;
684: const PetscInt dof = ctx->dof;
688: DMGetNumFields(dm, &Nf);
689: if (Nf) {DMGetField(dm, 0, NULL, (PetscObject *) &fem);}
690: DMGetCoordinatesLocal(dm, &coordsLocal);
691: DMGetCoordinateDM(dm, &dmCoord);
692: SNESCreate(PETSC_COMM_SELF, &snes);
693: SNESSetOptionsPrefix(snes, "quad_interp_");
694: VecCreate(PETSC_COMM_SELF, &r);
695: VecSetSizes(r, 2, 2);
696: VecSetType(r,dm->vectype);
697: VecDuplicate(r, &ref);
698: VecDuplicate(r, &real);
699: MatCreate(PETSC_COMM_SELF, &J);
700: MatSetSizes(J, 2, 2, 2, 2);
701: MatSetType(J, MATSEQDENSE);
702: MatSetUp(J);
703: SNESSetFunction(snes, r, QuadMap_Private, NULL);
704: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
705: SNESGetKSP(snes, &ksp);
706: KSPGetPC(ksp, &pc);
707: PCSetType(pc, PCLU);
708: SNESSetFromOptions(snes);
710: VecGetArrayRead(ctx->coords, &coords);
711: VecGetArray(v, &a);
712: PetscFECreateTabulation(fem, 1, 1, xir, 0, &T);
713: for (p = 0; p < ctx->n; ++p) {
714: PetscScalar *x = NULL, *vertices = NULL;
715: PetscScalar *xi;
716: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
718: /* Can make this do all points at once */
719: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
720: if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 4*2);
721: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
722: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
723: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
724: VecGetArray(real, &xi);
725: xi[0] = coords[p*ctx->dim+0];
726: xi[1] = coords[p*ctx->dim+1];
727: VecRestoreArray(real, &xi);
728: SNESSolve(snes, real, ref);
729: VecGetArray(ref, &xi);
730: xir[0] = PetscRealPart(xi[0]);
731: xir[1] = PetscRealPart(xi[1]);
732: if (4*dof != xSize) {
733: PetscInt d;
735: xir[0] = 2.0*xir[0] - 1.0; xir[1] = 2.0*xir[1] - 1.0;
736: PetscFEComputeTabulation(fem, 1, xir, 0, T);
737: for (comp = 0; comp < dof; ++comp) {
738: a[p*dof+comp] = 0.0;
739: for (d = 0; d < xSize/dof; ++d) {
740: a[p*dof+comp] += x[d*dof+comp]*T->T[0][d*dof+comp];
741: }
742: }
743: } else {
744: for (comp = 0; comp < dof; ++comp)
745: 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];
746: }
747: VecRestoreArray(ref, &xi);
748: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
749: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
750: }
751: PetscTabulationDestroy(&T);
752: VecRestoreArray(v, &a);
753: VecRestoreArrayRead(ctx->coords, &coords);
755: SNESDestroy(&snes);
756: VecDestroy(&r);
757: VecDestroy(&ref);
758: VecDestroy(&real);
759: MatDestroy(&J);
760: return(0);
761: }
763: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
764: {
765: const PetscScalar *vertices = (const PetscScalar*) ctx;
766: const PetscScalar x0 = vertices[0];
767: const PetscScalar y0 = vertices[1];
768: const PetscScalar z0 = vertices[2];
769: const PetscScalar x1 = vertices[9];
770: const PetscScalar y1 = vertices[10];
771: const PetscScalar z1 = vertices[11];
772: const PetscScalar x2 = vertices[6];
773: const PetscScalar y2 = vertices[7];
774: const PetscScalar z2 = vertices[8];
775: const PetscScalar x3 = vertices[3];
776: const PetscScalar y3 = vertices[4];
777: const PetscScalar z3 = vertices[5];
778: const PetscScalar x4 = vertices[12];
779: const PetscScalar y4 = vertices[13];
780: const PetscScalar z4 = vertices[14];
781: const PetscScalar x5 = vertices[15];
782: const PetscScalar y5 = vertices[16];
783: const PetscScalar z5 = vertices[17];
784: const PetscScalar x6 = vertices[18];
785: const PetscScalar y6 = vertices[19];
786: const PetscScalar z6 = vertices[20];
787: const PetscScalar x7 = vertices[21];
788: const PetscScalar y7 = vertices[22];
789: const PetscScalar z7 = vertices[23];
790: const PetscScalar f_1 = x1 - x0;
791: const PetscScalar g_1 = y1 - y0;
792: const PetscScalar h_1 = z1 - z0;
793: const PetscScalar f_3 = x3 - x0;
794: const PetscScalar g_3 = y3 - y0;
795: const PetscScalar h_3 = z3 - z0;
796: const PetscScalar f_4 = x4 - x0;
797: const PetscScalar g_4 = y4 - y0;
798: const PetscScalar h_4 = z4 - z0;
799: const PetscScalar f_01 = x2 - x1 - x3 + x0;
800: const PetscScalar g_01 = y2 - y1 - y3 + y0;
801: const PetscScalar h_01 = z2 - z1 - z3 + z0;
802: const PetscScalar f_12 = x7 - x3 - x4 + x0;
803: const PetscScalar g_12 = y7 - y3 - y4 + y0;
804: const PetscScalar h_12 = z7 - z3 - z4 + z0;
805: const PetscScalar f_02 = x5 - x1 - x4 + x0;
806: const PetscScalar g_02 = y5 - y1 - y4 + y0;
807: const PetscScalar h_02 = z5 - z1 - z4 + z0;
808: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
809: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
810: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
811: const PetscScalar *ref;
812: PetscScalar *real;
813: PetscErrorCode ierr;
816: VecGetArrayRead(Xref, &ref);
817: VecGetArray(Xreal, &real);
818: {
819: const PetscScalar p0 = ref[0];
820: const PetscScalar p1 = ref[1];
821: const PetscScalar p2 = ref[2];
823: 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;
824: 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;
825: 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;
826: }
827: PetscLogFlops(114);
828: VecRestoreArrayRead(Xref, &ref);
829: VecRestoreArray(Xreal, &real);
830: return(0);
831: }
833: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
834: {
835: const PetscScalar *vertices = (const PetscScalar*) ctx;
836: const PetscScalar x0 = vertices[0];
837: const PetscScalar y0 = vertices[1];
838: const PetscScalar z0 = vertices[2];
839: const PetscScalar x1 = vertices[9];
840: const PetscScalar y1 = vertices[10];
841: const PetscScalar z1 = vertices[11];
842: const PetscScalar x2 = vertices[6];
843: const PetscScalar y2 = vertices[7];
844: const PetscScalar z2 = vertices[8];
845: const PetscScalar x3 = vertices[3];
846: const PetscScalar y3 = vertices[4];
847: const PetscScalar z3 = vertices[5];
848: const PetscScalar x4 = vertices[12];
849: const PetscScalar y4 = vertices[13];
850: const PetscScalar z4 = vertices[14];
851: const PetscScalar x5 = vertices[15];
852: const PetscScalar y5 = vertices[16];
853: const PetscScalar z5 = vertices[17];
854: const PetscScalar x6 = vertices[18];
855: const PetscScalar y6 = vertices[19];
856: const PetscScalar z6 = vertices[20];
857: const PetscScalar x7 = vertices[21];
858: const PetscScalar y7 = vertices[22];
859: const PetscScalar z7 = vertices[23];
860: const PetscScalar f_xy = x2 - x1 - x3 + x0;
861: const PetscScalar g_xy = y2 - y1 - y3 + y0;
862: const PetscScalar h_xy = z2 - z1 - z3 + z0;
863: const PetscScalar f_yz = x7 - x3 - x4 + x0;
864: const PetscScalar g_yz = y7 - y3 - y4 + y0;
865: const PetscScalar h_yz = z7 - z3 - z4 + z0;
866: const PetscScalar f_xz = x5 - x1 - x4 + x0;
867: const PetscScalar g_xz = y5 - y1 - y4 + y0;
868: const PetscScalar h_xz = z5 - z1 - z4 + z0;
869: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
870: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
871: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
872: const PetscScalar *ref;
873: PetscErrorCode ierr;
876: VecGetArrayRead(Xref, &ref);
877: {
878: const PetscScalar x = ref[0];
879: const PetscScalar y = ref[1];
880: const PetscScalar z = ref[2];
881: const PetscInt rows[3] = {0, 1, 2};
882: PetscScalar values[9];
884: values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
885: values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
886: values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
887: values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
888: values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
889: values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
890: values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
891: values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
892: values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
894: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
895: }
896: PetscLogFlops(152);
897: VecRestoreArrayRead(Xref, &ref);
898: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
899: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
900: return(0);
901: }
903: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
904: {
905: DM dmCoord;
906: SNES snes;
907: KSP ksp;
908: PC pc;
909: Vec coordsLocal, r, ref, real;
910: Mat J;
911: const PetscScalar *coords;
912: PetscScalar *a;
913: PetscInt p;
917: DMGetCoordinatesLocal(dm, &coordsLocal);
918: DMGetCoordinateDM(dm, &dmCoord);
919: SNESCreate(PETSC_COMM_SELF, &snes);
920: SNESSetOptionsPrefix(snes, "hex_interp_");
921: VecCreate(PETSC_COMM_SELF, &r);
922: VecSetSizes(r, 3, 3);
923: VecSetType(r,dm->vectype);
924: VecDuplicate(r, &ref);
925: VecDuplicate(r, &real);
926: MatCreate(PETSC_COMM_SELF, &J);
927: MatSetSizes(J, 3, 3, 3, 3);
928: MatSetType(J, MATSEQDENSE);
929: MatSetUp(J);
930: SNESSetFunction(snes, r, HexMap_Private, NULL);
931: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
932: SNESGetKSP(snes, &ksp);
933: KSPGetPC(ksp, &pc);
934: PCSetType(pc, PCLU);
935: SNESSetFromOptions(snes);
937: VecGetArrayRead(ctx->coords, &coords);
938: VecGetArray(v, &a);
939: for (p = 0; p < ctx->n; ++p) {
940: PetscScalar *x = NULL, *vertices = NULL;
941: PetscScalar *xi;
942: PetscReal xir[3];
943: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
945: /* Can make this do all points at once */
946: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
947: if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %d", coordSize, 8*3);
948: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
949: if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %D should be %D", xSize, 8*ctx->dof);
950: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
951: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
952: VecGetArray(real, &xi);
953: xi[0] = coords[p*ctx->dim+0];
954: xi[1] = coords[p*ctx->dim+1];
955: xi[2] = coords[p*ctx->dim+2];
956: VecRestoreArray(real, &xi);
957: SNESSolve(snes, real, ref);
958: VecGetArray(ref, &xi);
959: xir[0] = PetscRealPart(xi[0]);
960: xir[1] = PetscRealPart(xi[1]);
961: xir[2] = PetscRealPart(xi[2]);
962: for (comp = 0; comp < ctx->dof; ++comp) {
963: a[p*ctx->dof+comp] =
964: x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
965: x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
966: x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
967: x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
968: x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
969: x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
970: x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
971: x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
972: }
973: VecRestoreArray(ref, &xi);
974: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
975: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
976: }
977: VecRestoreArray(v, &a);
978: VecRestoreArrayRead(ctx->coords, &coords);
980: SNESDestroy(&snes);
981: VecDestroy(&r);
982: VecDestroy(&ref);
983: VecDestroy(&real);
984: MatDestroy(&J);
985: return(0);
986: }
988: /*@C
989: DMInterpolationEvaluate - Using the input from dm and x, calculates interpolated field values at the interpolation points.
991: Input Parameters:
992: + ctx - The DMInterpolationInfo context
993: . dm - The DM
994: - x - The local vector containing the field to be interpolated
996: Output Parameters:
997: . v - The vector containing the interpolated values
999: Note: A suitable v can be obtained using DMInterpolationGetVector().
1001: Level: beginner
1003: .seealso: DMInterpolationGetVector(), DMInterpolationAddPoints(), DMInterpolationCreate()
1004: @*/
1005: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
1006: {
1007: PetscInt n;
1014: VecGetLocalSize(v, &n);
1015: if (n != ctx->n*ctx->dof) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid input vector size %D should be %D", n, ctx->n*ctx->dof);
1016: if (n) {
1017: PetscDS ds;
1018: DMPolytopeType ct;
1019: PetscBool done = PETSC_FALSE;
1021: DMGetDS(dm, &ds);
1022: if (ds) {
1023: const PetscScalar *coords;
1024: PetscScalar *interpolant;
1025: PetscInt cdim, d, p, Nf, field, c = 0;
1027: DMGetCoordinateDim(dm, &cdim);
1028: PetscDSGetNumFields(ds, &Nf);
1029: for (field = 0; field < Nf; ++field) {
1030: PetscTabulation T;
1031: PetscFE fe;
1032: PetscClassId id;
1033: PetscReal xi[3];
1034: PetscInt Nc, f, fc;
1036: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
1037: PetscObjectGetClassId((PetscObject) fe, &id);
1038: if (id != PETSCFE_CLASSID) break;
1039: PetscFEGetNumComponents(fe, &Nc);
1040: VecGetArrayRead(ctx->coords, &coords);
1041: VecGetArrayWrite(v, &interpolant);
1042: for (p = 0; p < ctx->n; ++p) {
1043: PetscScalar *xa = NULL;
1044: PetscReal pcoords[3];
1046: if (ctx->cells[p] < 0) continue;
1047: for (d = 0; d < cdim; ++d) pcoords[d] = PetscRealPart(coords[p*cdim+d]);
1048: DMPlexCoordinatesToReference(dm, ctx->cells[p], 1, pcoords, xi);
1049: DMPlexVecGetClosure(dm, NULL, x, ctx->cells[p], NULL, &xa);
1050: PetscFECreateTabulation(fe, 1, 1, xi, 0, &T);
1051: {
1052: const PetscReal *basis = T->T[0];
1053: const PetscInt Nb = T->Nb;
1054: const PetscInt Nc = T->Nc;
1055: for (fc = 0; fc < Nc; ++fc) {
1056: interpolant[p*ctx->dof+c+fc] = 0.0;
1057: for (f = 0; f < Nb; ++f) {
1058: interpolant[p*ctx->dof+c+fc] += xa[f]*basis[(0*Nb + f)*Nc + fc];
1059: }
1060: }
1061: }
1062: PetscTabulationDestroy(&T);
1063: DMPlexVecRestoreClosure(dm, NULL, x, ctx->cells[p], NULL, &xa);
1064: }
1065: VecRestoreArrayWrite(v, &interpolant);
1066: VecRestoreArrayRead(ctx->coords, &coords);
1067: c += Nc;
1068: }
1069: if (field == Nf) {
1070: done = PETSC_TRUE;
1071: if (c != ctx->dof) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Total components %D != %D dof specified for interpolation", c, ctx->dof);
1072: }
1073: }
1074: if (!done) {
1075: /* TODO Check each cell individually */
1076: DMPlexGetCellType(dm, ctx->cells[0], &ct);
1077: switch (ct) {
1078: case DM_POLYTOPE_TRIANGLE: DMInterpolate_Triangle_Private(ctx, dm, x, v);break;
1079: case DM_POLYTOPE_QUADRILATERAL: DMInterpolate_Quad_Private(ctx, dm, x, v);break;
1080: case DM_POLYTOPE_TETRAHEDRON: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);break;
1081: case DM_POLYTOPE_HEXAHEDRON: DMInterpolate_Hex_Private(ctx, dm, x, v);break;
1082: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support fpr cell type %s", DMPolytopeTypes[ct]);
1083: }
1084: }
1085: }
1086: return(0);
1087: }
1089: /*@C
1090: DMInterpolationDestroy - Destroys a DMInterpolationInfo context
1092: Collective on ctx
1094: Input Parameter:
1095: . ctx - the context
1097: Level: beginner
1099: .seealso: DMInterpolationEvaluate(), DMInterpolationAddPoints(), DMInterpolationCreate()
1100: @*/
1101: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
1102: {
1107: VecDestroy(&(*ctx)->coords);
1108: PetscFree((*ctx)->points);
1109: PetscFree((*ctx)->cells);
1110: PetscFree(*ctx);
1111: *ctx = NULL;
1112: return(0);
1113: }
1115: /*@C
1116: SNESMonitorFields - Monitors the residual for each field separately
1118: Collective on SNES
1120: Input Parameters:
1121: + snes - the SNES context
1122: . its - iteration number
1123: . fgnorm - 2-norm of residual
1124: - vf - PetscViewerAndFormat of type ASCII
1126: Notes:
1127: This routine prints the residual norm at each iteration.
1129: Level: intermediate
1131: .seealso: SNESMonitorSet(), SNESMonitorDefault()
1132: @*/
1133: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
1134: {
1135: PetscViewer viewer = vf->viewer;
1136: Vec res;
1137: DM dm;
1138: PetscSection s;
1139: const PetscScalar *r;
1140: PetscReal *lnorms, *norms;
1141: PetscInt numFields, f, pStart, pEnd, p;
1142: PetscErrorCode ierr;
1146: SNESGetFunction(snes, &res, NULL, NULL);
1147: SNESGetDM(snes, &dm);
1148: DMGetLocalSection(dm, &s);
1149: PetscSectionGetNumFields(s, &numFields);
1150: PetscSectionGetChart(s, &pStart, &pEnd);
1151: PetscCalloc2(numFields, &lnorms, numFields, &norms);
1152: VecGetArrayRead(res, &r);
1153: for (p = pStart; p < pEnd; ++p) {
1154: for (f = 0; f < numFields; ++f) {
1155: PetscInt fdof, foff, d;
1157: PetscSectionGetFieldDof(s, p, f, &fdof);
1158: PetscSectionGetFieldOffset(s, p, f, &foff);
1159: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
1160: }
1161: }
1162: VecRestoreArrayRead(res, &r);
1163: MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
1164: PetscViewerPushFormat(viewer,vf->format);
1165: PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
1166: PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
1167: for (f = 0; f < numFields; ++f) {
1168: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
1169: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
1170: }
1171: PetscViewerASCIIPrintf(viewer, "]\n");
1172: PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
1173: PetscViewerPopFormat(viewer);
1174: PetscFree2(lnorms, norms);
1175: return(0);
1176: }
1178: /********************* Residual Computation **************************/
1180: PetscErrorCode DMPlexGetAllCells_Internal(DM plex, IS *cellIS)
1181: {
1182: PetscInt depth;
1186: DMPlexGetDepth(plex, &depth);
1187: DMGetStratumIS(plex, "dim", depth, cellIS);
1188: if (!*cellIS) {DMGetStratumIS(plex, "depth", depth, cellIS);}
1189: return(0);
1190: }
1192: /*@
1193: DMPlexSNESComputeResidualFEM - Sums the local residual into vector F from the local input X using pointwise functions specified by the user
1195: Input Parameters:
1196: + dm - The mesh
1197: . X - Local solution
1198: - user - The user context
1200: Output Parameter:
1201: . F - Local output vector
1203: Notes:
1204: The residual is summed into F; the caller is responsible for using VecZeroEntries() or otherwise ensuring that any data in F is intentional.
1206: Level: developer
1208: .seealso: DMPlexComputeJacobianAction()
1209: @*/
1210: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1211: {
1212: DM plex;
1213: IS allcellIS;
1214: PetscInt Nds, s;
1218: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1219: DMPlexGetAllCells_Internal(plex, &allcellIS);
1220: DMGetNumDS(dm, &Nds);
1221: for (s = 0; s < Nds; ++s) {
1222: PetscDS ds;
1223: IS cellIS;
1224: PetscHashFormKey key;
1226: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1227: key.value = 0;
1228: key.field = 0;
1229: if (!key.label) {
1230: PetscObjectReference((PetscObject) allcellIS);
1231: cellIS = allcellIS;
1232: } else {
1233: IS pointIS;
1235: key.value = 1;
1236: DMLabelGetStratumIS(key.label, key.value, &pointIS);
1237: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1238: ISDestroy(&pointIS);
1239: }
1240: DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1241: ISDestroy(&cellIS);
1242: }
1243: ISDestroy(&allcellIS);
1244: DMDestroy(&plex);
1245: return(0);
1246: }
1248: PetscErrorCode DMSNESComputeResidual(DM dm, Vec X, Vec F, void *user)
1249: {
1250: DM plex;
1251: IS allcellIS;
1252: PetscInt Nds, s;
1256: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1257: DMPlexGetAllCells_Internal(plex, &allcellIS);
1258: DMGetNumDS(dm, &Nds);
1259: for (s = 0; s < Nds; ++s) {
1260: PetscDS ds;
1261: DMLabel label;
1262: IS cellIS;
1264: DMGetRegionNumDS(dm, s, &label, NULL, &ds);
1265: {
1266: PetscHMapForm resmap[2] = {ds->wf->f0, ds->wf->f1};
1267: PetscWeakForm wf;
1268: PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
1269: PetscHashFormKey *reskeys;
1271: /* Get unique residual keys */
1272: for (m = 0; m < Nm; ++m) {
1273: PetscInt Nkm;
1274: PetscHMapFormGetSize(resmap[m], &Nkm);
1275: Nk += Nkm;
1276: }
1277: PetscMalloc1(Nk, &reskeys);
1278: for (m = 0; m < Nm; ++m) {
1279: PetscHMapFormGetKeys(resmap[m], &off, reskeys);
1280: }
1281: if (off != Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %D should be %D", off, Nk);
1282: PetscHashFormKeySort(Nk, reskeys);
1283: for (k = 0, kp = 1; kp < Nk; ++kp) {
1284: if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
1285: ++k;
1286: if (kp != k) reskeys[k] = reskeys[kp];
1287: }
1288: }
1289: Nk = k;
1291: PetscDSGetWeakForm(ds, &wf);
1292: for (k = 0; k < Nk; ++k) {
1293: DMLabel label = reskeys[k].label;
1294: PetscInt val = reskeys[k].value;
1296: if (!label) {
1297: PetscObjectReference((PetscObject) allcellIS);
1298: cellIS = allcellIS;
1299: } else {
1300: IS pointIS;
1302: DMLabelGetStratumIS(label, val, &pointIS);
1303: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1304: ISDestroy(&pointIS);
1305: }
1306: DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user);
1307: ISDestroy(&cellIS);
1308: }
1309: PetscFree(reskeys);
1310: }
1311: }
1312: ISDestroy(&allcellIS);
1313: DMDestroy(&plex);
1314: return(0);
1315: }
1317: /*@
1318: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
1320: Input Parameters:
1321: + dm - The mesh
1322: - user - The user context
1324: Output Parameter:
1325: . X - Local solution
1327: Level: developer
1329: .seealso: DMPlexComputeJacobianAction()
1330: @*/
1331: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1332: {
1333: DM plex;
1337: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1338: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1339: DMDestroy(&plex);
1340: return(0);
1341: }
1343: /*@
1344: DMPlexComputeJacobianAction - Form the local portion of the Jacobian action Z = J(X) Y at the local solution X using pointwise functions specified by the user.
1346: Input Parameters:
1347: + dm - The mesh
1348: . cellIS - index set for the cells or NULL to use entire depth=dim stratum
1349: . t - The time
1350: . X_tShift - The multiplier for the Jacobian with respect to X_t
1351: . X - Local solution vector
1352: . X_t - Time-derivative of the local solution vector
1353: . Y - Local input vector
1354: - user - The user context
1356: Output Parameter:
1357: . Z - Local output vector
1359: Note:
1360: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1361: like a GPU, or vectorize on a multicore machine.
1363: Level: developer
1365: .seealso: DMPlexSNESComputeResidualFEM()
1366: @*/
1367: PetscErrorCode DMPlexComputeJacobianAction(DM dm, IS cellIS, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Vec Y, Vec Z, void *user)
1368: {
1369: DM_Plex *mesh = (DM_Plex *) dm->data;
1370: const char *name = "Jacobian";
1371: DM dmAux, plex, plexAux = NULL;
1372: DMEnclosureType encAux;
1373: Vec A;
1374: PetscDS prob, probAux = NULL;
1375: PetscQuadrature quad;
1376: PetscSection section, globalSection, sectionAux;
1377: PetscScalar *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z;
1378: PetscInt Nf, fieldI, fieldJ;
1379: PetscInt totDim, totDimAux = 0;
1380: const PetscInt *cells;
1381: PetscInt cStart, cEnd, numCells, c;
1382: PetscBool hasDyn;
1383: DMField coordField;
1384: PetscHashFormKey key;
1385: PetscErrorCode ierr;
1388: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1389: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1390: if (!cellIS) {
1391: PetscInt depth;
1393: DMPlexGetDepth(plex, &depth);
1394: DMGetStratumIS(plex, "dim", depth, &cellIS);
1395: if (!cellIS) {DMGetStratumIS(plex, "depth", depth, &cellIS);}
1396: } else {
1397: PetscObjectReference((PetscObject) cellIS);
1398: }
1399: key.label = NULL;
1400: key.value = 0;
1401: DMGetLocalSection(dm, §ion);
1402: DMGetGlobalSection(dm, &globalSection);
1403: ISGetLocalSize(cellIS, &numCells);
1404: ISGetPointRange(cellIS, &cStart, &cEnd, &cells);
1405: DMGetCellDS(dm, cells ? cells[cStart] : cStart, &prob);
1406: PetscDSGetNumFields(prob, &Nf);
1407: PetscDSGetTotalDimension(prob, &totDim);
1408: PetscDSHasDynamicJacobian(prob, &hasDyn);
1409: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
1410: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1411: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1412: if (dmAux) {
1413: DMGetEnclosureRelation(dmAux, dm, &encAux);
1414: DMConvert(dmAux, DMPLEX, &plexAux);
1415: DMGetLocalSection(plexAux, §ionAux);
1416: DMGetDS(dmAux, &probAux);
1417: PetscDSGetTotalDimension(probAux, &totDimAux);
1418: }
1419: VecSet(Z, 0.0);
1420: PetscMalloc6(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasDyn ? numCells*totDim*totDim : 0, &elemMatD,numCells*totDim,&y,totDim,&z);
1421: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1422: DMGetCoordinateField(dm, &coordField);
1423: for (c = cStart; c < cEnd; ++c) {
1424: const PetscInt cell = cells ? cells[c] : c;
1425: const PetscInt cind = c - cStart;
1426: PetscScalar *x = NULL, *x_t = NULL;
1427: PetscInt i;
1429: DMPlexVecGetClosure(dm, section, X, cell, NULL, &x);
1430: for (i = 0; i < totDim; ++i) u[cind*totDim+i] = x[i];
1431: DMPlexVecRestoreClosure(dm, section, X, cell, NULL, &x);
1432: if (X_t) {
1433: DMPlexVecGetClosure(dm, section, X_t, cell, NULL, &x_t);
1434: for (i = 0; i < totDim; ++i) u_t[cind*totDim+i] = x_t[i];
1435: DMPlexVecRestoreClosure(dm, section, X_t, cell, NULL, &x_t);
1436: }
1437: if (dmAux) {
1438: PetscInt subcell;
1439: DMGetEnclosurePoint(dmAux, dm, encAux, cell, &subcell);
1440: DMPlexVecGetClosure(plexAux, sectionAux, A, subcell, NULL, &x);
1441: for (i = 0; i < totDimAux; ++i) a[cind*totDimAux+i] = x[i];
1442: DMPlexVecRestoreClosure(plexAux, sectionAux, A, subcell, NULL, &x);
1443: }
1444: DMPlexVecGetClosure(dm, section, Y, cell, NULL, &x);
1445: for (i = 0; i < totDim; ++i) y[cind*totDim+i] = x[i];
1446: DMPlexVecRestoreClosure(dm, section, Y, cell, NULL, &x);
1447: }
1448: PetscArrayzero(elemMat, numCells*totDim*totDim);
1449: if (hasDyn) {PetscArrayzero(elemMatD, numCells*totDim*totDim);}
1450: for (fieldI = 0; fieldI < Nf; ++fieldI) {
1451: PetscFE fe;
1452: PetscInt Nb;
1453: /* Conforming batches */
1454: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1455: /* Remainder */
1456: PetscInt Nr, offset, Nq;
1457: PetscQuadrature qGeom = NULL;
1458: PetscInt maxDegree;
1459: PetscFEGeom *cgeomFEM, *chunkGeom = NULL, *remGeom = NULL;
1461: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
1462: PetscFEGetQuadrature(fe, &quad);
1463: PetscFEGetDimension(fe, &Nb);
1464: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1465: DMFieldGetDegree(coordField,cellIS,NULL,&maxDegree);
1466: if (maxDegree <= 1) {DMFieldCreateDefaultQuadrature(coordField,cellIS,&qGeom);}
1467: if (!qGeom) {
1468: PetscFEGetQuadrature(fe,&qGeom);
1469: PetscObjectReference((PetscObject)qGeom);
1470: }
1471: PetscQuadratureGetData(qGeom, NULL, NULL, &Nq, NULL, NULL);
1472: DMSNESGetFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
1473: blockSize = Nb;
1474: batchSize = numBlocks * blockSize;
1475: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1476: numChunks = numCells / (numBatches*batchSize);
1477: Ne = numChunks*numBatches*batchSize;
1478: Nr = numCells % (numBatches*batchSize);
1479: offset = numCells - Nr;
1480: PetscFEGeomGetChunk(cgeomFEM,0,offset,&chunkGeom);
1481: PetscFEGeomGetChunk(cgeomFEM,offset,numCells,&remGeom);
1482: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
1483: key.field = fieldI*Nf + fieldJ;
1484: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, key, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMat);
1485: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN, key, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
1486: if (hasDyn) {
1487: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, key, Ne, chunkGeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
1488: PetscFEIntegrateJacobian(prob, PETSCFE_JACOBIAN_DYN, key, Nr, remGeom, &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMatD[offset*totDim*totDim]);
1489: }
1490: }
1491: PetscFEGeomRestoreChunk(cgeomFEM,offset,numCells,&remGeom);
1492: PetscFEGeomRestoreChunk(cgeomFEM,0,offset,&chunkGeom);
1493: DMSNESRestoreFEGeom(coordField,cellIS,qGeom,PETSC_FALSE,&cgeomFEM);
1494: PetscQuadratureDestroy(&qGeom);
1495: }
1496: if (hasDyn) {
1497: for (c = 0; c < numCells*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
1498: }
1499: for (c = cStart; c < cEnd; ++c) {
1500: const PetscInt cell = cells ? cells[c] : c;
1501: const PetscInt cind = c - cStart;
1502: const PetscBLASInt M = totDim, one = 1;
1503: const PetscScalar a = 1.0, b = 0.0;
1505: PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[cind*totDim*totDim], &M, &y[cind*totDim], &one, &b, z, &one));
1506: if (mesh->printFEM > 1) {
1507: DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[cind*totDim*totDim]);
1508: DMPrintCellVector(c, "Y", totDim, &y[cind*totDim]);
1509: DMPrintCellVector(c, "Z", totDim, z);
1510: }
1511: DMPlexVecSetClosure(dm, section, Z, cell, z, ADD_VALUES);
1512: }
1513: PetscFree6(u,u_t,elemMat,elemMatD,y,z);
1514: if (mesh->printFEM) {
1515: PetscPrintf(PetscObjectComm((PetscObject)Z), "Z:\n");
1516: VecView(Z, NULL);
1517: }
1518: PetscFree(a);
1519: ISDestroy(&cellIS);
1520: DMDestroy(&plexAux);
1521: DMDestroy(&plex);
1522: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
1523: return(0);
1524: }
1526: /*@
1527: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
1529: Input Parameters:
1530: + dm - The mesh
1531: . X - Local input vector
1532: - user - The user context
1534: Output Parameter:
1535: . Jac - Jacobian matrix
1537: Note:
1538: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
1539: like a GPU, or vectorize on a multicore machine.
1541: Level: developer
1543: .seealso: FormFunctionLocal()
1544: @*/
1545: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
1546: {
1547: DM plex;
1548: IS allcellIS;
1549: PetscBool hasJac, hasPrec;
1550: PetscInt Nds, s;
1554: DMSNESConvertPlex(dm, &plex, PETSC_TRUE);
1555: DMPlexGetAllCells_Internal(plex, &allcellIS);
1556: DMGetNumDS(dm, &Nds);
1557: for (s = 0; s < Nds; ++s) {
1558: PetscDS ds;
1559: IS cellIS;
1560: PetscHashFormKey key;
1562: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
1563: key.value = 0;
1564: key.field = 0;
1565: if (!key.label) {
1566: PetscObjectReference((PetscObject) allcellIS);
1567: cellIS = allcellIS;
1568: } else {
1569: IS pointIS;
1571: key.value = 1;
1572: DMLabelGetStratumIS(key.label, key.value, &pointIS);
1573: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
1574: ISDestroy(&pointIS);
1575: }
1576: if (!s) {
1577: PetscDSHasJacobian(ds, &hasJac);
1578: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1579: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
1580: MatZeroEntries(JacP);
1581: }
1582: DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user);
1583: ISDestroy(&cellIS);
1584: }
1585: ISDestroy(&allcellIS);
1586: DMDestroy(&plex);
1587: return(0);
1588: }
1590: /*
1591: MatComputeNeumannOverlap - Computes an unassembled (Neumann) local overlapping Mat in nonlinear context.
1593: Input Parameters:
1594: + X - SNES linearization point
1595: . ovl - index set of overlapping subdomains
1597: Output Parameter:
1598: . J - unassembled (Neumann) local matrix
1600: Level: intermediate
1602: .seealso: DMCreateNeumannOverlap(), MATIS, PCHPDDMSetAuxiliaryMat()
1603: */
1604: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *ctx)
1605: {
1606: SNES snes;
1607: Mat pJ;
1608: DM ovldm,origdm;
1609: DMSNES sdm;
1610: PetscErrorCode (*bfun)(DM,Vec,void*);
1611: PetscErrorCode (*jfun)(DM,Vec,Mat,Mat,void*);
1612: void *bctx,*jctx;
1616: PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_MATIS",(PetscObject*)&pJ);
1617: if (!pJ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing overlapping Mat");
1618: PetscObjectQuery((PetscObject)ovl,"_DM_Original_HPDDM",(PetscObject*)&origdm);
1619: if (!origdm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing original DM");
1620: MatGetDM(pJ,&ovldm);
1621: DMSNESGetBoundaryLocal(origdm,&bfun,&bctx);
1622: DMSNESSetBoundaryLocal(ovldm,bfun,bctx);
1623: DMSNESGetJacobianLocal(origdm,&jfun,&jctx);
1624: DMSNESSetJacobianLocal(ovldm,jfun,jctx);
1625: PetscObjectQuery((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject*)&snes);
1626: if (!snes) {
1627: SNESCreate(PetscObjectComm((PetscObject)ovl),&snes);
1628: SNESSetDM(snes,ovldm);
1629: PetscObjectCompose((PetscObject)ovl,"_DM_Overlap_HPDDM_SNES",(PetscObject)snes);
1630: PetscObjectDereference((PetscObject)snes);
1631: }
1632: DMGetDMSNES(ovldm,&sdm);
1633: VecLockReadPush(X);
1634: PetscStackPush("SNES user Jacobian function");
1635: (*sdm->ops->computejacobian)(snes,X,pJ,pJ,sdm->jacobianctx);
1636: PetscStackPop;
1637: VecLockReadPop(X);
1638: /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
1639: {
1640: Mat locpJ;
1642: MatISGetLocalMat(pJ,&locpJ);
1643: MatCopy(locpJ,J,SAME_NONZERO_PATTERN);
1644: }
1645: return(0);
1646: }
1648: /*@
1649: DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.
1651: Input Parameters:
1652: + dm - The DM object
1653: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
1654: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
1655: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())
1657: Level: developer
1658: @*/
1659: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
1660: {
1664: DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
1665: DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
1666: DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
1667: PetscObjectComposeFunction((PetscObject)dm,"MatComputeNeumannOverlap_C",MatComputeNeumannOverlap_Plex);
1668: return(0);
1669: }
1671: /*@C
1672: DMSNESCheckDiscretization - Check the discretization error of the exact solution
1674: Input Parameters:
1675: + snes - the SNES object
1676: . dm - the DM
1677: . t - the time
1678: . u - a DM vector
1679: - tol - A tolerance for the check, or -1 to print the results instead
1681: Output Parameters:
1682: . error - An array which holds the discretization error in each field, or NULL
1684: Note: The user must call PetscDSSetExactSolution() beforehand
1686: Level: developer
1688: .seealso: DNSNESCheckFromOptions(), DMSNESCheckResidual(), DMSNESCheckJacobian(), PetscDSSetExactSolution()
1689: @*/
1690: PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
1691: {
1692: PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx);
1693: void **ectxs;
1694: PetscReal *err;
1695: MPI_Comm comm;
1696: PetscInt Nf, f;
1697: PetscErrorCode ierr;
1705: DMComputeExactSolution(dm, t, u, NULL);
1706: VecViewFromOptions(u, NULL, "-vec_view");
1708: PetscObjectGetComm((PetscObject) snes, &comm);
1709: DMGetNumFields(dm, &Nf);
1710: PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err);
1711: {
1712: PetscInt Nds, s;
1714: DMGetNumDS(dm, &Nds);
1715: for (s = 0; s < Nds; ++s) {
1716: PetscDS ds;
1717: DMLabel label;
1718: IS fieldIS;
1719: const PetscInt *fields;
1720: PetscInt dsNf, f;
1722: DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds);
1723: PetscDSGetNumFields(ds, &dsNf);
1724: ISGetIndices(fieldIS, &fields);
1725: for (f = 0; f < dsNf; ++f) {
1726: const PetscInt field = fields[f];
1727: PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]);
1728: }
1729: ISRestoreIndices(fieldIS, &fields);
1730: }
1731: }
1732: if (Nf > 1) {
1733: DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err);
1734: if (tol >= 0.0) {
1735: for (f = 0; f < Nf; ++f) {
1736: if (err[f] > tol) SETERRQ3(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %D exceeds tolerance %g", (double) err[f], f, (double) tol);
1737: }
1738: } else if (error) {
1739: for (f = 0; f < Nf; ++f) error[f] = err[f];
1740: } else {
1741: PetscPrintf(comm, "L_2 Error: [");
1742: for (f = 0; f < Nf; ++f) {
1743: if (f) {PetscPrintf(comm, ", ");}
1744: PetscPrintf(comm, "%g", (double)err[f]);
1745: }
1746: PetscPrintf(comm, "]\n");
1747: }
1748: } else {
1749: DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]);
1750: if (tol >= 0.0) {
1751: if (err[0] > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double) err[0], (double) tol);
1752: } else if (error) {
1753: error[0] = err[0];
1754: } else {
1755: PetscPrintf(comm, "L_2 Error: %g\n", (double) err[0]);
1756: }
1757: }
1758: PetscFree3(exacts, ectxs, err);
1759: return(0);
1760: }
1762: /*@C
1763: DMSNESCheckResidual - Check the residual of the exact solution
1765: Input Parameters:
1766: + snes - the SNES object
1767: . dm - the DM
1768: . u - a DM vector
1769: - tol - A tolerance for the check, or -1 to print the results instead
1771: Output Parameters:
1772: . residual - The residual norm of the exact solution, or NULL
1774: Level: developer
1776: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
1777: @*/
1778: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
1779: {
1780: MPI_Comm comm;
1781: Vec r;
1782: PetscReal res;
1790: PetscObjectGetComm((PetscObject) snes, &comm);
1791: DMComputeExactSolution(dm, 0.0, u, NULL);
1792: VecDuplicate(u, &r);
1793: SNESComputeFunction(snes, u, r);
1794: VecNorm(r, NORM_2, &res);
1795: if (tol >= 0.0) {
1796: if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
1797: } else if (residual) {
1798: *residual = res;
1799: } else {
1800: PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
1801: VecChop(r, 1.0e-10);
1802: PetscObjectSetName((PetscObject) r, "Initial Residual");
1803: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
1804: VecViewFromOptions(r, NULL, "-vec_view");
1805: }
1806: VecDestroy(&r);
1807: return(0);
1808: }
1810: /*@C
1811: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
1813: Input Parameters:
1814: + snes - the SNES object
1815: . dm - the DM
1816: . u - a DM vector
1817: - tol - A tolerance for the check, or -1 to print the results instead
1819: Output Parameters:
1820: + isLinear - Flag indicaing that the function looks linear, or NULL
1821: - convRate - The rate of convergence of the linear model, or NULL
1823: Level: developer
1825: .seealso: DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
1826: @*/
1827: PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
1828: {
1829: MPI_Comm comm;
1830: PetscDS ds;
1831: Mat J, M;
1832: MatNullSpace nullspace;
1833: PetscReal slope, intercept;
1834: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
1843: PetscObjectGetComm((PetscObject) snes, &comm);
1844: DMComputeExactSolution(dm, 0.0, u, NULL);
1845: /* Create and view matrices */
1846: DMCreateMatrix(dm, &J);
1847: DMGetDS(dm, &ds);
1848: PetscDSHasJacobian(ds, &hasJac);
1849: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
1850: if (hasJac && hasPrec) {
1851: DMCreateMatrix(dm, &M);
1852: SNESComputeJacobian(snes, u, J, M);
1853: PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
1854: PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
1855: MatViewFromOptions(M, NULL, "-mat_view");
1856: MatDestroy(&M);
1857: } else {
1858: SNESComputeJacobian(snes, u, J, J);
1859: }
1860: PetscObjectSetName((PetscObject) J, "Jacobian");
1861: PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
1862: MatViewFromOptions(J, NULL, "-mat_view");
1863: /* Check nullspace */
1864: MatGetNullSpace(J, &nullspace);
1865: if (nullspace) {
1866: PetscBool isNull;
1867: MatNullSpaceTest(nullspace, J, &isNull);
1868: if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
1869: }
1870: /* Taylor test */
1871: {
1872: PetscRandom rand;
1873: Vec du, uhat, r, rhat, df;
1874: PetscReal h;
1875: PetscReal *es, *hs, *errors;
1876: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
1877: PetscInt Nv, v;
1879: /* Choose a perturbation direction */
1880: PetscRandomCreate(comm, &rand);
1881: VecDuplicate(u, &du);
1882: VecSetRandom(du, rand);
1883: PetscRandomDestroy(&rand);
1884: VecDuplicate(u, &df);
1885: MatMult(J, du, df);
1886: /* Evaluate residual at u, F(u), save in vector r */
1887: VecDuplicate(u, &r);
1888: SNESComputeFunction(snes, u, r);
1889: /* Look at the convergence of our Taylor approximation as we approach u */
1890: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
1891: PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
1892: VecDuplicate(u, &uhat);
1893: VecDuplicate(u, &rhat);
1894: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
1895: VecWAXPY(uhat, h, du, u);
1896: /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
1897: SNESComputeFunction(snes, uhat, rhat);
1898: VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
1899: VecNorm(rhat, NORM_2, &errors[Nv]);
1901: es[Nv] = PetscLog10Real(errors[Nv]);
1902: hs[Nv] = PetscLog10Real(h);
1903: }
1904: VecDestroy(&uhat);
1905: VecDestroy(&rhat);
1906: VecDestroy(&df);
1907: VecDestroy(&r);
1908: VecDestroy(&du);
1909: for (v = 0; v < Nv; ++v) {
1910: if ((tol >= 0) && (errors[v] > tol)) break;
1911: else if (errors[v] > PETSC_SMALL) break;
1912: }
1913: if (v == Nv) isLin = PETSC_TRUE;
1914: PetscLinearRegression(Nv, hs, es, &slope, &intercept);
1915: PetscFree3(es, hs, errors);
1916: /* Slope should be about 2 */
1917: if (tol >= 0) {
1918: if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
1919: } else if (isLinear || convRate) {
1920: if (isLinear) *isLinear = isLin;
1921: if (convRate) *convRate = slope;
1922: } else {
1923: if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
1924: else {PetscPrintf(comm, "Function appears to be linear\n");}
1925: }
1926: }
1927: MatDestroy(&J);
1928: return(0);
1929: }
1931: PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1932: {
1936: DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL);
1937: DMSNESCheckResidual(snes, dm, u, -1.0, NULL);
1938: DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL);
1939: return(0);
1940: }
1942: /*@C
1943: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1945: Input Parameters:
1946: + snes - the SNES object
1947: - u - representative SNES vector
1949: Note: The user must call PetscDSSetExactSolution() beforehand
1951: Level: developer
1952: @*/
1953: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1954: {
1955: DM dm;
1956: Vec sol;
1957: PetscBool check;
1961: PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
1962: if (!check) return(0);
1963: SNESGetDM(snes, &dm);
1964: VecDuplicate(u, &sol);
1965: SNESSetSolution(snes, sol);
1966: DMSNESCheck_Internal(snes, dm, sol);
1967: VecDestroy(&sol);
1968: return(0);
1969: }