Actual source code: dmplexsnes.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petscblaslapack.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/petscfeimpl.h>
8: /************************** Interpolation *******************************/
10: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
11: {
12: PetscBool isPlex;
16: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
17: if (isPlex) {
18: *plex = dm;
19: PetscObjectReference((PetscObject) dm);
20: } else {
21: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
22: if (!*plex) {
23: DMConvert(dm,DMPLEX,plex);
24: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
25: if (copy) {
26: PetscInt i;
27: PetscObject obj;
28: const char *comps[3] = {"A","dmAux","dmCh"};
30: DMCopyDMSNES(dm, *plex);
31: for (i = 0; i < 3; i++) {
32: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
33: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
34: }
35: }
36: } else {
37: PetscObjectReference((PetscObject) *plex);
38: }
39: }
40: return(0);
41: }
43: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
44: {
49: PetscNew(ctx);
51: (*ctx)->comm = comm;
52: (*ctx)->dim = -1;
53: (*ctx)->nInput = 0;
54: (*ctx)->points = NULL;
55: (*ctx)->cells = NULL;
56: (*ctx)->n = -1;
57: (*ctx)->coords = NULL;
58: return(0);
59: }
61: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
62: {
64: if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
65: ctx->dim = dim;
66: return(0);
67: }
69: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
70: {
73: *dim = ctx->dim;
74: return(0);
75: }
77: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
78: {
80: if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
81: ctx->dof = dof;
82: return(0);
83: }
85: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
86: {
89: *dof = ctx->dof;
90: return(0);
91: }
93: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
94: {
98: if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
99: if (ctx->points) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
100: ctx->nInput = n;
102: PetscMalloc1(n*ctx->dim, &ctx->points);
103: PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
104: return(0);
105: }
107: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
108: {
109: MPI_Comm comm = ctx->comm;
110: PetscScalar *a;
111: PetscInt p, q, i;
112: PetscMPIInt rank, size;
113: PetscErrorCode ierr;
114: Vec pointVec;
115: PetscSF cellSF;
116: PetscLayout layout;
117: PetscReal *globalPoints;
118: PetscScalar *globalPointsScalar;
119: const PetscInt *ranges;
120: PetscMPIInt *counts, *displs;
121: const PetscSFNode *foundCells;
122: const PetscInt *foundPoints;
123: PetscMPIInt *foundProcs, *globalProcs;
124: PetscInt n, N, numFound;
128: MPI_Comm_size(comm, &size);
129: MPI_Comm_rank(comm, &rank);
130: if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
131: /* Locate points */
132: n = ctx->nInput;
133: if (!redundantPoints) {
134: PetscLayoutCreate(comm, &layout);
135: PetscLayoutSetBlockSize(layout, 1);
136: PetscLayoutSetLocalSize(layout, n);
137: PetscLayoutSetUp(layout);
138: PetscLayoutGetSize(layout, &N);
139: /* Communicate all points to all processes */
140: PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
141: PetscLayoutGetRanges(layout, &ranges);
142: for (p = 0; p < size; ++p) {
143: counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
144: displs[p] = ranges[p]*ctx->dim;
145: }
146: MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
147: } else {
148: N = n;
149: globalPoints = ctx->points;
150: counts = displs = NULL;
151: layout = NULL;
152: }
153: #if 0
154: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
155: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
156: #else
157: #if defined(PETSC_USE_COMPLEX)
158: PetscMalloc1(N*ctx->dim,&globalPointsScalar);
159: for (i=0; i<N*ctx->dim; i++) globalPointsScalar[i] = globalPoints[i];
160: #else
161: globalPointsScalar = globalPoints;
162: #endif
163: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
164: PetscMalloc2(N,&foundProcs,N,&globalProcs);
165: for (p = 0; p < N; ++p) {foundProcs[p] = size;}
166: cellSF = NULL;
167: DMLocatePoints(dm, pointVec, DM_POINTLOCATION_REMOVE, &cellSF);
168: PetscSFGetGraph(cellSF,NULL,&numFound,&foundPoints,&foundCells);
169: #endif
170: for (p = 0; p < numFound; ++p) {
171: if (foundCells[p].index >= 0) foundProcs[foundPoints ? foundPoints[p] : p] = rank;
172: }
173: /* Let the lowest rank process own each point */
174: MPIU_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
175: ctx->n = 0;
176: for (p = 0; p < N; ++p) {
177: if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, (double)globalPoints[p*ctx->dim+0], (double)(ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0), (double)(ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0));
178: else if (globalProcs[p] == rank) ctx->n++;
179: }
180: /* Create coordinates vector and array of owned cells */
181: PetscMalloc1(ctx->n, &ctx->cells);
182: VecCreate(comm, &ctx->coords);
183: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
184: VecSetBlockSize(ctx->coords, ctx->dim);
185: VecSetType(ctx->coords,VECSTANDARD);
186: VecGetArray(ctx->coords, &a);
187: for (p = 0, q = 0, i = 0; p < N; ++p) {
188: if (globalProcs[p] == rank) {
189: PetscInt d;
191: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
192: ctx->cells[q] = foundCells[q].index;
193: ++q;
194: }
195: }
196: VecRestoreArray(ctx->coords, &a);
197: #if 0
198: PetscFree3(foundCells,foundProcs,globalProcs);
199: #else
200: PetscFree2(foundProcs,globalProcs);
201: PetscSFDestroy(&cellSF);
202: VecDestroy(&pointVec);
203: #endif
204: if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
205: if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
206: PetscLayoutDestroy(&layout);
207: return(0);
208: }
210: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
211: {
214: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
215: *coordinates = ctx->coords;
216: return(0);
217: }
219: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
220: {
225: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
226: VecCreate(ctx->comm, v);
227: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
228: VecSetBlockSize(*v, ctx->dof);
229: VecSetType(*v,VECSTANDARD);
230: return(0);
231: }
233: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
234: {
239: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
240: VecDestroy(v);
241: return(0);
242: }
244: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
245: {
246: PetscReal *v0, *J, *invJ, detJ;
247: const PetscScalar *coords;
248: PetscScalar *a;
249: PetscInt p;
253: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
254: VecGetArrayRead(ctx->coords, &coords);
255: VecGetArray(v, &a);
256: for (p = 0; p < ctx->n; ++p) {
257: PetscInt c = ctx->cells[p];
258: PetscScalar *x = NULL;
259: PetscReal xi[4];
260: PetscInt d, f, comp;
262: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
263: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", (double)detJ, c);
264: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
265: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
267: for (d = 0; d < ctx->dim; ++d) {
268: xi[d] = 0.0;
269: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
270: 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];
271: }
272: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
273: }
274: VecRestoreArray(v, &a);
275: VecRestoreArrayRead(ctx->coords, &coords);
276: PetscFree3(v0, J, invJ);
277: return(0);
278: }
280: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
281: {
282: PetscReal *v0, *J, *invJ, detJ;
283: const PetscScalar *coords;
284: PetscScalar *a;
285: PetscInt p;
289: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
290: VecGetArrayRead(ctx->coords, &coords);
291: VecGetArray(v, &a);
292: for (p = 0; p < ctx->n; ++p) {
293: PetscInt c = ctx->cells[p];
294: const PetscInt order[3] = {2, 1, 3};
295: PetscScalar *x = NULL;
296: PetscReal xi[4];
297: PetscInt d, f, comp;
299: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
300: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", (double)detJ, c);
301: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
302: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
304: for (d = 0; d < ctx->dim; ++d) {
305: xi[d] = 0.0;
306: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
307: 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];
308: }
309: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
310: }
311: VecRestoreArray(v, &a);
312: VecRestoreArrayRead(ctx->coords, &coords);
313: PetscFree3(v0, J, invJ);
314: return(0);
315: }
317: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
318: {
319: const PetscScalar *vertices = (const PetscScalar*) ctx;
320: const PetscScalar x0 = vertices[0];
321: const PetscScalar y0 = vertices[1];
322: const PetscScalar x1 = vertices[2];
323: const PetscScalar y1 = vertices[3];
324: const PetscScalar x2 = vertices[4];
325: const PetscScalar y2 = vertices[5];
326: const PetscScalar x3 = vertices[6];
327: const PetscScalar y3 = vertices[7];
328: const PetscScalar f_1 = x1 - x0;
329: const PetscScalar g_1 = y1 - y0;
330: const PetscScalar f_3 = x3 - x0;
331: const PetscScalar g_3 = y3 - y0;
332: const PetscScalar f_01 = x2 - x1 - x3 + x0;
333: const PetscScalar g_01 = y2 - y1 - y3 + y0;
334: const PetscScalar *ref;
335: PetscScalar *real;
336: PetscErrorCode ierr;
339: VecGetArrayRead(Xref, &ref);
340: VecGetArray(Xreal, &real);
341: {
342: const PetscScalar p0 = ref[0];
343: const PetscScalar p1 = ref[1];
345: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
346: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
347: }
348: PetscLogFlops(28);
349: VecRestoreArrayRead(Xref, &ref);
350: VecRestoreArray(Xreal, &real);
351: return(0);
352: }
354: #include <petsc/private/dmimpl.h>
355: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
356: {
357: const PetscScalar *vertices = (const PetscScalar*) ctx;
358: const PetscScalar x0 = vertices[0];
359: const PetscScalar y0 = vertices[1];
360: const PetscScalar x1 = vertices[2];
361: const PetscScalar y1 = vertices[3];
362: const PetscScalar x2 = vertices[4];
363: const PetscScalar y2 = vertices[5];
364: const PetscScalar x3 = vertices[6];
365: const PetscScalar y3 = vertices[7];
366: const PetscScalar f_01 = x2 - x1 - x3 + x0;
367: const PetscScalar g_01 = y2 - y1 - y3 + y0;
368: const PetscScalar *ref;
369: PetscErrorCode ierr;
372: VecGetArrayRead(Xref, &ref);
373: {
374: const PetscScalar x = ref[0];
375: const PetscScalar y = ref[1];
376: const PetscInt rows[2] = {0, 1};
377: PetscScalar values[4];
379: values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
380: values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
381: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
382: }
383: PetscLogFlops(30);
384: VecRestoreArrayRead(Xref, &ref);
385: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
386: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
387: return(0);
388: }
390: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
391: {
392: DM dmCoord;
393: PetscFE fem = NULL;
394: SNES snes;
395: KSP ksp;
396: PC pc;
397: Vec coordsLocal, r, ref, real;
398: Mat J;
399: const PetscScalar *coords;
400: PetscScalar *a;
401: PetscInt Nf, p;
402: const PetscInt dof = ctx->dof;
406: DMGetNumFields(dm, &Nf);
407: if (Nf) {DMGetField(dm, 0, (PetscObject *) &fem);}
408: DMGetCoordinatesLocal(dm, &coordsLocal);
409: DMGetCoordinateDM(dm, &dmCoord);
410: SNESCreate(PETSC_COMM_SELF, &snes);
411: SNESSetOptionsPrefix(snes, "quad_interp_");
412: VecCreate(PETSC_COMM_SELF, &r);
413: VecSetSizes(r, 2, 2);
414: VecSetType(r,dm->vectype);
415: VecDuplicate(r, &ref);
416: VecDuplicate(r, &real);
417: MatCreate(PETSC_COMM_SELF, &J);
418: MatSetSizes(J, 2, 2, 2, 2);
419: MatSetType(J, MATSEQDENSE);
420: MatSetUp(J);
421: SNESSetFunction(snes, r, QuadMap_Private, NULL);
422: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
423: SNESGetKSP(snes, &ksp);
424: KSPGetPC(ksp, &pc);
425: PCSetType(pc, PCLU);
426: SNESSetFromOptions(snes);
428: VecGetArrayRead(ctx->coords, &coords);
429: VecGetArray(v, &a);
430: for (p = 0; p < ctx->n; ++p) {
431: PetscScalar *x = NULL, *vertices = NULL;
432: PetscScalar *xi;
433: PetscReal xir[2];
434: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
436: /* Can make this do all points at once */
437: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
438: if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
439: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
440: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
441: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
442: VecGetArray(real, &xi);
443: xi[0] = coords[p*ctx->dim+0];
444: xi[1] = coords[p*ctx->dim+1];
445: VecRestoreArray(real, &xi);
446: SNESSolve(snes, real, ref);
447: VecGetArray(ref, &xi);
448: xir[0] = PetscRealPart(xi[0]);
449: xir[1] = PetscRealPart(xi[1]);
450: if (4*dof != xSize) {
451: PetscReal *B;
452: PetscInt d;
454: xir[0] = 2.0*xir[0] - 1.0; xir[1] = 2.0*xir[1] - 1.0;
455: PetscFEGetTabulation(fem, 1, xir, &B, NULL, NULL);
456: for (comp = 0; comp < dof; ++comp) {
457: a[p*dof+comp] = 0.0;
458: for (d = 0; d < xSize/dof; ++d) {
459: a[p*dof+comp] += x[d*dof+comp]*B[d*dof+comp];
460: }
461: }
462: PetscFERestoreTabulation(fem, 1, xir, &B, NULL, NULL);
463: } else {
464: for (comp = 0; comp < dof; ++comp)
465: a[p*dof+comp] = x[0*dof+comp]*(1 - xir[0])*(1 - xir[1]) + x[1*dof+comp]*xir[0]*(1 - xir[1]) + x[2*dof+comp]*xir[0]*xir[1] + x[3*dof+comp]*(1 - xir[0])*xir[1];
466: }
467: VecRestoreArray(ref, &xi);
468: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
469: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
470: }
471: VecRestoreArray(v, &a);
472: VecRestoreArrayRead(ctx->coords, &coords);
474: SNESDestroy(&snes);
475: VecDestroy(&r);
476: VecDestroy(&ref);
477: VecDestroy(&real);
478: MatDestroy(&J);
479: return(0);
480: }
482: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
483: {
484: const PetscScalar *vertices = (const PetscScalar*) ctx;
485: const PetscScalar x0 = vertices[0];
486: const PetscScalar y0 = vertices[1];
487: const PetscScalar z0 = vertices[2];
488: const PetscScalar x1 = vertices[9];
489: const PetscScalar y1 = vertices[10];
490: const PetscScalar z1 = vertices[11];
491: const PetscScalar x2 = vertices[6];
492: const PetscScalar y2 = vertices[7];
493: const PetscScalar z2 = vertices[8];
494: const PetscScalar x3 = vertices[3];
495: const PetscScalar y3 = vertices[4];
496: const PetscScalar z3 = vertices[5];
497: const PetscScalar x4 = vertices[12];
498: const PetscScalar y4 = vertices[13];
499: const PetscScalar z4 = vertices[14];
500: const PetscScalar x5 = vertices[15];
501: const PetscScalar y5 = vertices[16];
502: const PetscScalar z5 = vertices[17];
503: const PetscScalar x6 = vertices[18];
504: const PetscScalar y6 = vertices[19];
505: const PetscScalar z6 = vertices[20];
506: const PetscScalar x7 = vertices[21];
507: const PetscScalar y7 = vertices[22];
508: const PetscScalar z7 = vertices[23];
509: const PetscScalar f_1 = x1 - x0;
510: const PetscScalar g_1 = y1 - y0;
511: const PetscScalar h_1 = z1 - z0;
512: const PetscScalar f_3 = x3 - x0;
513: const PetscScalar g_3 = y3 - y0;
514: const PetscScalar h_3 = z3 - z0;
515: const PetscScalar f_4 = x4 - x0;
516: const PetscScalar g_4 = y4 - y0;
517: const PetscScalar h_4 = z4 - z0;
518: const PetscScalar f_01 = x2 - x1 - x3 + x0;
519: const PetscScalar g_01 = y2 - y1 - y3 + y0;
520: const PetscScalar h_01 = z2 - z1 - z3 + z0;
521: const PetscScalar f_12 = x7 - x3 - x4 + x0;
522: const PetscScalar g_12 = y7 - y3 - y4 + y0;
523: const PetscScalar h_12 = z7 - z3 - z4 + z0;
524: const PetscScalar f_02 = x5 - x1 - x4 + x0;
525: const PetscScalar g_02 = y5 - y1 - y4 + y0;
526: const PetscScalar h_02 = z5 - z1 - z4 + z0;
527: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
528: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
529: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
530: const PetscScalar *ref;
531: PetscScalar *real;
532: PetscErrorCode ierr;
535: VecGetArrayRead(Xref, &ref);
536: VecGetArray(Xreal, &real);
537: {
538: const PetscScalar p0 = ref[0];
539: const PetscScalar p1 = ref[1];
540: const PetscScalar p2 = ref[2];
542: 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;
543: 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;
544: 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;
545: }
546: PetscLogFlops(114);
547: VecRestoreArrayRead(Xref, &ref);
548: VecRestoreArray(Xreal, &real);
549: return(0);
550: }
552: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
553: {
554: const PetscScalar *vertices = (const PetscScalar*) ctx;
555: const PetscScalar x0 = vertices[0];
556: const PetscScalar y0 = vertices[1];
557: const PetscScalar z0 = vertices[2];
558: const PetscScalar x1 = vertices[9];
559: const PetscScalar y1 = vertices[10];
560: const PetscScalar z1 = vertices[11];
561: const PetscScalar x2 = vertices[6];
562: const PetscScalar y2 = vertices[7];
563: const PetscScalar z2 = vertices[8];
564: const PetscScalar x3 = vertices[3];
565: const PetscScalar y3 = vertices[4];
566: const PetscScalar z3 = vertices[5];
567: const PetscScalar x4 = vertices[12];
568: const PetscScalar y4 = vertices[13];
569: const PetscScalar z4 = vertices[14];
570: const PetscScalar x5 = vertices[15];
571: const PetscScalar y5 = vertices[16];
572: const PetscScalar z5 = vertices[17];
573: const PetscScalar x6 = vertices[18];
574: const PetscScalar y6 = vertices[19];
575: const PetscScalar z6 = vertices[20];
576: const PetscScalar x7 = vertices[21];
577: const PetscScalar y7 = vertices[22];
578: const PetscScalar z7 = vertices[23];
579: const PetscScalar f_xy = x2 - x1 - x3 + x0;
580: const PetscScalar g_xy = y2 - y1 - y3 + y0;
581: const PetscScalar h_xy = z2 - z1 - z3 + z0;
582: const PetscScalar f_yz = x7 - x3 - x4 + x0;
583: const PetscScalar g_yz = y7 - y3 - y4 + y0;
584: const PetscScalar h_yz = z7 - z3 - z4 + z0;
585: const PetscScalar f_xz = x5 - x1 - x4 + x0;
586: const PetscScalar g_xz = y5 - y1 - y4 + y0;
587: const PetscScalar h_xz = z5 - z1 - z4 + z0;
588: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
589: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
590: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
591: const PetscScalar *ref;
592: PetscErrorCode ierr;
595: VecGetArrayRead(Xref, &ref);
596: {
597: const PetscScalar x = ref[0];
598: const PetscScalar y = ref[1];
599: const PetscScalar z = ref[2];
600: const PetscInt rows[3] = {0, 1, 2};
601: PetscScalar values[9];
603: values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
604: values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
605: values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
606: values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
607: values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
608: values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
609: values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
610: values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
611: values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
613: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
614: }
615: PetscLogFlops(152);
616: VecRestoreArrayRead(Xref, &ref);
617: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
618: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
619: return(0);
620: }
622: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
623: {
624: DM dmCoord;
625: SNES snes;
626: KSP ksp;
627: PC pc;
628: Vec coordsLocal, r, ref, real;
629: Mat J;
630: const PetscScalar *coords;
631: PetscScalar *a;
632: PetscInt p;
636: DMGetCoordinatesLocal(dm, &coordsLocal);
637: DMGetCoordinateDM(dm, &dmCoord);
638: SNESCreate(PETSC_COMM_SELF, &snes);
639: SNESSetOptionsPrefix(snes, "hex_interp_");
640: VecCreate(PETSC_COMM_SELF, &r);
641: VecSetSizes(r, 3, 3);
642: VecSetType(r,dm->vectype);
643: VecDuplicate(r, &ref);
644: VecDuplicate(r, &real);
645: MatCreate(PETSC_COMM_SELF, &J);
646: MatSetSizes(J, 3, 3, 3, 3);
647: MatSetType(J, MATSEQDENSE);
648: MatSetUp(J);
649: SNESSetFunction(snes, r, HexMap_Private, NULL);
650: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
651: SNESGetKSP(snes, &ksp);
652: KSPGetPC(ksp, &pc);
653: PCSetType(pc, PCLU);
654: SNESSetFromOptions(snes);
656: VecGetArrayRead(ctx->coords, &coords);
657: VecGetArray(v, &a);
658: for (p = 0; p < ctx->n; ++p) {
659: PetscScalar *x = NULL, *vertices = NULL;
660: PetscScalar *xi;
661: PetscReal xir[3];
662: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
664: /* Can make this do all points at once */
665: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
666: if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
667: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
668: if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
669: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
670: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
671: VecGetArray(real, &xi);
672: xi[0] = coords[p*ctx->dim+0];
673: xi[1] = coords[p*ctx->dim+1];
674: xi[2] = coords[p*ctx->dim+2];
675: VecRestoreArray(real, &xi);
676: SNESSolve(snes, real, ref);
677: VecGetArray(ref, &xi);
678: xir[0] = PetscRealPart(xi[0]);
679: xir[1] = PetscRealPart(xi[1]);
680: xir[2] = PetscRealPart(xi[2]);
681: for (comp = 0; comp < ctx->dof; ++comp) {
682: a[p*ctx->dof+comp] =
683: x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
684: x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
685: x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
686: x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
687: x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
688: x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
689: x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
690: x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
691: }
692: VecRestoreArray(ref, &xi);
693: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
694: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
695: }
696: VecRestoreArray(v, &a);
697: VecRestoreArrayRead(ctx->coords, &coords);
699: SNESDestroy(&snes);
700: VecDestroy(&r);
701: VecDestroy(&ref);
702: VecDestroy(&real);
703: MatDestroy(&J);
704: return(0);
705: }
707: /*
708: Input Parameters:
709: + ctx - The DMInterpolationInfo context
710: . dm - The DM
711: - x - The local vector containing the field to be interpolated
713: Output Parameters:
714: . v - The vector containing the interpolated values
715: */
716: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
717: {
718: PetscInt dim, coneSize, n;
725: VecGetLocalSize(v, &n);
726: 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);
727: if (n) {
728: DMGetDimension(dm, &dim);
729: DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
730: if (dim == 2) {
731: if (coneSize == 3) {
732: DMInterpolate_Triangle_Private(ctx, dm, x, v);
733: } else if (coneSize == 4) {
734: DMInterpolate_Quad_Private(ctx, dm, x, v);
735: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
736: } else if (dim == 3) {
737: if (coneSize == 4) {
738: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
739: } else {
740: DMInterpolate_Hex_Private(ctx, dm, x, v);
741: }
742: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
743: }
744: return(0);
745: }
747: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
748: {
753: VecDestroy(&(*ctx)->coords);
754: PetscFree((*ctx)->points);
755: PetscFree((*ctx)->cells);
756: PetscFree(*ctx);
757: *ctx = NULL;
758: return(0);
759: }
761: /*@C
762: SNESMonitorFields - Monitors the residual for each field separately
764: Collective on SNES
766: Input Parameters:
767: + snes - the SNES context
768: . its - iteration number
769: . fgnorm - 2-norm of residual
770: - vf - PetscViewerAndFormat of type ASCII
772: Notes:
773: This routine prints the residual norm at each iteration.
775: Level: intermediate
777: .keywords: SNES, nonlinear, default, monitor, norm
778: .seealso: SNESMonitorSet(), SNESMonitorDefault()
779: @*/
780: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
781: {
782: PetscViewer viewer = vf->viewer;
783: Vec res;
784: DM dm;
785: PetscSection s;
786: const PetscScalar *r;
787: PetscReal *lnorms, *norms;
788: PetscInt numFields, f, pStart, pEnd, p;
789: PetscErrorCode ierr;
793: SNESGetFunction(snes, &res, 0, 0);
794: SNESGetDM(snes, &dm);
795: DMGetDefaultSection(dm, &s);
796: PetscSectionGetNumFields(s, &numFields);
797: PetscSectionGetChart(s, &pStart, &pEnd);
798: PetscCalloc2(numFields, &lnorms, numFields, &norms);
799: VecGetArrayRead(res, &r);
800: for (p = pStart; p < pEnd; ++p) {
801: for (f = 0; f < numFields; ++f) {
802: PetscInt fdof, foff, d;
804: PetscSectionGetFieldDof(s, p, f, &fdof);
805: PetscSectionGetFieldOffset(s, p, f, &foff);
806: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
807: }
808: }
809: VecRestoreArrayRead(res, &r);
810: MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
811: PetscViewerPushFormat(viewer,vf->format);
812: PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
813: PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
814: for (f = 0; f < numFields; ++f) {
815: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
816: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
817: }
818: PetscViewerASCIIPrintf(viewer, "]\n");
819: PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
820: PetscViewerPopFormat(viewer);
821: PetscFree2(lnorms, norms);
822: return(0);
823: }
825: /********************* Residual Computation **************************/
827: /*@
828: DMPlexSNESGetGeometryFEM - Return precomputed geometric data
830: Input Parameter:
831: . dm - The DM
833: Output Parameters:
834: . cellgeom - The values precomputed from cell geometry
836: Level: developer
838: .seealso: DMPlexSNESSetFunctionLocal()
839: @*/
840: PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
841: {
842: DMSNES dmsnes;
843: PetscObject obj;
848: DMGetDMSNES(dm, &dmsnes);
849: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);
850: if (!obj) {
851: Vec cellgeom;
853: DMPlexComputeGeometryFEM(dm, &cellgeom);
854: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);
855: VecDestroy(&cellgeom);
856: }
858: return(0);
859: }
861: /*@
862: DMPlexSNESGetGeometryFVM - Return precomputed geometric data
864: Input Parameter:
865: . dm - The DM
867: Output Parameters:
868: + facegeom - The values precomputed from face geometry
869: . cellgeom - The values precomputed from cell geometry
870: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
872: Level: developer
874: .seealso: DMPlexTSSetRHSFunctionLocal()
875: @*/
876: PetscErrorCode DMPlexSNESGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
877: {
878: DM plex;
883: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
884: DMPlexGetDataFVM(plex, NULL, cellgeom, facegeom, NULL);
885: if (minRadius) {DMPlexGetMinRadius(plex, minRadius);}
886: DMDestroy(&plex);
887: return(0);
888: }
890: /*@
891: DMPlexSNESGetGradientDM - Return gradient data layout
893: Input Parameters:
894: + dm - The DM
895: - fv - The PetscFV
897: Output Parameter:
898: . dmGrad - The layout for gradient values
900: Level: developer
902: .seealso: DMPlexSNESGetGeometryFVM()
903: @*/
904: PetscErrorCode DMPlexSNESGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
905: {
906: DM plex;
907: PetscBool computeGradients;
914: PetscFVGetComputeGradients(fv, &computeGradients);
915: if (!computeGradients) {*dmGrad = NULL; return(0);}
916: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
917: DMPlexGetDataFVM(plex, fv, NULL, NULL, dmGrad);
918: DMDestroy(&plex);
919: return(0);
920: }
922: /*@C
923: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
925: Input Parameters:
926: + dm - The DM
927: . cStart - The first cell to include
928: . cEnd - The first cell to exclude
929: . locX - A local vector with the solution fields
930: . locX_t - A local vector with solution field time derivatives, or NULL
931: - locA - A local vector with auxiliary fields, or NULL
933: Output Parameters:
934: + u - The field coefficients
935: . u_t - The fields derivative coefficients
936: - a - The auxiliary field coefficients
938: Level: developer
940: .seealso: DMPlexGetFaceFields()
941: @*/
942: PetscErrorCode DMPlexGetCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
943: {
944: DM dmAux;
945: PetscSection section, sectionAux;
946: PetscDS prob;
947: PetscInt numCells = cEnd - cStart, totDim, totDimAux, c;
958: DMGetDefaultSection(dm, §ion);
959: DMGetDS(dm, &prob);
960: PetscDSGetTotalDimension(prob, &totDim);
961: if (locA) {
962: PetscDS probAux;
964: VecGetDM(locA, &dmAux);
965: DMGetDefaultSection(dmAux, §ionAux);
966: DMGetDS(dmAux, &probAux);
967: PetscDSGetTotalDimension(probAux, &totDimAux);
968: }
969: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u);
970: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, u_t);} else {*u_t = NULL;}
971: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, MPIU_SCALAR, a);} else {*a = NULL;}
972: for (c = cStart; c < cEnd; ++c) {
973: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
974: PetscInt i;
976: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
977: for (i = 0; i < totDim; ++i) ul[(c-cStart)*totDim+i] = x[i];
978: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
979: if (locX_t) {
980: DMPlexVecGetClosure(dm, section, locX_t, c, NULL, &x_t);
981: for (i = 0; i < totDim; ++i) ul_t[(c-cStart)*totDim+i] = x_t[i];
982: DMPlexVecRestoreClosure(dm, section, locX_t, c, NULL, &x_t);
983: }
984: if (locA) {
985: DM dmAuxPlex;
987: DMSNESConvertPlex(dmAux, &dmAuxPlex, PETSC_FALSE);
988: DMPlexVecGetClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
989: for (i = 0; i < totDimAux; ++i) al[(c-cStart)*totDimAux+i] = x[i];
990: DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, locA, c, NULL, &x);
991: DMDestroy(&dmAuxPlex);
992: }
993: }
994: return(0);
995: }
997: /*@C
998: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
1000: Input Parameters:
1001: + dm - The DM
1002: . cStart - The first cell to include
1003: . cEnd - The first cell to exclude
1004: . locX - A local vector with the solution fields
1005: . locX_t - A local vector with solution field time derivatives, or NULL
1006: - locA - A local vector with auxiliary fields, or NULL
1008: Output Parameters:
1009: + u - The field coefficients
1010: . u_t - The fields derivative coefficients
1011: - a - The auxiliary field coefficients
1013: Level: developer
1015: .seealso: DMPlexGetFaceFields()
1016: @*/
1017: PetscErrorCode DMPlexRestoreCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
1018: {
1022: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u);
1023: if (locX_t) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, u_t);}
1024: if (locA) {DMRestoreWorkArray(dm, 0, MPIU_SCALAR, a);}
1025: return(0);
1026: }
1028: /*@C
1029: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
1031: Input Parameters:
1032: + dm - The DM
1033: . fStart - The first face to include
1034: . fEnd - The first face to exclude
1035: . locX - A local vector with the solution fields
1036: . locX_t - A local vector with solution field time derivatives, or NULL
1037: . faceGeometry - A local vector with face geometry
1038: . cellGeometry - A local vector with cell geometry
1039: - locaGrad - A local vector with field gradients, or NULL
1041: Output Parameters:
1042: + Nface - The number of faces with field values
1043: . uL - The field values at the left side of the face
1044: - uR - The field values at the right side of the face
1046: Level: developer
1048: .seealso: DMPlexGetCellFields()
1049: @*/
1050: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
1051: {
1052: DM dmFace, dmCell, dmGrad = NULL;
1053: PetscSection section;
1054: PetscDS prob;
1055: DMLabel ghostLabel;
1056: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
1057: PetscBool *isFE;
1058: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
1059: PetscErrorCode ierr;
1070: DMGetDimension(dm, &dim);
1071: DMGetDS(dm, &prob);
1072: DMGetDefaultSection(dm, §ion);
1073: PetscDSGetNumFields(prob, &Nf);
1074: PetscDSGetTotalComponents(prob, &Nc);
1075: PetscMalloc1(Nf, &isFE);
1076: for (f = 0; f < Nf; ++f) {
1077: PetscObject obj;
1078: PetscClassId id;
1080: PetscDSGetDiscretization(prob, f, &obj);
1081: PetscObjectGetClassId(obj, &id);
1082: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
1083: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
1084: else {isFE[f] = PETSC_FALSE;}
1085: }
1086: DMGetLabel(dm, "ghost", &ghostLabel);
1087: VecGetArrayRead(locX, &x);
1088: VecGetDM(faceGeometry, &dmFace);
1089: VecGetArrayRead(faceGeometry, &facegeom);
1090: VecGetDM(cellGeometry, &dmCell);
1091: VecGetArrayRead(cellGeometry, &cellgeom);
1092: if (locGrad) {
1093: VecGetDM(locGrad, &dmGrad);
1094: VecGetArrayRead(locGrad, &lgrad);
1095: }
1096: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uL);
1097: DMGetWorkArray(dm, numFaces*Nc, MPIU_SCALAR, uR);
1098: /* Right now just eat the extra work for FE (could make a cell loop) */
1099: for (face = fStart, iface = 0; face < fEnd; ++face) {
1100: const PetscInt *cells;
1101: PetscFVFaceGeom *fg;
1102: PetscFVCellGeom *cgL, *cgR;
1103: PetscScalar *xL, *xR, *gL, *gR;
1104: PetscScalar *uLl = *uL, *uRl = *uR;
1105: PetscInt ghost, nsupp, nchild;
1107: DMLabelGetValue(ghostLabel, face, &ghost);
1108: DMPlexGetSupportSize(dm, face, &nsupp);
1109: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
1110: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
1111: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1112: DMPlexGetSupport(dm, face, &cells);
1113: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1114: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1115: for (f = 0; f < Nf; ++f) {
1116: PetscInt off;
1118: PetscDSGetComponentOffset(prob, f, &off);
1119: if (isFE[f]) {
1120: const PetscInt *cone;
1121: PetscInt comp, coneSizeL, coneSizeR, faceLocL, faceLocR, ldof, rdof, d;
1123: xL = xR = NULL;
1124: PetscSectionGetFieldComponents(section, f, &comp);
1125: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1126: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1127: DMPlexGetCone(dm, cells[0], &cone);
1128: DMPlexGetConeSize(dm, cells[0], &coneSizeL);
1129: for (faceLocL = 0; faceLocL < coneSizeL; ++faceLocL) if (cone[faceLocL] == face) break;
1130: DMPlexGetCone(dm, cells[1], &cone);
1131: DMPlexGetConeSize(dm, cells[1], &coneSizeR);
1132: for (faceLocR = 0; faceLocR < coneSizeR; ++faceLocR) if (cone[faceLocR] == face) break;
1133: if (faceLocL == coneSizeL && faceLocR == coneSizeR) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d or cell %d", face, cells[0], cells[1]);
1134: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
1135: /* TODO: this is a hack that might not be right for nonconforming */
1136: if (faceLocL < coneSizeL) {
1137: EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
1138: if (rdof == ldof && faceLocR < coneSizeR) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
1139: else {for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
1140: }
1141: else {
1142: EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);
1143: PetscSectionGetFieldComponents(section, f, &comp);
1144: for(d = 0; d < comp; ++d) uLl[iface*Nc+off+d] = uRl[iface*Nc+off+d];
1145: }
1146: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1147: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1148: } else {
1149: PetscFV fv;
1150: PetscInt numComp, c;
1152: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
1153: PetscFVGetNumComponents(fv, &numComp);
1154: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
1155: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
1156: if (dmGrad) {
1157: PetscReal dxL[3], dxR[3];
1159: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
1160: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
1161: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
1162: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
1163: for (c = 0; c < numComp; ++c) {
1164: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
1165: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
1166: }
1167: } else {
1168: for (c = 0; c < numComp; ++c) {
1169: uLl[iface*Nc+off+c] = xL[c];
1170: uRl[iface*Nc+off+c] = xR[c];
1171: }
1172: }
1173: }
1174: }
1175: ++iface;
1176: }
1177: *Nface = iface;
1178: VecRestoreArrayRead(locX, &x);
1179: VecRestoreArrayRead(faceGeometry, &facegeom);
1180: VecRestoreArrayRead(cellGeometry, &cellgeom);
1181: if (locGrad) {
1182: VecRestoreArrayRead(locGrad, &lgrad);
1183: }
1184: PetscFree(isFE);
1185: return(0);
1186: }
1188: /*@C
1189: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
1191: Input Parameters:
1192: + dm - The DM
1193: . fStart - The first face to include
1194: . fEnd - The first face to exclude
1195: . locX - A local vector with the solution fields
1196: . locX_t - A local vector with solution field time derivatives, or NULL
1197: . faceGeometry - A local vector with face geometry
1198: . cellGeometry - A local vector with cell geometry
1199: - locaGrad - A local vector with field gradients, or NULL
1201: Output Parameters:
1202: + Nface - The number of faces with field values
1203: . uL - The field values at the left side of the face
1204: - uR - The field values at the right side of the face
1206: Level: developer
1208: .seealso: DMPlexGetFaceFields()
1209: @*/
1210: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscInt *Nface, PetscScalar **uL, PetscScalar **uR)
1211: {
1215: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uL);
1216: DMRestoreWorkArray(dm, 0, MPIU_SCALAR, uR);
1217: return(0);
1218: }
1220: /*@C
1221: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
1223: Input Parameters:
1224: + dm - The DM
1225: . fStart - The first face to include
1226: . fEnd - The first face to exclude
1227: . faceGeometry - A local vector with face geometry
1228: - cellGeometry - A local vector with cell geometry
1230: Output Parameters:
1231: + Nface - The number of faces with field values
1232: . fgeom - The extract the face centroid and normal
1233: - vol - The cell volume
1235: Level: developer
1237: .seealso: DMPlexGetCellFields()
1238: @*/
1239: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
1240: {
1241: DM dmFace, dmCell;
1242: DMLabel ghostLabel;
1243: const PetscScalar *facegeom, *cellgeom;
1244: PetscInt dim, numFaces = fEnd - fStart, iface, face;
1245: PetscErrorCode ierr;
1253: DMGetDimension(dm, &dim);
1254: DMGetLabel(dm, "ghost", &ghostLabel);
1255: VecGetDM(faceGeometry, &dmFace);
1256: VecGetArrayRead(faceGeometry, &facegeom);
1257: VecGetDM(cellGeometry, &dmCell);
1258: VecGetArrayRead(cellGeometry, &cellgeom);
1259: PetscMalloc1(numFaces, fgeom);
1260: DMGetWorkArray(dm, numFaces*2, MPIU_SCALAR, vol);
1261: for (face = fStart, iface = 0; face < fEnd; ++face) {
1262: const PetscInt *cells;
1263: PetscFVFaceGeom *fg;
1264: PetscFVCellGeom *cgL, *cgR;
1265: PetscFVFaceGeom *fgeoml = *fgeom;
1266: PetscReal *voll = *vol;
1267: PetscInt ghost, d, nchild, nsupp;
1269: DMLabelGetValue(ghostLabel, face, &ghost);
1270: DMPlexGetSupportSize(dm, face, &nsupp);
1271: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
1272: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
1273: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1274: DMPlexGetSupport(dm, face, &cells);
1275: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1276: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1277: for (d = 0; d < dim; ++d) {
1278: fgeoml[iface].centroid[d] = fg->centroid[d];
1279: fgeoml[iface].normal[d] = fg->normal[d];
1280: }
1281: voll[iface*2+0] = cgL->volume;
1282: voll[iface*2+1] = cgR->volume;
1283: ++iface;
1284: }
1285: *Nface = iface;
1286: VecRestoreArrayRead(faceGeometry, &facegeom);
1287: VecRestoreArrayRead(cellGeometry, &cellgeom);
1288: return(0);
1289: }
1291: /*@C
1292: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
1294: Input Parameters:
1295: + dm - The DM
1296: . fStart - The first face to include
1297: . fEnd - The first face to exclude
1298: . faceGeometry - A local vector with face geometry
1299: - cellGeometry - A local vector with cell geometry
1301: Output Parameters:
1302: + Nface - The number of faces with field values
1303: . fgeom - The extract the face centroid and normal
1304: - vol - The cell volume
1306: Level: developer
1308: .seealso: DMPlexGetFaceFields()
1309: @*/
1310: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscInt *Nface, PetscFVFaceGeom **fgeom, PetscReal **vol)
1311: {
1315: PetscFree(*fgeom);
1316: DMRestoreWorkArray(dm, 0, MPIU_REAL, vol);
1317: return(0);
1318: }
1320: PetscErrorCode DMPlexComputeBdResidual_Single_Internal(DM dm, PetscReal t, DMLabel label, PetscInt numValues, const PetscInt values[], PetscInt field, Vec locX, Vec locX_t, Vec locF)
1321: {
1322: DM_Plex *mesh = (DM_Plex *) dm->data;
1323: DM plex = NULL;
1324: DMLabel depth;
1325: PetscDS prob, probAux = NULL;
1326: PetscSection section, sectionAux = NULL;
1327: Vec locA = NULL;
1328: PetscFEFaceGeom *fgeom;
1329: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemVec = NULL;
1330: PetscInt v;
1331: PetscInt dim, totDim, totDimAux = 0;
1332: PetscErrorCode ierr;
1335: DMGetDimension(dm, &dim);
1336: DMPlexGetDepthLabel(dm, &depth);
1337: DMGetDefaultSection(dm, §ion);
1338: DMGetDS(dm, &prob);
1339: PetscDSGetTotalDimension(prob, &totDim);
1340: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1341: if (locA) {
1342: DM dmAux;
1344: VecGetDM(locA, &dmAux);
1345: DMGetDS(dmAux, &probAux);
1346: PetscDSGetTotalDimension(probAux, &totDimAux);
1347: DMConvert(dmAux, DMPLEX, &plex);
1348: DMGetDefaultSection(plex, §ionAux);
1349: }
1350: for (v = 0; v < numValues; ++v) {
1351: IS pointIS;
1352: const PetscInt *points;
1353: PetscInt numPoints, p, dep, numFaces, face;
1355: DMLabelGetStratumSize(label, values[v], &numPoints);
1356: DMLabelGetStratumIS(label, values[v], &pointIS);
1357: if (!pointIS) continue; /* No points with that id on this process */
1358: ISGetIndices(pointIS, &points);
1359: for (p = 0, numFaces = 0; p < numPoints; ++p) {
1360: DMLabelGetValue(depth, points[p], &dep);
1361: if (dep == dim-1) ++numFaces;
1362: }
1363: PetscMalloc5(numFaces*totDim, &u, locX_t ? numFaces*totDim : 0, &u_t, numFaces, &fgeom, numFaces*totDim, &elemVec, locA ? numFaces*totDimAux : 0, &a);
1364: for (p = 0, face = 0; p < numPoints; ++p) {
1365: const PetscInt point = points[p], *support, *cone;
1366: PetscScalar *x = NULL;
1367: PetscReal dummyJ[9], dummyDetJ;
1368: PetscInt i, coneSize, faceLoc;
1370: DMLabelGetValue(depth, points[p], &dep);
1371: if (dep != dim-1) continue;
1372: DMPlexGetSupport(dm, point, &support);
1373: DMPlexComputeCellGeometryFEM(dm, support[0], NULL, NULL, dummyJ, fgeom[face].invJ[0], &dummyDetJ);
1374: DMPlexComputeCellGeometryFEM(dm, point, NULL, fgeom[face].v0, fgeom[face].J, NULL, &fgeom[face].detJ);
1375: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, fgeom[face].n);
1376: if (fgeom[face].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %D", (double) fgeom[face].detJ, point);
1377: DMPlexGetConeSize(dm, support[0], &coneSize);
1378: DMPlexGetCone(dm, support[0], &cone);
1379: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
1380: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %D in cone of support[0] %D", point, support[0]);
1381: fgeom[face].face[0] = faceLoc;
1382: DMPlexVecGetClosure(dm, section, locX, support[0], NULL, &x);
1383: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
1384: DMPlexVecRestoreClosure(dm, section, locX, support[0], NULL, &x);
1385: if (locX_t) {
1386: DMPlexVecGetClosure(dm, section, locX_t, support[0], NULL, &x);
1387: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
1388: DMPlexVecRestoreClosure(dm, section, locX_t, support[0], NULL, &x);
1389: }
1390: if (locA) {
1391: DMPlexVecGetClosure(plex, sectionAux, locA, support[0], NULL, &x);
1392: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
1393: DMPlexVecRestoreClosure(plex, sectionAux, locA, support[0], NULL, &x);
1394: }
1395: ++face;
1396: }
1397: if (plex) {DMDestroy(&plex);}
1398: PetscMemzero(elemVec, numFaces*totDim * sizeof(PetscScalar));
1399: {
1400: PetscFE fe;
1401: PetscQuadrature q;
1402: PetscInt numQuadPoints, Nb;
1403: /* Conforming batches */
1404: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1405: /* Remainder */
1406: PetscInt Nr, offset;
1408: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1409: PetscFEGetFaceQuadrature(fe, &q);
1410: PetscFEGetDimension(fe, &Nb);
1411: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1412: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
1413: blockSize = Nb*numQuadPoints;
1414: batchSize = numBlocks * blockSize;
1415: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1416: numChunks = numFaces / (numBatches*batchSize);
1417: Ne = numChunks*numBatches*batchSize;
1418: Nr = numFaces % (numBatches*batchSize);
1419: offset = numFaces - Nr;
1420: PetscFEIntegrateBdResidual(fe, prob, field, Ne, fgeom, u, u_t, probAux, a, t, elemVec);
1421: PetscFEIntegrateBdResidual(fe, prob, field, Nr, &fgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, a ? &a[offset*totDimAux] : NULL, t, &elemVec[offset*totDim]);
1422: }
1423: for (p = 0, face = 0; p < numPoints; ++p) {
1424: const PetscInt point = points[p], *support;
1426: DMLabelGetValue(depth, point, &dep);
1427: if (dep != dim-1) continue;
1428: if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDim, &elemVec[face*totDim]);}
1429: DMPlexGetSupport(dm, point, &support);
1430: DMPlexVecSetClosure(dm, NULL, locF, support[0], &elemVec[face*totDim], ADD_ALL_VALUES);
1431: ++face;
1432: }
1433: ISRestoreIndices(pointIS, &points);
1434: ISDestroy(&pointIS);
1435: PetscFree5(u, u_t, fgeom, elemVec, a);
1436: }
1437: return(0);
1438: }
1440: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
1441: {
1442: PetscDS prob;
1443: PetscInt numBd, bd;
1447: DMGetDS(dm, &prob);
1448: PetscDSGetNumBoundary(prob, &numBd);
1449: for (bd = 0; bd < numBd; ++bd) {
1450: DMBoundaryConditionType type;
1451: const char *bdLabel;
1452: DMLabel label;
1453: const PetscInt *values;
1454: PetscInt field, numValues;
1455: PetscObject obj;
1456: PetscClassId id;
1458: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1459: PetscDSGetDiscretization(prob, field, &obj);
1460: PetscObjectGetClassId(obj, &id);
1461: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
1462: DMGetLabel(dm, bdLabel, &label);
1463: DMPlexComputeBdResidual_Single_Internal(dm, t, label, numValues, values, field, locX, locX_t, locF);
1464: }
1465: return(0);
1466: }
1468: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
1469: {
1470: DM_Plex *mesh = (DM_Plex *) dm->data;
1471: const char *name = "Residual";
1472: DM dmAux = NULL;
1473: DM dmGrad = NULL;
1474: DMLabel ghostLabel = NULL;
1475: PetscDS prob = NULL;
1476: PetscDS probAux = NULL;
1477: PetscSection section = NULL;
1478: PetscBool useFEM = PETSC_FALSE;
1479: PetscBool useFVM = PETSC_FALSE;
1480: PetscBool isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1481: PetscFV fvm = NULL;
1482: PetscFECellGeom *cgeomFEM = NULL;
1483: PetscScalar *cgeomScal;
1484: PetscFVCellGeom *cgeomFVM = NULL;
1485: PetscFVFaceGeom *fgeomFVM = NULL;
1486: Vec locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1487: PetscScalar *u = NULL, *u_t, *a, *uL, *uR;
1488: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1489: PetscErrorCode ierr;
1492: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1493: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1494: /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1495: /* FEM+FVM */
1496: /* 1: Get sizes from dm and dmAux */
1497: DMGetDefaultSection(dm, §ion);
1498: DMGetLabel(dm, "ghost", &ghostLabel);
1499: DMGetDS(dm, &prob);
1500: PetscDSGetNumFields(prob, &Nf);
1501: PetscDSGetTotalDimension(prob, &totDim);
1502: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1503: if (locA) {
1504: VecGetDM(locA, &dmAux);
1505: DMGetDS(dmAux, &probAux);
1506: PetscDSGetTotalDimension(probAux, &totDimAux);
1507: }
1508: /* 2: Get geometric data */
1509: for (f = 0; f < Nf; ++f) {
1510: PetscObject obj;
1511: PetscClassId id;
1512: PetscBool fimp;
1514: PetscDSGetImplicit(prob, f, &fimp);
1515: if (isImplicit != fimp) continue;
1516: PetscDSGetDiscretization(prob, f, &obj);
1517: PetscObjectGetClassId(obj, &id);
1518: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1519: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1520: }
1521: if (useFEM) {
1522: DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1523: VecGetArray(cellGeometryFEM, &cgeomScal);
1524: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1525: DM dmCell;
1526: PetscInt c;
1528: VecGetDM(cellGeometryFEM,&dmCell);
1529: PetscMalloc1(cEnd-cStart,&cgeomFEM);
1530: for (c = 0; c < cEnd - cStart; c++) {
1531: PetscScalar *thisgeom;
1533: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1534: cgeomFEM[c] = *((PetscFECellGeom *) thisgeom);
1535: }
1536: }
1537: else {
1538: cgeomFEM = (PetscFECellGeom *) cgeomScal;
1539: }
1540: }
1541: if (useFVM) {
1542: DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1543: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1544: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1545: /* Reconstruct and limit cell gradients */
1546: DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1547: if (dmGrad) {
1548: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1549: DMGetGlobalVector(dmGrad, &grad);
1550: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1551: /* Communicate gradient values */
1552: DMGetLocalVector(dmGrad, &locGrad);
1553: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1554: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1555: DMRestoreGlobalVector(dmGrad, &grad);
1556: }
1557: /* Handle non-essential (e.g. outflow) boundary values */
1558: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1559: }
1560: /* Loop over chunks */
1561: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1562: numChunks = 1;
1563: cellChunkSize = (cEnd - cStart)/numChunks;
1564: faceChunkSize = (fEnd - fStart)/numChunks;
1565: for (chunk = 0; chunk < numChunks; ++chunk) {
1566: PetscScalar *elemVec, *fluxL, *fluxR;
1567: PetscReal *vol;
1568: PetscFVFaceGeom *fgeom;
1569: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1570: PetscInt fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = 0, face;
1572: /* Extract field coefficients */
1573: if (useFEM) {
1574: DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1575: DMGetWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
1576: PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1577: }
1578: if (useFVM) {
1579: DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
1580: DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
1581: DMGetWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxL);
1582: DMGetWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxR);
1583: PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1584: PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1585: }
1586: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1587: /* Loop over fields */
1588: for (f = 0; f < Nf; ++f) {
1589: PetscObject obj;
1590: PetscClassId id;
1591: PetscBool fimp;
1592: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1594: PetscDSGetImplicit(prob, f, &fimp);
1595: if (isImplicit != fimp) continue;
1596: PetscDSGetDiscretization(prob, f, &obj);
1597: PetscObjectGetClassId(obj, &id);
1598: if (id == PETSCFE_CLASSID) {
1599: PetscFE fe = (PetscFE) obj;
1600: PetscQuadrature q;
1601: PetscInt Nq, Nb;
1603: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1605: PetscFEGetQuadrature(fe, &q);
1606: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1607: PetscFEGetDimension(fe, &Nb);
1608: blockSize = Nb*Nq;
1609: batchSize = numBlocks * blockSize;
1610: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1611: numChunks = numCells / (numBatches*batchSize);
1612: Ne = numChunks*numBatches*batchSize;
1613: Nr = numCells % (numBatches*batchSize);
1614: offset = numCells - Nr;
1615: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1616: /* 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) */
1617: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, t, elemVec);
1618: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
1619: } else if (id == PETSCFV_CLASSID) {
1620: PetscFV fv = (PetscFV) obj;
1622: Ne = numFaces;
1623: /* Riemann solve over faces (need fields at face centroids) */
1624: /* We need to evaluate FE fields at those coordinates */
1625: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1626: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1627: }
1628: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1629: PetscFree(cgeomFEM);
1630: }
1631: else {
1632: cgeomFEM = NULL;
1633: }
1634: if (cellGeometryFEM) {VecRestoreArray(cellGeometryFEM, &cgeomScal);}
1635: /* Loop over domain */
1636: if (useFEM) {
1637: /* Add elemVec to locX */
1638: for (cell = cS; cell < cE; ++cell) {
1639: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1640: if (ghostLabel) {
1641: PetscInt ghostVal;
1643: DMLabelGetValue(ghostLabel,cell,&ghostVal);
1644: if (ghostVal > 0) continue;
1645: }
1646: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_ALL_VALUES);
1647: }
1648: }
1649: if (useFVM) {
1650: PetscScalar *fa;
1651: PetscInt iface;
1653: VecGetArray(locF, &fa);
1654: for (f = 0; f < Nf; ++f) {
1655: PetscFV fv;
1656: PetscObject obj;
1657: PetscClassId id;
1658: PetscInt foff, pdim;
1660: PetscDSGetDiscretization(prob, f, &obj);
1661: PetscDSGetFieldOffset(prob, f, &foff);
1662: PetscObjectGetClassId(obj, &id);
1663: if (id != PETSCFV_CLASSID) continue;
1664: fv = (PetscFV) obj;
1665: PetscFVGetNumComponents(fv, &pdim);
1666: /* Accumulate fluxes to cells */
1667: for (face = fS, iface = 0; face < fE; ++face) {
1668: const PetscInt *cells;
1669: PetscScalar *fL = NULL, *fR = NULL;
1670: PetscInt ghost, d, nsupp, nchild;
1672: DMLabelGetValue(ghostLabel, face, &ghost);
1673: DMPlexGetSupportSize(dm, face, &nsupp);
1674: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
1675: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
1676: DMPlexGetSupport(dm, face, &cells);
1677: DMLabelGetValue(ghostLabel,cells[0],&ghost);
1678: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, cells[0], f, fa, &fL);}
1679: DMLabelGetValue(ghostLabel,cells[1],&ghost);
1680: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, cells[1], f, fa, &fR);}
1681: for (d = 0; d < pdim; ++d) {
1682: if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1683: if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1684: }
1685: ++iface;
1686: }
1687: }
1688: VecRestoreArray(locF, &fa);
1689: }
1690: /* Handle time derivative */
1691: if (locX_t) {
1692: PetscScalar *x_t, *fa;
1694: VecGetArray(locF, &fa);
1695: VecGetArray(locX_t, &x_t);
1696: for (f = 0; f < Nf; ++f) {
1697: PetscFV fv;
1698: PetscObject obj;
1699: PetscClassId id;
1700: PetscInt pdim, d;
1702: PetscDSGetDiscretization(prob, f, &obj);
1703: PetscObjectGetClassId(obj, &id);
1704: if (id != PETSCFV_CLASSID) continue;
1705: fv = (PetscFV) obj;
1706: PetscFVGetNumComponents(fv, &pdim);
1707: for (cell = cS; cell < cE; ++cell) {
1708: PetscScalar *u_t, *r;
1710: if (ghostLabel) {
1711: PetscInt ghostVal;
1713: DMLabelGetValue(ghostLabel,cell,&ghostVal);
1714: if (ghostVal > 0) continue;
1715: }
1716: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1717: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1718: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1719: }
1720: }
1721: VecRestoreArray(locX_t, &x_t);
1722: VecRestoreArray(locF, &fa);
1723: }
1724: if (useFEM) {
1725: DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1726: DMRestoreWorkArray(dm, numCells*totDim, MPIU_SCALAR, &elemVec);
1727: }
1728: if (useFVM) {
1729: DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
1730: DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
1731: DMRestoreWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxL);
1732: DMRestoreWorkArray(dm, numFaces*totDim, MPIU_SCALAR, &fluxR);
1733: if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1734: }
1735: }
1737: if (useFEM) {DMPlexComputeBdResidual_Internal(dm, locX, locX_t, t, locF, user);}
1739: /* FEM */
1740: /* 1: Get sizes from dm and dmAux */
1741: /* 2: Get geometric data */
1742: /* 3: Handle boundary values */
1743: /* 4: Loop over domain */
1744: /* Extract coefficients */
1745: /* Loop over fields */
1746: /* Set tiling for FE*/
1747: /* Integrate FE residual to get elemVec */
1748: /* Loop over subdomain */
1749: /* Loop over quad points */
1750: /* Transform coords to real space */
1751: /* Evaluate field and aux fields at point */
1752: /* Evaluate residual at point */
1753: /* Transform residual to real space */
1754: /* Add residual to elemVec */
1755: /* Loop over domain */
1756: /* Add elemVec to locX */
1758: /* FVM */
1759: /* Get geometric data */
1760: /* If using gradients */
1761: /* Compute gradient data */
1762: /* Loop over domain faces */
1763: /* Count computational faces */
1764: /* Reconstruct cell gradient */
1765: /* Loop over domain cells */
1766: /* Limit cell gradients */
1767: /* Handle boundary values */
1768: /* Loop over domain faces */
1769: /* Read out field, centroid, normal, volume for each side of face */
1770: /* Riemann solve over faces */
1771: /* Loop over domain faces */
1772: /* Accumulate fluxes to cells */
1773: /* TODO Change printFEM to printDisc here */
1774: if (mesh->printFEM) {
1775: Vec locFbc;
1776: PetscInt pStart, pEnd, p, maxDof;
1777: PetscScalar *zeroes;
1779: VecDuplicate(locF,&locFbc);
1780: VecCopy(locF,locFbc);
1781: PetscSectionGetChart(section,&pStart,&pEnd);
1782: PetscSectionGetMaxDof(section,&maxDof);
1783: PetscCalloc1(maxDof,&zeroes);
1784: for (p = pStart; p < pEnd; p++) {
1785: VecSetValuesSection(locFbc,section,p,zeroes,INSERT_BC_VALUES);
1786: }
1787: PetscFree(zeroes);
1788: DMPrintLocalVec(dm, name, mesh->printTol, locFbc);
1789: VecDestroy(&locFbc);
1790: }
1791: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1792: return(0);
1793: }
1795: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, PetscReal t, Vec F, void *user)
1796: {
1797: DM dmCh, dmAux;
1798: Vec A, cellgeom;
1799: PetscDS prob, probCh, probAux = NULL;
1800: PetscQuadrature q;
1801: PetscSection section, sectionAux;
1802: PetscFECellGeom *cgeom = NULL;
1803: PetscScalar *cgeomScal;
1804: PetscScalar *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1805: PetscInt dim, Nf, f, numCells, cStart, cEnd, c;
1806: PetscInt totDim, totDimAux = 0, diffCell = 0;
1807: PetscErrorCode ierr;
1810: DMGetDimension(dm, &dim);
1811: DMGetDefaultSection(dm, §ion);
1812: DMGetDS(dm, &prob);
1813: PetscDSGetTotalDimension(prob, &totDim);
1814: PetscSectionGetNumFields(section, &Nf);
1815: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1816: numCells = cEnd - cStart;
1817: PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
1818: DMGetDS(dmCh, &probCh);
1819: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1820: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1821: if (dmAux) {
1822: DMGetDefaultSection(dmAux, §ionAux);
1823: DMGetDS(dmAux, &probAux);
1824: PetscDSGetTotalDimension(probAux, &totDimAux);
1825: }
1826: VecSet(F, 0.0);
1827: PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
1828: PetscMalloc1(numCells*totDim,&elemVecCh);
1829: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1830: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1831: VecGetArray(cellgeom, &cgeomScal);
1832: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1833: DM dmCell;
1835: VecGetDM(cellgeom,&dmCell);
1836: PetscMalloc1(cEnd-cStart,&cgeom);
1837: for (c = 0; c < cEnd - cStart; c++) {
1838: PetscScalar *thisgeom;
1840: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1841: cgeom[c] = *((PetscFECellGeom *) thisgeom);
1842: }
1843: }
1844: else {
1845: cgeom = (PetscFECellGeom *) cgeomScal;
1846: }
1847: for (c = cStart; c < cEnd; ++c) {
1848: PetscScalar *x = NULL, *x_t = NULL;
1849: PetscInt i;
1851: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1852: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1853: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1854: if (X_t) {
1855: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1856: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1857: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1858: }
1859: if (dmAux) {
1860: DM dmAuxPlex;
1862: DMSNESConvertPlex(dmAux,&dmAuxPlex, PETSC_FALSE);
1863: DMPlexVecGetClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
1864: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1865: DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
1866: DMDestroy(&dmAuxPlex);
1867: }
1868: }
1869: for (f = 0; f < Nf; ++f) {
1870: PetscFE fe, feCh;
1871: PetscInt numQuadPoints, Nb;
1872: /* Conforming batches */
1873: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1874: /* Remainder */
1875: PetscInt Nr, offset;
1877: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1878: PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
1879: PetscFEGetQuadrature(fe, &q);
1880: PetscFEGetDimension(fe, &Nb);
1881: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1882: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
1883: blockSize = Nb*numQuadPoints;
1884: batchSize = numBlocks * blockSize;
1885: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1886: numChunks = numCells / (numBatches*batchSize);
1887: Ne = numChunks*numBatches*batchSize;
1888: Nr = numCells % (numBatches*batchSize);
1889: offset = numCells - Nr;
1890: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, t, elemVec);
1891: PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, t, elemVecCh);
1892: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVec[offset*totDim]);
1893: PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, &elemVecCh[offset*totDim]);
1894: }
1895: for (c = cStart; c < cEnd; ++c) {
1896: PetscBool diff = PETSC_FALSE;
1897: PetscInt d;
1899: for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1900: if (diff) {
1901: PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
1902: DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
1903: DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
1904: ++diffCell;
1905: }
1906: if (diffCell > 9) break;
1907: DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_ALL_VALUES);
1908: }
1909: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1910: PetscFree(cgeom);
1911: }
1912: else {
1913: cgeom = NULL;
1914: }
1915: VecRestoreArray(cellgeom, &cgeomScal);
1916: PetscFree3(u,u_t,elemVec);
1917: PetscFree(elemVecCh);
1918: if (dmAux) {PetscFree(a);}
1919: return(0);
1920: }
1922: /*@
1923: DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1925: Input Parameters:
1926: + dm - The mesh
1927: . X - Local solution
1928: - user - The user context
1930: Output Parameter:
1931: . F - Local output vector
1933: Level: developer
1935: .seealso: DMPlexComputeJacobianActionFEM()
1936: @*/
1937: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1938: {
1939: PetscObject check;
1940: PetscInt cStart, cEnd, cEndInterior;
1941: DM plex;
1945: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1946: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
1947: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
1948: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1949: /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1950: PetscObjectQuery((PetscObject) plex, "dmCh", &check);
1951: if (check) {DMPlexComputeResidualFEM_Check_Internal(plex, X, NULL, 0.0, F, user);}
1952: else {DMPlexComputeResidual_Internal(plex, cStart, cEnd, PETSC_MIN_REAL, X, NULL, 0.0, F, user);}
1953: DMDestroy(&plex);
1954: return(0);
1955: }
1957: /*@
1958: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
1960: Input Parameters:
1961: + dm - The mesh
1962: - user - The user context
1964: Output Parameter:
1965: . X - Local solution
1967: Level: developer
1969: .seealso: DMPlexComputeJacobianActionFEM()
1970: @*/
1971: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1972: {
1973: DM plex;
1977: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1978: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1979: DMDestroy(&plex);
1980: return(0);
1981: }
1983: PetscErrorCode DMPlexComputeBdJacobian_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
1984: {
1985: DM_Plex *mesh = (DM_Plex *) dm->data;
1986: DM dmAux = NULL, plex = NULL;
1987: PetscSection section, globalSection, subSection, sectionAux = NULL;
1988: PetscDS prob, probAux = NULL;
1989: DMLabel depth;
1990: Vec locA = NULL;
1991: PetscFEFaceGeom *fgeom;
1992: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL;
1993: PetscInt dim, totDim, totDimAux, numBd, bd, Nf;
1994: PetscBool isMatISP;
1995: PetscErrorCode ierr;
1998: DMGetDimension(dm, &dim);
1999: DMGetDefaultSection(dm, §ion);
2000: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
2001: DMGetDefaultGlobalSection(dm, &globalSection);
2002: if (isMatISP) {
2003: DMPlexGetSubdomainSection(dm, &subSection);
2004: }
2005: DMPlexGetDepthLabel(dm, &depth);
2006: DMGetDS(dm, &prob);
2007: PetscDSGetNumFields(prob, &Nf);
2008: PetscDSGetTotalDimension(prob, &totDim);
2009: PetscDSGetNumBoundary(prob, &numBd);
2010: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
2011: if (locA) {
2012: VecGetDM(locA, &dmAux);
2013: DMConvert(dmAux, DMPLEX, &plex);
2014: DMGetDefaultSection(plex, §ionAux);
2015: DMGetDS(dmAux, &probAux);
2016: PetscDSGetTotalDimension(probAux, &totDimAux);
2017: }
2018: for (bd = 0; bd < numBd; ++bd) {
2019: DMBoundaryConditionType type;
2020: const char *bdLabel;
2021: DMLabel label;
2022: IS pointIS;
2023: const PetscInt *points;
2024: const PetscInt *values;
2025: PetscInt fieldI, fieldJ, numValues, v, numPoints, p, dep, numFaces, face;
2026: PetscObject obj;
2027: PetscClassId id;
2029: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &fieldI, NULL, NULL, NULL, &numValues, &values, NULL);
2030: PetscDSGetDiscretization(prob, fieldI, &obj);
2031: PetscObjectGetClassId(obj, &id);
2032: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
2033: DMGetLabel(dm, bdLabel, &label);
2034: for (v = 0; v < numValues; ++v) {
2035: DMLabelGetStratumSize(label, values[v], &numPoints);
2036: DMLabelGetStratumIS(label, values[v], &pointIS);
2037: if (!pointIS) continue; /* No points with that id on this process */
2038: ISGetIndices(pointIS, &points);
2039: for (p = 0, numFaces = 0; p < numPoints; ++p) {
2040: DMLabelGetValue(depth, points[p], &dep);
2041: if (dep == dim-1) ++numFaces;
2042: }
2043: PetscMalloc4(numFaces*totDim,&u,locX_t ? numFaces*totDim : 0,&u_t,numFaces,&fgeom,numFaces*totDim*totDim,&elemMat);
2044: if (locA) {PetscMalloc1(numFaces*totDimAux,&a);}
2045: for (p = 0, face = 0; p < numPoints; ++p) {
2046: const PetscInt point = points[p], *support, *cone;
2047: PetscScalar *x = NULL;
2048: PetscReal dummyJ[9], dummyDetJ;
2049: PetscInt i, coneSize, faceLoc;
2051: DMLabelGetValue(depth, points[p], &dep);
2052: if (dep != dim-1) continue;
2053: DMPlexGetSupport(dm, point, &support);
2054: DMPlexComputeCellGeometryFEM(dm, support[0], NULL, NULL, dummyJ, fgeom[face].invJ[0], &dummyDetJ);
2055: DMPlexComputeCellGeometryFEM(dm, point, NULL, fgeom[face].v0, fgeom[face].J, NULL, &fgeom[face].detJ);
2056: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, fgeom[face].n);
2057: if (fgeom[face].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", (double)fgeom[face].detJ, point);
2058: DMPlexGetConeSize(dm, support[0], &coneSize);
2059: DMPlexGetCone(dm, support[0], &cone);
2060: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
2061: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of support[0] %d", point, support[0]);
2062: fgeom[face].face[0] = faceLoc;
2063: DMPlexVecGetClosure(dm, section, locX, support[0], NULL, &x);
2064: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2065: DMPlexVecRestoreClosure(dm, section, locX, support[0], NULL, &x);
2066: if (locX_t) {
2067: DMPlexVecGetClosure(dm, section, locX_t, support[0], NULL, &x);
2068: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
2069: DMPlexVecRestoreClosure(dm, section, locX_t, support[0], NULL, &x);
2070: }
2071: if (locA) {
2072: DMPlexVecGetClosure(plex, sectionAux, locA, support[0], NULL, &x);
2073: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
2074: DMPlexVecRestoreClosure(plex, sectionAux, locA, support[0], NULL, &x);
2075: }
2076: ++face;
2077: }
2078: PetscMemzero(elemMat, numFaces*totDim*totDim * sizeof(PetscScalar));
2079: {
2080: PetscFE fe;
2081: PetscQuadrature q;
2082: PetscInt numQuadPoints, Nb;
2083: /* Conforming batches */
2084: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2085: /* Remainder */
2086: PetscInt Nr, offset;
2088: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2089: PetscFEGetFaceQuadrature(fe, &q);
2090: PetscFEGetDimension(fe, &Nb);
2091: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2092: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
2093: blockSize = Nb*numQuadPoints;
2094: batchSize = numBlocks * blockSize;
2095: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2096: numChunks = numFaces / (numBatches*batchSize);
2097: Ne = numChunks*numBatches*batchSize;
2098: Nr = numFaces % (numBatches*batchSize);
2099: offset = numFaces - Nr;
2100: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2101: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, fgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2102: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &fgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, a ? &a[offset*totDimAux] : NULL, t, X_tShift, &elemMat[offset*totDim*totDim]);
2103: }
2104: }
2105: for (p = 0, face = 0; p < numPoints; ++p) {
2106: const PetscInt point = points[p], *support;
2108: DMLabelGetValue(depth, point, &dep);
2109: if (dep != dim-1) continue;
2110: if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDim, totDim, &elemMat[face*totDim*totDim]);}
2111: DMPlexGetSupport(dm, point, &support);
2112: if (!isMatISP) {
2113: DMPlexMatSetClosure(dm, section, globalSection, JacP, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
2114: } else {
2115: Mat lJ;
2117: MatISGetLocalMat(JacP,&lJ);
2118: DMPlexMatSetClosure(dm, section, subSection, lJ, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
2119: }
2120: ++face;
2121: }
2122: ISRestoreIndices(pointIS, &points);
2123: ISDestroy(&pointIS);
2124: PetscFree4(u,u_t,fgeom,elemMat);
2125: if (locA) {PetscFree(a);}
2126: }
2127: }
2128: if (plex) {DMDestroy(&plex);}
2129: return(0);
2130: }
2132: 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)
2133: {
2134: DM_Plex *mesh = (DM_Plex *) dm->data;
2135: const char *name = "Jacobian";
2136: DM dmAux, plex;
2137: Vec A, cellgeom;
2138: PetscDS prob, probAux = NULL;
2139: PetscSection section, globalSection, subSection, sectionAux;
2140: PetscFECellGeom *cgeom = NULL;
2141: PetscScalar *cgeomScal;
2142: PetscScalar *elemMat, *elemMatP, *elemMatD, *u, *u_t, *a = NULL;
2143: PetscInt dim, Nf, fieldI, fieldJ, numCells, c;
2144: PetscInt totDim, totDimAux;
2145: PetscBool isMatIS, isMatISP, isShell, hasJac, hasPrec, hasDyn, hasFV = PETSC_FALSE;
2146: PetscErrorCode ierr;
2149: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2150: DMGetDimension(dm, &dim);
2151: DMGetDefaultSection(dm, §ion);
2152: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
2153: DMGetDefaultGlobalSection(dm, &globalSection);
2154: if (isMatISP) {
2155: DMPlexGetSubdomainSection(dm, &subSection);
2156: }
2157: DMGetDS(dm, &prob);
2158: PetscDSGetTotalDimension(prob, &totDim);
2159: PetscDSHasJacobian(prob, &hasJac);
2160: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2161: PetscDSHasDynamicJacobian(prob, &hasDyn);
2162: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2163: PetscSectionGetNumFields(section, &Nf);
2164: numCells = cEnd - cStart;
2165: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2166: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2167: if (dmAux) {
2168: DMConvert(dmAux, DMPLEX, &plex);
2169: DMGetDefaultSection(plex, §ionAux);
2170: DMGetDS(dmAux, &probAux);
2171: PetscDSGetTotalDimension(probAux, &totDimAux);
2172: }
2173: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
2174: MatZeroEntries(JacP);
2175: PetscMalloc5(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,hasJac ? numCells*totDim*totDim : 0,&elemMat,hasPrec ? numCells*totDim*totDim : 0, &elemMatP,hasDyn ? numCells*totDim*totDim : 0, &elemMatD);
2176: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2177: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2178: VecGetArray(cellgeom, &cgeomScal);
2179: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2180: DM dmCell;
2182: VecGetDM(cellgeom,&dmCell);
2183: PetscMalloc1(cEnd-cStart,&cgeom);
2184: for (c = 0; c < cEnd - cStart; c++) {
2185: PetscScalar *thisgeom;
2187: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2188: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2189: }
2190: } else {
2191: cgeom = (PetscFECellGeom *) cgeomScal;
2192: }
2193: for (c = cStart; c < cEnd; ++c) {
2194: PetscScalar *x = NULL, *x_t = NULL;
2195: PetscInt i;
2197: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2198: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2199: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2200: if (X_t) {
2201: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2202: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2203: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2204: }
2205: if (dmAux) {
2206: DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2207: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2208: DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2209: }
2210: }
2211: if (hasJac) {PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));}
2212: if (hasPrec) {PetscMemzero(elemMatP, numCells*totDim*totDim * sizeof(PetscScalar));}
2213: if (hasDyn) {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2214: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2215: PetscClassId id;
2216: PetscFE fe;
2217: PetscQuadrature q;
2218: PetscInt numQuadPoints, Nb;
2219: /* Conforming batches */
2220: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2221: /* Remainder */
2222: PetscInt Nr, offset;
2224: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2225: PetscObjectGetClassId((PetscObject) fe, &id);
2226: if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; continue;}
2227: PetscFEGetQuadrature(fe, &q);
2228: PetscFEGetDimension(fe, &Nb);
2229: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2230: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
2231: blockSize = Nb*numQuadPoints;
2232: batchSize = numBlocks * blockSize;
2233: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2234: numChunks = numCells / (numBatches*batchSize);
2235: Ne = numChunks*numBatches*batchSize;
2236: Nr = numCells % (numBatches*batchSize);
2237: offset = numCells - Nr;
2238: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2239: if (hasJac) {
2240: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2241: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
2242: }
2243: if (hasPrec) {
2244: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatP);
2245: 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], t, X_tShift, &elemMatP[offset*totDim*totDim]);
2246: }
2247: if (hasDyn) {
2248: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
2249: 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], t, X_tShift, &elemMatD[offset*totDim*totDim]);
2250: }
2251: }
2252: }
2253: if (hasDyn) {
2254: for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2255: }
2256: if (hasFV) {
2257: PetscClassId id;
2258: PetscFV fv;
2259: PetscInt offsetI, NcI, NbI = 1, fc, f;
2261: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2262: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
2263: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
2264: PetscObjectGetClassId((PetscObject) fv, &id);
2265: if (id != PETSCFV_CLASSID) continue;
2266: /* Put in the identity */
2267: PetscFVGetNumComponents(fv, &NcI);
2268: for (c = cStart; c < cEnd; ++c) {
2269: const PetscInt eOffset = (c-cStart)*totDim*totDim;
2270: for (fc = 0; fc < NcI; ++fc) {
2271: for (f = 0; f < NbI; ++f) {
2272: const PetscInt i = offsetI + f*NcI+fc;
2273: if (hasPrec) {
2274: if (hasJac) {elemMat[eOffset+i*totDim+i] = 1.0;}
2275: elemMatP[eOffset+i*totDim+i] = 1.0;
2276: } else {elemMat[eOffset+i*totDim+i] = 1.0;}
2277: }
2278: }
2279: }
2280: }
2281: /* No allocated space for FV stuff, so ignore the zero entries */
2282: MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
2283: }
2284: isMatIS = PETSC_FALSE;
2285: if (hasPrec && hasJac) {
2286: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatIS);
2287: }
2288: if (isMatIS && !subSection) {
2289: DMPlexGetSubdomainSection(dm, &subSection);
2290: }
2291: for (c = cStart; c < cEnd; ++c) {
2292: if (hasPrec) {
2293: if (hasJac) {
2294: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2295: if (!isMatIS) {
2296: DMPlexMatSetClosure(dm, section, globalSection, Jac, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2297: } else {
2298: Mat lJ;
2300: MatISGetLocalMat(Jac,&lJ);
2301: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2302: }
2303: }
2304: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
2305: if (!isMatISP) {
2306: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2307: } else {
2308: Mat lJ;
2310: MatISGetLocalMat(JacP,&lJ);
2311: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2312: }
2313: } else {
2314: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2315: if (!isMatISP) {
2316: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2317: } else {
2318: Mat lJ;
2320: MatISGetLocalMat(JacP,&lJ);
2321: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2322: }
2323: }
2324: }
2325: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
2326: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2327: PetscFree(cgeom);
2328: } else {
2329: cgeom = NULL;
2330: }
2331: VecRestoreArray(cellgeom, &cgeomScal);
2332: PetscFree5(u,u_t,elemMat,elemMatP,elemMatD);
2333: if (dmAux) {
2334: PetscFree(a);
2335: DMDestroy(&plex);
2336: }
2337: DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, user);
2338: if (hasJac && hasPrec) {
2339: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
2340: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
2341: }
2342: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2343: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2344: if (mesh->printFEM) {
2345: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2346: MatChop(JacP, 1.0e-10);
2347: MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2348: }
2349: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2350: PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2351: if (isShell) {
2352: JacActionCtx *jctx;
2354: MatShellGetContext(Jac, &jctx);
2355: VecCopy(X, jctx->u);
2356: }
2357: return(0);
2358: }
2361: PetscErrorCode DMPlexComputeJacobianAction_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal t, PetscReal X_tShift, Vec X, Vec X_t, Vec Y, Vec Z, void *user)
2362: {
2363: DM_Plex *mesh = (DM_Plex *) dm->data;
2364: const char *name = "Jacobian";
2365: DM dmAux, plex;
2366: Vec A, cellgeom;
2367: PetscDS prob, probAux = NULL;
2368: PetscQuadrature quad;
2369: PetscSection section, globalSection, sectionAux;
2370: PetscFECellGeom *cgeom = NULL;
2371: PetscScalar *cgeomScal;
2372: PetscScalar *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z;
2373: PetscInt dim, Nf, fieldI, fieldJ, numCells, c;
2374: PetscInt totDim, totDimAux = 0;
2375: PetscBool hasDyn;
2376: PetscErrorCode ierr;
2379: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2380: DMGetDimension(dm, &dim);
2381: DMGetDefaultSection(dm, §ion);
2382: DMGetDefaultGlobalSection(dm, &globalSection);
2383: DMGetDS(dm, &prob);
2384: PetscDSGetTotalDimension(prob, &totDim);
2385: PetscDSHasDynamicJacobian(prob, &hasDyn);
2386: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2387: PetscSectionGetNumFields(section, &Nf);
2388: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2389: numCells = cEnd - cStart;
2390: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2391: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2392: if (dmAux) {
2393: DMConvert(dmAux, DMPLEX, &plex);
2394: DMGetDefaultSection(plex, §ionAux);
2395: DMGetDS(dmAux, &probAux);
2396: PetscDSGetTotalDimension(probAux, &totDimAux);
2397: }
2398: VecSet(Z, 0.0);
2399: PetscMalloc6(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat,hasDyn ? numCells*totDim*totDim : 0, &elemMatD,numCells*totDim,&y,totDim,&z);
2400: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2401: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2402: VecGetArray(cellgeom, &cgeomScal);
2403: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2404: DM dmCell;
2406: VecGetDM(cellgeom,&dmCell);
2407: PetscMalloc1(cEnd-cStart,&cgeom);
2408: for (c = 0; c < cEnd - cStart; c++) {
2409: PetscScalar *thisgeom;
2411: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2412: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2413: }
2414: } else {
2415: cgeom = (PetscFECellGeom *) cgeomScal;
2416: }
2417: for (c = cStart; c < cEnd; ++c) {
2418: PetscScalar *x = NULL, *x_t = NULL;
2419: PetscInt i;
2421: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2422: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2423: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2424: if (X_t) {
2425: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2426: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2427: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2428: }
2429: if (dmAux) {
2430: DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2431: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2432: DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2433: }
2434: DMPlexVecGetClosure(dm, section, Y, c, NULL, &x);
2435: for (i = 0; i < totDim; ++i) y[(c-cStart)*totDim+i] = x[i];
2436: DMPlexVecRestoreClosure(dm, section, Y, c, NULL, &x);
2437: }
2438: PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2439: if (hasDyn) {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2440: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2441: PetscFE fe;
2442: PetscInt numQuadPoints, Nb;
2443: /* Conforming batches */
2444: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2445: /* Remainder */
2446: PetscInt Nr, offset;
2448: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2449: PetscFEGetQuadrature(fe, &quad);
2450: PetscFEGetDimension(fe, &Nb);
2451: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2452: PetscQuadratureGetData(quad, NULL, NULL, &numQuadPoints, NULL, NULL);
2453: blockSize = Nb*numQuadPoints;
2454: batchSize = numBlocks * blockSize;
2455: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2456: numChunks = numCells / (numBatches*batchSize);
2457: Ne = numChunks*numBatches*batchSize;
2458: Nr = numCells % (numBatches*batchSize);
2459: offset = numCells - Nr;
2460: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2461: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2462: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], t, X_tShift, &elemMat[offset*totDim*totDim]);
2463: if (hasDyn) {
2464: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
2465: 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], t, X_tShift, &elemMatD[offset*totDim*totDim]);
2466: }
2467: }
2468: }
2469: if (hasDyn) {
2470: for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2471: }
2472: for (c = cStart; c < cEnd; ++c) {
2473: const PetscBLASInt M = totDim, one = 1;
2474: const PetscScalar a = 1.0, b = 0.0;
2476: PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[(c-cStart)*totDim*totDim], &M, &y[(c-cStart)*totDim], &one, &b, z, &one));
2477: if (mesh->printFEM > 1) {
2478: DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);
2479: DMPrintCellVector(c, "Y", totDim, &y[(c-cStart)*totDim]);
2480: DMPrintCellVector(c, "Z", totDim, z);
2481: }
2482: DMPlexVecSetClosure(dm, section, Z, c, z, ADD_VALUES);
2483: }
2484: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {PetscFree(cgeom);}
2485: else {cgeom = NULL;}
2486: VecRestoreArray(cellgeom, &cgeomScal);
2487: PetscFree6(u,u_t,elemMat,elemMatD,y,z);
2488: if (dmAux) {
2489: PetscFree(a);
2490: DMDestroy(&plex);
2491: }
2492: if (mesh->printFEM) {
2493: PetscPrintf(PETSC_COMM_WORLD, "Z:\n");
2494: VecView(Z, PETSC_VIEWER_STDOUT_WORLD);
2495: }
2496: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2497: return(0);
2498: }
2500: /*@
2501: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
2503: Input Parameters:
2504: + dm - The mesh
2505: . X - Local input vector
2506: - user - The user context
2508: Output Parameter:
2509: . Jac - Jacobian matrix
2511: Note:
2512: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2513: like a GPU, or vectorize on a multicore machine.
2515: Level: developer
2517: .seealso: FormFunctionLocal()
2518: @*/
2519: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2520: {
2521: PetscInt cStart, cEnd, cEndInterior;
2522: DM plex;
2526: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2527: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2528: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2529: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2530: DMPlexComputeJacobian_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2531: DMDestroy(&plex);
2532: return(0);
2533: }
2535: /*@
2536: DMPlexSNESComputeJacobianActionFEM - Form the local portion of the Jacobian action Z = J(X) Y at the local solution X using pointwise functions specified by the user.
2538: Input Parameters:
2539: + dm - The mesh
2540: . X - Local solution vector
2541: . Y - Local input vector
2542: - user - The user context
2544: Output Parameter:
2545: . Z - Local output vector
2547: Note:
2548: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2549: like a GPU, or vectorize on a multicore machine.
2551: Level: developer
2553: .seealso: FormFunctionLocal()
2554: @*/
2555: PetscErrorCode DMPlexSNESComputeJacobianActionFEM(DM dm, Vec X, Vec Y, Vec Z, void *user)
2556: {
2557: PetscInt cStart, cEnd, cEndInterior;
2558: DM plex;
2562: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2563: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2564: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2565: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2566: DMPlexComputeJacobianAction_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Y, Z, user);
2567: DMDestroy(&plex);
2568: return(0);
2569: }
2571: /*@
2572: DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.
2574: Input Parameters:
2575: + dm - The DM object
2576: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
2577: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
2578: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())
2580: Level: developer
2581: @*/
2582: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
2583: {
2587: DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
2588: DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
2589: DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
2590: return(0);
2591: }
2593: 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)
2594: {
2595: PetscDS prob;
2596: Mat J, M;
2597: Vec r, b;
2598: MatNullSpace nullSpace;
2599: PetscReal *error, res = 0.0;
2600: PetscInt numFields;
2601: PetscBool hasJac, hasPrec;
2605: VecDuplicate(u, &r);
2606: DMCreateMatrix(dm, &J);
2607: /* TODO Null space for J */
2608: /* Check discretization error */
2609: DMGetNumFields(dm, &numFields);
2610: PetscMalloc1(PetscMax(1, numFields), &error);
2611: if (numFields > 1) {
2612: PetscInt f;
2614: DMComputeL2FieldDiff(dm, 0.0, exactFuncs, ctxs, u, error);
2615: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2616: for (f = 0; f < numFields; ++f) {
2617: if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2618: if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", (double)error[f]);}
2619: else {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2620: }
2621: PetscPrintf(PETSC_COMM_WORLD, "]\n");
2622: } else {
2623: DMComputeL2Diff(dm, 0.0, exactFuncs, ctxs, u, &error[0]);
2624: if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)error[0]);}
2625: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2626: }
2627: PetscFree(error);
2628: /* Check residual */
2629: SNESComputeFunction(snes, u, r);
2630: VecNorm(r, NORM_2, &res);
2631: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
2632: VecChop(r, 1.0e-10);
2633: PetscObjectSetName((PetscObject) r, "Initial Residual");
2634: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2635: VecViewFromOptions(r, NULL, "-vec_view");
2636: /* Check Jacobian */
2637: DMGetDS(dm, &prob);
2638: PetscDSHasJacobian(prob, &hasJac);
2639: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2640: if (hasJac && hasPrec) {
2641: DMCreateMatrix(dm, &M);
2642: SNESComputeJacobian(snes, u, J, M);
2643: PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
2644: MatViewFromOptions(M, NULL, "-mat_view");
2645: MatDestroy(&M);
2646: } else {
2647: SNESComputeJacobian(snes, u, J, J);
2648: }
2649: PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
2650: MatViewFromOptions(J, NULL, "-mat_view");
2651: MatGetNullSpace(J, &nullSpace);
2652: if (nullSpace) {
2653: PetscBool isNull;
2654: MatNullSpaceTest(nullSpace, J, &isNull);
2655: if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2656: }
2657: VecDuplicate(u, &b);
2658: VecSet(r, 0.0);
2659: SNESComputeFunction(snes, r, b);
2660: MatMult(J, u, r);
2661: VecAXPY(r, 1.0, b);
2662: VecDestroy(&b);
2663: VecNorm(r, NORM_2, &res);
2664: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
2665: VecChop(r, 1.0e-10);
2666: PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2667: PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2668: VecViewFromOptions(r, NULL, "-vec_view");
2669: VecDestroy(&r);
2670: MatNullSpaceDestroy(&nullSpace);
2671: MatDestroy(&J);
2672: return(0);
2673: }
2675: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2676: {
2677: DM dm;
2678: Vec sol;
2679: PetscBool check;
2683: PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
2684: if (!check) return(0);
2685: SNESGetDM(snes, &dm);
2686: VecDuplicate(u, &sol);
2687: SNESSetSolution(snes, sol);
2688: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2689: VecDestroy(&sol);
2690: return(0);
2691: }