Actual source code: dmplexsnes.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3: #include <petscds.h>
4: #include <petsc/private/petscimpl.h>
5: #include <petsc/private/petscfeimpl.h>
7: /************************** Interpolation *******************************/
11: PetscErrorCode DMInterpolationCreate(MPI_Comm comm, DMInterpolationInfo *ctx)
12: {
17: PetscMalloc(sizeof(struct _DMInterpolationInfo), ctx);
19: (*ctx)->comm = comm;
20: (*ctx)->dim = -1;
21: (*ctx)->nInput = 0;
22: (*ctx)->points = NULL;
23: (*ctx)->cells = NULL;
24: (*ctx)->n = -1;
25: (*ctx)->coords = NULL;
26: return(0);
27: }
31: PetscErrorCode DMInterpolationSetDim(DMInterpolationInfo ctx, PetscInt dim)
32: {
34: if ((dim < 1) || (dim > 3)) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension for points: %d", dim);
35: ctx->dim = dim;
36: return(0);
37: }
41: PetscErrorCode DMInterpolationGetDim(DMInterpolationInfo ctx, PetscInt *dim)
42: {
45: *dim = ctx->dim;
46: return(0);
47: }
51: PetscErrorCode DMInterpolationSetDof(DMInterpolationInfo ctx, PetscInt dof)
52: {
54: if (dof < 1) SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid number of components: %d", dof);
55: ctx->dof = dof;
56: return(0);
57: }
61: PetscErrorCode DMInterpolationGetDof(DMInterpolationInfo ctx, PetscInt *dof)
62: {
65: *dof = ctx->dof;
66: return(0);
67: }
71: PetscErrorCode DMInterpolationAddPoints(DMInterpolationInfo ctx, PetscInt n, PetscReal points[])
72: {
76: if (ctx->dim < 0) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
77: if (ctx->points) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "Cannot add points multiple times yet");
78: ctx->nInput = n;
80: PetscMalloc1(n*ctx->dim, &ctx->points);
81: PetscMemcpy(ctx->points, points, n*ctx->dim * sizeof(PetscReal));
82: return(0);
83: }
87: PetscErrorCode DMInterpolationSetUp(DMInterpolationInfo ctx, DM dm, PetscBool redundantPoints)
88: {
89: MPI_Comm comm = ctx->comm;
90: PetscScalar *a;
91: PetscInt p, q, i;
92: PetscMPIInt rank, size;
94: Vec pointVec;
95: IS cellIS;
96: PetscLayout layout;
97: PetscReal *globalPoints;
98: PetscScalar *globalPointsScalar;
99: const PetscInt *ranges;
100: PetscMPIInt *counts, *displs;
101: const PetscInt *foundCells;
102: PetscMPIInt *foundProcs, *globalProcs;
103: PetscInt n, N;
107: MPI_Comm_size(comm, &size);
108: MPI_Comm_rank(comm, &rank);
109: if (ctx->dim < 0) SETERRQ(comm, PETSC_ERR_ARG_WRONGSTATE, "The spatial dimension has not been set");
110: /* Locate points */
111: n = ctx->nInput;
112: if (!redundantPoints) {
113: PetscLayoutCreate(comm, &layout);
114: PetscLayoutSetBlockSize(layout, 1);
115: PetscLayoutSetLocalSize(layout, n);
116: PetscLayoutSetUp(layout);
117: PetscLayoutGetSize(layout, &N);
118: /* Communicate all points to all processes */
119: PetscMalloc3(N*ctx->dim,&globalPoints,size,&counts,size,&displs);
120: PetscLayoutGetRanges(layout, &ranges);
121: for (p = 0; p < size; ++p) {
122: counts[p] = (ranges[p+1] - ranges[p])*ctx->dim;
123: displs[p] = ranges[p]*ctx->dim;
124: }
125: MPI_Allgatherv(ctx->points, n*ctx->dim, MPIU_REAL, globalPoints, counts, displs, MPIU_REAL, comm);
126: } else {
127: N = n;
128: globalPoints = ctx->points;
129: counts = displs = NULL;
130: layout = NULL;
131: }
132: #if 0
133: PetscMalloc3(N,&foundCells,N,&foundProcs,N,&globalProcs);
134: /* foundCells[p] = m->locatePoint(&globalPoints[p*ctx->dim]); */
135: #else
136: #if defined(PETSC_USE_COMPLEX)
137: PetscMalloc1(N,&globalPointsScalar);
138: for (i=0; i<N; i++) globalPointsScalar[i] = globalPoints[i];
139: #else
140: globalPointsScalar = globalPoints;
141: #endif
142: VecCreateSeqWithArray(PETSC_COMM_SELF, ctx->dim, N*ctx->dim, globalPointsScalar, &pointVec);
143: PetscMalloc2(N,&foundProcs,N,&globalProcs);
144: DMLocatePoints(dm, pointVec, &cellIS);
145: ISGetIndices(cellIS, &foundCells);
146: #endif
147: for (p = 0; p < N; ++p) {
148: if (foundCells[p] >= 0) foundProcs[p] = rank;
149: else foundProcs[p] = size;
150: }
151: /* Let the lowest rank process own each point */
152: MPI_Allreduce(foundProcs, globalProcs, N, MPI_INT, MPI_MIN, comm);
153: ctx->n = 0;
154: for (p = 0; p < N; ++p) {
155: if (globalProcs[p] == size) SETERRQ4(comm, PETSC_ERR_PLIB, "Point %d: %g %g %g not located in mesh", p, globalPoints[p*ctx->dim+0], ctx->dim > 1 ? globalPoints[p*ctx->dim+1] : 0.0, ctx->dim > 2 ? globalPoints[p*ctx->dim+2] : 0.0);
156: else if (globalProcs[p] == rank) ctx->n++;
157: }
158: /* Create coordinates vector and array of owned cells */
159: PetscMalloc1(ctx->n, &ctx->cells);
160: VecCreate(comm, &ctx->coords);
161: VecSetSizes(ctx->coords, ctx->n*ctx->dim, PETSC_DECIDE);
162: VecSetBlockSize(ctx->coords, ctx->dim);
163: VecSetType(ctx->coords,VECSTANDARD);
164: VecGetArray(ctx->coords, &a);
165: for (p = 0, q = 0, i = 0; p < N; ++p) {
166: if (globalProcs[p] == rank) {
167: PetscInt d;
169: for (d = 0; d < ctx->dim; ++d, ++i) a[i] = globalPoints[p*ctx->dim+d];
170: ctx->cells[q++] = foundCells[p];
171: }
172: }
173: VecRestoreArray(ctx->coords, &a);
174: #if 0
175: PetscFree3(foundCells,foundProcs,globalProcs);
176: #else
177: PetscFree2(foundProcs,globalProcs);
178: ISRestoreIndices(cellIS, &foundCells);
179: ISDestroy(&cellIS);
180: VecDestroy(&pointVec);
181: #endif
182: if ((void*)globalPointsScalar != (void*)globalPoints) {PetscFree(globalPointsScalar);}
183: if (!redundantPoints) {PetscFree3(globalPoints,counts,displs);}
184: PetscLayoutDestroy(&layout);
185: return(0);
186: }
190: PetscErrorCode DMInterpolationGetCoordinates(DMInterpolationInfo ctx, Vec *coordinates)
191: {
194: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
195: *coordinates = ctx->coords;
196: return(0);
197: }
201: PetscErrorCode DMInterpolationGetVector(DMInterpolationInfo ctx, Vec *v)
202: {
207: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
208: VecCreate(ctx->comm, v);
209: VecSetSizes(*v, ctx->n*ctx->dof, PETSC_DECIDE);
210: VecSetBlockSize(*v, ctx->dof);
211: VecSetType(*v,VECSTANDARD);
212: return(0);
213: }
217: PetscErrorCode DMInterpolationRestoreVector(DMInterpolationInfo ctx, Vec *v)
218: {
223: if (!ctx->coords) SETERRQ(ctx->comm, PETSC_ERR_ARG_WRONGSTATE, "The interpolation context has not been setup.");
224: VecDestroy(v);
225: return(0);
226: }
230: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Triangle_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
231: {
232: PetscReal *v0, *J, *invJ, detJ;
233: const PetscScalar *coords;
234: PetscScalar *a;
235: PetscInt p;
239: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
240: VecGetArrayRead(ctx->coords, &coords);
241: VecGetArray(v, &a);
242: for (p = 0; p < ctx->n; ++p) {
243: PetscInt c = ctx->cells[p];
244: PetscScalar *x = NULL;
245: PetscReal xi[4];
246: PetscInt d, f, comp;
248: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
249: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
250: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
251: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
253: for (d = 0; d < ctx->dim; ++d) {
254: xi[d] = 0.0;
255: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
256: 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];
257: }
258: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
259: }
260: VecRestoreArray(v, &a);
261: VecRestoreArrayRead(ctx->coords, &coords);
262: PetscFree3(v0, J, invJ);
263: return(0);
264: }
268: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Tetrahedron_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
269: {
270: PetscReal *v0, *J, *invJ, detJ;
271: const PetscScalar *coords;
272: PetscScalar *a;
273: PetscInt p;
277: PetscMalloc3(ctx->dim,&v0,ctx->dim*ctx->dim,&J,ctx->dim*ctx->dim,&invJ);
278: VecGetArrayRead(ctx->coords, &coords);
279: VecGetArray(v, &a);
280: for (p = 0; p < ctx->n; ++p) {
281: PetscInt c = ctx->cells[p];
282: const PetscInt order[3] = {2, 1, 3};
283: PetscScalar *x = NULL;
284: PetscReal xi[4];
285: PetscInt d, f, comp;
287: DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
288: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
289: DMPlexVecGetClosure(dm, NULL, xLocal, c, NULL, &x);
290: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp];
292: for (d = 0; d < ctx->dim; ++d) {
293: xi[d] = 0.0;
294: for (f = 0; f < ctx->dim; ++f) xi[d] += invJ[d*ctx->dim+f]*0.5*PetscRealPart(coords[p*ctx->dim+f] - v0[f]);
295: 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];
296: }
297: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, NULL, &x);
298: }
299: VecRestoreArray(v, &a);
300: VecRestoreArrayRead(ctx->coords, &coords);
301: PetscFree3(v0, J, invJ);
302: return(0);
303: }
307: PETSC_STATIC_INLINE PetscErrorCode QuadMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
308: {
309: const PetscScalar *vertices = (const PetscScalar*) ctx;
310: const PetscScalar x0 = vertices[0];
311: const PetscScalar y0 = vertices[1];
312: const PetscScalar x1 = vertices[2];
313: const PetscScalar y1 = vertices[3];
314: const PetscScalar x2 = vertices[4];
315: const PetscScalar y2 = vertices[5];
316: const PetscScalar x3 = vertices[6];
317: const PetscScalar y3 = vertices[7];
318: const PetscScalar f_1 = x1 - x0;
319: const PetscScalar g_1 = y1 - y0;
320: const PetscScalar f_3 = x3 - x0;
321: const PetscScalar g_3 = y3 - y0;
322: const PetscScalar f_01 = x2 - x1 - x3 + x0;
323: const PetscScalar g_01 = y2 - y1 - y3 + y0;
324: const PetscScalar *ref;
325: PetscScalar *real;
326: PetscErrorCode ierr;
329: VecGetArrayRead(Xref, &ref);
330: VecGetArray(Xreal, &real);
331: {
332: const PetscScalar p0 = ref[0];
333: const PetscScalar p1 = ref[1];
335: real[0] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
336: real[1] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
337: }
338: PetscLogFlops(28);
339: VecRestoreArrayRead(Xref, &ref);
340: VecRestoreArray(Xreal, &real);
341: return(0);
342: }
344: #include <petsc/private/dmimpl.h>
347: PETSC_STATIC_INLINE PetscErrorCode QuadJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
348: {
349: const PetscScalar *vertices = (const PetscScalar*) ctx;
350: const PetscScalar x0 = vertices[0];
351: const PetscScalar y0 = vertices[1];
352: const PetscScalar x1 = vertices[2];
353: const PetscScalar y1 = vertices[3];
354: const PetscScalar x2 = vertices[4];
355: const PetscScalar y2 = vertices[5];
356: const PetscScalar x3 = vertices[6];
357: const PetscScalar y3 = vertices[7];
358: const PetscScalar f_01 = x2 - x1 - x3 + x0;
359: const PetscScalar g_01 = y2 - y1 - y3 + y0;
360: const PetscScalar *ref;
361: PetscErrorCode ierr;
364: VecGetArrayRead(Xref, &ref);
365: {
366: const PetscScalar x = ref[0];
367: const PetscScalar y = ref[1];
368: const PetscInt rows[2] = {0, 1};
369: PetscScalar values[4];
371: values[0] = (x1 - x0 + f_01*y) * 0.5; values[1] = (x3 - x0 + f_01*x) * 0.5;
372: values[2] = (y1 - y0 + g_01*y) * 0.5; values[3] = (y3 - y0 + g_01*x) * 0.5;
373: MatSetValues(J, 2, rows, 2, rows, values, INSERT_VALUES);
374: }
375: PetscLogFlops(30);
376: VecRestoreArrayRead(Xref, &ref);
377: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
378: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
379: return(0);
380: }
384: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Quad_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
385: {
386: DM dmCoord;
387: SNES snes;
388: KSP ksp;
389: PC pc;
390: Vec coordsLocal, r, ref, real;
391: Mat J;
392: const PetscScalar *coords;
393: PetscScalar *a;
394: PetscInt p;
398: DMGetCoordinatesLocal(dm, &coordsLocal);
399: DMGetCoordinateDM(dm, &dmCoord);
400: SNESCreate(PETSC_COMM_SELF, &snes);
401: SNESSetOptionsPrefix(snes, "quad_interp_");
402: VecCreate(PETSC_COMM_SELF, &r);
403: VecSetSizes(r, 2, 2);
404: VecSetType(r,dm->vectype);
405: VecDuplicate(r, &ref);
406: VecDuplicate(r, &real);
407: MatCreate(PETSC_COMM_SELF, &J);
408: MatSetSizes(J, 2, 2, 2, 2);
409: MatSetType(J, MATSEQDENSE);
410: MatSetUp(J);
411: SNESSetFunction(snes, r, QuadMap_Private, NULL);
412: SNESSetJacobian(snes, J, J, QuadJacobian_Private, NULL);
413: SNESGetKSP(snes, &ksp);
414: KSPGetPC(ksp, &pc);
415: PCSetType(pc, PCLU);
416: SNESSetFromOptions(snes);
418: VecGetArrayRead(ctx->coords, &coords);
419: VecGetArray(v, &a);
420: for (p = 0; p < ctx->n; ++p) {
421: PetscScalar *x = NULL, *vertices = NULL;
422: PetscScalar *xi;
423: PetscReal xir[2];
424: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
426: /* Can make this do all points at once */
427: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
428: if (4*2 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 4*2);
429: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
430: if (4*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 4*ctx->dof);
431: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
432: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
433: VecGetArray(real, &xi);
434: xi[0] = coords[p*ctx->dim+0];
435: xi[1] = coords[p*ctx->dim+1];
436: VecRestoreArray(real, &xi);
437: SNESSolve(snes, real, ref);
438: VecGetArray(ref, &xi);
439: xir[0] = PetscRealPart(xi[0]);
440: xir[1] = PetscRealPart(xi[1]);
441: for (comp = 0; comp < ctx->dof; ++comp) a[p*ctx->dof+comp] = x[0*ctx->dof+comp]*(1 - xir[0])*(1 - xir[1]) + x[1*ctx->dof+comp]*xir[0]*(1 - xir[1]) + x[2*ctx->dof+comp]*xir[0]*xir[1] + x[3*ctx->dof+comp]*(1 - xir[0])*xir[1];
443: VecRestoreArray(ref, &xi);
444: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
445: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
446: }
447: VecRestoreArray(v, &a);
448: VecRestoreArrayRead(ctx->coords, &coords);
450: SNESDestroy(&snes);
451: VecDestroy(&r);
452: VecDestroy(&ref);
453: VecDestroy(&real);
454: MatDestroy(&J);
455: return(0);
456: }
460: PETSC_STATIC_INLINE PetscErrorCode HexMap_Private(SNES snes, Vec Xref, Vec Xreal, void *ctx)
461: {
462: const PetscScalar *vertices = (const PetscScalar*) ctx;
463: const PetscScalar x0 = vertices[0];
464: const PetscScalar y0 = vertices[1];
465: const PetscScalar z0 = vertices[2];
466: const PetscScalar x1 = vertices[9];
467: const PetscScalar y1 = vertices[10];
468: const PetscScalar z1 = vertices[11];
469: const PetscScalar x2 = vertices[6];
470: const PetscScalar y2 = vertices[7];
471: const PetscScalar z2 = vertices[8];
472: const PetscScalar x3 = vertices[3];
473: const PetscScalar y3 = vertices[4];
474: const PetscScalar z3 = vertices[5];
475: const PetscScalar x4 = vertices[12];
476: const PetscScalar y4 = vertices[13];
477: const PetscScalar z4 = vertices[14];
478: const PetscScalar x5 = vertices[15];
479: const PetscScalar y5 = vertices[16];
480: const PetscScalar z5 = vertices[17];
481: const PetscScalar x6 = vertices[18];
482: const PetscScalar y6 = vertices[19];
483: const PetscScalar z6 = vertices[20];
484: const PetscScalar x7 = vertices[21];
485: const PetscScalar y7 = vertices[22];
486: const PetscScalar z7 = vertices[23];
487: const PetscScalar f_1 = x1 - x0;
488: const PetscScalar g_1 = y1 - y0;
489: const PetscScalar h_1 = z1 - z0;
490: const PetscScalar f_3 = x3 - x0;
491: const PetscScalar g_3 = y3 - y0;
492: const PetscScalar h_3 = z3 - z0;
493: const PetscScalar f_4 = x4 - x0;
494: const PetscScalar g_4 = y4 - y0;
495: const PetscScalar h_4 = z4 - z0;
496: const PetscScalar f_01 = x2 - x1 - x3 + x0;
497: const PetscScalar g_01 = y2 - y1 - y3 + y0;
498: const PetscScalar h_01 = z2 - z1 - z3 + z0;
499: const PetscScalar f_12 = x7 - x3 - x4 + x0;
500: const PetscScalar g_12 = y7 - y3 - y4 + y0;
501: const PetscScalar h_12 = z7 - z3 - z4 + z0;
502: const PetscScalar f_02 = x5 - x1 - x4 + x0;
503: const PetscScalar g_02 = y5 - y1 - y4 + y0;
504: const PetscScalar h_02 = z5 - z1 - z4 + z0;
505: const PetscScalar f_012 = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
506: const PetscScalar g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
507: const PetscScalar h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
508: const PetscScalar *ref;
509: PetscScalar *real;
510: PetscErrorCode ierr;
513: VecGetArrayRead(Xref, &ref);
514: VecGetArray(Xreal, &real);
515: {
516: const PetscScalar p0 = ref[0];
517: const PetscScalar p1 = ref[1];
518: const PetscScalar p2 = ref[2];
520: 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;
521: 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;
522: 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;
523: }
524: PetscLogFlops(114);
525: VecRestoreArrayRead(Xref, &ref);
526: VecRestoreArray(Xreal, &real);
527: return(0);
528: }
532: PETSC_STATIC_INLINE PetscErrorCode HexJacobian_Private(SNES snes, Vec Xref, Mat J, Mat M, void *ctx)
533: {
534: const PetscScalar *vertices = (const PetscScalar*) ctx;
535: const PetscScalar x0 = vertices[0];
536: const PetscScalar y0 = vertices[1];
537: const PetscScalar z0 = vertices[2];
538: const PetscScalar x1 = vertices[9];
539: const PetscScalar y1 = vertices[10];
540: const PetscScalar z1 = vertices[11];
541: const PetscScalar x2 = vertices[6];
542: const PetscScalar y2 = vertices[7];
543: const PetscScalar z2 = vertices[8];
544: const PetscScalar x3 = vertices[3];
545: const PetscScalar y3 = vertices[4];
546: const PetscScalar z3 = vertices[5];
547: const PetscScalar x4 = vertices[12];
548: const PetscScalar y4 = vertices[13];
549: const PetscScalar z4 = vertices[14];
550: const PetscScalar x5 = vertices[15];
551: const PetscScalar y5 = vertices[16];
552: const PetscScalar z5 = vertices[17];
553: const PetscScalar x6 = vertices[18];
554: const PetscScalar y6 = vertices[19];
555: const PetscScalar z6 = vertices[20];
556: const PetscScalar x7 = vertices[21];
557: const PetscScalar y7 = vertices[22];
558: const PetscScalar z7 = vertices[23];
559: const PetscScalar f_xy = x2 - x1 - x3 + x0;
560: const PetscScalar g_xy = y2 - y1 - y3 + y0;
561: const PetscScalar h_xy = z2 - z1 - z3 + z0;
562: const PetscScalar f_yz = x7 - x3 - x4 + x0;
563: const PetscScalar g_yz = y7 - y3 - y4 + y0;
564: const PetscScalar h_yz = z7 - z3 - z4 + z0;
565: const PetscScalar f_xz = x5 - x1 - x4 + x0;
566: const PetscScalar g_xz = y5 - y1 - y4 + y0;
567: const PetscScalar h_xz = z5 - z1 - z4 + z0;
568: const PetscScalar f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
569: const PetscScalar g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
570: const PetscScalar h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
571: const PetscScalar *ref;
572: PetscErrorCode ierr;
575: VecGetArrayRead(Xref, &ref);
576: {
577: const PetscScalar x = ref[0];
578: const PetscScalar y = ref[1];
579: const PetscScalar z = ref[2];
580: const PetscInt rows[3] = {0, 1, 2};
581: PetscScalar values[9];
583: values[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
584: values[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
585: values[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
586: values[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
587: values[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
588: values[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
589: values[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
590: values[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
591: values[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
593: MatSetValues(J, 3, rows, 3, rows, values, INSERT_VALUES);
594: }
595: PetscLogFlops(152);
596: VecRestoreArrayRead(Xref, &ref);
597: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
598: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
599: return(0);
600: }
604: PETSC_STATIC_INLINE PetscErrorCode DMInterpolate_Hex_Private(DMInterpolationInfo ctx, DM dm, Vec xLocal, Vec v)
605: {
606: DM dmCoord;
607: SNES snes;
608: KSP ksp;
609: PC pc;
610: Vec coordsLocal, r, ref, real;
611: Mat J;
612: const PetscScalar *coords;
613: PetscScalar *a;
614: PetscInt p;
618: DMGetCoordinatesLocal(dm, &coordsLocal);
619: DMGetCoordinateDM(dm, &dmCoord);
620: SNESCreate(PETSC_COMM_SELF, &snes);
621: SNESSetOptionsPrefix(snes, "hex_interp_");
622: VecCreate(PETSC_COMM_SELF, &r);
623: VecSetSizes(r, 3, 3);
624: VecSetType(r,dm->vectype);
625: VecDuplicate(r, &ref);
626: VecDuplicate(r, &real);
627: MatCreate(PETSC_COMM_SELF, &J);
628: MatSetSizes(J, 3, 3, 3, 3);
629: MatSetType(J, MATSEQDENSE);
630: MatSetUp(J);
631: SNESSetFunction(snes, r, HexMap_Private, NULL);
632: SNESSetJacobian(snes, J, J, HexJacobian_Private, NULL);
633: SNESGetKSP(snes, &ksp);
634: KSPGetPC(ksp, &pc);
635: PCSetType(pc, PCLU);
636: SNESSetFromOptions(snes);
638: VecGetArrayRead(ctx->coords, &coords);
639: VecGetArray(v, &a);
640: for (p = 0; p < ctx->n; ++p) {
641: PetscScalar *x = NULL, *vertices = NULL;
642: PetscScalar *xi;
643: PetscReal xir[3];
644: PetscInt c = ctx->cells[p], comp, coordSize, xSize;
646: /* Can make this do all points at once */
647: DMPlexVecGetClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
648: if (8*3 != coordSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", coordSize, 8*3);
649: DMPlexVecGetClosure(dm, NULL, xLocal, c, &xSize, &x);
650: if (8*ctx->dof != xSize) SETERRQ2(ctx->comm, PETSC_ERR_ARG_SIZ, "Invalid closure size %d should be %d", xSize, 8*ctx->dof);
651: SNESSetFunction(snes, NULL, NULL, (void*) vertices);
652: SNESSetJacobian(snes, NULL, NULL, NULL, (void*) vertices);
653: VecGetArray(real, &xi);
654: xi[0] = coords[p*ctx->dim+0];
655: xi[1] = coords[p*ctx->dim+1];
656: xi[2] = coords[p*ctx->dim+2];
657: VecRestoreArray(real, &xi);
658: SNESSolve(snes, real, ref);
659: VecGetArray(ref, &xi);
660: xir[0] = PetscRealPart(xi[0]);
661: xir[1] = PetscRealPart(xi[1]);
662: xir[2] = PetscRealPart(xi[2]);
663: for (comp = 0; comp < ctx->dof; ++comp) {
664: a[p*ctx->dof+comp] =
665: x[0*ctx->dof+comp]*(1-xir[0])*(1-xir[1])*(1-xir[2]) +
666: x[3*ctx->dof+comp]* xir[0]*(1-xir[1])*(1-xir[2]) +
667: x[2*ctx->dof+comp]* xir[0]* xir[1]*(1-xir[2]) +
668: x[1*ctx->dof+comp]*(1-xir[0])* xir[1]*(1-xir[2]) +
669: x[4*ctx->dof+comp]*(1-xir[0])*(1-xir[1])* xir[2] +
670: x[5*ctx->dof+comp]* xir[0]*(1-xir[1])* xir[2] +
671: x[6*ctx->dof+comp]* xir[0]* xir[1]* xir[2] +
672: x[7*ctx->dof+comp]*(1-xir[0])* xir[1]* xir[2];
673: }
674: VecRestoreArray(ref, &xi);
675: DMPlexVecRestoreClosure(dmCoord, NULL, coordsLocal, c, &coordSize, &vertices);
676: DMPlexVecRestoreClosure(dm, NULL, xLocal, c, &xSize, &x);
677: }
678: VecRestoreArray(v, &a);
679: VecRestoreArrayRead(ctx->coords, &coords);
681: SNESDestroy(&snes);
682: VecDestroy(&r);
683: VecDestroy(&ref);
684: VecDestroy(&real);
685: MatDestroy(&J);
686: return(0);
687: }
691: /*
692: Input Parameters:
693: + ctx - The DMInterpolationInfo context
694: . dm - The DM
695: - x - The local vector containing the field to be interpolated
697: Output Parameters:
698: . v - The vector containing the interpolated values
699: */
700: PetscErrorCode DMInterpolationEvaluate(DMInterpolationInfo ctx, DM dm, Vec x, Vec v)
701: {
702: PetscInt dim, coneSize, n;
709: VecGetLocalSize(v, &n);
710: 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);
711: if (n) {
712: DMGetDimension(dm, &dim);
713: DMPlexGetConeSize(dm, ctx->cells[0], &coneSize);
714: if (dim == 2) {
715: if (coneSize == 3) {
716: DMInterpolate_Triangle_Private(ctx, dm, x, v);
717: } else if (coneSize == 4) {
718: DMInterpolate_Quad_Private(ctx, dm, x, v);
719: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
720: } else if (dim == 3) {
721: if (coneSize == 4) {
722: DMInterpolate_Tetrahedron_Private(ctx, dm, x, v);
723: } else {
724: DMInterpolate_Hex_Private(ctx, dm, x, v);
725: }
726: } else SETERRQ1(ctx->comm, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported dimension %d for point interpolation", dim);
727: }
728: return(0);
729: }
733: PetscErrorCode DMInterpolationDestroy(DMInterpolationInfo *ctx)
734: {
739: VecDestroy(&(*ctx)->coords);
740: PetscFree((*ctx)->points);
741: PetscFree((*ctx)->cells);
742: PetscFree(*ctx);
743: *ctx = NULL;
744: return(0);
745: }
749: /*@C
750: SNESMonitorFields - Monitors the residual for each field separately
752: Collective on SNES
754: Input Parameters:
755: + snes - the SNES context
756: . its - iteration number
757: . fgnorm - 2-norm of residual
758: - dummy - unused context
760: Notes:
761: This routine prints the residual norm at each iteration.
763: Level: intermediate
765: .keywords: SNES, nonlinear, default, monitor, norm
766: .seealso: SNESMonitorSet(), SNESMonitorDefault()
767: @*/
768: PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
769: {
770: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) snes));
771: Vec res;
772: DM dm;
773: PetscSection s;
774: const PetscScalar *r;
775: PetscReal *lnorms, *norms;
776: PetscInt numFields, f, pStart, pEnd, p;
777: PetscErrorCode ierr;
780: SNESGetFunction(snes, &res, 0, 0);
781: SNESGetDM(snes, &dm);
782: DMGetDefaultSection(dm, &s);
783: PetscSectionGetNumFields(s, &numFields);
784: PetscSectionGetChart(s, &pStart, &pEnd);
785: PetscCalloc2(numFields, &lnorms, numFields, &norms);
786: VecGetArrayRead(res, &r);
787: for (p = pStart; p < pEnd; ++p) {
788: for (f = 0; f < numFields; ++f) {
789: PetscInt fdof, foff, d;
791: PetscSectionGetFieldDof(s, p, f, &fdof);
792: PetscSectionGetFieldOffset(s, p, f, &foff);
793: for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff+d]));
794: }
795: }
796: VecRestoreArrayRead(res, &r);
797: MPI_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));
798: PetscViewerASCIIAddTab(viewer, ((PetscObject) snes)->tablevel);
799: PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e [", its, (double) fgnorm);
800: for (f = 0; f < numFields; ++f) {
801: if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
802: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) PetscSqrtReal(norms[f]));
803: }
804: PetscViewerASCIIPrintf(viewer, "]\n");
805: PetscViewerASCIISubtractTab(viewer, ((PetscObject) snes)->tablevel);
806: PetscFree2(lnorms, norms);
807: return(0);
808: }
810: /********************* Residual Computation **************************/
814: /*@
815: DMPlexSNESGetGeometryFEM - Return precomputed geometric data
817: Input Parameter:
818: . dm - The DM
820: Output Parameters:
821: . cellgeom - The values precomputed from cell geometry
823: Level: developer
825: .seealso: DMPlexSNESSetFunctionLocal()
826: @*/
827: PetscErrorCode DMPlexSNESGetGeometryFEM(DM dm, Vec *cellgeom)
828: {
829: DMSNES dmsnes;
830: PetscObject obj;
835: DMGetDMSNES(dm, &dmsnes);
836: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", &obj);
837: if (!obj) {
838: Vec cellgeom;
840: DMPlexComputeGeometryFEM(dm, &cellgeom);
841: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fem", (PetscObject) cellgeom);
842: VecDestroy(&cellgeom);
843: }
845: return(0);
846: }
850: /*@
851: DMPlexSNESGetGeometryFVM - Return precomputed geometric data
853: Input Parameter:
854: . dm - The DM
856: Output Parameters:
857: + facegeom - The values precomputed from face geometry
858: . cellgeom - The values precomputed from cell geometry
859: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
861: Level: developer
863: .seealso: DMPlexTSSetRHSFunctionLocal()
864: @*/
865: PetscErrorCode DMPlexSNESGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
866: {
867: DMSNES dmsnes;
868: PetscObject obj;
873: DMGetDMSNES(dm, &dmsnes);
874: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", &obj);
875: if (!obj) {
876: Vec cellgeom, facegeom;
878: DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);
879: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_facegeom_fvm", (PetscObject) facegeom);
880: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_cellgeom_fvm", (PetscObject) cellgeom);
881: VecDestroy(&facegeom);
882: VecDestroy(&cellgeom);
883: }
886: if (minRadius) {DMPlexGetMinRadius(dm, minRadius);}
887: return(0);
888: }
892: /*@
893: DMPlexSNESGetGradientDM - Return gradient data layout
895: Input Parameters:
896: + dm - The DM
897: - fv - The PetscFV
899: Output Parameter:
900: . dmGrad - The layout for gradient values
902: Level: developer
904: .seealso: DMPlexSNESGetGeometryFVM()
905: @*/
906: PetscErrorCode DMPlexSNESGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
907: {
908: DMSNES dmsnes;
909: PetscObject obj;
910: PetscBool computeGradients;
917: PetscFVGetComputeGradients(fv, &computeGradients);
918: if (!computeGradients) {*dmGrad = NULL; return(0);}
919: DMGetDMSNES(dm, &dmsnes);
920: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", &obj);
921: if (!obj) {
922: DM dmGrad;
923: Vec faceGeometry, cellGeometry;
925: DMPlexSNESGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);
926: DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);
927: PetscObjectCompose((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject) dmGrad);
928: DMDestroy(&dmGrad);
929: }
930: PetscObjectQuery((PetscObject) dmsnes, "DMPlexSNES_dmgrad_fvm", (PetscObject *) dmGrad);
931: return(0);
932: }
936: /*@C
937: DMPlexGetCellFields - Retrieve the field values values for a chunk of cells
939: Input Parameters:
940: + dm - The DM
941: . cStart - The first cell to include
942: . cEnd - The first cell to exclude
943: . locX - A local vector with the solution fields
944: . locX_t - A local vector with solution field time derivatives, or NULL
945: - locA - A local vector with auxiliary fields, or NULL
947: Output Parameters:
948: + u - The field coefficients
949: . u_t - The fields derivative coefficients
950: - a - The auxiliary field coefficients
952: Level: developer
954: .seealso: DMPlexGetFaceFields()
955: @*/
956: PetscErrorCode DMPlexGetCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
957: {
958: DM dmAux;
959: PetscSection section, sectionAux;
960: PetscDS prob;
961: PetscInt numCells = cEnd - cStart, totDim, totDimAux, c;
972: DMGetDefaultSection(dm, §ion);
973: DMGetDS(dm, &prob);
974: PetscDSGetTotalDimension(prob, &totDim);
975: if (locA) {
976: PetscDS probAux;
978: VecGetDM(locA, &dmAux);
979: DMGetDefaultSection(dmAux, §ionAux);
980: DMGetDS(dmAux, &probAux);
981: PetscDSGetTotalDimension(probAux, &totDimAux);
982: }
983: DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u);
984: if (locX_t) {DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, u_t);} else {*u_t = NULL;}
985: if (locA) {DMGetWorkArray(dm, numCells*totDimAux, PETSC_SCALAR, a);} else {*a = NULL;}
986: for (c = cStart; c < cEnd; ++c) {
987: PetscScalar *x = NULL, *x_t = NULL, *ul = *u, *ul_t = *u_t, *al = *a;
988: PetscInt i;
990: DMPlexVecGetClosure(dm, section, locX, c, NULL, &x);
991: for (i = 0; i < totDim; ++i) ul[(c-cStart)*totDim+i] = x[i];
992: DMPlexVecRestoreClosure(dm, section, locX, c, NULL, &x);
993: if (locX_t) {
994: DMPlexVecGetClosure(dm, section, locX_t, c, NULL, &x_t);
995: for (i = 0; i < totDim; ++i) ul_t[(c-cStart)*totDim+i] = x_t[i];
996: DMPlexVecRestoreClosure(dm, section, locX_t, c, NULL, &x_t);
997: }
998: if (locA) {
999: DMPlexVecGetClosure(dmAux, sectionAux, locA, c, NULL, &x);
1000: for (i = 0; i < totDimAux; ++i) al[(c-cStart)*totDimAux+i] = x[i];
1001: DMPlexVecRestoreClosure(dmAux, sectionAux, locA, c, NULL, &x);
1002: }
1003: }
1004: return(0);
1005: }
1009: /*@C
1010: DMPlexRestoreCellFields - Restore the field values values for a chunk of cells
1012: Input Parameters:
1013: + dm - The DM
1014: . cStart - The first cell to include
1015: . cEnd - The first cell to exclude
1016: . locX - A local vector with the solution fields
1017: . locX_t - A local vector with solution field time derivatives, or NULL
1018: - locA - A local vector with auxiliary fields, or NULL
1020: Output Parameters:
1021: + u - The field coefficients
1022: . u_t - The fields derivative coefficients
1023: - a - The auxiliary field coefficients
1025: Level: developer
1027: .seealso: DMPlexGetFaceFields()
1028: @*/
1029: PetscErrorCode DMPlexRestoreCellFields(DM dm, PetscInt cStart, PetscInt cEnd, Vec locX, Vec locX_t, Vec locA, PetscScalar **u, PetscScalar **u_t, PetscScalar **a)
1030: {
1034: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u);
1035: if (*u_t) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, u_t);}
1036: if (*a) {DMRestoreWorkArray(dm, 0, PETSC_SCALAR, a);}
1037: return(0);
1038: }
1042: /*@C
1043: DMPlexGetFaceFields - Retrieve the field values values for a chunk of faces
1045: Input Parameters:
1046: + dm - The DM
1047: . fStart - The first face to include
1048: . fEnd - The first face to exclude
1049: . locX - A local vector with the solution fields
1050: . locX_t - A local vector with solution field time derivatives, or NULL
1051: . faceGeometry - A local vector with face geometry
1052: . cellGeometry - A local vector with cell geometry
1053: - locaGrad - A local vector with field gradients, or NULL
1055: Output Parameters:
1056: + uL - The field values at the left side of the face
1057: - uR - The field values at the right side of the face
1059: Level: developer
1061: .seealso: DMPlexGetCellFields()
1062: @*/
1063: PetscErrorCode DMPlexGetFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1064: {
1065: DM dmFace, dmCell, dmGrad = NULL;
1066: PetscSection section;
1067: PetscDS prob;
1068: DMLabel ghostLabel;
1069: const PetscScalar *facegeom, *cellgeom, *x, *lgrad;
1070: PetscBool *isFE;
1071: PetscInt dim, Nf, f, Nc, numFaces = fEnd - fStart, iface, face;
1072: PetscErrorCode ierr;
1083: DMGetDimension(dm, &dim);
1084: DMGetDS(dm, &prob);
1085: DMGetDefaultSection(dm, §ion);
1086: PetscDSGetNumFields(prob, &Nf);
1087: PetscDSGetTotalComponents(prob, &Nc);
1088: PetscMalloc1(Nf, &isFE);
1089: for (f = 0; f < Nf; ++f) {
1090: PetscObject obj;
1091: PetscClassId id;
1093: DMGetField(dm, f, &obj);
1094: PetscObjectGetClassId(obj, &id);
1095: if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;}
1096: else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;}
1097: else {isFE[f] = PETSC_FALSE;}
1098: }
1099: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1100: VecGetArrayRead(locX, &x);
1101: VecGetDM(faceGeometry, &dmFace);
1102: VecGetArrayRead(faceGeometry, &facegeom);
1103: VecGetDM(cellGeometry, &dmCell);
1104: VecGetArrayRead(cellGeometry, &cellgeom);
1105: if (locGrad) {
1106: VecGetDM(locGrad, &dmGrad);
1107: VecGetArrayRead(locGrad, &lgrad);
1108: }
1109: DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uL);
1110: DMGetWorkArray(dm, numFaces*Nc, PETSC_SCALAR, uR);
1111: /* Right now just eat the extra work for FE (could make a cell loop) */
1112: for (face = fStart, iface = 0; face < fEnd; ++face) {
1113: const PetscInt *cells;
1114: const PetscFVFaceGeom *fg;
1115: const PetscFVCellGeom *cgL, *cgR;
1116: const PetscScalar *xL, *xR, *gL, *gR;
1117: PetscScalar *uLl = *uL, *uRl = *uR;
1118: PetscInt ghost;
1120: DMLabelGetValue(ghostLabel, face, &ghost);
1121: if (ghost >= 0) continue;
1122: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1123: DMPlexGetSupport(dm, face, &cells);
1124: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1125: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1126: for (f = 0; f < Nf; ++f) {
1127: PetscInt off;
1129: PetscDSGetComponentOffset(prob, f, &off);
1130: if (isFE[f]) {
1131: const PetscInt *cone;
1132: PetscInt comp, coneSize, faceLocL, faceLocR, ldof, rdof, d;
1134: xL = xR = NULL;
1135: DMPlexVecGetClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1136: DMPlexVecGetClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1137: DMPlexGetCone(dm, cells[0], &cone);
1138: DMPlexGetConeSize(dm, cells[0], &coneSize);
1139: for (faceLocL = 0; faceLocL < coneSize; ++faceLocL) if (cone[faceLocL] == face) break;
1140: if (faceLocL == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[0]);
1141: DMPlexGetCone(dm, cells[1], &cone);
1142: DMPlexGetConeSize(dm, cells[1], &coneSize);
1143: for (faceLocR = 0; faceLocR < coneSize; ++faceLocR) if (cone[faceLocR] == face) break;
1144: if (faceLocR == coneSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not find face %d in cone of cell %d", face, cells[1]);
1145: /* Check that FEM field has values in the right cell (sometimes its an FV ghost cell) */
1146: EvaluateFaceFields(prob, f, faceLocL, xL, &uLl[iface*Nc+off]);
1147: if (rdof == ldof) {EvaluateFaceFields(prob, f, faceLocR, xR, &uRl[iface*Nc+off]);}
1148: else {PetscSectionGetFieldComponents(section, f, &comp); for(d = 0; d < comp; ++d) uRl[iface*Nc+off+d] = uLl[iface*Nc+off+d];}
1149: DMPlexVecRestoreClosure(dm, section, locX, cells[0], &ldof, (PetscScalar **) &xL);
1150: DMPlexVecRestoreClosure(dm, section, locX, cells[1], &rdof, (PetscScalar **) &xR);
1151: } else {
1152: PetscFV fv;
1153: PetscInt numComp, c;
1155: PetscDSGetDiscretization(prob, f, (PetscObject *) &fv);
1156: PetscFVGetNumComponents(fv, &numComp);
1157: DMPlexPointLocalFieldRead(dm, cells[0], f, x, &xL);
1158: DMPlexPointLocalFieldRead(dm, cells[1], f, x, &xR);
1159: if (dmGrad) {
1160: PetscReal dxL[3], dxR[3];
1162: DMPlexPointLocalRead(dmGrad, cells[0], lgrad, &gL);
1163: DMPlexPointLocalRead(dmGrad, cells[1], lgrad, &gR);
1164: DMPlex_WaxpyD_Internal(dim, -1, cgL->centroid, fg->centroid, dxL);
1165: DMPlex_WaxpyD_Internal(dim, -1, cgR->centroid, fg->centroid, dxR);
1166: for (c = 0; c < numComp; ++c) {
1167: uLl[iface*Nc+off+c] = xL[c] + DMPlex_DotD_Internal(dim, &gL[c*dim], dxL);
1168: uRl[iface*Nc+off+c] = xR[c] + DMPlex_DotD_Internal(dim, &gR[c*dim], dxR);
1169: }
1170: } else {
1171: for (c = 0; c < numComp; ++c) {
1172: uLl[iface*Nc+off+c] = xL[c];
1173: uRl[iface*Nc+off+c] = xR[c];
1174: }
1175: }
1176: }
1177: }
1178: ++iface;
1179: }
1180: VecRestoreArrayRead(locX, &x);
1181: VecRestoreArrayRead(faceGeometry, &facegeom);
1182: VecRestoreArrayRead(cellGeometry, &cellgeom);
1183: if (locGrad) {
1184: VecRestoreArrayRead(locGrad, &lgrad);
1185: }
1186: PetscFree(isFE);
1187: return(0);
1188: }
1192: /*@C
1193: DMPlexRestoreFaceFields - Restore the field values values for a chunk of faces
1195: Input Parameters:
1196: + dm - The DM
1197: . fStart - The first face to include
1198: . fEnd - The first face to exclude
1199: . locX - A local vector with the solution fields
1200: . locX_t - A local vector with solution field time derivatives, or NULL
1201: . faceGeometry - A local vector with face geometry
1202: . cellGeometry - A local vector with cell geometry
1203: - locaGrad - A local vector with field gradients, or NULL
1205: Output Parameters:
1206: + uL - The field values at the left side of the face
1207: - uR - The field values at the right side of the face
1209: Level: developer
1211: .seealso: DMPlexGetFaceFields()
1212: @*/
1213: PetscErrorCode DMPlexRestoreFaceFields(DM dm, PetscInt fStart, PetscInt fEnd, Vec locX, Vec locX_t, Vec faceGeometry, Vec cellGeometry, Vec locGrad, PetscScalar **uL, PetscScalar **uR)
1214: {
1218: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uL);
1219: DMRestoreWorkArray(dm, 0, PETSC_SCALAR, uR);
1220: return(0);
1221: }
1225: /*@C
1226: DMPlexGetFaceGeometry - Retrieve the geometric values for a chunk of faces
1228: Input Parameters:
1229: + dm - The DM
1230: . fStart - The first face to include
1231: . fEnd - The first face to exclude
1232: . faceGeometry - A local vector with face geometry
1233: - cellGeometry - A local vector with cell geometry
1235: Output Parameters:
1236: + fgeom - The extract the face centroid and normal
1237: - vol - The cell volume
1239: Level: developer
1241: .seealso: DMPlexGetCellFields()
1242: @*/
1243: PetscErrorCode DMPlexGetFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1244: {
1245: DM dmFace, dmCell;
1246: DMLabel ghostLabel;
1247: const PetscScalar *facegeom, *cellgeom;
1248: PetscInt dim, numFaces = fEnd - fStart, iface, face;
1249: PetscErrorCode ierr;
1257: DMGetDimension(dm, &dim);
1258: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1259: VecGetDM(faceGeometry, &dmFace);
1260: VecGetArrayRead(faceGeometry, &facegeom);
1261: VecGetDM(cellGeometry, &dmCell);
1262: VecGetArrayRead(cellGeometry, &cellgeom);
1263: PetscMalloc1(numFaces, fgeom);
1264: DMGetWorkArray(dm, numFaces*2, PETSC_SCALAR, vol);
1265: for (face = fStart, iface = 0; face < fEnd; ++face) {
1266: const PetscInt *cells;
1267: const PetscFVFaceGeom *fg;
1268: const PetscFVCellGeom *cgL, *cgR;
1269: PetscFVFaceGeom *fgeoml = *fgeom;
1270: PetscReal *voll = *vol;
1271: PetscInt ghost, d;
1273: DMLabelGetValue(ghostLabel, face, &ghost);
1274: if (ghost >= 0) continue;
1275: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1276: DMPlexGetSupport(dm, face, &cells);
1277: DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);
1278: DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);
1279: for (d = 0; d < dim; ++d) {
1280: fgeoml[iface].centroid[d] = fg->centroid[d];
1281: fgeoml[iface].normal[d] = fg->normal[d];
1282: }
1283: voll[iface*2+0] = cgL->volume;
1284: voll[iface*2+1] = cgR->volume;
1285: ++iface;
1286: }
1287: VecRestoreArrayRead(faceGeometry, &facegeom);
1288: VecRestoreArrayRead(cellGeometry, &cellgeom);
1289: return(0);
1290: }
1294: /*@C
1295: DMPlexRestoreFaceGeometry - Restore the field values values for a chunk of faces
1297: Input Parameters:
1298: + dm - The DM
1299: . fStart - The first face to include
1300: . fEnd - The first face to exclude
1301: . faceGeometry - A local vector with face geometry
1302: - cellGeometry - A local vector with cell geometry
1304: Output Parameters:
1305: + fgeom - The extract the face centroid and normal
1306: - vol - The cell volume
1308: Level: developer
1310: .seealso: DMPlexGetFaceFields()
1311: @*/
1312: PetscErrorCode DMPlexRestoreFaceGeometry(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, PetscFVFaceGeom **fgeom, PetscReal **vol)
1313: {
1317: PetscFree(*fgeom);
1318: DMRestoreWorkArray(dm, 0, PETSC_REAL, vol);
1319: return(0);
1320: }
1324: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
1325: {
1326: DM dmFace, dmCell, dmGrad;
1327: DMLabel ghostLabel;
1328: PetscDS prob;
1329: PetscFV fvm;
1330: PetscLimiter lim;
1331: const PetscScalar *facegeom, *cellgeom, *x;
1332: PetscScalar *gr;
1333: PetscReal *cellPhi;
1334: PetscInt dim, face, cell, totDim, cStart, cEnd, cEndInterior;
1335: PetscErrorCode ierr;
1338: DMGetDimension(dm, &dim);
1339: DMGetDS(dm, &prob);
1340: PetscDSGetTotalDimension(prob, &totDim);
1341: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1342: PetscDSGetDiscretization(prob, 0, (PetscObject *) &fvm);
1343: PetscFVGetLimiter(fvm, &lim);
1344: VecGetDM(faceGeometry, &dmFace);
1345: VecGetArrayRead(faceGeometry, &facegeom);
1346: VecGetDM(cellGeometry, &dmCell);
1347: VecGetArrayRead(cellGeometry, &cellgeom);
1348: VecGetArrayRead(locX, &x);
1349: VecGetDM(grad, &dmGrad);
1350: VecZeroEntries(grad);
1351: VecGetArray(grad, &gr);
1352: /* Reconstruct gradients */
1353: for (face = fStart; face < fEnd; ++face) {
1354: const PetscInt *cells;
1355: const PetscFVFaceGeom *fg;
1356: const PetscScalar *cx[2];
1357: PetscScalar *cgrad[2];
1358: PetscBool boundary;
1359: PetscInt ghost, c, pd, d;
1361: DMLabelGetValue(ghostLabel, face, &ghost);
1362: if (ghost >= 0) continue;
1363: DMPlexIsBoundaryPoint(dm, face, &boundary);
1364: if (boundary) continue;
1365: DMPlexGetSupport(dm, face, &cells);
1366: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
1367: for (c = 0; c < 2; ++c) {
1368: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
1369: DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
1370: }
1371: for (pd = 0; pd < totDim; ++pd) {
1372: PetscScalar delta = cx[1][pd] - cx[0][pd];
1374: for (d = 0; d < dim; ++d) {
1375: if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta;
1376: if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta;
1377: }
1378: }
1379: }
1380: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
1381: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1382: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1383: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
1384: DMGetWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1385: for (cell = dmGrad && lim ? cStart : cEnd; cell < cEndInterior; ++cell) {
1386: const PetscInt *faces;
1387: const PetscScalar *cx;
1388: const PetscFVCellGeom *cg;
1389: PetscScalar *cgrad;
1390: PetscInt coneSize, f, pd, d;
1392: DMPlexGetConeSize(dm, cell, &coneSize);
1393: DMPlexGetCone(dm, cell, &faces);
1394: DMPlexPointLocalRead(dm, cell, x, &cx);
1395: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
1396: DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
1397: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
1398: /* Limiter will be minimum value over all neighbors */
1399: for (d = 0; d < totDim; ++d) cellPhi[d] = PETSC_MAX_REAL;
1400: for (f = 0; f < coneSize; ++f) {
1401: const PetscScalar *ncx;
1402: const PetscFVCellGeom *ncg;
1403: const PetscInt *fcells;
1404: PetscInt face = faces[f], ncell, ghost;
1405: PetscReal v[3];
1406: PetscBool boundary;
1408: DMLabelGetValue(ghostLabel, face, &ghost);
1409: DMPlexIsBoundaryPoint(dm, face, &boundary);
1410: if ((ghost >= 0) || boundary) continue;
1411: DMPlexGetSupport(dm, face, &fcells);
1412: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
1413: DMPlexPointLocalRead(dm, ncell, x, &ncx);
1414: DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
1415: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
1416: for (d = 0; d < totDim; ++d) {
1417: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
1418: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / DMPlex_DotD_Internal(dim, &cgrad[d*dim], v);
1420: PetscLimiterLimit(lim, flim, &phi);
1421: cellPhi[d] = PetscMin(cellPhi[d], phi);
1422: }
1423: }
1424: /* Apply limiter to gradient */
1425: for (pd = 0; pd < totDim; ++pd)
1426: /* Scalar limiter applied to each component separately */
1427: for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
1428: }
1429: DMRestoreWorkArray(dm, totDim, PETSC_REAL, &cellPhi);
1430: VecRestoreArrayRead(faceGeometry, &facegeom);
1431: VecRestoreArrayRead(cellGeometry, &cellgeom);
1432: VecRestoreArrayRead(locX, &x);
1433: VecRestoreArray(grad, &gr);
1434: return(0);
1435: }
1439: PetscErrorCode DMPlexComputeBdResidual_Internal(DM dm, Vec locX, Vec locX_t, Vec locF, void *user)
1440: {
1441: DM_Plex *mesh = (DM_Plex *) dm->data;
1442: PetscSection section;
1443: PetscDS prob;
1444: DMLabel depth;
1445: PetscFECellGeom *cgeom;
1446: PetscScalar *u = NULL, *u_t = NULL, *elemVec = NULL;
1447: PetscInt dim, Nf, f, totDimBd, numBd, bd;
1448: PetscErrorCode ierr;
1451: DMGetDimension(dm, &dim);
1452: DMGetDefaultSection(dm, §ion);
1453: DMGetDS(dm, &prob);
1454: PetscDSGetNumFields(prob, &Nf);
1455: PetscDSGetTotalBdDimension(prob, &totDimBd);
1456: DMPlexGetDepthLabel(dm, &depth);
1457: DMPlexGetNumBoundary(dm, &numBd);
1458: for (bd = 0; bd < numBd; ++bd) {
1459: const char *bdLabel;
1460: DMLabel label;
1461: IS pointIS;
1462: const PetscInt *points;
1463: const PetscInt *values;
1464: PetscInt field, numValues, numPoints, p, dep, numFaces;
1465: PetscBool isEssential;
1467: DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
1468: if (isEssential) continue;
1469: if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
1470: DMPlexGetLabel(dm, bdLabel, &label);
1471: DMLabelGetStratumSize(label, 1, &numPoints);
1472: DMLabelGetStratumIS(label, 1, &pointIS);
1473: ISGetIndices(pointIS, &points);
1474: for (p = 0, numFaces = 0; p < numPoints; ++p) {
1475: DMLabelGetValue(depth, points[p], &dep);
1476: if (dep == dim-1) ++numFaces;
1477: }
1478: PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd,&elemVec);
1479: if (locX_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
1480: for (p = 0, f = 0; p < numPoints; ++p) {
1481: const PetscInt point = points[p];
1482: PetscScalar *x = NULL;
1483: PetscInt i;
1485: DMLabelGetValue(depth, points[p], &dep);
1486: if (dep != dim-1) continue;
1487: DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
1488: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
1489: if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
1490: DMPlexVecGetClosure(dm, section, locX, point, NULL, &x);
1491: for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
1492: DMPlexVecRestoreClosure(dm, section, locX, point, NULL, &x);
1493: if (locX_t) {
1494: DMPlexVecGetClosure(dm, section, locX_t, point, NULL, &x);
1495: for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
1496: DMPlexVecRestoreClosure(dm, section, locX_t, point, NULL, &x);
1497: }
1498: ++f;
1499: }
1500: for (f = 0; f < Nf; ++f) {
1501: PetscFE fe;
1502: PetscQuadrature q;
1503: PetscInt numQuadPoints, Nb;
1504: /* Conforming batches */
1505: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1506: /* Remainder */
1507: PetscInt Nr, offset;
1509: PetscDSGetBdDiscretization(prob, f, (PetscObject *) &fe);
1510: PetscFEGetQuadrature(fe, &q);
1511: PetscFEGetDimension(fe, &Nb);
1512: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1513: PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1514: blockSize = Nb*numQuadPoints;
1515: batchSize = numBlocks * blockSize;
1516: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1517: numChunks = numFaces / (numBatches*batchSize);
1518: Ne = numChunks*numBatches*batchSize;
1519: Nr = numFaces % (numBatches*batchSize);
1520: offset = numFaces - Nr;
1521: PetscFEIntegrateBdResidual(fe, prob, f, Ne, cgeom, u, u_t, NULL, NULL, elemVec);
1522: PetscFEIntegrateBdResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemVec[offset*totDimBd]);
1523: }
1524: for (p = 0, f = 0; p < numPoints; ++p) {
1525: const PetscInt point = points[p];
1527: DMLabelGetValue(depth, point, &dep);
1528: if (dep != dim-1) continue;
1529: if (mesh->printFEM > 1) {DMPrintCellVector(point, "BdResidual", totDimBd, &elemVec[f*totDimBd]);}
1530: DMPlexVecSetClosure(dm, NULL, locF, point, &elemVec[f*totDimBd], ADD_VALUES);
1531: ++f;
1532: }
1533: ISRestoreIndices(pointIS, &points);
1534: ISDestroy(&pointIS);
1535: PetscFree3(u,cgeom,elemVec);
1536: if (locX_t) {PetscFree(u_t);}
1537: }
1538: return(0);
1539: }
1543: PetscErrorCode DMPlexComputeResidual_Internal(DM dm, PetscInt cStart, PetscInt cEnd, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
1544: {
1545: DM_Plex *mesh = (DM_Plex *) dm->data;
1546: const char *name = "Residual";
1547: DM dmAux = NULL;
1548: DM dmGrad = NULL;
1549: DMLabel ghostLabel = NULL;
1550: PetscDS prob = NULL;
1551: PetscDS probAux = NULL;
1552: PetscSection section = NULL;
1553: PetscBool useFEM = PETSC_FALSE;
1554: PetscBool useFVM = PETSC_FALSE;
1555: PetscBool isImplicit = (locX_t || time == PETSC_MIN_REAL) ? PETSC_TRUE : PETSC_FALSE;
1556: PetscFV fvm = NULL;
1557: PetscFECellGeom *cgeomFEM = NULL;
1558: PetscFVCellGeom *cgeomFVM = NULL;
1559: PetscFVFaceGeom *fgeomFVM = NULL;
1560: Vec locA, cellGeometryFEM = NULL, cellGeometryFVM = NULL, faceGeometryFVM = NULL, grad, locGrad = NULL;
1561: PetscScalar *u, *u_t, *a, *uL, *uR;
1562: PetscInt Nf, f, totDim, totDimAux, numChunks, cellChunkSize, faceChunkSize, chunk, fStart, fEnd;
1563: PetscErrorCode ierr;
1566: PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);
1567: /* TODO The places where we have to use isFE are probably the member functions for the PetscDisc class */
1568: /* TODO The FVM geometry is over-manipulated. Make the precalc functions return exactly what we need */
1569: /* FEM+FVM */
1570: /* 1: Get sizes from dm and dmAux */
1571: DMGetDefaultSection(dm, §ion);
1572: DMPlexGetLabel(dm, "ghost", &ghostLabel);
1573: DMGetDS(dm, &prob);
1574: PetscDSGetNumFields(prob, &Nf);
1575: PetscDSGetTotalDimension(prob, &totDim);
1576: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &locA);
1577: if (locA) {
1578: VecGetDM(locA, &dmAux);
1579: DMGetDS(dmAux, &probAux);
1580: PetscDSGetTotalDimension(probAux, &totDimAux);
1581: }
1582: /* 2: Get geometric data */
1583: for (f = 0; f < Nf; ++f) {
1584: PetscObject obj;
1585: PetscClassId id;
1586: PetscBool fimp;
1588: PetscDSGetImplicit(prob, f, &fimp);
1589: if (isImplicit != fimp) continue;
1590: PetscDSGetDiscretization(prob, f, &obj);
1591: PetscObjectGetClassId(obj, &id);
1592: if (id == PETSCFE_CLASSID) {useFEM = PETSC_TRUE;}
1593: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
1594: }
1595: if (useFEM) {
1596: DMPlexSNESGetGeometryFEM(dm, &cellGeometryFEM);
1597: VecGetArrayRead(cellGeometryFEM, (const PetscScalar **) &cgeomFEM);
1598: }
1599: if (useFVM) {
1600: DMPlexSNESGetGeometryFVM(dm, &faceGeometryFVM, &cellGeometryFVM, NULL);
1601: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
1602: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
1603: /* Reconstruct and limit cell gradients */
1604: DMPlexSNESGetGradientDM(dm, fvm, &dmGrad);
1605: if (dmGrad) {
1606: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1607: DMGetGlobalVector(dmGrad, &grad);
1608: DMPlexReconstructGradients_Internal(dm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
1609: /* Communicate gradient values */
1610: DMGetLocalVector(dmGrad, &locGrad);
1611: DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);
1612: DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);
1613: DMRestoreGlobalVector(dmGrad, &grad);
1614: }
1615: }
1616: /* Handle boundary values */
1617: DMPlexInsertBoundaryValues(dm, locX, time, faceGeometryFVM, cellGeometryFVM, locGrad);
1618: /* Loop over chunks */
1619: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1620: numChunks = 1;
1621: cellChunkSize = (cEnd - cStart)/numChunks;
1622: faceChunkSize = (fEnd - fStart)/numChunks;
1623: for (chunk = 0; chunk < numChunks; ++chunk) {
1624: PetscScalar *elemVec, *fluxL, *fluxR;
1625: PetscReal *vol;
1626: PetscFVFaceGeom *fgeom;
1627: PetscInt cS = cStart+chunk*cellChunkSize, cE = PetscMin(cS+cellChunkSize, cEnd), numCells = cE - cS, cell;
1628: PetscInt fS = fStart+chunk*faceChunkSize, fE = PetscMin(fS+faceChunkSize, fEnd), numFaces = fE - fS, face;
1630: /* Extract field coefficients */
1631: if (useFEM) {
1632: DMPlexGetCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1633: DMGetWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1634: PetscMemzero(elemVec, numCells*totDim * sizeof(PetscScalar));
1635: }
1636: if (useFVM) {
1637: DMPlexGetFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1638: DMPlexGetFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1639: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1640: DMGetWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1641: PetscMemzero(fluxL, numFaces*totDim * sizeof(PetscScalar));
1642: PetscMemzero(fluxR, numFaces*totDim * sizeof(PetscScalar));
1643: }
1644: /* TODO We will interlace both our field coefficients (u, u_t, uL, uR, etc.) and our output (elemVec, fL, fR). I think this works */
1645: /* Loop over fields */
1646: for (f = 0; f < Nf; ++f) {
1647: PetscObject obj;
1648: PetscClassId id;
1649: PetscBool fimp;
1650: PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset;
1652: PetscDSGetImplicit(prob, f, &fimp);
1653: if (isImplicit != fimp) continue;
1654: PetscDSGetDiscretization(prob, f, &obj);
1655: PetscObjectGetClassId(obj, &id);
1656: if (id == PETSCFE_CLASSID) {
1657: PetscFE fe = (PetscFE) obj;
1658: PetscQuadrature q;
1659: PetscInt Nq, Nb;
1661: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1663: PetscFEGetQuadrature(fe, &q);
1664: PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);
1665: PetscFEGetDimension(fe, &Nb);
1666: blockSize = Nb*Nq;
1667: batchSize = numBlocks * blockSize;
1668: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1669: numChunks = numCells / (numBatches*batchSize);
1670: Ne = numChunks*numBatches*batchSize;
1671: Nr = numCells % (numBatches*batchSize);
1672: offset = numCells - Nr;
1673: /* Integrate FE residual to get elemVec (need fields at quadrature points) */
1674: /* 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) */
1675: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeomFEM, u, u_t, probAux, a, elemVec);
1676: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1677: } else if (id == PETSCFV_CLASSID) {
1678: PetscFV fv = (PetscFV) obj;
1680: Ne = numFaces;
1681: Nr = 0;
1682: /* Riemann solve over faces (need fields at face centroids) */
1683: /* We need to evaluate FE fields at those coordinates */
1684: PetscFVIntegrateRHSFunction(fv, prob, f, Ne, fgeom, vol, uL, uR, fluxL, fluxR);
1685: } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
1686: }
1687: /* Loop over domain */
1688: if (useFEM) {
1689: /* Add elemVec to locX */
1690: for (cell = cS; cell < cE; ++cell) {
1691: if (mesh->printFEM > 1) {DMPrintCellVector(cell, name, totDim, &elemVec[cell*totDim]);}
1692: DMPlexVecSetClosure(dm, section, locF, cell, &elemVec[cell*totDim], ADD_VALUES);
1693: }
1694: }
1695: if (useFVM) {
1696: PetscScalar *fa;
1697: PetscInt iface;
1699: VecGetArray(locF, &fa);
1700: for (f = 0; f < Nf; ++f) {
1701: PetscFV fv;
1702: PetscObject obj;
1703: PetscClassId id;
1704: PetscInt foff, pdim;
1706: PetscDSGetDiscretization(prob, f, &obj);
1707: PetscDSGetFieldOffset(prob, f, &foff);
1708: PetscObjectGetClassId(obj, &id);
1709: if (id != PETSCFV_CLASSID) continue;
1710: fv = (PetscFV) obj;
1711: PetscFVGetNumComponents(fv, &pdim);
1712: /* Accumulate fluxes to cells */
1713: for (face = fS, iface = 0; face < fE; ++face) {
1714: const PetscInt *cells;
1715: PetscScalar *fL, *fR;
1716: PetscInt ghost, d;
1718: DMLabelGetValue(ghostLabel, face, &ghost);
1719: if (ghost >= 0) continue;
1720: DMPlexGetSupport(dm, face, &cells);
1721: DMPlexPointGlobalFieldRef(dm, cells[0], f, fa, &fL);
1722: DMPlexPointGlobalFieldRef(dm, cells[1], f, fa, &fR);
1723: for (d = 0; d < pdim; ++d) {
1724: if (fL) fL[d] -= fluxL[iface*totDim+foff+d];
1725: if (fR) fR[d] += fluxR[iface*totDim+foff+d];
1726: }
1727: ++iface;
1728: }
1729: }
1730: VecRestoreArray(locF, &fa);
1731: }
1732: /* Handle time derivative */
1733: if (locX_t) {
1734: PetscScalar *x_t, *fa;
1736: VecGetArray(locF, &fa);
1737: VecGetArray(locX_t, &x_t);
1738: for (f = 0; f < Nf; ++f) {
1739: PetscFV fv;
1740: PetscObject obj;
1741: PetscClassId id;
1742: PetscInt pdim, d;
1744: PetscDSGetDiscretization(prob, f, &obj);
1745: PetscObjectGetClassId(obj, &id);
1746: if (id != PETSCFV_CLASSID) continue;
1747: fv = (PetscFV) obj;
1748: PetscFVGetNumComponents(fv, &pdim);
1749: for (cell = cS; cell < cE; ++cell) {
1750: PetscScalar *u_t, *r;
1752: DMPlexPointLocalFieldRead(dm, cell, f, x_t, &u_t);
1753: DMPlexPointLocalFieldRef(dm, cell, f, fa, &r);
1754: for (d = 0; d < pdim; ++d) r[d] += u_t[d];
1755: }
1756: }
1757: VecRestoreArray(locX_t, &x_t);
1758: VecRestoreArray(locF, &fa);
1759: }
1760: if (useFEM) {
1761: DMPlexRestoreCellFields(dm, cS, cE, locX, locX_t, locA, &u, &u_t, &a);
1762: DMRestoreWorkArray(dm, numCells*totDim, PETSC_SCALAR, &elemVec);
1763: }
1764: if (useFVM) {
1765: DMPlexRestoreFaceFields(dm, fS, fE, locX, locX_t, faceGeometryFVM, cellGeometryFVM, locGrad, &uL, &uR);
1766: DMPlexRestoreFaceGeometry(dm, fS, fE, faceGeometryFVM, cellGeometryFVM, &fgeom, &vol);
1767: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxL);
1768: DMRestoreWorkArray(dm, numFaces*totDim, PETSC_SCALAR, &fluxR);
1769: if (dmGrad) {DMRestoreLocalVector(dmGrad, &locGrad);}
1770: }
1771: }
1773: if (useFEM) {DMPlexComputeBdResidual_Internal(dm, locX, locX_t, locF, user);}
1775: /* FEM */
1776: /* 1: Get sizes from dm and dmAux */
1777: /* 2: Get geometric data */
1778: /* 3: Handle boundary values */
1779: /* 4: Loop over domain */
1780: /* Extract coefficients */
1781: /* Loop over fields */
1782: /* Set tiling for FE*/
1783: /* Integrate FE residual to get elemVec */
1784: /* Loop over subdomain */
1785: /* Loop over quad points */
1786: /* Transform coords to real space */
1787: /* Evaluate field and aux fields at point */
1788: /* Evaluate residual at point */
1789: /* Transform residual to real space */
1790: /* Add residual to elemVec */
1791: /* Loop over domain */
1792: /* Add elemVec to locX */
1794: /* FVM */
1795: /* Get geometric data */
1796: /* If using gradients */
1797: /* Compute gradient data */
1798: /* Loop over domain faces */
1799: /* Count computational faces */
1800: /* Reconstruct cell gradient */
1801: /* Loop over domain cells */
1802: /* Limit cell gradients */
1803: /* Handle boundary values */
1804: /* Loop over domain faces */
1805: /* Read out field, centroid, normal, volume for each side of face */
1806: /* Riemann solve over faces */
1807: /* Loop over domain faces */
1808: /* Accumulate fluxes to cells */
1809: /* TODO Change printFEM to printDisc here */
1810: if (mesh->printFEM) {DMPrintLocalVec(dm, name, mesh->printTol, locF);}
1811: PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);
1812: return(0);
1813: }
1817: static PetscErrorCode DMPlexComputeResidualFEM_Check_Internal(DM dm, Vec X, Vec X_t, Vec F, void *user)
1818: {
1819: DM dmCh, dmAux;
1820: Vec A, cellgeom;
1821: PetscDS prob, probCh, probAux = NULL;
1822: PetscQuadrature q;
1823: PetscSection section, sectionAux;
1824: PetscFECellGeom *cgeom;
1825: PetscScalar *elemVec, *elemVecCh, *u, *u_t, *a = NULL;
1826: PetscInt dim, Nf, f, numCells, cStart, cEnd, c;
1827: PetscInt totDim, totDimAux, diffCell = 0;
1828: PetscErrorCode ierr;
1831: DMGetDimension(dm, &dim);
1832: DMGetDefaultSection(dm, §ion);
1833: DMGetDS(dm, &prob);
1834: PetscDSGetTotalDimension(prob, &totDim);
1835: PetscSectionGetNumFields(section, &Nf);
1836: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1837: numCells = cEnd - cStart;
1838: PetscObjectQuery((PetscObject) dm, "dmCh", (PetscObject *) &dmCh);
1839: DMGetDS(dmCh, &probCh);
1840: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1841: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1842: if (dmAux) {
1843: DMGetDefaultSection(dmAux, §ionAux);
1844: DMGetDS(dmAux, &probAux);
1845: PetscDSGetTotalDimension(probAux, &totDimAux);
1846: }
1847: DMPlexInsertBoundaryValues(dm, X, 0.0, NULL, NULL, NULL);
1848: VecSet(F, 0.0);
1849: PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim,&elemVec);
1850: PetscMalloc1(numCells*totDim,&elemVecCh);
1851: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1852: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1853: VecGetArray(cellgeom, (PetscScalar **) &cgeom);
1854: for (c = cStart; c < cEnd; ++c) {
1855: PetscScalar *x = NULL, *x_t = NULL;
1856: PetscInt i;
1858: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
1859: for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i];
1860: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
1861: if (X_t) {
1862: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
1863: for (i = 0; i < totDim; ++i) u_t[c*totDim+i] = x_t[i];
1864: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
1865: }
1866: if (dmAux) {
1867: DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
1868: for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i];
1869: DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
1870: }
1871: }
1872: for (f = 0; f < Nf; ++f) {
1873: PetscFE fe, feCh;
1874: PetscInt numQuadPoints, Nb;
1875: /* Conforming batches */
1876: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
1877: /* Remainder */
1878: PetscInt Nr, offset;
1880: PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);
1881: PetscDSGetDiscretization(probCh, f, (PetscObject *) &feCh);
1882: PetscFEGetQuadrature(fe, &q);
1883: PetscFEGetDimension(fe, &Nb);
1884: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
1885: PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);
1886: blockSize = Nb*numQuadPoints;
1887: batchSize = numBlocks * blockSize;
1888: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
1889: numChunks = numCells / (numBatches*batchSize);
1890: Ne = numChunks*numBatches*batchSize;
1891: Nr = numCells % (numBatches*batchSize);
1892: offset = numCells - Nr;
1893: PetscFEIntegrateResidual(fe, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVec);
1894: PetscFEIntegrateResidual(feCh, prob, f, Ne, cgeom, u, u_t, probAux, a, elemVecCh);
1895: PetscFEIntegrateResidual(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVec[offset*totDim]);
1896: PetscFEIntegrateResidual(feCh, prob, f, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemVecCh[offset*totDim]);
1897: }
1898: for (c = cStart; c < cEnd; ++c) {
1899: PetscBool diff = PETSC_FALSE;
1900: PetscInt d;
1902: for (d = 0; d < totDim; ++d) if (PetscAbsScalar(elemVec[c*totDim+d] - elemVecCh[c*totDim+d]) > 1.0e-7) {diff = PETSC_TRUE;break;}
1903: if (diff) {
1904: PetscPrintf(PetscObjectComm((PetscObject) dm), "Different cell %d\n", c);
1905: DMPrintCellVector(c, "Residual", totDim, &elemVec[c*totDim]);
1906: DMPrintCellVector(c, "Check Residual", totDim, &elemVecCh[c*totDim]);
1907: ++diffCell;
1908: }
1909: if (diffCell > 9) break;
1910: DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*totDim], ADD_VALUES);
1911: }
1912: VecRestoreArray(cellgeom, (PetscScalar **) &cgeom);
1913: PetscFree3(u,u_t,elemVec);
1914: PetscFree(elemVecCh);
1915: if (dmAux) {PetscFree(a);}
1916: return(0);
1917: }
1921: /*@
1922: DMPlexSNESComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
1924: Input Parameters:
1925: + dm - The mesh
1926: . X - Local solution
1927: - user - The user context
1929: Output Parameter:
1930: . F - Local output vector
1932: Level: developer
1934: .seealso: DMPlexComputeJacobianActionFEM()
1935: @*/
1936: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
1937: {
1938: PetscObject check;
1939: PetscInt cStart, cEnd, cEndInterior;
1943: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1944: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1945: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
1946: /* The dmCh is used to check two mathematically equivalent discretizations for computational equivalence */
1947: PetscObjectQuery((PetscObject) dm, "dmCh", &check);
1948: if (check) {DMPlexComputeResidualFEM_Check_Internal(dm, X, NULL, F, user);}
1949: else {DMPlexComputeResidual_Internal(dm, cStart, cEnd, PETSC_MIN_REAL, X, NULL, F, user);}
1950: return(0);
1951: }
1955: 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)
1956: {
1957: DM_Plex *mesh = (DM_Plex *) dm->data;
1958: const char *name = "Jacobian";
1959: DM dmAux;
1960: DMLabel depth;
1961: Vec A, cellgeom;
1962: PetscDS prob, probAux = NULL;
1963: PetscQuadrature quad;
1964: PetscSection section, globalSection, sectionAux;
1965: PetscFECellGeom *cgeom;
1966: PetscScalar *elemMat, *u, *u_t, *a = NULL;
1967: PetscInt dim, Nf, f, fieldI, fieldJ, numCells, c;
1968: PetscInt totDim, totDimBd, totDimAux, numBd, bd;
1969: PetscBool isShell;
1970: PetscErrorCode ierr;
1973: PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);
1974: DMGetDimension(dm, &dim);
1975: DMGetDefaultSection(dm, §ion);
1976: DMGetDefaultGlobalSection(dm, &globalSection);
1977: DMGetDS(dm, &prob);
1978: PetscDSGetTotalDimension(prob, &totDim);
1979: PetscDSGetTotalBdDimension(prob, &totDimBd);
1980: PetscSectionGetNumFields(section, &Nf);
1981: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1982: numCells = cEnd - cStart;
1983: PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);
1984: PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);
1985: if (dmAux) {
1986: DMGetDefaultSection(dmAux, §ionAux);
1987: DMGetDS(dmAux, &probAux);
1988: PetscDSGetTotalDimension(probAux, &totDimAux);
1989: }
1990: DMPlexInsertBoundaryValues(dm, X, 0.0, NULL, NULL, NULL);
1991: MatZeroEntries(JacP);
1992: PetscMalloc3(numCells*totDim,&u,X_t ? numCells*totDim : 0,&u_t,numCells*totDim*totDim,&elemMat);
1993: if (dmAux) {PetscMalloc1(numCells*totDimAux, &a);}
1994: DMPlexSNESGetGeometryFEM(dm, &cellgeom);
1995: VecGetArray(cellgeom, (PetscScalar **) &cgeom);
1996: for (c = cStart; c < cEnd; ++c) {
1997: PetscScalar *x = NULL, *x_t = NULL;
1998: PetscInt i;
2000: DMPlexVecGetClosure(dm, section, X, c, NULL, &x);
2001: for (i = 0; i < totDim; ++i) u[(c-cStart)*totDim+i] = x[i];
2002: DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);
2003: if (X_t) {
2004: DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);
2005: for (i = 0; i < totDim; ++i) u_t[(c-cStart)*totDim+i] = x_t[i];
2006: DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);
2007: }
2008: if (dmAux) {
2009: DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);
2010: for (i = 0; i < totDimAux; ++i) a[(c-cStart)*totDimAux+i] = x[i];
2011: DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);
2012: }
2013: }
2014: PetscMemzero(elemMat, numCells*totDim*totDim * sizeof(PetscScalar));
2015: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2016: PetscFE fe;
2017: PetscInt numQuadPoints, Nb;
2018: /* Conforming batches */
2019: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2020: /* Remainder */
2021: PetscInt Nr, offset;
2023: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
2024: PetscFEGetQuadrature(fe, &quad);
2025: PetscFEGetDimension(fe, &Nb);
2026: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2027: PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2028: blockSize = Nb*numQuadPoints;
2029: batchSize = numBlocks * blockSize;
2030: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2031: numChunks = numCells / (numBatches*batchSize);
2032: Ne = numChunks*numBatches*batchSize;
2033: Nr = numCells % (numBatches*batchSize);
2034: offset = numCells - Nr;
2035: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2036: PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, probAux, a, elemMat);
2037: PetscFEIntegrateJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDim], u_t ? &u_t[offset*totDim] : NULL, probAux, &a[offset*totDimAux], &elemMat[offset*totDim*totDim]);
2038: }
2039: }
2040: for (c = cStart; c < cEnd; ++c) {
2041: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, name, totDim, totDim, &elemMat[(c-cStart)*totDim*totDim]);}
2042: DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[(c-cStart)*totDim*totDim], ADD_VALUES);
2043: }
2044: VecGetArray(cellgeom, (PetscScalar **) &cgeom);
2045: PetscFree3(u,u_t,elemMat);
2046: if (dmAux) {PetscFree(a);}
2047: DMPlexGetDepthLabel(dm, &depth);
2048: DMPlexGetNumBoundary(dm, &numBd);
2049: DMPlexGetDepthLabel(dm, &depth);
2050: DMPlexGetNumBoundary(dm, &numBd);
2051: for (bd = 0; bd < numBd; ++bd) {
2052: const char *bdLabel;
2053: DMLabel label;
2054: IS pointIS;
2055: const PetscInt *points;
2056: const PetscInt *values;
2057: PetscInt field, numValues, numPoints, p, dep, numFaces;
2058: PetscBool isEssential;
2060: DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, NULL, NULL, &numValues, &values, NULL);
2061: if (isEssential) continue;
2062: if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this");
2063: DMPlexGetLabel(dm, bdLabel, &label);
2064: DMLabelGetStratumSize(label, 1, &numPoints);
2065: DMLabelGetStratumIS(label, 1, &pointIS);
2066: ISGetIndices(pointIS, &points);
2067: for (p = 0, numFaces = 0; p < numPoints; ++p) {
2068: DMLabelGetValue(depth, points[p], &dep);
2069: if (dep == dim-1) ++numFaces;
2070: }
2071: PetscMalloc3(numFaces*totDimBd,&u,numFaces,&cgeom,numFaces*totDimBd*totDimBd,&elemMat);
2072: if (X_t) {PetscMalloc1(numFaces*totDimBd,&u_t);}
2073: for (p = 0, f = 0; p < numPoints; ++p) {
2074: const PetscInt point = points[p];
2075: PetscScalar *x = NULL;
2076: PetscInt i;
2078: DMLabelGetValue(depth, points[p], &dep);
2079: if (dep != dim-1) continue;
2080: DMPlexComputeCellGeometryFEM(dm, point, NULL, cgeom[f].v0, cgeom[f].J, cgeom[f].invJ, &cgeom[f].detJ);
2081: DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, cgeom[f].n);
2082: if (cgeom[f].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", cgeom[f].detJ, point);
2083: DMPlexVecGetClosure(dm, section, X, point, NULL, &x);
2084: for (i = 0; i < totDimBd; ++i) u[f*totDimBd+i] = x[i];
2085: DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);
2086: if (X_t) {
2087: DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);
2088: for (i = 0; i < totDimBd; ++i) u_t[f*totDimBd+i] = x[i];
2089: DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);
2090: }
2091: ++f;
2092: }
2093: PetscMemzero(elemMat, numFaces*totDimBd*totDimBd * sizeof(PetscScalar));
2094: for (fieldI = 0; fieldI < Nf; ++fieldI) {
2095: PetscFE fe;
2096: PetscInt numQuadPoints, Nb;
2097: /* Conforming batches */
2098: PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize;
2099: /* Remainder */
2100: PetscInt Nr, offset;
2102: PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
2103: PetscFEGetQuadrature(fe, &quad);
2104: PetscFEGetDimension(fe, &Nb);
2105: PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);
2106: PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);
2107: blockSize = Nb*numQuadPoints;
2108: batchSize = numBlocks * blockSize;
2109: PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);
2110: numChunks = numFaces / (numBatches*batchSize);
2111: Ne = numChunks*numBatches*batchSize;
2112: Nr = numFaces % (numBatches*batchSize);
2113: offset = numFaces - Nr;
2114: for (fieldJ = 0; fieldJ < Nf; ++fieldJ) {
2115: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Ne, cgeom, u, u_t, NULL, NULL, elemMat);
2116: PetscFEIntegrateBdJacobian(fe, prob, fieldI, fieldJ, Nr, &cgeom[offset], &u[offset*totDimBd], u_t ? &u_t[offset*totDimBd] : NULL, NULL, NULL, &elemMat[offset*totDimBd*totDimBd]);
2117: }
2118: }
2119: for (p = 0, f = 0; p < numPoints; ++p) {
2120: const PetscInt point = points[p];
2122: DMLabelGetValue(depth, point, &dep);
2123: if (dep != dim-1) continue;
2124: if (mesh->printFEM > 1) {DMPrintCellMatrix(point, "BdJacobian", totDimBd, totDimBd, &elemMat[f*totDimBd*totDimBd]);}
2125: DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*totDimBd*totDimBd], ADD_VALUES);
2126: ++f;
2127: }
2128: ISRestoreIndices(pointIS, &points);
2129: ISDestroy(&pointIS);
2130: PetscFree3(u,cgeom,elemMat);
2131: if (X_t) {PetscFree(u_t);}
2132: }
2133: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
2134: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
2135: if (mesh->printFEM) {
2136: PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);
2137: MatChop(JacP, 1.0e-10);
2138: MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
2139: }
2140: PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);
2141: PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);
2142: if (isShell) {
2143: JacActionCtx *jctx;
2145: MatShellGetContext(Jac, &jctx);
2146: VecCopy(X, jctx->u);
2147: }
2148: return(0);
2149: }
2153: /*@
2154: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
2156: Input Parameters:
2157: + dm - The mesh
2158: . X - Local input vector
2159: - user - The user context
2161: Output Parameter:
2162: . Jac - Jacobian matrix
2164: Note:
2165: The first member of the user context must be an FEMContext.
2167: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
2168: like a GPU, or vectorize on a multicore machine.
2170: Level: developer
2172: .seealso: FormFunctionLocal()
2173: @*/
2174: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user)
2175: {
2176: PetscInt cStart, cEnd, cEndInterior;
2180: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2181: DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
2182: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
2183: DMPlexComputeJacobian_Internal(dm, cStart, cEnd, 0.0, 0.0, X, NULL, Jac, JacP, user);
2184: return(0);
2185: }
2189: PetscErrorCode DMSNESCheckFromOptions_Internal(SNES snes, DM dm, Vec u, Vec sol, PetscErrorCode (**exactFuncs)(PetscInt, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx), void **ctxs)
2190: {
2191: Mat J, M;
2192: Vec r, b;
2193: MatNullSpace nullSpace;
2194: PetscReal *error, res = 0.0;
2195: PetscInt numFields;
2199: VecDuplicate(u, &r);
2200: DMCreateMatrix(dm, &J);
2201: M = J;
2202: /* TODO Null space for J */
2203: /* Check discretization error */
2204: DMGetNumFields(dm, &numFields);
2205: PetscMalloc1(PetscMax(1, numFields), &error);
2206: if (numFields > 1) {
2207: PetscInt f;
2209: DMPlexComputeL2FieldDiff(dm, exactFuncs, ctxs, u, error);
2210: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: [");
2211: for (f = 0; f < numFields; ++f) {
2212: if (f) {PetscPrintf(PETSC_COMM_WORLD, ", ");}
2213: if (error[f] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "%g", error[f]);}
2214: else {PetscPrintf(PETSC_COMM_WORLD, "< 1.0e-11");}
2215: }
2216: PetscPrintf(PETSC_COMM_WORLD, "]\n");
2217: } else {
2218: DMPlexComputeL2Diff(dm, exactFuncs, ctxs, u, &error[0]);
2219: if (error[0] >= 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error[0]);}
2220: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
2221: }
2222: PetscFree(error);
2223: /* Check residual */
2224: SNESComputeFunction(snes, u, r);
2225: VecNorm(r, NORM_2, &res);
2226: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
2227: VecChop(r, 1.0e-10);
2228: PetscObjectSetName((PetscObject) r, "Initial Residual");
2229: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
2230: VecViewFromOptions(r, NULL, "-vec_view");
2231: /* Check Jacobian */
2232: SNESComputeJacobian(snes, u, M, M);
2233: MatGetNullSpace(J, &nullSpace);
2234: if (nullSpace) {
2235: PetscBool isNull;
2236: MatNullSpaceTest(nullSpace, J, &isNull);
2237: if (!isNull) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
2238: }
2239: VecDuplicate(u, &b);
2240: VecSet(r, 0.0);
2241: SNESComputeFunction(snes, r, b);
2242: MatMult(M, u, r);
2243: VecAXPY(r, 1.0, b);
2244: VecDestroy(&b);
2245: VecNorm(r, NORM_2, &res);
2246: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", res);
2247: VecChop(r, 1.0e-10);
2248: PetscObjectSetName((PetscObject) r, "Au - b = Au + F(0)");
2249: PetscObjectSetOptionsPrefix((PetscObject)r,"linear_res_");
2250: VecViewFromOptions(r, NULL, "-vec_view");
2251: VecDestroy(&r);
2252: MatNullSpaceDestroy(&nullSpace);
2253: MatDestroy(&J);
2254: return(0);
2255: }
2259: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
2260: {
2261: DM dm;
2262: Vec sol;
2263: PetscBool check;
2267: PetscOptionsHasName(snes->hdr.prefix, "-dmsnes_check", &check);
2268: if (!check) return(0);
2269: SNESGetDM(snes, &dm);
2270: VecDuplicate(u, &sol);
2271: SNESSetSolution(snes, sol);
2272: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
2273: VecDestroy(&sol);
2274: return(0);
2275: }