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