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,
8: PetscReal *cellPhi, const PetscScalar *x, const PetscScalar *cellgeom, const PetscFVCellGeom *cg, const PetscScalar *cx, const PetscScalar *cgrad)
9: {
10: const PetscInt *children;
11: PetscInt numChildren;
13: 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) {
21: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,dof,cell,field,childFace,fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
22: }
23: }
24: } else {
25: PetscScalar *ncx;
26: PetscFVCellGeom *ncg;
27: const PetscInt *fcells;
28: PetscInt ncell, d;
29: PetscReal v[3];
31: DMPlexGetSupport(dm, face, &fcells);
32: ncell = cell == fcells[0] ? fcells[1] : fcells[0];
33: if (field >= 0) {
34: DMPlexPointLocalFieldRead(dm, ncell, field, x, &ncx);
35: } else {
36: DMPlexPointLocalRead(dm, ncell, x, &ncx);
37: }
38: 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 phi, flim = 0.5 * PetscRealPart(ncx[d] - cx[d]) / denom;
45: PetscLimiterLimit(lim, flim, &phi);
46: cellPhi[d] = PetscMin(cellPhi[d], phi);
47: }
48: }
49: return 0;
50: }
52: PetscErrorCode DMPlexReconstructGradients_Internal(DM dm, PetscFV fvm, PetscInt fStart, PetscInt fEnd, Vec faceGeometry, Vec cellGeometry, Vec locX, Vec grad)
53: {
54: DM dmFace, dmCell, dmGrad;
55: DMLabel ghostLabel;
56: PetscDS prob;
57: PetscLimiter lim;
58: const PetscScalar *facegeom, *cellgeom, *x;
59: PetscScalar *gr;
60: PetscReal *cellPhi;
61: PetscInt dim, face, cell, field, dof, cStart, cEnd, nFields;
63: DMGetDimension(dm, &dim);
64: DMGetDS(dm, &prob);
65: PetscDSGetNumFields(prob, &nFields);
66: PetscDSGetFieldIndex(prob, (PetscObject) fvm, &field);
67: PetscDSGetFieldSize(prob, field, &dof);
68: DMGetLabel(dm, "ghost", &ghostLabel);
69: PetscFVGetLimiter(fvm, &lim);
70: VecGetDM(faceGeometry, &dmFace);
71: VecGetArrayRead(faceGeometry, &facegeom);
72: VecGetDM(cellGeometry, &dmCell);
73: VecGetArrayRead(cellGeometry, &cellgeom);
74: VecGetArrayRead(locX, &x);
75: VecGetDM(grad, &dmGrad);
76: VecZeroEntries(grad);
77: VecGetArray(grad, &gr);
78: /* Reconstruct gradients */
79: for (face = fStart; face < fEnd; ++face) {
80: const PetscInt *cells;
81: PetscFVFaceGeom *fg;
82: PetscScalar *cx[2];
83: PetscScalar *cgrad[2];
84: PetscBool boundary;
85: PetscInt ghost, c, pd, d, numChildren, numCells;
87: DMLabelGetValue(ghostLabel, face, &ghost);
88: DMIsBoundaryPoint(dm, face, &boundary);
89: DMPlexGetTreeChildren(dm, face, &numChildren, NULL);
90: if (ghost >= 0 || boundary || numChildren) continue;
91: DMPlexGetSupportSize(dm, face, &numCells);
93: DMPlexGetSupport(dm, face, &cells);
94: DMPlexPointLocalRead(dmFace, face, facegeom, &fg);
95: for (c = 0; c < 2; ++c) {
96: if (nFields > 1) {
97: DMPlexPointLocalFieldRead(dm, cells[c], field, x, &cx[c]);
98: } else {
99: DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);
100: }
101: 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: DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
114: DMGetWorkArray(dm, dof, MPIU_REAL, &cellPhi);
115: for (cell = (dmGrad && lim) ? cStart : cEnd; cell < cEnd; ++cell) {
116: const PetscInt *faces;
117: PetscScalar *cx;
118: PetscFVCellGeom *cg;
119: PetscScalar *cgrad;
120: PetscInt coneSize, f, pd, d;
122: DMPlexGetConeSize(dm, cell, &coneSize);
123: DMPlexGetCone(dm, cell, &faces);
124: if (nFields > 1) {
125: DMPlexPointLocalFieldRead(dm, cell, field, x, &cx);
126: }
127: else {
128: DMPlexPointLocalRead(dm, cell, x, &cx);
129: }
130: DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);
131: DMPlexPointGlobalRef(dmGrad, cell, gr, &cgrad);
132: if (!cgrad) continue; /* Unowned overlap cell, we do not compute */
133: /* Limiter will be minimum value over all neighbors */
134: for (d = 0; d < dof; ++d) cellPhi[d] = PETSC_MAX_REAL;
135: for (f = 0; f < coneSize; ++f) {
136: DMPlexApplyLimiter_Internal(dm,dmCell,lim,dim,dof,cell,nFields > 1 ? field : -1,faces[f],fStart,fEnd,cellPhi,x,cellgeom,cg,cx,cgrad);
137: }
138: /* Apply limiter to gradient */
139: for (pd = 0; pd < dof; ++pd)
140: /* Scalar limiter applied to each component separately */
141: for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd];
142: }
143: DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi);
144: VecRestoreArrayRead(faceGeometry, &facegeom);
145: VecRestoreArrayRead(cellGeometry, &cellgeom);
146: VecRestoreArrayRead(locX, &x);
147: VecRestoreArray(grad, &gr);
148: return 0;
149: }
151: /*@
152: DMPlexReconstructGradientsFVM - reconstruct the gradient of a vector using a finite volume method.
154: Input Parameters:
155: + dm - the mesh
156: - locX - the local representation of the vector
158: Output Parameter:
159: . grad - the global representation of the gradient
161: Level: developer
163: .seealso: DMPlexGetGradientDM()
164: @*/
165: PetscErrorCode DMPlexReconstructGradientsFVM(DM dm, Vec locX, Vec grad)
166: {
167: PetscDS prob;
168: PetscInt Nf, f, fStart, fEnd;
169: PetscBool useFVM = PETSC_FALSE;
170: PetscFV fvm = NULL;
171: Vec faceGeometryFVM, cellGeometryFVM;
172: PetscFVCellGeom *cgeomFVM = NULL;
173: PetscFVFaceGeom *fgeomFVM = NULL;
174: DM dmGrad = NULL;
176: DMGetDS(dm, &prob);
177: PetscDSGetNumFields(prob, &Nf);
178: for (f = 0; f < Nf; ++f) {
179: PetscObject obj;
180: PetscClassId id;
182: PetscDSGetDiscretization(prob, f, &obj);
183: PetscObjectGetClassId(obj, &id);
184: if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;}
185: }
187: DMPlexGetDataFVM(dm, fvm, &cellGeometryFVM, &faceGeometryFVM, &dmGrad);
189: VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);
190: VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);
191: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
192: DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, locX, grad);
193: return 0;
194: }