Actual source code: dmplexsnes.c
petsc-3.8.4 2018-03-24
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, PETSC_SCALAR, u);
970: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u_t);} else {*u_t = NULL;}
971: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, PETSC_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, PETSC_SCALAR, u);
1023: if (*u_t) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u_t);}
1024: if (*a) {DMRestoreWorkArray(dm, 0, PETSC_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, PETSC_SCALAR, uL);
1097: DMGetWorkArray(dm, numFaces*Nc, PETSC_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, PETSC_SCALAR, uL);
1216: DMRestoreWorkArray(dm, 0, PETSC_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, PETSC_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, PETSC_REAL, vol);
1317: return(0);
1318: }
1320: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
1321: {
1322: DM_Plex *mesh = (DM_Plex *) dm->data;
1323: DM dmAux = NULL, plex = NULL;
1324: PetscSection section, sectionAux = NULL;
1325: PetscDS prob, probAux = NULL;
1326: DMLabel depth;
1327: Vec locA = NULL;
1328: PetscFEFaceGeom *fgeom;
1329: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemVec = NULL;
1330: PetscInt dim, totDim, totDimAux, numBd, bd;
1331: PetscErrorCode ierr;
1334: DMGetDimension(dm, &dim);
1335: DMGetDefaultSection(dm, §ion);
1336: DMPlexGetDepthLabel(dm, &depth);
1337: DMGetDS(dm, &prob);
1338: PetscDSGetTotalDimension(prob, &totDim);
1339: PetscDSGetNumBoundary(prob, &numBd);
1340: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1341: if (locA) {
1342: VecGetDM(locA, &dmAux);
1343: DMConvert(dmAux, DMPLEX, &plex);
1344: DMGetDefaultSection(plex, §ionAux);
1345: DMGetDS(dmAux, &probAux);
1346: PetscDSGetTotalDimension(probAux, &totDimAux);
1347: }
1348: for (bd = 0; bd < numBd; ++bd) {
1349: DMBoundaryConditionType type;
1350: const char *bdLabel;
1351: DMLabel label;
1352: IS pointIS;
1353: const PetscInt *points;
1354: const PetscInt *values;
1355: PetscInt field, numValues, v, numPoints, p, dep, numFaces, face;
1356: PetscObject obj;
1357: PetscClassId id;
1359: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1360: PetscDSGetDiscretization(prob, field, &obj);
1361: PetscObjectGetClassId(obj, &id);
1362: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
1363: DMGetLabel(dm, bdLabel, &label);
1364: for (v = 0; v < numValues; ++v) {
1365: DMLabelGetStratumSize(label, values[v], &numPoints);
1366: DMLabelGetStratumIS(label, values[v], &pointIS);
1367: if (!pointIS) continue; /* No points with that id on this process */
1368: ISGetIndices(pointIS, &points);
1369: for (p = 0, numFaces = 0; p < numPoints; ++p) {
1370: DMLabelGetValue(depth, points[p], &dep);
1371: if (dep == dim-1) ++numFaces;
1372: }
1373: PetscMalloc4(numFaces*totDim,&u,locX_t ? numFaces*totDim : 0, &u_t, numFaces,&fgeom, numFaces*totDim,&elemVec);
1374: if (locA) {PetscMalloc1(numFaces*totDimAux,&a);}
1375: for (p = 0, face = 0; p < numPoints; ++p) {
1376: const PetscInt point = points[p], *support, *cone;
1377: PetscScalar *x = NULL;
1378: PetscReal dummyJ[9], dummyDetJ;
1379: PetscInt i, coneSize, faceLoc;
1381: DMLabelGetValue(depth, points[p], &dep);
1382: if (dep != dim-1) continue;
1383: DMPlexGetSupport(dm, point, &support);
1384: DMPlexComputeCellGeometryFEM(dm, support[0], NULL, NULL, dummyJ, fgeom[face].invJ[0], &dummyDetJ);
1385: DMPlexComputeCellGeometryFEM(dm, point, NULL, fgeom[face].v0, fgeom[face].J, NULL, &fgeom[face].detJ);
1386: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, fgeom[face].n);
1387: 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);
1388: DMPlexGetConeSize(dm, support[0], &coneSize);
1389: DMPlexGetCone(dm, support[0], &cone);
1390: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
1391: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of support[0] %d", point, support[0]);
1392: fgeom[face].face[0] = faceLoc;
1393: DMPlexVecGetClosure(dm, section, locX, support[0], NULL, &x);
1394: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
1395: DMPlexVecRestoreClosure(dm, section, locX, support[0], NULL, &x);
1396: if (locX_t) {
1397: DMPlexVecGetClosure(dm, section, locX_t, support[0], NULL, &x);
1398: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
1399: DMPlexVecRestoreClosure(dm, section, locX_t, support[0], NULL, &x);
1400: }
1401: if (locA) {
1402: DMPlexVecGetClosure(plex, sectionAux, locA, support[0], NULL, &x);
1403: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
1404: DMPlexVecRestoreClosure(plex, sectionAux, locA, support[0], NULL, &x);
1405: }
1406: ++face;
1407: }
1408: PetscMemzero(elemVec, numFaces* totDim * sizeof(PetscScalar));
1409: {
1410: PetscFE fe;
1411: PetscQuadrature q;
1412: PetscInt numQuadPoints, Nb;
1413: /* Conforming batches */
1414: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1415: /* Remainder */
1416: PetscInt Nr, offset;
1418: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
1419: PetscFEGetFaceQuadrature(fe, &q);
1420: PetscFEGetDimension(fe, &Nb);
1421: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1422: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
1423: blockSize = Nb*numQuadPoints;
1424: batchSize = numBlocks * blockSize;
1425: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1426: numChunks = numFaces / (numBatches*batchSize);
1427: Ne = numChunks*numBatches*batchSize;
1428: Nr = numFaces % (numBatches*batchSize);
1429: offset = numFaces - Nr;
1430: PetscFEIntegrateBdResidual(fe, prob, field, Ne, fgeom, u, u_t, probAux, a, t, elemVec);
1431: 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]);
1432: }
1433: for (p = 0, face = 0; p < numPoints; ++p) {
1434: const PetscInt point = points[p], *support;
1436: DMLabelGetValue(depth, point, &dep);
1437: if (dep != dim-1) continue;
1438: if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDim, &elemVec[face*totDim]);}
1439: DMPlexGetSupport(dm, point, &support);
1440: DMPlexVecSetClosure(dm, NULL, locF, support[0], &elemVec[face*totDim], ADD_ALL_VALUES);
1441: ++face;
1442: }
1443: ISRestoreIndices(pointIS, &points);
1444: ISDestroy(&pointIS);
1445: PetscFree4(u,u_t,fgeom,elemVec);
1446: if (locA) {PetscFree(a);}
1447: }
1448: }
1449: if (plex) {DMDestroy(&plex);}
1450: return(0);
1451: }
1453: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *user)
1454: {
1455: DM_Plex *mesh = (DM_Plex *) dm->data;
1456: const char *name = "Residual";
1457: DM dmAux = NULL;
1458: DM dmGrad = NULL;
1459: DMLabel ghostLabel = NULL;
1460: PetscDS prob = NULL;
1461: PetscDS probAux = NULL;
1462: PetscSection section = NULL;
1463: PetscBool useFEM = PETSC_FALSE;
1464: PetscBool useFVM = PETSC_FALSE;
1465: PetscBool isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1466: PetscFV fvm = NULL;
1467: PetscFECellGeom *cgeomFEM = NULL;
1468: PetscScalar *cgeomScal;
1469: PetscFVCellGeom *cgeomFVM = NULL;
1470: PetscFVFaceGeom *fgeomFVM = NULL;
1471: Vec locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1472: PetscScalar *u = NULL, *u_t, *a, *uL, *uR;
1473: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1474: PetscErrorCode ierr;
1477: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1478: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1479: /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1480: /* FEM+FVM */
1481: /* 1: Get sizes from dm and dmAux */
1482: DMGetDefaultSection(dm, §ion);
1483: DMGetLabel(dm, "ghost", &ghostLabel);
1484: DMGetDS(dm, &prob);
1485: PetscDSGetNumFields(prob, &Nf);
1486: PetscDSGetTotalDimension(prob, &totDim);
1487: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1488: if (locA) {
1489: VecGetDM(locA, &dmAux);
1490: DMGetDS(dmAux, &probAux);
1491: PetscDSGetTotalDimension(probAux, &totDimAux);
1492: }
1493: /* 2: Get geometric data */
1494: for (f = 0; f < Nf; ++f) {
1495: PetscObject obj;
1496: PetscClassId id;
1497: PetscBool fimp;
1499: PetscDSGetImplicit(prob, f, &fimp);
1500: if (isImplicit != fimp) continue;
1501: PetscDSGetDiscretization(prob, f, &obj);
1502: PetscObjectGetClassId(obj, &id);
1503: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1504: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1505: }
1506: if (useFEM) {
1507: DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1508: VecGetArray(cellGeometryFEM, &cgeomScal);
1509: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1510: DM dmCell;
1511: PetscInt c;
1513: VecGetDM(cellGeometryFEM,&dmCell);
1514: PetscMalloc1(cEnd-cStart,&cgeomFEM);
1515: for (c = 0; c < cEnd - cStart; c++) {
1516: PetscScalar *thisgeom;
1518: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1519: cgeomFEM[c] = *((PetscFECellGeom *) thisgeom);
1520: }
1521: }
1522: else {
1523: cgeomFEM = (PetscFECellGeom *) cgeomScal;
1524: }
1525: }
1526: if (useFVM) {
1527: DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1528: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1529: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1530: /* Reconstruct and limit cell gradients */
1531: DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1532: if (dmGrad) {
1533: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1534: DMGetGlobalVector(dmGrad, &grad);
1535: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1536: /* Communicate gradient values */
1537: DMGetLocalVector(dmGrad, &locGrad);
1538: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1539: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1540: DMRestoreGlobalVector(dmGrad, &grad);
1541: }
1542: /* Handle non-essential (e.g. outflow) boundary values */
1543: DMPlexInsertBoundaryValues(dm, PETSC_FALSE, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1544: }
1545: /* Loop over chunks */
1546: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1547: numChunks = 1;
1548: cellChunkSize = (cEnd - cStart)/numChunks;
1549: faceChunkSize = (fEnd - fStart)/numChunks;
1550: for (chunk = 0; chunk < numChunks; ++chunk) {
1551: PetscScalar *elemVec, *fluxL, *fluxR;
1552: PetscReal *vol;
1553: PetscFVFaceGeom *fgeom;
1554: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1555: PetscInt fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = 0, face;
1557: /* Extract field coefficients */
1558: if (useFEM) {
1559: DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1560: DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1561: PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1562: }
1563: if (useFVM) {
1564: DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
1565: DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
1566: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1567: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1568: PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1569: PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1570: }
1571: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1572: /* Loop over fields */
1573: for (f = 0; f < Nf; ++f) {
1574: PetscObject obj;
1575: PetscClassId id;
1576: PetscBool fimp;
1577: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1579: PetscDSGetImplicit(prob, f, &fimp);
1580: if (isImplicit != fimp) continue;
1581: PetscDSGetDiscretization(prob, f, &obj);
1582: PetscObjectGetClassId(obj, &id);
1583: if (id == PETSCFE_CLASSID) {
1584: PetscFE fe = (PetscFE) obj;
1585: PetscQuadrature q;
1586: PetscInt Nq, Nb;
1588: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1590: PetscFEGetQuadrature(fe, &q);
1591: PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);
1592: PetscFEGetDimension(fe, &Nb);
1593: blockSize = Nb*Nq;
1594: batchSize = numBlocks * blockSize;
1595: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1596: numChunks = numCells / (numBatches*batchSize);
1597: Ne = numChunks*numBatches*batchSize;
1598: Nr = numCells % (numBatches*batchSize);
1599: offset = numCells - Nr;
1600: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1601: /* 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) */
1602: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, t, elemVec);
1603: 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]);
1604: } else if (id == PETSCFV_CLASSID) {
1605: PetscFV fv = (PetscFV) obj;
1607: Ne = numFaces;
1608: /* Riemann solve over faces (need fields at face centroids) */
1609: /* We need to evaluate FE fields at those coordinates */
1610: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1611: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1612: }
1613: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1614: PetscFree(cgeomFEM);
1615: }
1616: else {
1617: cgeomFEM = NULL;
1618: }
1619: if (cellGeometryFEM) {VecRestoreArray(cellGeometryFEM, &cgeomScal);}
1620: /* Loop over domain */
1621: if (useFEM) {
1622: /* Add elemVec to locX */
1623: for (cell = cS; cell < cE; ++cell) {
1624: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1625: if (ghostLabel) {
1626: PetscInt ghostVal;
1628: DMLabelGetValue(ghostLabel,cell,&ghostVal);
1629: if (ghostVal > 0) continue;
1630: }
1631: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_ALL_VALUES);
1632: }
1633: }
1634: if (useFVM) {
1635: PetscScalar *fa;
1636: PetscInt iface;
1638: VecGetArray(locF, &fa);
1639: for (f = 0; f < Nf; ++f) {
1640: PetscFV fv;
1641: PetscObject obj;
1642: PetscClassId id;
1643: PetscInt foff, pdim;
1645: PetscDSGetDiscretization(prob, f, &obj);
1646: PetscDSGetFieldOffset(prob, f, &foff);
1647: PetscObjectGetClassId(obj, &id);
1648: if (id != PETSCFV_CLASSID) continue;
1649: fv = (PetscFV) obj;
1650: PetscFVGetNumComponents(fv, &pdim);
1651: /* Accumulate fluxes to cells */
1652: for (face = fS, iface = 0; face < fE; ++face) {
1653: const PetscInt *cells;
1654: PetscScalar *fL = NULL, *fR = NULL;
1655: PetscInt ghost, d, nsupp, nchild;
1657: DMLabelGetValue(ghostLabel, face, &ghost);
1658: DMPlexGetSupportSize(dm, face, &nsupp);
1659: DMPlexGetTreeChildren(dm, face, &nchild, NULL);
1660: if (ghost >= 0 || nsupp > 2 || nchild > 0) continue;
1661: DMPlexGetSupport(dm, face, &cells);
1662: DMLabelGetValue(ghostLabel,cells[0],&ghost);
1663: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, cells[0], f, fa, &fL);}
1664: DMLabelGetValue(ghostLabel,cells[1],&ghost);
1665: if (ghost <= 0) {DMPlexPointLocalFieldRef(dm, cells[1], f, fa, &fR);}
1666: for (d = 0; d < pdim; ++d) {
1667: if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1668: if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1669: }
1670: ++iface;
1671: }
1672: }
1673: VecRestoreArray(locF, &fa);
1674: }
1675: /* Handle time derivative */
1676: if (locX_t) {
1677: PetscScalar *x_t, *fa;
1679: VecGetArray(locF, &fa);
1680: VecGetArray(locX_t, &x_t);
1681: for (f = 0; f < Nf; ++f) {
1682: PetscFV fv;
1683: PetscObject obj;
1684: PetscClassId id;
1685: PetscInt pdim, d;
1687: PetscDSGetDiscretization(prob, f, &obj);
1688: PetscObjectGetClassId(obj, &id);
1689: if (id != PETSCFV_CLASSID) continue;
1690: fv = (PetscFV) obj;
1691: PetscFVGetNumComponents(fv, &pdim);
1692: for (cell = cS; cell < cE; ++cell) {
1693: PetscScalar *u_t, *r;
1695: if (ghostLabel) {
1696: PetscInt ghostVal;
1698: DMLabelGetValue(ghostLabel,cell,&ghostVal);
1699: if (ghostVal > 0) continue;
1700: }
1701: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1702: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1703: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1704: }
1705: }
1706: VecRestoreArray(locX_t, &x_t);
1707: VecRestoreArray(locF, &fa);
1708: }
1709: if (useFEM) {
1710: DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1711: DMRestoreWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1712: }
1713: if (useFVM) {
1714: DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &numFaces, &uL, &uR);
1715: DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &numFaces, &fgeom, &vol);
1716: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1717: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1718: if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1719: }
1720: }
1722: if (useFEM) {DMPlexComputeBdResidual_Internal(dm, locX, locX_t, t, locF, user);}
1724: /* FEM */
1725: /* 1: Get sizes from dm and dmAux */
1726: /* 2: Get geometric data */
1727: /* 3: Handle boundary values */
1728: /* 4: Loop over domain */
1729: /* Extract coefficients */
1730: /* Loop over fields */
1731: /* Set tiling for FE*/
1732: /* Integrate FE residual to get elemVec */
1733: /* Loop over subdomain */
1734: /* Loop over quad points */
1735: /* Transform coords to real space */
1736: /* Evaluate field and aux fields at point */
1737: /* Evaluate residual at point */
1738: /* Transform residual to real space */
1739: /* Add residual to elemVec */
1740: /* Loop over domain */
1741: /* Add elemVec to locX */
1743: /* FVM */
1744: /* Get geometric data */
1745: /* If using gradients */
1746: /* Compute gradient data */
1747: /* Loop over domain faces */
1748: /* Count computational faces */
1749: /* Reconstruct cell gradient */
1750: /* Loop over domain cells */
1751: /* Limit cell gradients */
1752: /* Handle boundary values */
1753: /* Loop over domain faces */
1754: /* Read out field, centroid, normal, volume for each side of face */
1755: /* Riemann solve over faces */
1756: /* Loop over domain faces */
1757: /* Accumulate fluxes to cells */
1758: /* TODO Change printFEM to printDisc here */
1759: if (mesh->printFEM) {
1760: Vec locFbc;
1761: PetscInt pStart, pEnd, p, maxDof;
1762: PetscScalar *zeroes;
1764: VecDuplicate(locF,&locFbc);
1765: VecCopy(locF,locFbc);
1766: PetscSectionGetChart(section,&pStart,&pEnd);
1767: PetscSectionGetMaxDof(section,&maxDof);
1768: PetscCalloc1(maxDof,&zeroes);
1769: for (p = pStart; p < pEnd; p++) {
1770: VecSetValuesSection(locFbc,section,p,zeroes,INSERT_BC_VALUES);
1771: }
1772: PetscFree(zeroes);
1773: DMPrintLocalVec(dm, name, mesh->printTol, locFbc);
1774: VecDestroy(&locFbc);
1775: }
1776: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1777: return(0);
1778: }
1780: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, PetscReal t, Vec F, void *user)
1781: {
1782: DM dmCh, dmAux;
1783: Vec A, cellgeom;
1784: PetscDS prob, probCh, probAux = NULL;
1785: PetscQuadrature q;
1786: PetscSection section, sectionAux;
1787: PetscFECellGeom *cgeom = NULL;
1788: PetscScalar *cgeomScal;
1789: PetscScalar *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1790: PetscInt dim, Nf, f, numCells, cStart, cEnd, c;
1791: PetscInt totDim, totDimAux = 0, diffCell = 0;
1792: PetscErrorCode ierr;
1795: DMGetDimension(dm, &dim);
1796: DMGetDefaultSection(dm, §ion);
1797: DMGetDS(dm, &prob);
1798: PetscDSGetTotalDimension(prob, &totDim);
1799: PetscSectionGetNumFields(section, &Nf);
1800: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1801: numCells = cEnd - cStart;
1802: PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
1803: DMGetDS(dmCh, &probCh);
1804: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1805: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1806: if (dmAux) {
1807: DMGetDefaultSection(dmAux, §ionAux);
1808: DMGetDS(dmAux, &probAux);
1809: PetscDSGetTotalDimension(probAux, &totDimAux);
1810: }
1811: VecSet(F, 0.0);
1812: PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
1813: PetscMalloc1(numCells*totDim,&elemVecCh);
1814: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1815: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1816: VecGetArray(cellgeom, &cgeomScal);
1817: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1818: DM dmCell;
1820: VecGetDM(cellgeom,&dmCell);
1821: PetscMalloc1(cEnd-cStart,&cgeom);
1822: for (c = 0; c < cEnd - cStart; c++) {
1823: PetscScalar *thisgeom;
1825: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
1826: cgeom[c] = *((PetscFECellGeom *) thisgeom);
1827: }
1828: }
1829: else {
1830: cgeom = (PetscFECellGeom *) cgeomScal;
1831: }
1832: for (c = cStart; c < cEnd; ++c) {
1833: PetscScalar *x = NULL, *x_t = NULL;
1834: PetscInt i;
1836: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1837: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1838: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1839: if (X_t) {
1840: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1841: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1842: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1843: }
1844: if (dmAux) {
1845: DM dmAuxPlex;
1847: DMSNESConvertPlex(dmAux,&dmAuxPlex, PETSC_FALSE);
1848: DMPlexVecGetClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
1849: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1850: DMPlexVecRestoreClosure(dmAuxPlex, sectionAux, A, c, NULL, &x);
1851: DMDestroy(&dmAuxPlex);
1852: }
1853: }
1854: for (f = 0; f < Nf; ++f) {
1855: PetscFE fe, feCh;
1856: PetscInt numQuadPoints, Nb;
1857: /* Conforming batches */
1858: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1859: /* Remainder */
1860: PetscInt Nr, offset;
1862: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1863: PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
1864: PetscFEGetQuadrature(fe, &q);
1865: PetscFEGetDimension(fe, &Nb);
1866: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1867: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
1868: blockSize = Nb*numQuadPoints;
1869: batchSize = numBlocks * blockSize;
1870: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1871: numChunks = numCells / (numBatches*batchSize);
1872: Ne = numChunks*numBatches*batchSize;
1873: Nr = numCells % (numBatches*batchSize);
1874: offset = numCells - Nr;
1875: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, t, elemVec);
1876: PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, t, elemVecCh);
1877: 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]);
1878: 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]);
1879: }
1880: for (c = cStart; c < cEnd; ++c) {
1881: PetscBool diff = PETSC_FALSE;
1882: PetscInt d;
1884: for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1885: if (diff) {
1886: PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
1887: DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
1888: DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
1889: ++diffCell;
1890: }
1891: if (diffCell > 9) break;
1892: DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_ALL_VALUES);
1893: }
1894: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
1895: PetscFree(cgeom);
1896: }
1897: else {
1898: cgeom = NULL;
1899: }
1900: VecRestoreArray(cellgeom, &cgeomScal);
1901: PetscFree3(u,u_t,elemVec);
1902: PetscFree(elemVecCh);
1903: if (dmAux) {PetscFree(a);}
1904: return(0);
1905: }
1907: /*@
1908: DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1910: Input Parameters:
1911: + dm - The mesh
1912: . X - Local solution
1913: - user - The user context
1915: Output Parameter:
1916: . F - Local output vector
1918: Level: developer
1920: .seealso: DMPlexComputeJacobianActionFEM()
1921: @*/
1922: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1923: {
1924: PetscObject check;
1925: PetscInt cStart, cEnd, cEndInterior;
1926: DM plex;
1930: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1931: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
1932: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
1933: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1934: /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1935: PetscObjectQuery((PetscObject) plex, "dmCh", &check);
1936: if (check) {DMPlexComputeResidualFEM_Check_Internal(plex, X, NULL, 0.0, F, user);}
1937: else {DMPlexComputeResidual_Internal(plex, cStart, cEnd, PETSC_MIN_REAL, X, NULL, 0.0, F, user);}
1938: DMDestroy(&plex);
1939: return(0);
1940: }
1942: /*@
1943: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input X
1945: Input Parameters:
1946: + dm - The mesh
1947: - user - The user context
1949: Output Parameter:
1950: . X - Local solution
1952: Level: developer
1954: .seealso: DMPlexComputeJacobianActionFEM()
1955: @*/
1956: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
1957: {
1958: DM plex;
1962: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
1963: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL);
1964: DMDestroy(&plex);
1965: return(0);
1966: }
1968: PetscErrorCode DMPlexComputeBdJacobian_Internal(DM dm, Vec locX, Vec locX_t, PetscReal t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
1969: {
1970: DM_Plex *mesh = (DM_Plex *) dm->data;
1971: DM dmAux = NULL, plex = NULL;
1972: PetscSection section, globalSection, subSection, sectionAux = NULL;
1973: PetscDS prob, probAux = NULL;
1974: DMLabel depth;
1975: Vec locA = NULL;
1976: PetscFEFaceGeom *fgeom;
1977: PetscScalar *u = NULL, *u_t = NULL, *a = NULL, *elemMat = NULL;
1978: PetscInt dim, totDim, totDimAux, numBd, bd, Nf;
1979: PetscBool isMatISP;
1980: PetscErrorCode ierr;
1983: DMGetDimension(dm, &dim);
1984: DMGetDefaultSection(dm, §ion);
1985: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
1986: DMGetDefaultGlobalSection(dm, &globalSection);
1987: if (isMatISP) {
1988: DMPlexGetSubdomainSection(dm, &subSection);
1989: }
1990: DMPlexGetDepthLabel(dm, &depth);
1991: DMGetDS(dm, &prob);
1992: PetscDSGetNumFields(prob, &Nf);
1993: PetscDSGetTotalDimension(prob, &totDim);
1994: PetscDSGetNumBoundary(prob, &numBd);
1995: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1996: if (locA) {
1997: VecGetDM(locA, &dmAux);
1998: DMConvert(dmAux, DMPLEX, &plex);
1999: DMGetDefaultSection(plex, §ionAux);
2000: DMGetDS(dmAux, &probAux);
2001: PetscDSGetTotalDimension(probAux, &totDimAux);
2002: }
2003: for (bd = 0; bd < numBd; ++bd) {
2004: DMBoundaryConditionType type;
2005: const char *bdLabel;
2006: DMLabel label;
2007: IS pointIS;
2008: const PetscInt *points;
2009: const PetscInt *values;
2010: PetscInt fieldI, fieldJ, numValues, v, numPoints, p, dep, numFaces, face;
2011: PetscObject obj;
2012: PetscClassId id;
2014: PetscDSGetBoundary(prob, bd, &type, NULL, &bdLabel, &fieldI, NULL, NULL, NULL, &numValues, &values, NULL);
2015: PetscDSGetDiscretization(prob, fieldI, &obj);
2016: PetscObjectGetClassId(obj, &id);
2017: if ((id != PETSCFE_CLASSID) || (type & DM_BC_ESSENTIAL)) continue;
2018: DMGetLabel(dm, bdLabel, &label);
2019: for (v = 0; v < numValues; ++v) {
2020: DMLabelGetStratumSize(label, values[v], &numPoints);
2021: DMLabelGetStratumIS(label, values[v], &pointIS);
2022: if (!pointIS) continue; /* No points with that id on this process */
2023: ISGetIndices(pointIS, &points);
2024: for (p = 0, numFaces = 0; p < numPoints; ++p) {
2025: DMLabelGetValue(depth, points[p], &dep);
2026: if (dep == dim-1) ++numFaces;
2027: }
2028: PetscMalloc4(numFaces*totDim,&u,locX_t ? numFaces*totDim : 0,&u_t,numFaces,&fgeom,numFaces*totDim*totDim,&elemMat);
2029: if (locA) {PetscMalloc1(numFaces*totDimAux,&a);}
2030: for (p = 0, face = 0; p < numPoints; ++p) {
2031: const PetscInt point = points[p], *support, *cone;
2032: PetscScalar *x = NULL;
2033: PetscReal dummyJ[9], dummyDetJ;
2034: PetscInt i, coneSize, faceLoc;
2036: DMLabelGetValue(depth, points[p], &dep);
2037: if (dep != dim-1) continue;
2038: DMPlexGetSupport(dm, point, &support);
2039: DMPlexComputeCellGeometryFEM(dm, support[0], NULL, NULL, dummyJ, fgeom[face].invJ[0], &dummyDetJ);
2040: DMPlexComputeCellGeometryFEM(dm, point, NULL, fgeom[face].v0, fgeom[face].J, NULL, &fgeom[face].detJ);
2041: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, fgeom[face].n);
2042: 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);
2043: DMPlexGetConeSize(dm, support[0], &coneSize);
2044: DMPlexGetCone(dm, support[0], &cone);
2045: for (faceLoc = 0; faceLoc < coneSize; ++faceLoc) if (cone[faceLoc] == point) break;
2046: if (faceLoc == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of support[0] %d", point, support[0]);
2047: fgeom[face].face[0] = faceLoc;
2048: DMPlexVecGetClosure(dm, section, locX, support[0], NULL, &x);
2049: for (i = 0; i < totDim; ++i) u[face*totDim+i] = x[i];
2050: DMPlexVecRestoreClosure(dm, section, locX, support[0], NULL, &x);
2051: if (locX_t) {
2052: DMPlexVecGetClosure(dm, section, locX_t, support[0], NULL, &x);
2053: for (i = 0; i < totDim; ++i) u_t[face*totDim+i] = x[i];
2054: DMPlexVecRestoreClosure(dm, section, locX_t, support[0], NULL, &x);
2055: }
2056: if (locA) {
2057: DMPlexVecGetClosure(plex, sectionAux, locA, support[0], NULL, &x);
2058: for (i = 0; i < totDimAux; ++i) a[face*totDimAux+i] = x[i];
2059: DMPlexVecRestoreClosure(plex, sectionAux, locA, support[0], NULL, &x);
2060: }
2061: ++face;
2062: }
2063: PetscMemzero(elemMat, numFaces*totDim*totDim * sizeof(PetscScalar));
2064: {
2065: PetscFE fe;
2066: PetscQuadrature q;
2067: PetscInt numQuadPoints, Nb;
2068: /* Conforming batches */
2069: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2070: /* Remainder */
2071: PetscInt Nr, offset;
2073: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2074: PetscFEGetFaceQuadrature(fe, &q);
2075: PetscFEGetDimension(fe, &Nb);
2076: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2077: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
2078: blockSize = Nb*numQuadPoints;
2079: batchSize = numBlocks * blockSize;
2080: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2081: numChunks = numFaces / (numBatches*batchSize);
2082: Ne = numChunks*numBatches*batchSize;
2083: Nr = numFaces % (numBatches*batchSize);
2084: offset = numFaces - Nr;
2085: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2086: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, fgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2087: 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]);
2088: }
2089: }
2090: for (p = 0, face = 0; p < numPoints; ++p) {
2091: const PetscInt point = points[p], *support;
2093: DMLabelGetValue(depth, point, &dep);
2094: if (dep != dim-1) continue;
2095: if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDim, totDim, &elemMat[face*totDim*totDim]);}
2096: DMPlexGetSupport(dm, point, &support);
2097: if (!isMatISP) {
2098: DMPlexMatSetClosure(dm, section, globalSection, JacP, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
2099: } else {
2100: Mat lJ;
2102: MatISGetLocalMat(JacP,&lJ);
2103: DMPlexMatSetClosure(dm, section, subSection, lJ, support[0], &elemMat[face*totDim*totDim], ADD_VALUES);
2104: }
2105: ++face;
2106: }
2107: ISRestoreIndices(pointIS, &points);
2108: ISDestroy(&pointIS);
2109: PetscFree4(u,u_t,fgeom,elemMat);
2110: if (locA) {PetscFree(a);}
2111: }
2112: }
2113: if (plex) {DMDestroy(&plex);}
2114: return(0);
2115: }
2117: 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)
2118: {
2119: DM_Plex *mesh = (DM_Plex *) dm->data;
2120: const char *name = "Jacobian";
2121: DM dmAux, plex;
2122: Vec A, cellgeom;
2123: PetscDS prob, probAux = NULL;
2124: PetscSection section, globalSection, subSection, sectionAux;
2125: PetscFECellGeom *cgeom = NULL;
2126: PetscScalar *cgeomScal;
2127: PetscScalar *elemMat, *elemMatP, *elemMatD, *u, *u_t, *a = NULL;
2128: PetscInt dim, Nf, fieldI, fieldJ, numCells, c;
2129: PetscInt totDim, totDimAux;
2130: PetscBool isMatIS, isMatISP, isShell, hasJac, hasPrec, hasDyn, hasFV = PETSC_FALSE;
2131: PetscErrorCode ierr;
2134: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2135: DMGetDimension(dm, &dim);
2136: DMGetDefaultSection(dm, §ion);
2137: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatISP);
2138: DMGetDefaultGlobalSection(dm, &globalSection);
2139: if (isMatISP) {
2140: DMPlexGetSubdomainSection(dm, &subSection);
2141: }
2142: DMGetDS(dm, &prob);
2143: PetscDSGetTotalDimension(prob, &totDim);
2144: PetscDSHasJacobian(prob, &hasJac);
2145: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2146: PetscDSHasDynamicJacobian(prob, &hasDyn);
2147: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2148: PetscSectionGetNumFields(section, &Nf);
2149: numCells = cEnd - cStart;
2150: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2151: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2152: if (dmAux) {
2153: DMConvert(dmAux, DMPLEX, &plex);
2154: DMGetDefaultSection(plex, §ionAux);
2155: DMGetDS(dmAux, &probAux);
2156: PetscDSGetTotalDimension(probAux, &totDimAux);
2157: }
2158: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
2159: MatZeroEntries(JacP);
2160: 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);
2161: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2162: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2163: VecGetArray(cellgeom, &cgeomScal);
2164: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2165: DM dmCell;
2167: VecGetDM(cellgeom,&dmCell);
2168: PetscMalloc1(cEnd-cStart,&cgeom);
2169: for (c = 0; c < cEnd - cStart; c++) {
2170: PetscScalar *thisgeom;
2172: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2173: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2174: }
2175: } else {
2176: cgeom = (PetscFECellGeom *) cgeomScal;
2177: }
2178: for (c = cStart; c < cEnd; ++c) {
2179: PetscScalar *x = NULL, *x_t = NULL;
2180: PetscInt i;
2182: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2183: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2184: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2185: if (X_t) {
2186: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2187: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2188: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2189: }
2190: if (dmAux) {
2191: DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2192: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2193: DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2194: }
2195: }
2196: if (hasJac) {PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));}
2197: if (hasPrec) {PetscMemzero(elemMatP, numCells*totDim*totDim * sizeof(PetscScalar));}
2198: if (hasDyn) {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2199: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2200: PetscClassId id;
2201: PetscFE fe;
2202: PetscQuadrature q;
2203: PetscInt numQuadPoints, Nb;
2204: /* Conforming batches */
2205: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2206: /* Remainder */
2207: PetscInt Nr, offset;
2209: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2210: PetscObjectGetClassId((PetscObject) fe, &id);
2211: if (id == PETSCFV_CLASSID) {hasFV = PETSC_TRUE; continue;}
2212: PetscFEGetQuadrature(fe, &q);
2213: PetscFEGetDimension(fe, &Nb);
2214: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2215: PetscQuadratureGetData(q, NULL, NULL, &numQuadPoints, NULL, NULL);
2216: blockSize = Nb*numQuadPoints;
2217: batchSize = numBlocks * blockSize;
2218: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2219: numChunks = numCells / (numBatches*batchSize);
2220: Ne = numChunks*numBatches*batchSize;
2221: Nr = numCells % (numBatches*batchSize);
2222: offset = numCells - Nr;
2223: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2224: if (hasJac) {
2225: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2226: 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]);
2227: }
2228: if (hasPrec) {
2229: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_PRE, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatP);
2230: 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]);
2231: }
2232: if (hasDyn) {
2233: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
2234: 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]);
2235: }
2236: }
2237: }
2238: if (hasDyn) {
2239: for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2240: }
2241: if (hasFV) {
2242: PetscClassId id;
2243: PetscFV fv;
2244: PetscInt offsetI, NcI, NbI = 1, fc, f;
2246: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2247: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fv);
2248: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
2249: PetscObjectGetClassId((PetscObject) fv, &id);
2250: if (id != PETSCFV_CLASSID) continue;
2251: /* Put in the identity */
2252: PetscFVGetNumComponents(fv, &NcI);
2253: for (c = cStart; c < cEnd; ++c) {
2254: const PetscInt eOffset = (c-cStart)*totDim*totDim;
2255: for (fc = 0; fc < NcI; ++fc) {
2256: for (f = 0; f < NbI; ++f) {
2257: const PetscInt i = offsetI + f*NcI+fc;
2258: if (hasPrec) {
2259: if (hasJac) {elemMat[eOffset+i*totDim+i] = 1.0;}
2260: elemMatP[eOffset+i*totDim+i] = 1.0;
2261: } else {elemMat[eOffset+i*totDim+i] = 1.0;}
2262: }
2263: }
2264: }
2265: }
2266: /* No allocated space for FV stuff, so ignore the zero entries */
2267: MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
2268: }
2269: isMatIS = PETSC_FALSE;
2270: if (hasPrec && hasJac) {
2271: PetscObjectTypeCompare((PetscObject) JacP, MATIS, &isMatIS);
2272: }
2273: if (isMatIS && !subSection) {
2274: DMPlexGetSubdomainSection(dm, &subSection);
2275: }
2276: for (c = cStart; c < cEnd; ++c) {
2277: if (hasPrec) {
2278: if (hasJac) {
2279: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2280: if (!isMatIS) {
2281: DMPlexMatSetClosure(dm, section, globalSection, Jac, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2282: } else {
2283: Mat lJ;
2285: MatISGetLocalMat(Jac,&lJ);
2286: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2287: }
2288: }
2289: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMatP[(c-cStart)*totDim*totDim]);}
2290: if (!isMatISP) {
2291: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2292: } else {
2293: Mat lJ;
2295: MatISGetLocalMat(JacP,&lJ);
2296: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMatP[(c-cStart)*totDim*totDim], ADD_VALUES);
2297: }
2298: } else {
2299: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2300: if (!isMatISP) {
2301: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2302: } else {
2303: Mat lJ;
2305: MatISGetLocalMat(JacP,&lJ);
2306: DMPlexMatSetClosure(dm, section, subSection, lJ, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2307: }
2308: }
2309: }
2310: if (hasFV) {MatSetOption(JacP, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);}
2311: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2312: PetscFree(cgeom);
2313: } else {
2314: cgeom = NULL;
2315: }
2316: VecRestoreArray(cellgeom, &cgeomScal);
2317: PetscFree5(u,u_t,elemMat,elemMatP,elemMatD);
2318: if (dmAux) {
2319: PetscFree(a);
2320: DMDestroy(&plex);
2321: }
2322: DMPlexComputeBdJacobian_Internal(dm, X, X_t, t, X_tShift, Jac, JacP, user);
2323: if (hasJac && hasPrec) {
2324: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
2325: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
2326: }
2327: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2328: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2329: if (mesh->printFEM) {
2330: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2331: MatChop(JacP, 1.0e-10);
2332: MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2333: }
2334: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2335: PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2336: if (isShell) {
2337: JacActionCtx *jctx;
2339: MatShellGetContext(Jac, &jctx);
2340: VecCopy(X, jctx->u);
2341: }
2342: return(0);
2343: }
2346: 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)
2347: {
2348: DM_Plex *mesh = (DM_Plex *) dm->data;
2349: const char *name = "Jacobian";
2350: DM dmAux, plex;
2351: Vec A, cellgeom;
2352: PetscDS prob, probAux = NULL;
2353: PetscQuadrature quad;
2354: PetscSection section, globalSection, sectionAux;
2355: PetscFECellGeom *cgeom = NULL;
2356: PetscScalar *cgeomScal;
2357: PetscScalar *elemMat, *elemMatD, *u, *u_t, *a = NULL, *y, *z;
2358: PetscInt dim, Nf, fieldI, fieldJ, numCells, c;
2359: PetscInt totDim, totDimAux = 0;
2360: PetscBool hasDyn;
2361: PetscErrorCode ierr;
2364: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
2365: DMGetDimension(dm, &dim);
2366: DMGetDefaultSection(dm, §ion);
2367: DMGetDefaultGlobalSection(dm, &globalSection);
2368: DMGetDS(dm, &prob);
2369: PetscDSGetTotalDimension(prob, &totDim);
2370: PetscDSHasDynamicJacobian(prob, &hasDyn);
2371: hasDyn = hasDyn && (X_tShift != 0.0) ? PETSC_TRUE : PETSC_FALSE;
2372: PetscSectionGetNumFields(section, &Nf);
2373: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2374: numCells = cEnd - cStart;
2375: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
2376: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
2377: if (dmAux) {
2378: DMConvert(dmAux, DMPLEX, &plex);
2379: DMGetDefaultSection(plex, §ionAux);
2380: DMGetDS(dmAux, &probAux);
2381: PetscDSGetTotalDimension(probAux, &totDimAux);
2382: }
2383: VecSet(Z, 0.0);
2384: 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);
2385: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
2386: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
2387: VecGetArray(cellgeom, &cgeomScal);
2388: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {
2389: DM dmCell;
2391: VecGetDM(cellgeom,&dmCell);
2392: PetscMalloc1(cEnd-cStart,&cgeom);
2393: for (c = 0; c < cEnd - cStart; c++) {
2394: PetscScalar *thisgeom;
2396: DMPlexPointLocalRef(dmCell, c + cStart, cgeomScal, &thisgeom);
2397: cgeom[c] = *((PetscFECellGeom *) thisgeom);
2398: }
2399: } else {
2400: cgeom = (PetscFECellGeom *) cgeomScal;
2401: }
2402: for (c = cStart; c < cEnd; ++c) {
2403: PetscScalar *x = NULL, *x_t = NULL;
2404: PetscInt i;
2406: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2407: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2408: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2409: if (X_t) {
2410: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2411: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2412: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2413: }
2414: if (dmAux) {
2415: DMPlexVecGetClosure(plex, sectionAux, A, c, NULL, &x);
2416: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2417: DMPlexVecRestoreClosure(plex, sectionAux, A, c, NULL, &x);
2418: }
2419: DMPlexVecGetClosure(dm, section, Y, c, NULL, &x);
2420: for (i = 0; i < totDim; ++i) y[(c-cStart)*totDim+i] = x[i];
2421: DMPlexVecRestoreClosure(dm, section, Y, c, NULL, &x);
2422: }
2423: PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2424: if (hasDyn) {PetscMemzero(elemMatD, numCells*totDim*totDim * sizeof(PetscScalar));}
2425: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2426: PetscFE fe;
2427: PetscInt numQuadPoints, Nb;
2428: /* Conforming batches */
2429: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2430: /* Remainder */
2431: PetscInt Nr, offset;
2433: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2434: PetscFEGetQuadrature(fe, &quad);
2435: PetscFEGetDimension(fe, &Nb);
2436: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2437: PetscQuadratureGetData(quad, NULL, NULL, &numQuadPoints, NULL, NULL);
2438: blockSize = Nb*numQuadPoints;
2439: batchSize = numBlocks * blockSize;
2440: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2441: numChunks = numCells / (numBatches*batchSize);
2442: Ne = numChunks*numBatches*batchSize;
2443: Nr = numCells % (numBatches*batchSize);
2444: offset = numCells - Nr;
2445: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2446: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMat);
2447: 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]);
2448: if (hasDyn) {
2449: PetscFEIntegrateJacobian(fe, prob, PETSCFE_JACOBIAN_DYN, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, t, X_tShift, elemMatD);
2450: 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]);
2451: }
2452: }
2453: }
2454: if (hasDyn) {
2455: for (c = 0; c < (cEnd - cStart)*totDim*totDim; ++c) elemMat[c] += X_tShift*elemMatD[c];
2456: }
2457: for (c = cStart; c < cEnd; ++c) {
2458: const PetscBLASInt M = totDim, one = 1;
2459: const PetscScalar a = 1.0, b = 0.0;
2461: PetscStackCallBLAS("BLASgemv", BLASgemv_("N", &M, &M, &a, &elemMat[(c-cStart)*totDim*totDim], &M, &y[(c-cStart)*totDim], &one, &b, z, &one));
2462: if (mesh->printFEM > 1) {
2463: DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);
2464: DMPrintCellVector(c, "Y", totDim, &y[(c-cStart)*totDim]);
2465: DMPrintCellVector(c, "Z", totDim, z);
2466: }
2467: DMPlexVecSetClosure(dm, section, Z, c, z, ADD_VALUES);
2468: }
2469: if (sizeof(PetscFECellGeom) % sizeof(PetscScalar)) {PetscFree(cgeom);}
2470: else {cgeom = NULL;}
2471: VecRestoreArray(cellgeom, &cgeomScal);
2472: PetscFree6(u,u_t,elemMat,elemMatD,y,z);
2473: if (dmAux) {
2474: PetscFree(a);
2475: DMDestroy(&plex);
2476: }
2477: if (mesh->printFEM) {
2478: PetscPrintf(PETSC_COMM_WORLD, "Z:\n");
2479: VecView(Z, PETSC_VIEWER_STDOUT_WORLD);
2480: }
2481: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2482: return(0);
2483: }
2485: /*@
2486: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
2488: Input Parameters:
2489: + dm - The mesh
2490: . X - Local input vector
2491: - user - The user context
2493: Output Parameter:
2494: . Jac - Jacobian matrix
2496: Note:
2497: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2498: like a GPU, or vectorize on a multicore machine.
2500: Level: developer
2502: .seealso: FormFunctionLocal()
2503: @*/
2504: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2505: {
2506: PetscInt cStart, cEnd, cEndInterior;
2507: DM plex;
2511: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2512: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2513: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2514: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2515: DMPlexComputeJacobian_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2516: DMDestroy(&plex);
2517: return(0);
2518: }
2520: /*@
2521: 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.
2523: Input Parameters:
2524: + dm - The mesh
2525: . X - Local solution vector
2526: . Y - Local input vector
2527: - user - The user context
2529: Output Parameter:
2530: . Z - Local output vector
2532: Note:
2533: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2534: like a GPU, or vectorize on a multicore machine.
2536: Level: developer
2538: .seealso: FormFunctionLocal()
2539: @*/
2540: PetscErrorCode DMPlexSNESComputeJacobianActionFEM(DM dm, Vec X, Vec Y, Vec Z, void *user)
2541: {
2542: PetscInt cStart, cEnd, cEndInterior;
2543: DM plex;
2547: DMSNESConvertPlex(dm,&plex,PETSC_TRUE);
2548: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
2549: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
2550: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2551: DMPlexComputeJacobianAction_Internal(plex, cStart, cEnd, 0.0, 0.0, X, NULL, Y, Z, user);
2552: DMDestroy(&plex);
2553: return(0);
2554: }
2556: /*@
2557: DMPlexSetSNESLocalFEM - Use DMPlex's internal FEM routines to compute SNES boundary values, residual, and Jacobian.
2559: Input Parameters:
2560: + dm - The DM object
2561: . boundaryctx - the user context that will be passed to pointwise evaluation of boundary values (see PetscDSAddBoundary())
2562: . residualctx - the user context that will be passed to pointwise evaluation of finite element residual computations (see PetscDSSetResidual())
2563: - jacobianctx - the user context that will be passed to pointwise evaluation of finite element Jacobian construction (see PetscDSSetJacobian())
2565: Level: developer
2566: @*/
2567: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, void *boundaryctx, void *residualctx, void *jacobianctx)
2568: {
2572: DMSNESSetBoundaryLocal(dm,DMPlexSNESComputeBoundaryFEM,boundaryctx);
2573: DMSNESSetFunctionLocal(dm,DMPlexSNESComputeResidualFEM,residualctx);
2574: DMSNESSetJacobianLocal(dm,DMPlexSNESComputeJacobianFEM,jacobianctx);
2575: return(0);
2576: }
2578: 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)
2579: {
2580: PetscDS prob;
2581: Mat J, M;
2582: Vec r, b;
2583: MatNullSpace nullSpace;
2584: PetscReal *error, res = 0.0;
2585: PetscInt numFields;
2586: PetscBool hasJac, hasPrec;
2590: VecDuplicate(u, &r);
2591: DMCreateMatrix(dm, &J);
2592: /* TODO Null space for J */
2593: /* Check discretization error */
2594: DMGetNumFields(dm, &numFields);
2595: PetscMalloc1(PetscMax(1, numFields), &error);
2596: if (numFields > 1) {
2597: PetscInt f;
2599: DMComputeL2FieldDiff(dm, 0.0, exactFuncs, ctxs, u, error);
2600: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2601: for (f = 0; f < numFields; ++f) {
2602: if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2603: if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", (double)error[f]);}
2604: else {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2605: }
2606: PetscPrintf(PETSC_COMM_WORLD, "]\n");
2607: } else {
2608: DMComputeL2Diff(dm, 0.0, exactFuncs, ctxs, u, &error[0]);
2609: if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)error[0]);}
2610: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2611: }
2612: PetscFree(error);
2613: /* Check residual */
2614: SNESComputeFunction(snes, u, r);
2615: VecNorm(r, NORM_2, &res);
2616: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
2617: VecChop(r, 1.0e-10);
2618: PetscObjectSetName((PetscObject) r, "Initial Residual");
2619: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2620: VecViewFromOptions(r, NULL, "-vec_view");
2621: /* Check Jacobian */
2622: DMGetDS(dm, &prob);
2623: PetscDSHasJacobian(prob, &hasJac);
2624: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
2625: if (hasJac && hasPrec) {
2626: DMCreateMatrix(dm, &M);
2627: SNESComputeJacobian(snes, u, J, M);
2628: PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
2629: MatViewFromOptions(M, NULL, "-mat_view");
2630: MatDestroy(&M);
2631: } else {
2632: SNESComputeJacobian(snes, u, J, J);
2633: }
2634: PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
2635: MatViewFromOptions(J, NULL, "-mat_view");
2636: MatGetNullSpace(J, &nullSpace);
2637: if (nullSpace) {
2638: PetscBool isNull;
2639: MatNullSpaceTest(nullSpace, J, &isNull);
2640: if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2641: }
2642: VecDuplicate(u, &b);
2643: VecSet(r, 0.0);
2644: SNESComputeFunction(snes, r, b);
2645: MatMult(J, u, r);
2646: VecAXPY(r, 1.0, b);
2647: VecDestroy(&b);
2648: VecNorm(r, NORM_2, &res);
2649: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
2650: VecChop(r, 1.0e-10);
2651: PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2652: PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2653: VecViewFromOptions(r, NULL, "-vec_view");
2654: VecDestroy(&r);
2655: MatNullSpaceDestroy(&nullSpace);
2656: MatDestroy(&J);
2657: return(0);
2658: }
2660: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2661: {
2662: DM dm;
2663: Vec sol;
2664: PetscBool check;
2668: PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix, "-dmsnes_check", &check);
2669: if (!check) return(0);
2670: SNESGetDM(snes, &dm);
2671: VecDuplicate(u, &sol);
2672: SNESSetSolution(snes, sol);
2673: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2674: VecDestroy(&sol);
2675: return(0);
2676: }