Actual source code: plexfvm.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petscsf.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petsc/private/petscfvimpl.h>
7: static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter lim, PetscInt dim, PetscInt dof, PetscInt cell, PetscInt field, PetscInt face, PetscInt fStart, PetscInt fEnd, PetscReal *cellPhi, const PetscScalar *x, const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
8: {
9: const PetscInt *children;
10: PetscInt numChildren;
12: DMPlexGetTreeChildren(dm, face, &numChildren, &children);
13: if (numChildren) {
14: PetscInt c;
16: for (c = 0; c < numChildren; c++) {
17: PetscInt childFace = children[c];
19: if (childFace >= fStart && childFace < fEnd) DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field, childFace, fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad);
20: }
21: } else {
22: PetscScalar *ncx;
23: PetscFVCellGeom *ncg;
24: const PetscInt *fcells;
25: PetscInt ncell, d;
26: PetscReal v[3];
28: DMPlexGetSupport(dm, face, &fcells);
29: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
30: if (field >= 0) {
31: DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx);
32: } else {
33: DMPlexPointLocalRead(dm, ncell, x, &ncx);
34: }
35: DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);
36: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
37: for (d = 0; d < dof; ++d) {
38: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
39: PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v);
40: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom;
42: PetscLimiterLimit(lim, flim, &phi);
43: cellPhi[d] = PetscMin(cellPhi[d], phi);
44: }
45: }
46: return 0;
47: }
49: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscFV fvm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
50: {
51: DM dmFace, dmCell, dmGrad;
52: DMLabel ghostLabel;
53: PetscDS prob;
54: PetscLimiter lim;
55: const PetscScalar *facegeom, *cellgeom, *x;
56: PetscScalar *gr;
57: PetscReal *cellPhi;
58: PetscInt dim, face, cell, field, dof, cStart, cEnd, nFields;
60: DMGetDimension(dm, &dim);
61: DMGetDS(dm, &prob);
62: PetscDSGetNumFields(prob, &nFields);
63: PetscDSGetFieldIndex(prob, (PetscObject)fvm, &field);
64: PetscDSGetFieldSize(prob, field, &dof);
65: DMGetLabel(dm, "ghost", &ghostLabel);
66: PetscFVGetLimiter(fvm, &lim);
67: VecGetDM(faceGeometry, &dmFace);
68: VecGetArrayRead(faceGeometry, &facegeom);
69: VecGetDM(cellGeometry, &dmCell);
70: VecGetArrayRead(cellGeometry, &cellgeom);
71: VecGetArrayRead(locX, &x);
72: VecGetDM(grad, &dmGrad);
73: VecZeroEntries(grad);
74: VecGetArray(grad, &gr);
75: /* Reconstruct gradients */
76: for (face = fStart; face < fEnd; ++face) {
77: const PetscInt *cells;
78: PetscFVFaceGeom *fg;
79: PetscScalar *cx[2];
80: PetscScalar *cgrad[2];
81: PetscBool boundary;
82: PetscInt ghost, c, pd, d, numChildren, numCells;
84: DMLabelGetValue(ghostLabel, face, &ghost);
85: DMIsBoundaryPoint(dm, face, &boundary);
86: DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
87: if (ghost >= 0 || boundary || numChildren) continue;
88: DMPlexGetSupportSize(dm, face, &numCells);
90: DMPlexGetSupport(dm, face, &cells);
91: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
92: for (c = 0; c < 2; ++c) {
93: if (nFields > 1) {
94: DMPlexPointLocalFieldRead(dm, cells[c], field, x, &cx[c]);
95: } else {
96: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
97: }
98: DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]);
99: }
100: for (pd = 0; pd < dof; ++pd) {
101: PetscScalar delta = cx[1][pd] - cx[0][pd];
103: for (d = 0; d < dim; ++d) {
104: if (cgrad[0]) cgrad[0][pd * dim + d] += fg->grad[0][d] * delta;
105: if (cgrad[1]) cgrad[1][pd * dim + d] -= fg->grad[1][d] * delta;
106: }
107: }
108: }
109: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
110: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
111: DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi);
112: for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEnd; ++cell) {
113: const PetscInt *faces;
114: PetscScalar *cx;
115: PetscFVCellGeom *cg;
116: PetscScalar *cgrad;
117: PetscInt coneSize, f, pd, d;
119: DMPlexGetConeSize(dm, cell, &coneSize);
120: DMPlexGetCone(dm, cell, &faces);
121: if (nFields > 1) {
122: DMPlexPointLocalFieldRead(dm, cell, field, x, &cx);
123: } else {
124: DMPlexPointLocalRead(dm, cell, x, &cx);
125: }
126: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
127: DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
128: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
129: /* Limiter will be minimum value over all neighbors */
130: for (d = 0; d < dof; ++d) cellPhi[d] = PETSC_MAX_REAL;
131: for (f = 0; f < coneSize; ++f) DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, nFields > 1 ? field : -1, faces[f], fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad);
132: /* Apply limiter to gradient */
133: for (pd = 0; pd < dof; ++pd) /* Scalar limiter applied to each component separately */
134: for (d = 0; d < dim; ++d) cgrad[pd * dim + d] *= cellPhi[pd];
135: }
136: DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi);
137: VecRestoreArrayRead(faceGeometry, &facegeom);
138: VecRestoreArrayRead(cellGeometry, &cellgeom);
139: VecRestoreArrayRead(locX, &x);
140: VecRestoreArray(grad, &gr);
141: return 0;
142: }
144: /*@
145: DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.
147: Input Parameters:
148: + dm - the mesh
149: - locX - the local representation of the vector
151: Output Parameter:
152: . grad - the global representation of the gradient
154: Level: developer
156: .seealso: [](chapter_unstructured), `DM`, `Vec`, `DMPlexGetGradientDM()`
157: @*/
158: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
159: {
160: PetscDS prob;
161: PetscInt Nf, f, fStart, fEnd;
162: PetscBool useFVM = PETSC_FALSE;
163: PetscFV fvm = NULL;
164: Vec faceGeometryFVM, cellGeometryFVM;
165: PetscFVCellGeom *cgeomFVM = NULL;
166: PetscFVFaceGeom *fgeomFVM = NULL;
167: DM dmGrad = NULL;
169: DMGetDS(dm, &prob);
170: PetscDSGetNumFields(prob, &Nf);
171: for (f = 0; f < Nf; ++f) {
172: PetscObject obj;
173: PetscClassId id;
175: PetscDSGetDiscretization(prob, f, &obj);
176: PetscObjectGetClassId(obj, &id);
177: if (id == PETSCFV_CLASSID) {
178: useFVM = PETSC_TRUE;
179: fvm = (PetscFV)obj;
180: }
181: }
183: DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad);
185: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **)&fgeomFVM);
186: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **)&cgeomFVM);
187: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
188: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
189: return 0;
190: }