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