Actual source code: dmplexsnes.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3: #include <petscds.h>
4: #include <petsc/private/petscimpl.h>
5: #include <petsc/private/petscfeimpl.h>
7: /************************** Interpolation *******************************/
11: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
12: {
13: PetscBool isPlex;
17: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
18: if (isPlex) {
19: *plex = dm;
20: PetscObjectReference((PetscObject) dm);
21: } else {
22: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
23: if (!*plex) {
24: DMConvert(dm,DMPLEX,plex);
25: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
26: if (copy) {
27: PetscInt i;
28: PetscObject obj;
29: const char *comps[3] = {"A","dmAux","dmCh"};
31: DMCopyDMSNES(dm, *plex);
32: for (i = 0; i < 3; i++) {
33: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
34: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
35: }
36: }
37: } else {
38: PetscObjectReference((PetscObject) *plex);
39: }
40: }
41: return(0);
42: }
46: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
47: {
52: PetscMalloc(sizeof(struct _DMInterpolationInfo), ctx);
54: (*ctx)->comm = comm;
55: (*ctx)->dim = -1;
56: (*ctx)->nInput = 0;
57: (*ctx)->points = NULL;
58: (*ctx)->cells = NULL;
59: (*ctx)->n = -1;
60: (*ctx)->coords = NULL;
61: return(0);
62: }
66: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
67: {
69: if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
70: ctx->dim = dim;
71: return(0);
72: }
76: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
77: {
80: *dim = ctx->dim;
81: return(0);
82: }
86: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
87: {
89: if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
90: ctx->dof = dof;
91: return(0);
92: }
96: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
97: {
100: *dof = ctx->dof;
101: return(0);
102: }
106: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
107: {
111: if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
112: if (ctx->points) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
113: ctx->nInput = n;
115: PetscMalloc1(n*ctx->dim, &ctx->points);
116: PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
117: return(0);
118: }
122: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
123: {
124: MPI_Comm comm = ctx->comm;
125: PetscScalar *a;
126: PetscInt p, q, i;
127: PetscMPIInt rank, size;
128: PetscErrorCode ierr;
129: Vec pointVec;
130: PetscSF cellSF;
131: PetscLayout layout;
132: PetscReal *globalPoints;
133: PetscScalar *globalPointsScalar;
134: const PetscInt *ranges;
135: PetscMPIInt *counts, *displs;
136: const PetscSFNode *foundCells;
137: const PetscInt *foundPoints;
138: PetscMPIInt *foundProcs, *globalProcs;
139: PetscInt n, N, numFound;
143: MPI_Comm_size(comm, &size);
144: MPI_Comm_rank(comm, &rank);
145: if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
146: /* Locate points */
147: n = ctx->nInput;
148: if (!redundantPoints) {
149: PetscLayoutCreate(comm, &layout);
150: PetscLayoutSetBlockSize(layout, 1);
151: PetscLayoutSetLocalSize(layout, n);
152: PetscLayoutSetUp(layout);
153: PetscLayoutGetSize(layout, &N);
154: /* Communicate all points to all processes */
155: PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
156: PetscLayoutGetRanges(layout, &ranges);
157: for (p = 0; p < size; ++p) {
158: counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
159: displs[p] = ranges[p]*ctx->dim;
160: }
161: MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
162: } else {
163: N = n;
164: globalPoints = ctx->points;
165: counts = displs = NULL;
166: layout = NULL;
167: }
168: #if 0
169: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
170: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
171: #else
172: #if defined(PETSC_USE_COMPLEX)
173: PetscMalloc1(N,&globalPointsScalar);
174: for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
175: #else
176: globalPointsScalar = globalPoints;
177: #endif
178: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
179: PetscMalloc2(N,&foundProcs,N,&globalProcs);
180: cellSF = NULL;
181: DMLocatePoints(dm, pointVec, &cellSF);
182: PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
183: #endif
184: for (p = 0; p < numFound; ++p) {
185: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
186: else foundProcs[foundPoints ? foundPoints[p] : p] = size;
187: }
188: /* Let the lowest rank process own each point */
189: MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
190: ctx->n = 0;
191: for (p = 0; p < N; ++p) {
192: if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, globalPoints[p*ctx->dim+0], ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0, ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0);
193: else if (globalProcs[p] == rank) ctx->n++;
194: }
195: /* Create coordinates vector and array of owned cells */
196: PetscMalloc1(ctx->n, &ctx->cells);
197: VecCreate(comm, &ctx->coords);
198: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
199: VecSetBlockSize(ctx->coords, ctx->dim);
200: VecSetType(ctx->coords,VECSTANDARD);
201: VecGetArray(ctx->coords, &a);
202: for (p = 0, q = 0, i = 0; p < N; ++p) {
203: if (globalProcs[p] == rank) {
204: PetscInt d;
206: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
207: ctx->cells[q++] = foundCells[p].index;
208: }
209: }
210: VecRestoreArray(ctx->coords, &a);
211: #if 0
212: PetscFree3(foundCells,foundProcs,globalProcs);
213: #else
214: PetscFree2(foundProcs,globalProcs);
215: PetscSFDestroy(&cellSF);
216: VecDestroy(&pointVec);
217: #endif
218: if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
219: if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
220: PetscLayoutDestroy(&layout);
221: return(0);
222: }
226: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
227: {
230: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
231: *coordinates = ctx->coords;
232: return(0);
233: }
237: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
238: {
243: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
244: VecCreate(ctx->comm, v);
245: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
246: VecSetBlockSize(*v, ctx->dof);
247: VecSetType(*v,VECSTANDARD);
248: return(0);
249: }
253: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
254: {
259: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
260: VecDestroy(v);
261: return(0);
262: }
266: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
267: {
268: PetscReal *v0, *J, *invJ, detJ;
269: const PetscScalar *coords;
270: PetscScalar *a;
271: PetscInt p;
275: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
276: VecGetArrayRead(ctx->coords, &coords);
277: VecGetArray(v, &a);
278: for (p = 0; p < ctx->n; ++p) {
279: PetscInt c = ctx->cells[p];
280: PetscScalar *x = NULL;
281: PetscReal xi[4];
282: PetscInt d, f, comp;
284: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
285: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
286: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
287: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
289: for (d = 0; d < ctx->dim; ++d) {
290: xi[d] = 0.0;
291: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
292: 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];
293: }
294: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
295: }
296: VecRestoreArray(v, &a);
297: VecRestoreArrayRead(ctx->coords, &coords);
298: PetscFree3(v0, J, invJ);
299: return(0);
300: }
304: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
305: {
306: PetscReal *v0, *J, *invJ, detJ;
307: const PetscScalar *coords;
308: PetscScalar *a;
309: PetscInt p;
313: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
314: VecGetArrayRead(ctx->coords, &coords);
315: VecGetArray(v, &a);
316: for (p = 0; p < ctx->n; ++p) {
317: PetscInt c = ctx->cells[p];
318: const PetscInt order[3] = {2, 1, 3};
319: PetscScalar *x = NULL;
320: PetscReal xi[4];
321: PetscInt d, f, comp;
323: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
324: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
325: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
326: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
328: for (d = 0; d < ctx->dim; ++d) {
329: xi[d] = 0.0;
330: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
331: 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];
332: }
333: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
334: }
335: VecRestoreArray(v, &a);
336: VecRestoreArrayRead(ctx->coords, &coords);
337: PetscFree3(v0, J, invJ);
338: return(0);
339: }
343: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
344: {
345: const PetscScalar *vertices = (const PetscScalar*) ctx;
346: const PetscScalar x0 = vertices[0];
347: const PetscScalar y0 = vertices[1];
348: const PetscScalar x1 = vertices[2];
349: const PetscScalar y1 = vertices[3];
350: const PetscScalar x2 = vertices[4];
351: const PetscScalar y2 = vertices[5];
352: const PetscScalar x3 = vertices[6];
353: const PetscScalar y3 = vertices[7];
354: const PetscScalar f_1 = x1 - x0;
355: const PetscScalar g_1 = y1 - y0;
356: const PetscScalar f_3 = x3 - x0;
357: const PetscScalar g_3 = y3 - y0;
358: const PetscScalar f_01 = x2 - x1 - x3 + x0;
359: const PetscScalar g_01 = y2 - y1 - y3 + y0;
360: const PetscScalar *ref;
361: PetscScalar *real;
362: PetscErrorCode ierr;
365: VecGetArrayRead(Xref, &ref);
366: VecGetArray(Xreal, &real);
367: {
368: const PetscScalar p0 = ref[0];
369: const PetscScalar p1 = ref[1];
371: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
372: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
373: }
374: PetscLogFlops(28);
375: VecRestoreArrayRead(Xref, &ref);
376: VecRestoreArray(Xreal, &real);
377: return(0);
378: }
380: #include <petsc/private/dmimpl.h>
383: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
384: {
385: const PetscScalar *vertices = (const PetscScalar*) ctx;
386: const PetscScalar x0 = vertices[0];
387: const PetscScalar y0 = vertices[1];
388: const PetscScalar x1 = vertices[2];
389: const PetscScalar y1 = vertices[3];
390: const PetscScalar x2 = vertices[4];
391: const PetscScalar y2 = vertices[5];
392: const PetscScalar x3 = vertices[6];
393: const PetscScalar y3 = vertices[7];
394: const PetscScalar f_01 = x2 - x1 - x3 + x0;
395: const PetscScalar g_01 = y2 - y1 - y3 + y0;
396: const PetscScalar *ref;
397: PetscErrorCode ierr;
400: VecGetArrayRead(Xref, &ref);
401: {
402: const PetscScalar x = ref[0];
403: const PetscScalar y = ref[1];
404: const PetscInt rows[2] = {0, 1};
405: PetscScalar values[4];
407: values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
408: values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
409: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
410: }
411: PetscLogFlops(30);
412: VecRestoreArrayRead(Xref, &ref);
413: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
414: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
415: return(0);
416: }
420: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
421: {
422: DM dmCoord;
423: SNES snes;
424: KSP ksp;
425: PC pc;
426: Vec coordsLocal, r, ref, real;
427: Mat J;
428: const PetscScalar *coords;
429: PetscScalar *a;
430: PetscInt p;
434: DMGetCoordinatesLocal(dm, &coordsLocal);
435: DMGetCoordinateDM(dm, &dmCoord);
436: SNESCreate(PETSC_COMM_SELF, &snes);
437: SNESSetOptionsPrefix(snes, "quad_interp_");
438: VecCreate(PETSC_COMM_SELF, &r);
439: VecSetSizes(r, 2, 2);
440: VecSetType(r,dm->vectype);
441: VecDuplicate(r, &ref);
442: VecDuplicate(r, &real);
443: MatCreate(PETSC_COMM_SELF, &J);
444: MatSetSizes(J, 2, 2, 2, 2);
445: MatSetType(J, MATSEQDENSE);
446: MatSetUp(J);
447: SNESSetFunction(snes, r, QuadMap_Private, NULL);
448: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
449: SNESGetKSP(snes, &ksp);
450: KSPGetPC(ksp, &pc);
451: PCSetType(pc, PCLU);
452: SNESSetFromOptions(snes);
454: VecGetArrayRead(ctx->coords, &coords);
455: VecGetArray(v, &a);
456: for (p = 0; p < ctx->n; ++p) {
457: PetscScalar *x = NULL, *vertices = NULL;
458: PetscScalar *xi;
459: PetscReal xir[2];
460: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
462: /* Can make this do all points at once */
463: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
464: if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
465: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
466: if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
467: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
468: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
469: VecGetArray(real, &xi);
470: xi[0] = coords[p*ctx->dim+0];
471: xi[1] = coords[p*ctx->dim+1];
472: VecRestoreArray(real, &xi);
473: SNESSolve(snes, real, ref);
474: VecGetArray(ref, &xi);
475: xir[0] = PetscRealPart(xi[0]);
476: xir[1] = PetscRealPart(xi[1]);
477: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp]*(1 - xir[0])*(1 - xir[1]) + x[1*ctx->dof+comp]*xir[0]*(1 - xir[1]) + x[2*ctx->dof+comp]*xir[0]*xir[1] + x[3*ctx->dof+comp]*(1 - xir[0])*xir[1];
479: VecRestoreArray(ref, &xi);
480: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
481: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
482: }
483: VecRestoreArray(v, &a);
484: VecRestoreArrayRead(ctx->coords, &coords);
486: SNESDestroy(&snes);
487: VecDestroy(&r);
488: VecDestroy(&ref);
489: VecDestroy(&real);
490: MatDestroy(&J);
491: return(0);
492: }
496: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
497: {
498: const PetscScalar *vertices = (const PetscScalar*) ctx;
499: const PetscScalar x0 = vertices[0];
500: const PetscScalar y0 = vertices[1];
501: const PetscScalar z0 = vertices[2];
502: const PetscScalar x1 = vertices[9];
503: const PetscScalar y1 = vertices[10];
504: const PetscScalar z1 = vertices[11];
505: const PetscScalar x2 = vertices[6];
506: const PetscScalar y2 = vertices[7];
507: const PetscScalar z2 = vertices[8];
508: const PetscScalar x3 = vertices[3];
509: const PetscScalar y3 = vertices[4];
510: const PetscScalar z3 = vertices[5];
511: const PetscScalar x4 = vertices[12];
512: const PetscScalar y4 = vertices[13];
513: const PetscScalar z4 = vertices[14];
514: const PetscScalar x5 = vertices[15];
515: const PetscScalar y5 = vertices[16];
516: const PetscScalar z5 = vertices[17];
517: const PetscScalar x6 = vertices[18];
518: const PetscScalar y6 = vertices[19];
519: const PetscScalar z6 = vertices[20];
520: const PetscScalar x7 = vertices[21];
521: const PetscScalar y7 = vertices[22];
522: const PetscScalar z7 = vertices[23];
523: const PetscScalar f_1 = x1 - x0;
524: const PetscScalar g_1 = y1 - y0;
525: const PetscScalar h_1 = z1 - z0;
526: const PetscScalar f_3 = x3 - x0;
527: const PetscScalar g_3 = y3 - y0;
528: const PetscScalar h_3 = z3 - z0;
529: const PetscScalar f_4 = x4 - x0;
530: const PetscScalar g_4 = y4 - y0;
531: const PetscScalar h_4 = z4 - z0;
532: const PetscScalar f_01 = x2 - x1 - x3 + x0;
533: const PetscScalar g_01 = y2 - y1 - y3 + y0;
534: const PetscScalar h_01 = z2 - z1 - z3 + z0;
535: const PetscScalar f_12 = x7 - x3 - x4 + x0;
536: const PetscScalar g_12 = y7 - y3 - y4 + y0;
537: const PetscScalar h_12 = z7 - z3 - z4 + z0;
538: const PetscScalar f_02 = x5 - x1 - x4 + x0;
539: const PetscScalar g_02 = y5 - y1 - y4 + y0;
540: const PetscScalar h_02 = z5 - z1 - z4 + z0;
541: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
542: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
543: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
544: const PetscScalar *ref;
545: PetscScalar *real;
546: PetscErrorCode ierr;
549: VecGetArrayRead(Xref, &ref);
550: VecGetArray(Xreal, &real);
551: {
552: const PetscScalar p0 = ref[0];
553: const PetscScalar p1 = ref[1];
554: const PetscScalar p2 = ref[2];
556: 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;
557: 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;
558: 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;
559: }
560: PetscLogFlops(114);
561: VecRestoreArrayRead(Xref, &ref);
562: VecRestoreArray(Xreal, &real);
563: return(0);
564: }
568: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
569: {
570: const PetscScalar *vertices = (const PetscScalar*) ctx;
571: const PetscScalar x0 = vertices[0];
572: const PetscScalar y0 = vertices[1];
573: const PetscScalar z0 = vertices[2];
574: const PetscScalar x1 = vertices[9];
575: const PetscScalar y1 = vertices[10];
576: const PetscScalar z1 = vertices[11];
577: const PetscScalar x2 = vertices[6];
578: const PetscScalar y2 = vertices[7];
579: const PetscScalar z2 = vertices[8];
580: const PetscScalar x3 = vertices[3];
581: const PetscScalar y3 = vertices[4];
582: const PetscScalar z3 = vertices[5];
583: const PetscScalar x4 = vertices[12];
584: const PetscScalar y4 = vertices[13];
585: const PetscScalar z4 = vertices[14];
586: const PetscScalar x5 = vertices[15];
587: const PetscScalar y5 = vertices[16];
588: const PetscScalar z5 = vertices[17];
589: const PetscScalar x6 = vertices[18];
590: const PetscScalar y6 = vertices[19];
591: const PetscScalar z6 = vertices[20];
592: const PetscScalar x7 = vertices[21];
593: const PetscScalar y7 = vertices[22];
594: const PetscScalar z7 = vertices[23];
595: const PetscScalar f_xy = x2 - x1 - x3 + x0;
596: const PetscScalar g_xy = y2 - y1 - y3 + y0;
597: const PetscScalar h_xy = z2 - z1 - z3 + z0;
598: const PetscScalar f_yz = x7 - x3 - x4 + x0;
599: const PetscScalar g_yz = y7 - y3 - y4 + y0;
600: const PetscScalar h_yz = z7 - z3 - z4 + z0;
601: const PetscScalar f_xz = x5 - x1 - x4 + x0;
602: const PetscScalar g_xz = y5 - y1 - y4 + y0;
603: const PetscScalar h_xz = z5 - z1 - z4 + z0;
604: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
605: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
606: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
607: const PetscScalar *ref;
608: PetscErrorCode ierr;
611: VecGetArrayRead(Xref, &ref);
612: {
613: const PetscScalar x = ref[0];
614: const PetscScalar y = ref[1];
615: const PetscScalar z = ref[2];
616: const PetscInt rows[3] = {0, 1, 2};
617: PetscScalar values[9];
619: values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
620: values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
621: values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
622: values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
623: values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
624: values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
625: values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
626: values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
627: values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
629: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
630: }
631: PetscLogFlops(152);
632: VecRestoreArrayRead(Xref, &ref);
633: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
634: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
635: return(0);
636: }
640: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
641: {
642: DM dmCoord;
643: SNES snes;
644: KSP ksp;
645: PC pc;
646: Vec coordsLocal, r, ref, real;
647: Mat J;
648: const PetscScalar *coords;
649: PetscScalar *a;
650: PetscInt p;
654: DMGetCoordinatesLocal(dm, &coordsLocal);
655: DMGetCoordinateDM(dm, &dmCoord);
656: SNESCreate(PETSC_COMM_SELF, &snes);
657: SNESSetOptionsPrefix(snes, "hex_interp_");
658: VecCreate(PETSC_COMM_SELF, &r);
659: VecSetSizes(r, 3, 3);
660: VecSetType(r,dm->vectype);
661: VecDuplicate(r, &ref);
662: VecDuplicate(r, &real);
663: MatCreate(PETSC_COMM_SELF, &J);
664: MatSetSizes(J, 3, 3, 3, 3);
665: MatSetType(J, MATSEQDENSE);
666: MatSetUp(J);
667: SNESSetFunction(snes, r, HexMap_Private, NULL);
668: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
669: SNESGetKSP(snes, &ksp);
670: KSPGetPC(ksp, &pc);
671: PCSetType(pc, PCLU);
672: SNESSetFromOptions(snes);
674: VecGetArrayRead(ctx->coords, &coords);
675: VecGetArray(v, &a);
676: for (p = 0; p < ctx->n; ++p) {
677: PetscScalar *x = NULL, *vertices = NULL;
678: PetscScalar *xi;
679: PetscReal xir[3];
680: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
682: /* Can make this do all points at once */
683: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
684: if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
685: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
686: if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
687: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
688: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
689: VecGetArray(real, &xi);
690: xi[0] = coords[p*ctx->dim+0];
691: xi[1] = coords[p*ctx->dim+1];
692: xi[2] = coords[p*ctx->dim+2];
693: VecRestoreArray(real, &xi);
694: SNESSolve(snes, real, ref);
695: VecGetArray(ref, &xi);
696: xir[0] = PetscRealPart(xi[0]);
697: xir[1] = PetscRealPart(xi[1]);
698: xir[2] = PetscRealPart(xi[2]);
699: for (comp = 0; comp < ctx->dof; ++comp) {
700: a[p*ctx->dof+comp] =
701: x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
702: x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
703: x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
704: x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
705: x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
706: x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
707: x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
708: x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
709: }
710: VecRestoreArray(ref, &xi);
711: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
712: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
713: }
714: VecRestoreArray(v, &a);
715: VecRestoreArrayRead(ctx->coords, &coords);
717: SNESDestroy(&snes);
718: VecDestroy(&r);
719: VecDestroy(&ref);
720: VecDestroy(&real);
721: MatDestroy(&J);
722: return(0);
723: }
727: /*
728: Input Parameters:
729: + ctx - The DMInterpolationInfo context
730: . dm - The DM
731: - x - The local vector containing the field to be interpolated
733: Output Parameters:
734: . v - The vector containing the interpolated values
735: */
736: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
737: {
738: PetscInt dim, coneSize, n;
745: VecGetLocalSize(v, &n);
746: 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);
747: if (n) {
748: DMGetDimension(dm, &dim);
749: DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
750: if (dim == 2) {
751: if (coneSize == 3) {
752: DMInterpolate_Triangle_Private(ctx, dm, x, v);
753: } else if (coneSize == 4) {
754: DMInterpolate_Quad_Private(ctx, dm, x, v);
755: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
756: } else if (dim == 3) {
757: if (coneSize == 4) {
758: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
759: } else {
760: DMInterpolate_Hex_Private(ctx, dm, x, v);
761: }
762: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
763: }
764: return(0);
765: }
769: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
770: {
775: VecDestroy(&(*ctx)->coords);
776: PetscFree((*ctx)->points);
777: PetscFree((*ctx)->cells);
778: PetscFree(*ctx);
779: *ctx = NULL;
780: return(0);
781: }
785: /*@C
786: SNESMonitorFields - Monitors the residual for each field separately
788: Collective on SNES
790: Input Parameters:
791: + snes - the SNES context
792: . its - iteration number
793: . fgnorm - 2-norm of residual
794: - vf - PetscViewerAndFormat of type ASCII
796: Notes:
797: This routine prints the residual norm at each iteration.
799: Level: intermediate
801: .keywords: SNES, nonlinear, default, monitor, norm
802: .seealso: SNESMonitorSet(), SNESMonitorDefault()
803: @*/
804: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
805: {
806: PetscViewer viewer = vf->viewer;
807: Vec res;
808: DM dm;
809: PetscSection s;
810: const PetscScalar *r;
811: PetscReal *lnorms, *norms;
812: PetscInt numFields, f, pStart, pEnd, p;
813: PetscErrorCode ierr;
817: SNESGetFunction(snes, &res, 0, 0);
818: SNESGetDM(snes, &dm);
819: DMGetDefaultSection(dm, &s);
820: PetscSectionGetNumFields(s, &numFields);
821: PetscSectionGetChart(s, &pStart, &pEnd);
822: PetscCalloc2(numFields, &lnorms, numFields, &norms);
823: VecGetArrayRead(res, &r);
824: for (p = pStart; p < pEnd; ++p) {
825: for (f = 0; f < numFields; ++f) {
826: PetscInt fdof, foff, d;
828: PetscSectionGetFieldDof(s, p, f, &fdof);
829: PetscSectionGetFieldOffset(s, p, f, &foff);
830: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
831: }
832: }
833: VecRestoreArrayRead(res, &r);
834: MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
835: PetscViewerPushFormat(viewer,vf->format);
836: PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
837: PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
838: for (f = 0; f < numFields; ++f) {
839: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
840: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
841: }
842: PetscViewerASCIIPrintf(viewer, "]\n");
843: PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
844: PetscViewerPopFormat(viewer);
845: PetscFree2(lnorms, norms);
846: return(0);
847: }
849: /********************* Residual Computation **************************/
853: /*@
854: DMPlexSNESGetGeometryFEM - Return precomputed geometric data
856: Input Parameter:
857: . dm - The DM
859: Output Parameters:
860: . cellgeom - The values precomputed from cell geometry
862: Level: developer
864: .seealso: DMPlexSNESSetFunctionLocal()
865: @*/
866: PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
867: {
868: DMSNES dmsnes;
869: PetscObject obj;
874: DMGetDMSNES(dm, &dmsnes);
875: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);
876: if (!obj) {
877: Vec cellgeom;
879: DMPlexComputeGeometryFEM(dm, &cellgeom);
880: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);
881: VecDestroy(&cellgeom);
882: }
884: return(0);
885: }
889: /*@
890: DMPlexSNESGetGeometryFVM - Return precomputed geometric data
892: Input Parameter:
893: . dm - The DM
895: Output Parameters:
896: + facegeom - The values precomputed from face geometry
897: . cellgeom - The values precomputed from cell geometry
898: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
900: Level: developer
902: .seealso: DMPlexTSSetRHSFunctionLocal()
903: @*/
904: PetscErrorCode DMPlexSNESGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
905: {
906: DMSNES dmsnes;
907: PetscObject obj;
912: DMGetDMSNES(dm, &dmsnes);
913: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", &obj);
914: if (!obj) {
915: Vec cellgeom, facegeom;
917: DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);
918: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", (PetscObject) facegeom);
919: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fvm", (PetscObject) cellgeom);
920: VecDestroy(&facegeom);
921: VecDestroy(&cellgeom);
922: }
925: if (minRadius) {DMPlexGetMinRadius(dm, minRadius);}
926: return(0);
927: }
931: /*@
932: DMPlexSNESGetGradientDM - Return gradient data layout
934: Input Parameters:
935: + dm - The DM
936: - fv - The PetscFV
938: Output Parameter:
939: . dmGrad - The layout for gradient values
941: Level: developer
943: .seealso: DMPlexSNESGetGeometryFVM()
944: @*/
945: PetscErrorCode DMPlexSNESGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
946: {
947: DMSNES dmsnes;
948: PetscObject obj;
949: PetscBool computeGradients;
956: PetscFVGetComputeGradients(fv, &computeGradients);
957: if (!computeGradients) {*dmGrad = NULL; return(0);}
958: DMGetDMSNES(dm, &dmsnes);
959: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", &obj);
960: if (!obj) {
961: DM dmGrad;
962: Vec faceGeometry, cellGeometry;
964: DMPlexSNESGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);
965: DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);
966: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject) dmGrad);
967: DMDestroy(&dmGrad);
968: }
969: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject *) dmGrad);
970: return(0);
971: }
975: /*@C
976: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
978: Input Parameters:
979: + dm - The DM
980: . cStart - The first cell to include
981: . cEnd - The first cell to exclude
982: . locX - A local vector with the solution fields
983: . locX_t - A local vector with solution field time derivatives, or NULL
984: - locA - A local vector with auxiliary fields, or NULL
986: Output Parameters:
987: + u - The field coefficients
988: . u_t - The fields derivative coefficients
989: - a - The auxiliary field coefficients
991: Level: developer
993: .seealso: DMPlexGetFaceFields()
994: @*/
995: PetscErrorCode DMPlexGetCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
996: {
997: DM dmAux;
998: PetscSection section, sectionAux;
999: PetscDS prob;
1000: PetscInt numCells = cEnd - cStart, totDim, totDimAux, c;
1011: DMGetDefaultSection(dm, §ion);
1012: DMGetDS(dm, &prob);
1013: PetscDSGetTotalDimension(prob, &totDim);
1014: if (locA) {
1015: PetscDS probAux;
1017: VecGetDM(locA, &dmAux);
1018: DMGetDefaultSection(dmAux, §ionAux);
1019: DMGetDS(dmAux, &probAux);
1020: PetscDSGetTotalDimension(probAux, &totDimAux);
1021: }
1022: DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u);
1023: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u_t);} else {*u_t = NULL;}
1024: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, PETSC_SCALAR, a);} else {*a = NULL;}
1025: for (c = cStart; c < cEnd; ++c) {
1026: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
1027: PetscInt i;
1029: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
1030: for (i = 0; i < totDim; ++i) ul[(c-cStart)*totDim+i] = x[i];
1031: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
1032: if (locX_t) {
1033: DMPlexVecGetClosure(dm, section, locX_t, c, NULL, &x_t);
1034: for (i = 0; i < totDim; ++i) ul_t[(c-cStart)*totDim+i] = x_t[i];
1035: DMPlexVecRestoreClosure(dm, section, locX_t, c, NULL, &x_t);
1036: }
1037: if (locA) {
1038: DM dmAuxPlex;
1040: DMSNESConvertPlex(dmAux, &dmAuxPlex, PETSC_FALSE);
1041: DMPlexVecGetClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
1042: for (i = 0; i < totDimAux; ++i) al[(c-cStart)*totDimAux+i] = x[i];
1043: DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
1044: DMDestroy(&dmAuxPlex);
1045: }
1046: }
1047: return(0);
1048: }
1052: /*@C
1053: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
1055: Input Parameters:
1056: + dm - The DM
1057: . cStart - The first cell to include
1058: . cEnd - The first cell to exclude
1059: . locX - A local vector with the solution fields
1060: . locX_t - A local vector with solution field time derivatives, or NULL
1061: - locA - A local vector with auxiliary fields, or NULL
1063: Output Parameters:
1064: + u - The field coefficients
1065: . u_t - The fields derivative coefficients
1066: - a - The auxiliary field coefficients
1068: Level: developer
1070: .seealso: DMPlexGetFaceFields()
1071: @*/
1072: PetscErrorCode DMPlexRestoreCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
1073: {
1077: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u);
1078: if (*u_t) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u_t);}
1079: if (*a) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, a);}
1080: return(0);
1081: }
1085: /*@C
1086: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
1088: Input Parameters:
1089: + dm - The DM
1090: . fStart - The first face to include
1091: . fEnd - The first face to exclude
1092: . locX - A local vector with the solution fields
1093: . locX_t - A local vector with solution field time derivatives, or NULL
1094: . faceGeometry - A local vector with face geometry
1095: . cellGeometry - A local vector with cell geometry
1096: - locaGrad - A local vector with field gradients, or NULL
1098: Output Parameters:
1099: + uL - The field values at the left side of the face
1100: - uR - The field values at the right side of the face
1102: Level: developer
1104: .seealso: DMPlexGetCellFields()
1105: @*/
1106: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1107: {
1108: DM dmFace, dmCell, dmGrad = NULL;
1109: PetscSection section;
1110: PetscDS prob;
1111: DMLabel ghostLabel;
1112: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
1113: PetscBool *isFE;
1114: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
1115: PetscErrorCode ierr;
1126: DMGetDimension(dm, &dim);
1127: DMGetDS(dm, &prob);
1128: DMGetDefaultSection(dm, §ion);
1129: PetscDSGetNumFields(prob, &Nf);
1130: PetscDSGetTotalComponents(prob, &Nc);
1131: PetscMalloc1(Nf, &isFE);
1132: for (f = 0; f < Nf; ++f) {
1133: PetscObject obj;
1134: PetscClassId id;
1136: DMGetField(dm, f, &obj);
1137: PetscObjectGetClassId(obj, &id);
1138: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
1139: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
1140: else {isFE[f] = PETSC_FALSE;}
1141: }
1142: DMGetLabel(dm, "ghost", &ghostLabel);
1143: VecGetArrayRead(locX, &x);
1144: VecGetDM(faceGeometry, &dmFace);
1145: VecGetArrayRead(faceGeometry, &facegeom);
1146: VecGetDM(cellGeometry, &dmCell);
1147: VecGetArrayRead(cellGeometry, &cellgeom);
1148: if (locGrad) {
1149: VecGetDM(locGrad, &dmGrad);
1150: VecGetArrayRead(locGrad, &lgrad);
1151: }
1152: DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uL);
1153: DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uR);
1154: /* Right now just eat the extra work for FE (could make a cell loop) */
1155: for (face = fStart, iface = 0; face < fEnd; ++face) {
1156: const PetscInt *cells;
1157: PetscFVFaceGeom *fg;
1158: PetscFVCellGeom *cgL, *cgR;
1159: PetscScalar *xL, *xR, *gL, *gR;
1160: PetscScalar *uLl = *uL, *uRl = *uR;
1161: PetscInt ghost, nsupp;
1163: DMLabelGetValue(ghostLabel, face, &ghost);
1164: DMPlexGetSupportSize(dm, face, &nsupp);
1165: if (ghost >= 0) continue;
1166: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1167: DMPlexGetSupport(dm, face, &cells);
1168: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1169: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1170: for (f = 0; f < Nf; ++f) {
1171: PetscInt off;
1173: PetscDSGetComponentOffset(prob, f, &off);
1174: if (isFE[f]) {
1175: const PetscInt *cone;
1176: PetscInt comp, coneSize, faceLocL, faceLocR, ldof, rdof, d;
1178: xL = xR = NULL;
1179: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1180: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1181: DMPlexGetCone(dm, cells[0], &cone);
1182: DMPlexGetConeSize(dm, cells[0], &coneSize);
1183: for (faceLocL = 0; faceLocL < coneSize; ++faceLocL) if (cone[faceLocL] == face) break;
1184: if (faceLocL == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[0]);
1185: DMPlexGetCone(dm, cells[1], &cone);
1186: DMPlexGetConeSize(dm, cells[1], &coneSize);
1187: for (faceLocR = 0; faceLocR < coneSize; ++faceLocR) if (cone[faceLocR] == face) break;
1188: if (faceLocR == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[1]);
1189: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
1190: EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
1191: if (rdof == ldof) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
1192: else {PetscSectionGetFieldComponents(section, f, &comp); for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
1193: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1194: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1195: } else {
1196: PetscFV fv;
1197: PetscInt numComp, c;
1199: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
1200: PetscFVGetNumComponents(fv, &numComp);
1201: if (nsupp > 2) {
1202: for (f = 0; f < Nf; ++f) {
1203: PetscInt off;
1205: PetscDSGetComponentOffset(prob, f, &off);
1206: PetscFVGetNumComponents(fv, &numComp);
1207: for (c = 0; c < numComp; ++c) {
1208: uLl[iface*Nc+off+c] = 0.;
1209: uRl[iface*Nc+off+c] = 0.;
1210: }
1211: }
1212: continue;
1213: }
1214: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
1215: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
1216: if (dmGrad) {
1217: PetscReal dxL[3], dxR[3];
1219: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
1220: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
1221: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
1222: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
1223: for (c = 0; c < numComp; ++c) {
1224: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
1225: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
1226: }
1227: } else {
1228: for (c = 0; c < numComp; ++c) {
1229: uLl[iface*Nc+off+c] = xL[c];
1230: uRl[iface*Nc+off+c] = xR[c];
1231: }
1232: }
1233: }
1234: }
1235: ++iface;
1236: }
1237: VecRestoreArrayRead(locX, &x);
1238: VecRestoreArrayRead(faceGeometry, &facegeom);
1239: VecRestoreArrayRead(cellGeometry, &cellgeom);
1240: if (locGrad) {
1241: VecRestoreArrayRead(locGrad, &lgrad);
1242: }
1243: PetscFree(isFE);
1244: return(0);
1245: }
1249: /*@C
1250: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
1252: Input Parameters:
1253: + dm - The DM
1254: . fStart - The first face to include
1255: . fEnd - The first face to exclude
1256: . locX - A local vector with the solution fields
1257: . locX_t - A local vector with solution field time derivatives, or NULL
1258: . faceGeometry - A local vector with face geometry
1259: . cellGeometry - A local vector with cell geometry
1260: - locaGrad - A local vector with field gradients, or NULL
1262: Output Parameters:
1263: + uL - The field values at the left side of the face
1264: - uR - The field values at the right side of the face
1266: Level: developer
1268: .seealso: DMPlexGetFaceFields()
1269: @*/
1270: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1271: {
1275: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uL);
1276: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uR);
1277: return(0);
1278: }
1282: /*@C
1283: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
1285: Input Parameters:
1286: + dm - The DM
1287: . fStart - The first face to include
1288: . fEnd - The first face to exclude
1289: . faceGeometry - A local vector with face geometry
1290: - cellGeometry - A local vector with cell geometry
1292: Output Parameters:
1293: + fgeom - The extract the face centroid and normal
1294: - vol - The cell volume
1296: Level: developer
1298: .seealso: DMPlexGetCellFields()
1299: @*/
1300: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1301: {
1302: DM dmFace, dmCell;
1303: DMLabel ghostLabel;
1304: const PetscScalar *facegeom, *cellgeom;
1305: PetscInt dim, numFaces = fEnd - fStart, iface, face;
1306: PetscErrorCode ierr;
1314: DMGetDimension(dm, &dim);
1315: DMGetLabel(dm, "ghost", &ghostLabel);
1316: VecGetDM(faceGeometry, &dmFace);
1317: VecGetArrayRead(faceGeometry, &facegeom);
1318: VecGetDM(cellGeometry, &dmCell);
1319: VecGetArrayRead(cellGeometry, &cellgeom);
1320: PetscMalloc1(numFaces, fgeom);
1321: DMGetWorkArray(dm, numFaces*2, PETSC_SCALAR, vol);
1322: for (face = fStart, iface = 0; face < fEnd; ++face) {
1323: const PetscInt *cells;
1324: PetscFVFaceGeom *fg;
1325: PetscFVCellGeom *cgL, *cgR;
1326: PetscFVFaceGeom *fgeoml = *fgeom;
1327: PetscReal *voll = *vol;
1328: PetscInt ghost, d;
1330: DMLabelGetValue(ghostLabel, face, &ghost);
1331: if (ghost >= 0) continue;
1332: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1333: DMPlexGetSupport(dm, face, &cells);
1334: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1335: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1336: for (d = 0; d < dim; ++d) {
1337: fgeoml[iface].centroid[d] = fg->centroid[d];
1338: fgeoml[iface].normal[d] = fg->normal[d];
1339: }
1340: voll[iface*2+0] = cgL->volume;
1341: voll[iface*2+1] = cgR->volume;
1342: ++iface;
1343: }
1344: VecRestoreArrayRead(faceGeometry, &facegeom);
1345: VecRestoreArrayRead(cellGeometry, &cellgeom);
1346: return(0);
1347: }
1351: /*@C
1352: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
1354: Input Parameters:
1355: + dm - The DM
1356: . fStart - The first face to include
1357: . fEnd - The first face to exclude
1358: . faceGeometry - A local vector with face geometry
1359: - cellGeometry - A local vector with cell geometry
1361: Output Parameters:
1362: + fgeom - The extract the face centroid and normal
1363: - vol - The cell volume
1365: Level: developer
1367: .seealso: DMPlexGetFaceFields()
1368: @*/
1369: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1370: {
1374: PetscFree(*fgeom);
1375: DMRestoreWorkArray(dm, 0, PETSC_REAL, vol);
1376: return(0);
1377: }
1381: static PetscErrorCode DMPlexApplyLimiter_Internal (DM dm, DM dmCell, PetscLimiter lim, PetscInt dim, PetscInt totDim, PetscInt cell, PetscInt face, PetscInt fStart, PetscInt fEnd, PetscReal *cellPhi, const PetscScalar *x,
1382: const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
1383: {
1384: const PetscInt *children;
1385: PetscInt numChildren;
1386: PetscErrorCode ierr;
1389: DMPlexGetTreeChildren(dm,face,&numChildren,&children);
1390: if (numChildren) {
1391: PetscInt c;
1393: for (c = 0; c < numChildren; c++) {
1394: PetscInt childFace = children[c];
1396: if (childFace >= fStart && childFace < fEnd) {
1397: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,totDim,cell,childFace,fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
1398: }
1399: }
1400: }
1401: else {
1402: PetscScalar *ncx;
1403: PetscFVCellGeom *ncg;
1404: const PetscInt *fcells;
1405: PetscInt ncell, d;
1406: PetscReal v[3];
1408: DMPlexGetSupport(dm, face, &fcells);
1409: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
1410: DMPlexPointLocalRead(dm, ncell, x, &ncx);
1411: DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
1412: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
1413: for (d = 0; d < totDim; ++d) {
1414: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
1415: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DMPlex_DotD_Internal(dim, &cgrad[d*dim], v);
1417: PetscLimiterLimit(lim, flim, &phi);
1418: cellPhi[d] = PetscMin(cellPhi[d], phi);
1419: }
1420: }
1421: return(0);
1422: }
1426: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
1427: {
1428: DM dmFace, dmCell, dmGrad;
1429: DMLabel ghostLabel;
1430: PetscDS prob;
1431: PetscFV fvm;
1432: PetscLimiter lim;
1433: const PetscScalar *facegeom, *cellgeom, *x;
1434: PetscScalar *gr;
1435: PetscReal *cellPhi;
1436: PetscInt dim, face, cell, totDim, cStart, cEnd, cEndInterior;
1437: PetscErrorCode ierr;
1440: DMGetDimension(dm, &dim);
1441: DMGetDS(dm, &prob);
1442: PetscDSGetTotalDimension(prob, &totDim);
1443: DMGetLabel(dm, "ghost", &ghostLabel);
1444: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fvm);
1445: PetscFVGetLimiter(fvm, &lim);
1446: VecGetDM(faceGeometry, &dmFace);
1447: VecGetArrayRead(faceGeometry, &facegeom);
1448: VecGetDM(cellGeometry, &dmCell);
1449: VecGetArrayRead(cellGeometry, &cellgeom);
1450: VecGetArrayRead(locX, &x);
1451: VecGetDM(grad, &dmGrad);
1452: VecZeroEntries(grad);
1453: VecGetArray(grad, &gr);
1454: /* Reconstruct gradients */
1455: for (face = fStart; face < fEnd; ++face) {
1456: const PetscInt *cells;
1457: PetscFVFaceGeom *fg;
1458: PetscScalar *cx[2];
1459: PetscScalar *cgrad[2];
1460: PetscBool boundary;
1461: PetscInt ghost, c, pd, d, numChildren, numCells;
1463: DMLabelGetValue(ghostLabel, face, &ghost);
1464: DMIsBoundaryPoint(dm, face, &boundary);
1465: DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
1466: if (ghost >= 0 || boundary || numChildren) continue;
1467: DMPlexGetSupportSize(dm, face, &numCells);
1468: if (numCells != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "facet %d has %d support points: expected 2",face,numCells);
1469: DMPlexGetSupport(dm, face, &cells);
1470: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1471: for (c = 0; c < 2; ++c) {
1472: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
1473: DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
1474: }
1475: for (pd = 0; pd < totDim; ++pd) {
1476: PetscScalar delta = cx[1][pd] - cx[0][pd];
1478: for (d = 0; d < dim; ++d) {
1479: if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
1480: if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
1481: }
1482: }
1483: }
1484: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
1485: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1486: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1487: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
1488: DMGetWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1489: for (cell = dmGrad && lim ? cStart : cEnd; cell < cEndInterior; ++cell) {
1490: const PetscInt *faces;
1491: PetscScalar *cx;
1492: PetscFVCellGeom *cg;
1493: PetscScalar *cgrad;
1494: PetscInt coneSize, f, pd, d;
1496: DMPlexGetConeSize(dm, cell, &coneSize);
1497: DMPlexGetCone(dm, cell, &faces);
1498: DMPlexPointLocalRead(dm, cell, x, &cx);
1499: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
1500: DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
1501: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
1502: /* Limiter will be minimum value over all neighbors */
1503: for (d = 0; d < totDim; ++d) cellPhi[d] = PETSC_MAX_REAL;
1504: for (f = 0; f < coneSize; ++f) {
1505: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,totDim,cell,faces[f],fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
1506: }
1507: /* Apply limiter to gradient */
1508: for (pd = 0; pd < totDim; ++pd)
1509: /* Scalar limiter applied to each component separately */
1510: for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
1511: }
1512: DMRestoreWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1513: VecRestoreArrayRead(faceGeometry, &facegeom);
1514: VecRestoreArrayRead(cellGeometry, &cellgeom);
1515: VecRestoreArrayRead(locX, &x);
1516: VecRestoreArray(grad, &gr);
1517: return(0);
1518: }
1522: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, Vec locF, void *user)
1523: {
1524: DM_Plex *mesh = (DM_Plex *) dm->data;
1525: PetscSection section;
1526: PetscDS prob;
1527: DMLabel depth;
1528: PetscFECellGeom *cgeom;
1529: PetscScalar *u = NULL, *u_t = NULL, *elemVec = NULL;
1530: PetscInt dim, Nf, f, totDimBd, numBd, bd;
1531: PetscErrorCode ierr;
1534: DMGetDimension(dm, &dim);
1535: DMGetDefaultSection(dm, §ion);
1536: DMGetDS(dm, &prob);
1537: PetscDSGetNumFields(prob, &Nf);
1538: PetscDSGetTotalBdDimension(prob, &totDimBd);
1539: DMPlexGetDepthLabel(dm, &depth);
1540: DMGetNumBoundary(dm, &numBd);
1541: for (bd = 0; bd < numBd; ++bd) {
1542: const char *bdLabel;
1543: DMLabel label;
1544: IS pointIS;
1545: const PetscInt *points;
1546: const PetscInt *values;
1547: PetscInt field, numValues, v, numPoints, p, dep, numFaces;
1548: PetscBool isEssential;
1549: PetscObject obj;
1550: PetscClassId id;
1552: DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1553: DMGetField(dm, field, &obj);
1554: PetscObjectGetClassId(obj, &id);
1555: if ((id != PETSCFE_CLASSID) || isEssential) continue;
1556: DMGetLabel(dm, bdLabel, &label);
1557: for (v = 0; v < numValues; ++v) {
1558: DMLabelGetStratumSize(label, values[v], &numPoints);
1559: DMLabelGetStratumIS(label, values[v], &pointIS);
1560: if (!pointIS) continue; /* No points with that id on this process */
1561: ISGetIndices(pointIS, &points);
1562: for (p = 0, numFaces = 0; p < numPoints; ++p) {
1563: DMLabelGetValue(depth, points[p], &dep);
1564: if (dep == dim-1) ++numFaces;
1565: }
1566: PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd,&elemVec);
1567: if (locX_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1568: for (p = 0, f = 0; p < numPoints; ++p) {
1569: const PetscInt point = points[p];
1570: PetscScalar *x = NULL;
1571: PetscInt i;
1573: DMLabelGetValue(depth, points[p], &dep);
1574: if (dep != dim-1) continue;
1575: DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
1576: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1577: if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1578: /* TODO: Matt, this is wrong if feBd does not match fe: i.e., if the order differs. */
1579: DMPlexVecGetClosure(dm, section, locX, point, NULL, &x);
1580: for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1581: DMPlexVecRestoreClosure(dm, section, locX, point, NULL, &x);
1582: if (locX_t) {
1583: DMPlexVecGetClosure(dm, section, locX_t, point, NULL, &x);
1584: for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1585: DMPlexVecRestoreClosure(dm, section, locX_t, point, NULL, &x);
1586: }
1587: ++f;
1588: }
1589: for (f = 0; f < Nf; ++f) {
1590: PetscFE fe;
1591: PetscQuadrature q;
1592: PetscInt numQuadPoints, Nb;
1593: /* Conforming batches */
1594: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1595: /* Remainder */
1596: PetscInt Nr, offset;
1598: PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);
1599: PetscFEGetQuadrature(fe, &q);
1600: PetscFEGetDimension(fe, &Nb);
1601: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1602: PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1603: blockSize = Nb*numQuadPoints;
1604: batchSize = numBlocks * blockSize;
1605: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1606: numChunks = numFaces / (numBatches*batchSize);
1607: Ne = numChunks*numBatches*batchSize;
1608: Nr = numFaces % (numBatches*batchSize);
1609: offset = numFaces - Nr;
1610: PetscFEIntegrateBdResidual(fe, prob, f, Ne, cgeom, u, u_t, NULL, NULL, elemVec);
1611: PetscFEIntegrateBdResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);
1612: }
1613: for (p = 0, f = 0; p < numPoints; ++p) {
1614: const PetscInt point = points[p];
1616: DMLabelGetValue(depth, point, &dep);
1617: if (dep != dim-1) continue;
1618: if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);}
1619: DMPlexVecSetClosure(dm, NULL, locF, point, &elemVec[f*totDimBd], ADD_ALL_VALUES);
1620: ++f;
1621: }
1622: ISRestoreIndices(pointIS, &points);
1623: ISDestroy(&pointIS);
1624: PetscFree3(u,cgeom,elemVec);
1625: if (locX_t) {PetscFree(u_t);}
1626: }
1627: }
1628: return(0);
1629: }
1633: /*@
1634: DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.
1636: Input Parameters:
1637: + dm - the mesh
1638: - locX - the local representation of the vector
1640: Output Parameter:
1641: . grad - the global representation of the gradient
1643: Level: developer
1645: .seealso: DMPlexSNESGetGradientDM()
1646: @*/
1647: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
1648: {
1649: PetscDS prob;
1650: PetscInt Nf, f, fStart, fEnd;
1651: PetscBool useFVM = PETSC_FALSE;
1652: PetscFV fvm = NULL;
1653: Vec faceGeometryFVM, cellGeometryFVM;
1654: PetscFVCellGeom *cgeomFVM = NULL;
1655: PetscFVFaceGeom *fgeomFVM = NULL;
1656: DM dmGrad = NULL;
1657: PetscErrorCode ierr;
1660: DMGetDS(dm, &prob);
1661: PetscDSGetNumFields(prob, &Nf);
1662: for (f = 0; f < Nf; ++f) {
1663: PetscObject obj;
1664: PetscClassId id;
1666: PetscDSGetDiscretization(prob, f, &obj);
1667: PetscObjectGetClassId(obj, &id);
1668: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1669: }
1670: if (!useFVM) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm does not have a finite volume discretization");
1671: DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1672: if (!dmGrad) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This dm's finite volume discretization does not reconstruct gradients");
1673: DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1674: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1675: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1676: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1677: DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1678: return(0);
1679: }
1683: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
1684: {
1685: DM_Plex *mesh = (DM_Plex *) dm->data;
1686: const char *name = "Residual";
1687: DM dmAux = NULL;
1688: DM dmGrad = NULL;
1689: DMLabel ghostLabel = NULL;
1690: PetscDS prob = NULL;
1691: PetscDS probAux = NULL;
1692: PetscSection section = NULL;
1693: PetscBool useFEM = PETSC_FALSE;
1694: PetscBool useFVM = PETSC_FALSE;
1695: PetscBool isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1696: PetscFV fvm = NULL;
1697: PetscFECellGeom *cgeomFEM = NULL;
1698: PetscScalar *cgeomScal;
1699: PetscFVCellGeom *cgeomFVM = NULL;
1700: PetscFVFaceGeom *fgeomFVM = NULL;
1701: Vec locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1702: PetscScalar *u, *u_t, *a, *uL, *uR;
1703: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1704: PetscErrorCode ierr;
1707: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1708: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1709: /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1710: /* FEM+FVM */
1711: /* 1: Get sizes from dm and dmAux */
1712: DMGetDefaultSection(dm, §ion);
1713: DMGetLabel(dm, "ghost", &ghostLabel);
1714: DMGetDS(dm, &prob);
1715: PetscDSGetNumFields(prob, &Nf);
1716: PetscDSGetTotalDimension(prob, &totDim);
1717: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1718: if (locA) {
1719: VecGetDM(locA, &dmAux);
1720: DMGetDS(dmAux, &probAux);
1721: PetscDSGetTotalDimension(probAux, &totDimAux);
1722: }
1723: /* 2: Get geometric data */
1724: for (f = 0; f < Nf; ++f) {
1725: PetscObject obj;
1726: PetscClassId id;
1727: PetscBool fimp;
1729: PetscDSGetImplicit(prob, f, &fimp);
1730: if (isImplicit != fimp) continue;
1731: PetscDSGetDiscretization(prob, f, &obj);
1732: PetscObjectGetClassId(obj, &id);
1733: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1734: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1735: }
1736: if (useFEM) {
1737: DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1738: VecGetArray(cellGeometryFEM, &cgeomScal);
1739: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1740: DM dmCell;
1741: PetscInt c;
1743: VecGetDM(cellGeometryFEM,&dmCell);
1744: PetscMalloc1(cEnd-cStart,&cgeomFEM);
1745: for (c = 0; c < cEnd - cStart; c++) {
1746: PetscScalar *thisgeom;
1748: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1749: cgeomFEM[c] = *((PetscFECellGeom *) thisgeom);
1750: }
1751: }
1752: else {
1753: cgeomFEM = (PetscFECellGeom *) cgeomScal;
1754: }
1755: }
1756: if (useFVM) {
1757: DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1758: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1759: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1760: /* Reconstruct and limit cell gradients */
1761: DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1762: if (dmGrad) {
1763: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1764: DMGetGlobalVector(dmGrad, &grad);
1765: DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1766: /* Communicate gradient values */
1767: DMGetLocalVector(dmGrad, &locGrad);
1768: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1769: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1770: DMRestoreGlobalVector(dmGrad, &grad);
1771: }
1772: /* Handle non-essential (e.g. outflow) boundary values */
1773: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1774: }
1775: /* Loop over chunks */
1776: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1777: numChunks = 1;
1778: cellChunkSize = (cEnd - cStart)/numChunks;
1779: faceChunkSize = (fEnd - fStart)/numChunks;
1780: for (chunk = 0; chunk < numChunks; ++chunk) {
1781: PetscScalar *elemVec, *fluxL, *fluxR;
1782: PetscReal *vol;
1783: PetscFVFaceGeom *fgeom;
1784: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1785: PetscInt fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = fE - fS, face;
1787: /* Extract field coefficients */
1788: if (useFEM) {
1789: DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1790: DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1791: PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1792: }
1793: if (useFVM) {
1794: DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1795: DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1796: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1797: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1798: PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1799: PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1800: }
1801: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1802: /* Loop over fields */
1803: for (f = 0; f < Nf; ++f) {
1804: PetscObject obj;
1805: PetscClassId id;
1806: PetscBool fimp;
1807: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1809: PetscDSGetImplicit(prob, f, &fimp);
1810: if (isImplicit != fimp) continue;
1811: PetscDSGetDiscretization(prob, f, &obj);
1812: PetscObjectGetClassId(obj, &id);
1813: if (id == PETSCFE_CLASSID) {
1814: PetscFE fe = (PetscFE) obj;
1815: PetscQuadrature q;
1816: PetscInt Nq, Nb;
1818: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1820: PetscFEGetQuadrature(fe, &q);
1821: PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1822: PetscFEGetDimension(fe, &Nb);
1823: blockSize = Nb*Nq;
1824: batchSize = numBlocks * blockSize;
1825: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1826: numChunks = numCells / (numBatches*batchSize);
1827: Ne = numChunks*numBatches*batchSize;
1828: Nr = numCells % (numBatches*batchSize);
1829: offset = numCells - Nr;
1830: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1831: /* For FV, I think we use a P0 basis and the cell coefficients (for subdivided cells, we can tweak the basis tabulation to be the indicator function) */
1832: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, elemVec);
1833: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1834: } else if (id == PETSCFV_CLASSID) {
1835: PetscFV fv = (PetscFV) obj;
1837: Ne = numFaces;
1838: /* Riemann solve over faces (need fields at face centroids) */
1839: /* We need to evaluate FE fields at those coordinates */
1840: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1841: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1842: }
1843: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1844: PetscFree(cgeomFEM);
1845: }
1846: else {
1847: cgeomFEM = NULL;
1848: }
1849: if (cellGeometryFEM) {VecRestoreArray(cellGeometryFEM, &cgeomScal);}
1850: /* Loop over domain */
1851: if (useFEM) {
1852: /* Add elemVec to locX */
1853: for (cell = cS; cell < cE; ++cell) {
1854: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1855: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_ALL_VALUES);
1856: }
1857: }
1858: if (useFVM) {
1859: PetscScalar *fa;
1860: PetscInt iface;
1862: VecGetArray(locF, &fa);
1863: for (f = 0; f < Nf; ++f) {
1864: PetscFV fv;
1865: PetscObject obj;
1866: PetscClassId id;
1867: PetscInt foff, pdim;
1869: PetscDSGetDiscretization(prob, f, &obj);
1870: PetscDSGetFieldOffset(prob, f, &foff);
1871: PetscObjectGetClassId(obj, &id);
1872: if (id != PETSCFV_CLASSID) continue;
1873: fv = (PetscFV) obj;
1874: PetscFVGetNumComponents(fv, &pdim);
1875: /* Accumulate fluxes to cells */
1876: for (face = fS, iface = 0; face < fE; ++face) {
1877: const PetscInt *cells;
1878: PetscScalar *fL, *fR;
1879: PetscInt ghost, d, nsupp;
1881: DMLabelGetValue(ghostLabel, face, &ghost);
1882: DMPlexGetSupportSize(dm, face, &nsupp);
1883: if (ghost >= 0) continue;
1884: if (nsupp > 2) { /* noop */
1885: ++iface;
1886: continue;
1887: }
1888: DMPlexGetSupport(dm, face, &cells);
1889: DMPlexPointGlobalFieldRef(dm, cells[0], f, fa, &fL);
1890: DMPlexPointGlobalFieldRef(dm, cells[1], f, fa, &fR);
1891: for (d = 0; d < pdim; ++d) {
1892: if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1893: if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1894: }
1895: ++iface;
1896: }
1897: }
1898: VecRestoreArray(locF, &fa);
1899: }
1900: /* Handle time derivative */
1901: if (locX_t) {
1902: PetscScalar *x_t, *fa;
1904: VecGetArray(locF, &fa);
1905: VecGetArray(locX_t, &x_t);
1906: for (f = 0; f < Nf; ++f) {
1907: PetscFV fv;
1908: PetscObject obj;
1909: PetscClassId id;
1910: PetscInt pdim, d;
1912: PetscDSGetDiscretization(prob, f, &obj);
1913: PetscObjectGetClassId(obj, &id);
1914: if (id != PETSCFV_CLASSID) continue;
1915: fv = (PetscFV) obj;
1916: PetscFVGetNumComponents(fv, &pdim);
1917: for (cell = cS; cell < cE; ++cell) {
1918: PetscScalar *u_t, *r;
1920: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1921: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1922: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1923: }
1924: }
1925: VecRestoreArray(locX_t, &x_t);
1926: VecRestoreArray(locF, &fa);
1927: }
1928: if (useFEM) {
1929: DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1930: DMRestoreWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1931: }
1932: if (useFVM) {
1933: DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1934: DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1935: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1936: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1937: if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1938: }
1939: }
1941: if (useFEM) {DMPlexComputeBdResidual_Internal(dm, locX, locX_t, locF, user);}
1943: /* FEM */
1944: /* 1: Get sizes from dm and dmAux */
1945: /* 2: Get geometric data */
1946: /* 3: Handle boundary values */
1947: /* 4: Loop over domain */
1948: /* Extract coefficients */
1949: /* Loop over fields */
1950: /* Set tiling for FE*/
1951: /* Integrate FE residual to get elemVec */
1952: /* Loop over subdomain */
1953: /* Loop over quad points */
1954: /* Transform coords to real space */
1955: /* Evaluate field and aux fields at point */
1956: /* Evaluate residual at point */
1957: /* Transform residual to real space */
1958: /* Add residual to elemVec */
1959: /* Loop over domain */
1960: /* Add elemVec to locX */
1962: /* FVM */
1963: /* Get geometric data */
1964: /* If using gradients */
1965: /* Compute gradient data */
1966: /* Loop over domain faces */
1967: /* Count computational faces */
1968: /* Reconstruct cell gradient */
1969: /* Loop over domain cells */
1970: /* Limit cell gradients */
1971: /* Handle boundary values */
1972: /* Loop over domain faces */
1973: /* Read out field, centroid, normal, volume for each side of face */
1974: /* Riemann solve over faces */
1975: /* Loop over domain faces */
1976: /* Accumulate fluxes to cells */
1977: /* TODO Change printFEM to printDisc here */
1978: if (mesh->printFEM) {
1979: Vec locFbc;
1980: PetscInt pStart, pEnd, p, maxDof;
1981: PetscScalar *zeroes;
1983: VecDuplicate(locF,&locFbc);
1984: VecCopy(locF,locFbc);
1985: PetscSectionGetChart(section,&pStart,&pEnd);
1986: PetscSectionGetMaxDof(section,&maxDof);
1987: PetscCalloc1(maxDof,&zeroes);
1988: for (p = pStart; p < pEnd; p++) {
1989: VecSetValuesSection(locFbc,section,p,zeroes,INSERT_BC_VALUES);
1990: }
1991: PetscFree(zeroes);
1992: DMPrintLocalVec(dm, name, mesh->printTol, locFbc);
1993: VecDestroy(&locFbc);
1994: }
1995: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1996: return(0);
1997: }
2001: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
2002: {
2003: DM dmCh, dmAux;
2004: Vec A, cellgeom;
2005: PetscDS prob, probCh, probAux = NULL;
2006: PetscQuadrature q;
2007: PetscSection section, sectionAux;
2008: PetscFECellGeom *cgeom = NULL;
2009: PetscScalar *cgeomScal;
2010: PetscScalar *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
2011: PetscInt dim, Nf, f, numCells, cStart, cEnd, c;
2012: PetscInt totDim, totDimAux, diffCell = 0;
2013: PetscErrorCode ierr;
2016: DMGetDimension(dm, &dim);
2017: DMGetDefaultSection(dm, §ion);
2018: DMGetDS(dm, &prob);
2019: PetscDSGetTotalDimension(prob, &totDim);
2020: PetscSectionGetNumFields(section, &Nf);
2021: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2022: numCells = cEnd - cStart;
2023: PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
2024: DMGetDS(dmCh, &probCh);
2025: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2026: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2027: if (dmAux) {
2028: DMGetDefaultSection(dmAux, §ionAux);
2029: DMGetDS(dmAux, &probAux);
2030: PetscDSGetTotalDimension(probAux, &totDimAux);
2031: }
2032: VecSet(F, 0.0);
2033: PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
2034: PetscMalloc1(numCells*totDim,&elemVecCh);
2035: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2036: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2037: VecGetArray(cellgeom, &cgeomScal);
2038: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2039: DM dmCell;
2041: VecGetDM(cellgeom,&dmCell);
2042: PetscMalloc1(cEnd-cStart,&cgeom);
2043: for (c = 0; c < cEnd - cStart; c++) {
2044: PetscScalar *thisgeom;
2046: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2047: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2048: }
2049: }
2050: else {
2051: cgeom = (PetscFECellGeom *) cgeomScal;
2052: }
2053: for (c = cStart; c < cEnd; ++c) {
2054: PetscScalar *x = NULL, *x_t = NULL;
2055: PetscInt i;
2057: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2058: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
2059: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2060: if (X_t) {
2061: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2062: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
2063: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2064: }
2065: if (dmAux) {
2066: DM dmAuxPlex;
2068: DMSNESConvertPlex(dmAux,&dmAuxPlex, PETSC_FALSE);
2069: DMPlexVecGetClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
2070: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
2071: DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
2072: DMDestroy(&dmAuxPlex);
2073: }
2074: }
2075: for (f = 0; f < Nf; ++f) {
2076: PetscFE fe, feCh;
2077: PetscInt numQuadPoints, Nb;
2078: /* Conforming batches */
2079: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2080: /* Remainder */
2081: PetscInt Nr, offset;
2083: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
2084: PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
2085: PetscFEGetQuadrature(fe, &q);
2086: PetscFEGetDimension(fe, &Nb);
2087: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2088: PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
2089: blockSize = Nb*numQuadPoints;
2090: batchSize = numBlocks * blockSize;
2091: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2092: numChunks = numCells / (numBatches*batchSize);
2093: Ne = numChunks*numBatches*batchSize;
2094: Nr = numCells % (numBatches*batchSize);
2095: offset = numCells - Nr;
2096: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);
2097: PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVecCh);
2098: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
2099: PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);
2100: }
2101: for (c = cStart; c < cEnd; ++c) {
2102: PetscBool diff = PETSC_FALSE;
2103: PetscInt d;
2105: for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
2106: if (diff) {
2107: PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
2108: DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
2109: DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
2110: ++diffCell;
2111: }
2112: if (diffCell > 9) break;
2113: DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_ALL_VALUES);
2114: }
2115: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2116: PetscFree(cgeom);
2117: }
2118: else {
2119: cgeom = NULL;
2120: }
2121: VecRestoreArray(cellgeom, &cgeomScal);
2122: PetscFree3(u,u_t,elemVec);
2123: PetscFree(elemVecCh);
2124: if (dmAux) {PetscFree(a);}
2125: return(0);
2126: }
2130: /*@
2131: DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
2133: Input Parameters:
2134: + dm - The mesh
2135: . X - Local solution
2136: - user - The user context
2138: Output Parameter:
2139: . F - Local output vector
2141: Level: developer
2143: .seealso: DMPlexComputeJacobianActionFEM()
2144: @*/
2145: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
2146: {
2147: PetscObject check;
2148: PetscInt cStart, cEnd, cEndInterior;
2149: DM plex;
2153: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2154: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2155: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2156: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2157: /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
2158: PetscObjectQuery((PetscObject) plex, "dmCh", &check);
2159: if (check) {DMPlexComputeResidualFEM_Check_Internal(plex, X, NULL, F, user);}
2160: else {DMPlexComputeResidual_Internal(plex, cStart, cEnd, PETSC_MIN_REAL, X, NULL, F, user);}
2161: DMDestroy(&plex);
2162: return(0);
2163: }
2167: /*@
2168: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
2170: Input Parameters:
2171: + dm - The mesh
2172: - user - The user context
2174: Output Parameter:
2175: . X - Local solution
2177: Level: developer
2179: .seealso: DMPlexComputeJacobianActionFEM()
2180: @*/
2181: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
2182: {
2183: DM plex;
2187: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2188: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
2189: DMDestroy(&plex);
2190: return(0);
2191: }
2195: PetscErrorCode DMPlexComputeJacobian_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Mat Jac, Mat JacP,void *user)
2196: {
2197: DM_Plex *mesh = (DM_Plex *) dm->data;
2198: const char *name = "Jacobian";
2199: DM dmAux, plex;
2200: DMLabel depth;
2201: Vec A, cellgeom;
2202: PetscDS prob, probAux = NULL;
2203: PetscQuadrature quad;
2204: PetscSection section, globalSection, sectionAux;
2205: PetscFECellGeom *cgeom = NULL;
2206: PetscScalar *cgeomScal;
2207: PetscScalar *elemMat, *elemMatP, *elemMatD, *u, *u_t, *a = NULL;
2208: PetscInt dim, Nf, f, fieldI, fieldJ, numCells, c;
2209: PetscInt totDim, totDimBd, totDimAux, numBd, bd;
2210: PetscBool isShell, hasPrec, hasDyn;
2211: PetscErrorCode ierr;
2214: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2215: DMGetDimension(dm, &dim);
2216: DMGetDefaultSection(dm, §ion);
2217: DMGetDefaultGlobalSection(dm, &globalSection);
2218: DMGetDS(dm, &prob);
2219: PetscDSGetTotalDimension(prob, &totDim);
2220: PetscDSGetTotalBdDimension(prob, &totDimBd);
2221: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2222: PetscDSHasDynamicJacobian(prob, &hasDyn);
2223: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2224: PetscSectionGetNumFields(section, &Nf);
2225: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2226: numCells = cEnd - cStart;
2227: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2228: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2229: if (dmAux) {
2230: DMConvert(dmAux, DMPLEX, &plex);
2231: DMGetDefaultSection(plex, §ionAux);
2232: DMGetDS(dmAux, &probAux);
2233: PetscDSGetTotalDimension(probAux, &totDimAux);
2234: }
2235: MatZeroEntries(JacP);
2236: PetscMalloc5(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasPrec ? numCells*totDim*totDim : 0, &elemMatP,hasDyn ? numCells*totDim*totDim : 0, &elemMatD);
2237: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2238: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2239: VecGetArray(cellgeom, &cgeomScal);
2240: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2241: DM dmCell;
2243: VecGetDM(cellgeom,&dmCell);
2244: PetscMalloc1(cEnd-cStart,&cgeom);
2245: for (c = 0; c < cEnd - cStart; c++) {
2246: PetscScalar *thisgeom;
2248: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2249: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2250: }
2251: }
2252: else {
2253: cgeom = (PetscFECellGeom *) cgeomScal;
2254: }
2255: for (c = cStart; c < cEnd; ++c) {
2256: PetscScalar *x = NULL, *x_t = NULL;
2257: PetscInt i;
2259: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2260: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2261: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2262: if (X_t) {
2263: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2264: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2265: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2266: }
2267: if (dmAux) {
2268: DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2269: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2270: DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2271: }
2272: }
2273: PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2274: if (hasPrec) {PetscMemzero(elemMatP, numCells*totDim*totDim * sizeof(PetscScalar));}
2275: if (hasDyn) {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2276: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2277: PetscFE fe;
2278: PetscInt numQuadPoints, Nb;
2279: /* Conforming batches */
2280: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2281: /* Remainder */
2282: PetscInt Nr, offset;
2284: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2285: PetscFEGetQuadrature(fe, &quad);
2286: PetscFEGetDimension(fe, &Nb);
2287: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2288: PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2289: blockSize = Nb*numQuadPoints;
2290: batchSize = numBlocks * blockSize;
2291: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2292: numChunks = numCells / (numBatches*batchSize);
2293: Ne = numChunks*numBatches*batchSize;
2294: Nr = numCells % (numBatches*batchSize);
2295: offset = numCells - Nr;
2296: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2297: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMat);
2298: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);
2299: if (hasPrec) {
2300: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMatP);
2301: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMatP[offset*totDim*totDim]);
2302: }
2303: if (hasDyn) {
2304: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMatD);
2305: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMatD[offset*totDim*totDim]);
2306: }
2307: }
2308: }
2309: if (hasDyn) {
2310: for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2311: }
2312: for (c = cStart; c < cEnd; ++c) {
2313: if (hasPrec) {
2314: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2315: DMPlexMatSetClosure(dm, section, globalSection, Jac, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2316: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
2317: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2318: } else {
2319: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2320: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2321: }
2322: }
2323: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2324: PetscFree(cgeom);
2325: }
2326: else {
2327: cgeom = NULL;
2328: }
2329: VecRestoreArray(cellgeom, &cgeomScal);
2330: PetscFree5(u,u_t,elemMat,elemMatP,elemMatD);
2331: if (dmAux) {
2332: PetscFree(a);
2333: DMDestroy(&plex);
2334: }
2335: DMPlexGetDepthLabel(dm, &depth);
2336: DMGetNumBoundary(dm, &numBd);
2337: DMPlexGetDepthLabel(dm, &depth);
2338: DMGetNumBoundary(dm, &numBd);
2339: for (bd = 0; bd < numBd; ++bd) {
2340: const char *bdLabel;
2341: DMLabel label;
2342: IS pointIS;
2343: const PetscInt *points;
2344: const PetscInt *values;
2345: PetscInt field, numValues, v, numPoints, p, dep, numFaces;
2346: PetscBool isEssential;
2347: PetscObject obj;
2348: PetscClassId id;
2350: DMGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
2351: DMGetField(dm, field, &obj);
2352: PetscObjectGetClassId(obj, &id);
2353: if ((id != PETSCFE_CLASSID) || isEssential) continue;
2354: DMGetLabel(dm, bdLabel, &label);
2355: for (v = 0; v < numValues; ++v) {
2356: DMLabelGetStratumSize(label, values[v], &numPoints);
2357: DMLabelGetStratumIS(label, values[v], &pointIS);
2358: if (!pointIS) continue; /* No points with that id on this process */
2359: ISGetIndices(pointIS, &points);
2360: for (p = 0, numFaces = 0; p < numPoints; ++p) {
2361: DMLabelGetValue(depth, points[p], &dep);
2362: if (dep == dim-1) ++numFaces;
2363: }
2364: PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd*totDimBd,&elemMat);
2365: if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
2366: for (p = 0, f = 0; p < numPoints; ++p) {
2367: const PetscInt point = points[p];
2368: PetscScalar *x = NULL;
2369: PetscInt i;
2371: DMLabelGetValue(depth, points[p], &dep);
2372: if (dep != dim-1) continue;
2373: DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
2374: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
2375: if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
2376: DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
2377: for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
2378: DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
2379: if (X_t) {
2380: DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
2381: for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
2382: DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
2383: }
2384: ++f;
2385: }
2386: PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));
2387: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2388: PetscFE fe;
2389: PetscInt numQuadPoints, Nb;
2390: /* Conforming batches */
2391: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2392: /* Remainder */
2393: PetscInt Nr, offset;
2395: PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
2396: PetscFEGetQuadrature(fe, &quad);
2397: PetscFEGetDimension(fe, &Nb);
2398: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2399: PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2400: blockSize = Nb*numQuadPoints;
2401: batchSize = numBlocks * blockSize;
2402: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2403: numChunks = numFaces / (numBatches*batchSize);
2404: Ne = numChunks*numBatches*batchSize;
2405: Nr = numFaces % (numBatches*batchSize);
2406: offset = numFaces - Nr;
2407: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2408: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, NULL, NULL, elemMat);
2409: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);
2410: }
2411: }
2412: for (p = 0, f = 0; p < numPoints; ++p) {
2413: const PetscInt point = points[p];
2415: DMLabelGetValue(depth, point, &dep);
2416: if (dep != dim-1) continue;
2417: if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);}
2418: DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);
2419: ++f;
2420: }
2421: ISRestoreIndices(pointIS, &points);
2422: ISDestroy(&pointIS);
2423: PetscFree3(u,cgeom,elemMat);
2424: if (X_t) {PetscFree(u_t);}
2425: }
2426: }
2427: if (hasPrec) {
2428: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
2429: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
2430: }
2431: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2432: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2433: if (mesh->printFEM) {
2434: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2435: MatChop(JacP, 1.0e-10);
2436: MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2437: }
2438: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2439: PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2440: if (isShell) {
2441: JacActionCtx *jctx;
2443: MatShellGetContext(Jac, &jctx);
2444: VecCopy(X, jctx->u);
2445: }
2446: return(0);
2447: }
2451: /*@
2452: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
2454: Input Parameters:
2455: + dm - The mesh
2456: . X - Local input vector
2457: - user - The user context
2459: Output Parameter:
2460: . Jac - Jacobian matrix
2462: Note:
2463: The first member of the user context must be an FEMContext.
2465: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2466: like a GPU, or vectorize on a multicore machine.
2468: Level: developer
2470: .seealso: FormFunctionLocal()
2471: @*/
2472: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2473: {
2474: PetscInt cStart, cEnd, cEndInterior;
2475: DM plex;
2479: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2480: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2481: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2482: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2483: DMPlexComputeJacobian_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2484: DMDestroy(&plex);
2485: return(0);
2486: }
2490: /*@
2491: DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.
2493: Input Parameters:
2494: + dm - The DM object
2495: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see DMAddBoundary())
2496: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
2497: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())
2499: Level: developer
2500: @*/
2501: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
2502: {
2506: DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
2507: DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
2508: DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
2509: return(0);
2510: }
2514: PetscErrorCode DMSNESCheckFromOptions_Internal(SNES snes, DM dm, Vec u, Vec sol, PetscErrorCode (**exactFuncs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx), void **ctxs)
2515: {
2516: Mat J, M;
2517: Vec r, b;
2518: MatNullSpace nullSpace;
2519: PetscReal *error, res = 0.0;
2520: PetscInt numFields;
2524: VecDuplicate(u, &r);
2525: DMCreateMatrix(dm, &J);
2526: M = J;
2527: /* TODO Null space for J */
2528: /* Check discretization error */
2529: DMGetNumFields(dm, &numFields);
2530: PetscMalloc1(PetscMax(1, numFields), &error);
2531: if (numFields > 1) {
2532: PetscInt f;
2534: DMComputeL2FieldDiff(dm, 0.0, exactFuncs, ctxs, u, error);
2535: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2536: for (f = 0; f < numFields; ++f) {
2537: if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2538: if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", error[f]);}
2539: else {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2540: }
2541: PetscPrintf(PETSC_COMM_WORLD, "]\n");
2542: } else {
2543: DMComputeL2Diff(dm, 0.0, exactFuncs, ctxs, u, &error[0]);
2544: if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error[0]);}
2545: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2546: }
2547: PetscFree(error);
2548: /* Check residual */
2549: SNESComputeFunction(snes, u, r);
2550: VecNorm(r, NORM_2, &res);
2551: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
2552: VecChop(r, 1.0e-10);
2553: PetscObjectSetName((PetscObject) r, "Initial Residual");
2554: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2555: VecViewFromOptions(r, NULL, "-vec_view");
2556: /* Check Jacobian */
2557: SNESComputeJacobian(snes, u, M, M);
2558: MatGetNullSpace(J, &nullSpace);
2559: if (nullSpace) {
2560: PetscBool isNull;
2561: MatNullSpaceTest(nullSpace, J, &isNull);
2562: if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2563: }
2564: VecDuplicate(u, &b);
2565: VecSet(r, 0.0);
2566: SNESComputeFunction(snes, r, b);
2567: MatMult(M, u, r);
2568: VecAXPY(r, 1.0, b);
2569: VecDestroy(&b);
2570: VecNorm(r, NORM_2, &res);
2571: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
2572: VecChop(r, 1.0e-10);
2573: PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2574: PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2575: VecViewFromOptions(r, NULL, "-vec_view");
2576: VecDestroy(&r);
2577: MatNullSpaceDestroy(&nullSpace);
2578: MatDestroy(&J);
2579: return(0);
2580: }
2584: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2585: {
2586: DM dm;
2587: Vec sol;
2588: PetscBool check;
2592: PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
2593: if (!check) return(0);
2594: SNESGetDM(snes, &dm);
2595: VecDuplicate(u, &sol);
2596: SNESSetSolution(snes, sol);
2597: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2598: VecDestroy(&sol);
2599: return(0);
2600: }