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: PetscFunctionBegin;
13: PetscCall(DMPlexGetTreeChildren(dm, face, &numChildren, &children));
14: if (numChildren) {
15: PetscInt c;
17: for (c = 0; c < numChildren; c++) {
18: PetscInt childFace = children[c];
20: if (childFace >= fStart && childFace < fEnd) PetscCall(DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, field, childFace, fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad));
21: }
22: } else {
23: PetscScalar *ncx;
24: PetscFVCellGeom *ncg;
25: const PetscInt *fcells;
26: PetscInt Ns, ncell, d;
27: PetscReal v[3];
29: PetscCall(DMPlexGetSupport(dm, face, &fcells));
30: PetscCall(DMPlexGetSupportSize(dm, face, &Ns));
31: if (Ns < 2) {
32: for (d = 0; d < dof; ++d) cellPhi[d] = 1.0;
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
35: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
36: if (field >= 0) PetscCall(DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx));
37: else PetscCall(DMPlexPointLocalRead(dm, ncell, x, &ncx));
38: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg));
39: DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, ncg->centroid, v);
40: for (d = 0; d < dof; ++d) {
41: /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */
42: PetscReal denom = DMPlex_DotD_Internal(dim, &cgrad[d * dim], v);
43: PetscReal fact = denom == 0 ? 1.0e+30 : 1 / denom;
44: PetscReal phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) * fact;
46: PetscCall(PetscLimiterLimit(lim, flim, &phi));
47: cellPhi[d] = PetscMin(cellPhi[d], phi);
48: }
49: }
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscFV fvm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
54: {
55: DM dmFace, dmCell, dmGrad;
56: DMLabel ghostLabel;
57: PetscDS prob;
58: PetscLimiter lim;
59: PetscLimiterType limType;
60: const PetscScalar *facegeom, *cellgeom, *x;
61: PetscScalar *gr;
62: PetscReal *cellPhi;
63: PetscInt dim, face, cell, field, dof, cStart, cEnd, nFields;
65: PetscFunctionBegin;
66: PetscCall(DMGetDimension(dm, &dim));
67: PetscCall(DMGetDS(dm, &prob));
68: PetscCall(PetscDSGetNumFields(prob, &nFields));
69: PetscCall(PetscDSGetFieldIndex(prob, (PetscObject)fvm, &field));
70: PetscCall(PetscDSGetFieldSize(prob, field, &dof));
71: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
72: PetscCall(PetscFVGetLimiter(fvm, &lim));
73: PetscCall(VecGetDM(faceGeometry, &dmFace));
74: PetscCall(VecGetArrayRead(faceGeometry, &facegeom));
75: PetscCall(VecGetDM(cellGeometry, &dmCell));
76: PetscCall(VecGetArrayRead(cellGeometry, &cellgeom));
77: PetscCall(VecGetArrayRead(locX, &x));
78: PetscCall(VecGetDM(grad, &dmGrad));
79: PetscCall(VecZeroEntries(grad));
80: PetscCall(VecGetArray(grad, &gr));
81: /* Reconstruct gradients */
82: for (face = fStart; face < fEnd; ++face) {
83: const PetscInt *cells;
84: PetscFVFaceGeom *fg;
85: PetscScalar *cx[2];
86: PetscScalar *cgrad[2];
87: PetscBool boundary;
88: PetscInt ghost, c, pd, d, numChildren, numCells;
90: PetscCall(DMLabelGetValue(ghostLabel, face, &ghost));
91: PetscCall(DMIsBoundaryPoint(dm, face, &boundary));
92: PetscCall(DMPlexGetTreeChildren(dm, face, &numChildren, NULL));
93: if (ghost >= 0 || boundary || numChildren) continue;
94: PetscCall(DMPlexGetSupportSize(dm, face, &numCells));
95: PetscCheck(numCells == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "facet %" PetscInt_FMT " has %" PetscInt_FMT " support points: expected 2", face, numCells);
96: PetscCall(DMPlexGetSupport(dm, face, &cells));
97: PetscCall(DMPlexPointLocalRead(dmFace, face, facegeom, &fg));
98: for (c = 0; c < 2; ++c) {
99: if (nFields > 1) PetscCall(DMPlexPointLocalFieldRead(dm, cells[c], field, x, &cx[c]));
100: else PetscCall(DMPlexPointLocalRead(dm, cells[c], x, &cx[c]));
101: PetscCall(DMPlexPointGlobalRef(dmGrad, cells[c], gr, &cgrad[c]));
102: }
103: for (pd = 0; pd < dof; ++pd) {
104: PetscScalar delta = cx[1][pd] - cx[0][pd];
106: for (d = 0; d < dim; ++d) {
107: if (cgrad[0]) cgrad[0][pd * dim + d] += fg->grad[0][d] * delta;
108: if (cgrad[1]) cgrad[1][pd * dim + d] -= fg->grad[1][d] * delta;
109: }
110: }
111: }
112: /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */
113: PetscCall(PetscLimiterGetType(lim, &limType));
114: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
115: PetscCall(DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi));
116: for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEnd; ++cell) {
117: const PetscInt *faces;
118: PetscScalar *cx;
119: PetscFVCellGeom *cg;
120: PetscScalar *cgrad;
121: PetscInt coneSize, f, pd, d;
123: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
124: PetscCall(DMPlexGetCone(dm, cell, &faces));
125: if (nFields > 1) PetscCall(DMPlexPointLocalFieldRead(dm, cell, field, x, &cx));
126: else PetscCall(DMPlexPointLocalRead(dm, cell, x, &cx));
127: PetscCall(DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg));
128: PetscCall(DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad));
129: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
130: if (limType) {
131: /* Limiter will be minimum value over all neighbors */
132: for (d = 0; d < dof; ++d) cellPhi[d] = PETSC_MAX_REAL;
133: for (f = 0; f < coneSize; ++f) PetscCall(DMPlexApplyLimiter_Internal(dm, dmCell, lim, dim, dof, cell, nFields > 1 ? field : -1, faces[f], fStart, fEnd, cellPhi, x, cellgeom, cg, cx, cgrad));
134: } else {
135: for (d = 0; d < dof; ++d) cellPhi[d] = 1.0;
136: }
137: /* Apply limiter to gradient */
138: for (pd = 0; pd < dof; ++pd) /* Scalar limiter applied to each component separately */
139: for (d = 0; d < dim; ++d) cgrad[pd * dim + d] *= cellPhi[pd];
140: }
141: PetscCall(DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi));
142: PetscCall(VecRestoreArrayRead(faceGeometry, &facegeom));
143: PetscCall(VecRestoreArrayRead(cellGeometry, &cellgeom));
144: PetscCall(VecRestoreArrayRead(locX, &x));
145: PetscCall(VecRestoreArray(grad, &gr));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.
152: Input Parameters:
153: + dm - the mesh
154: - locX - the local representation of the vector
156: Output Parameter:
157: . grad - the global representation of the gradient
159: Level: developer
161: .seealso: [](ch_unstructured), `DM`, `Vec`, `DMPlexGetGradientDM()`
162: @*/
163: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
164: {
165: PetscDS prob;
166: PetscInt Nf, f, fStart, fEnd;
167: PetscBool useFVM = PETSC_FALSE;
168: PetscFV fvm = NULL;
169: Vec faceGeometryFVM, cellGeometryFVM;
170: PetscFVCellGeom *cgeomFVM = NULL;
171: PetscFVFaceGeom *fgeomFVM = NULL;
172: DM dmGrad = NULL;
174: PetscFunctionBegin;
175: PetscCall(DMGetDS(dm, &prob));
176: PetscCall(PetscDSGetNumFields(prob, &Nf));
177: for (f = 0; f < Nf; ++f) {
178: PetscObject obj;
179: PetscClassId id;
181: PetscCall(PetscDSGetDiscretization(prob, f, &obj));
182: PetscCall(PetscObjectGetClassId(obj, &id));
183: if (id == PETSCFV_CLASSID) {
184: useFVM = PETSC_TRUE;
185: fvm = (PetscFV)obj;
186: }
187: }
188: PetscCheck(useFVM, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This dm does not have a finite volume discretization");
189: PetscCall(DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad));
190: PetscCheck(dmGrad, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This dm's finite volume discretization does not reconstruct gradients");
191: PetscCall(VecGetArrayRead(faceGeometryFVM, (const PetscScalar **)&fgeomFVM));
192: PetscCall(VecGetArrayRead(cellGeometryFVM, (const PetscScalar **)&cgeomFVM));
193: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
194: PetscCall(DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }