1 #include <petscdmplex.h> /*I "petscdmplex.h" I*/ 2 #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/ 3 #include <petscfv.h> 4 5 /* TODO Move LS stuff to dtfv.c */ 6 #include <petscblaslapack.h> 7 8 PETSC_STATIC_INLINE void WaxpyD(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w) {PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];} 9 PETSC_STATIC_INLINE PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y) {PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;} 10 PETSC_STATIC_INLINE PetscReal NormD(PetscInt dim, const PetscScalar *x) {return PetscSqrtReal(PetscAbsScalar(DotD(dim,x,x)));} 11 12 typedef struct { 13 PetscBool setupGeom; /* Flag for geometry setup */ 14 PetscBool setupGrad; /* Flag for gradient calculation setup */ 15 Vec facegeom; /* FaceGeom struct for each face */ 16 Vec cellgeom; /* CellGeom struct for each cell */ 17 DM dmGrad; /* Layout for the gradient data */ 18 PetscReal minradius; /* Minimum distance from centroid to face */ 19 void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *); 20 void *rhsfunctionlocalctx; 21 } DMTS_Plex; 22 23 #undef __FUNCT__ 24 #define __FUNCT__ "DMTSDestroy_Plex" 25 static PetscErrorCode DMTSDestroy_Plex(DMTS dmts) 26 { 27 DMTS_Plex *dmplexts = (DMTS_Plex *) dmts->data; 28 PetscErrorCode ierr; 29 30 PetscFunctionBegin; 31 ierr = DMDestroy(&dmplexts->dmGrad);CHKERRQ(ierr); 32 ierr = VecDestroy(&dmplexts->facegeom);CHKERRQ(ierr); 33 ierr = VecDestroy(&dmplexts->cellgeom);CHKERRQ(ierr); 34 ierr = PetscFree(dmts->data);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "DMTSDuplicate_Plex" 40 static PetscErrorCode DMTSDuplicate_Plex(DMTS olddmts, DMTS dmts) 41 { 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = PetscNewLog(dmts, (DMTS_Plex **) &dmts->data);CHKERRQ(ierr); 46 if (olddmts->data) {ierr = PetscMemcpy(dmts->data, olddmts->data, sizeof(DMTS_Plex));CHKERRQ(ierr);} 47 PetscFunctionReturn(0); 48 } 49 50 #undef __FUNCT__ 51 #define __FUNCT__ "DMPlexTSGetContext" 52 static PetscErrorCode DMPlexTSGetContext(DM dm, DMTS dmts, DMTS_Plex **dmplexts) 53 { 54 PetscErrorCode ierr; 55 56 PetscFunctionBegin; 57 *dmplexts = NULL; 58 if (!dmts->data) { 59 ierr = PetscNewLog(dm, (DMTS_Plex **) &dmts->data);CHKERRQ(ierr); 60 dmts->ops->destroy = DMTSDestroy_Plex; 61 dmts->ops->duplicate = DMTSDuplicate_Plex; 62 } 63 *dmplexts = (DMTS_Plex *) dmts->data; 64 PetscFunctionReturn(0); 65 } 66 67 #undef __FUNCT__ 68 #define __FUNCT__ "DMPlexTSSetupGeometry" 69 static PetscErrorCode DMPlexTSSetupGeometry(DM dm, PetscFV fvm, DMTS_Plex *dmplexts) 70 { 71 DM dmFace, dmCell; 72 DMLabel ghostLabel; 73 PetscSection sectionFace, sectionCell; 74 PetscSection coordSection; 75 Vec coordinates; 76 PetscReal minradius; 77 PetscScalar *fgeom, *cgeom; 78 PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f; 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 if (dmplexts->setupGeom) PetscFunctionReturn(0); 83 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 84 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 85 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 86 /* Make cell centroids and volumes */ 87 ierr = DMClone(dm, &dmCell);CHKERRQ(ierr); 88 ierr = DMSetCoordinateSection(dmCell, coordSection);CHKERRQ(ierr); 89 ierr = DMSetCoordinatesLocal(dmCell, coordinates);CHKERRQ(ierr); 90 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionCell);CHKERRQ(ierr); 91 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 92 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 93 ierr = PetscSectionSetChart(sectionCell, cStart, cEnd);CHKERRQ(ierr); 94 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionCell, c, sizeof(CellGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 95 ierr = PetscSectionSetUp(sectionCell);CHKERRQ(ierr); 96 ierr = DMSetDefaultSection(dmCell, sectionCell);CHKERRQ(ierr); 97 ierr = PetscSectionDestroy(§ionCell);CHKERRQ(ierr); 98 ierr = DMCreateLocalVector(dmCell, &dmplexts->cellgeom);CHKERRQ(ierr); 99 ierr = VecGetArray(dmplexts->cellgeom, &cgeom);CHKERRQ(ierr); 100 for (c = cStart; c < cEndInterior; ++c) { 101 CellGeom *cg; 102 103 ierr = DMPlexPointLocalRef(dmCell, c, cgeom, &cg);CHKERRQ(ierr); 104 ierr = PetscMemzero(cg, sizeof(*cg));CHKERRQ(ierr); 105 ierr = DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);CHKERRQ(ierr); 106 } 107 /* Compute face normals and minimum cell radius */ 108 ierr = DMClone(dm, &dmFace);CHKERRQ(ierr); 109 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionFace);CHKERRQ(ierr); 110 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 111 ierr = PetscSectionSetChart(sectionFace, fStart, fEnd);CHKERRQ(ierr); 112 for (f = fStart; f < fEnd; ++f) {ierr = PetscSectionSetDof(sectionFace, f, sizeof(FaceGeom)/sizeof(PetscScalar));CHKERRQ(ierr);} 113 ierr = PetscSectionSetUp(sectionFace);CHKERRQ(ierr); 114 ierr = DMSetDefaultSection(dmFace, sectionFace);CHKERRQ(ierr); 115 ierr = PetscSectionDestroy(§ionFace);CHKERRQ(ierr); 116 ierr = DMCreateLocalVector(dmFace, &dmplexts->facegeom);CHKERRQ(ierr); 117 ierr = VecGetArray(dmplexts->facegeom, &fgeom);CHKERRQ(ierr); 118 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 119 minradius = PETSC_MAX_REAL; 120 for (f = fStart; f < fEnd; ++f) { 121 FaceGeom *fg; 122 PetscReal area; 123 PetscInt ghost, d; 124 125 ierr = DMLabelGetValue(ghostLabel, f, &ghost);CHKERRQ(ierr); 126 if (ghost >= 0) continue; 127 ierr = DMPlexPointLocalRef(dmFace, f, fgeom, &fg);CHKERRQ(ierr); 128 ierr = DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);CHKERRQ(ierr); 129 for (d = 0; d < dim; ++d) fg->normal[d] *= area; 130 /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */ 131 { 132 CellGeom *cL, *cR; 133 const PetscInt *cells; 134 PetscReal *lcentroid, *rcentroid; 135 PetscScalar v[3]; 136 137 ierr = DMPlexGetSupport(dm, f, &cells);CHKERRQ(ierr); 138 ierr = DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);CHKERRQ(ierr); 139 ierr = DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);CHKERRQ(ierr); 140 lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid; 141 rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid; 142 WaxpyD(dim, -1, lcentroid, rcentroid, v); 143 if (DotD(dim, fg->normal, v) < 0) { 144 for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d]; 145 } 146 if (DotD(dim, fg->normal, v) <= 0) { 147 if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, fg->normal[0], fg->normal[1], v[0], v[1]); 148 if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, fg->normal[0], fg->normal[1], fg->normal[2], v[0], v[1], v[2]); 149 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f); 150 } 151 if (cells[0] < cEndInterior) { 152 WaxpyD(dim, -1, fg->centroid, cL->centroid, v); 153 minradius = PetscMin(minradius, NormD(dim, v)); 154 } 155 if (cells[1] < cEndInterior) { 156 WaxpyD(dim, -1, fg->centroid, cR->centroid, v); 157 minradius = PetscMin(minradius, NormD(dim, v)); 158 } 159 } 160 } 161 ierr = MPI_Allreduce(&minradius, &dmplexts->minradius, 1, MPIU_SCALAR, MPI_MIN, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 162 /* Compute centroids of ghost cells */ 163 for (c = cEndInterior; c < cEnd; ++c) { 164 FaceGeom *fg; 165 const PetscInt *cone, *support; 166 PetscInt coneSize, supportSize, s; 167 168 ierr = DMPlexGetConeSize(dmCell, c, &coneSize);CHKERRQ(ierr); 169 if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize); 170 ierr = DMPlexGetCone(dmCell, c, &cone);CHKERRQ(ierr); 171 ierr = DMPlexGetSupportSize(dmCell, cone[0], &supportSize);CHKERRQ(ierr); 172 if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 1", cone[0], supportSize); 173 ierr = DMPlexGetSupport(dmCell, cone[0], &support);CHKERRQ(ierr); 174 ierr = DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);CHKERRQ(ierr); 175 for (s = 0; s < 2; ++s) { 176 /* Reflect ghost centroid across plane of face */ 177 if (support[s] == c) { 178 const CellGeom *ci; 179 CellGeom *cg; 180 PetscScalar c2f[3], a; 181 182 ierr = DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);CHKERRQ(ierr); 183 WaxpyD(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */ 184 a = DotD(dim, c2f, fg->normal)/DotD(dim, fg->normal, fg->normal); 185 ierr = DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);CHKERRQ(ierr); 186 WaxpyD(dim, 2*a, fg->normal, ci->centroid, cg->centroid); 187 cg->volume = ci->volume; 188 } 189 } 190 } 191 ierr = VecRestoreArray(dmplexts->facegeom, &fgeom);CHKERRQ(ierr); 192 ierr = VecRestoreArray(dmplexts->cellgeom, &cgeom);CHKERRQ(ierr); 193 ierr = DMDestroy(&dmCell);CHKERRQ(ierr); 194 ierr = DMDestroy(&dmFace);CHKERRQ(ierr); 195 dmplexts->setupGeom = PETSC_TRUE; 196 PetscFunctionReturn(0); 197 } 198 199 #undef __FUNCT__ 200 #define __FUNCT__ "PseudoInverse" 201 /* Overwrites A. Can only handle full-rank problems with m>=n */ 202 static PetscErrorCode PseudoInverse(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 203 { 204 PetscBool debug = PETSC_FALSE; 205 PetscErrorCode ierr; 206 PetscBLASInt M,N,K,lda,ldb,ldwork,info; 207 PetscScalar *R,*Q,*Aback,Alpha; 208 209 PetscFunctionBegin; 210 if (debug) { 211 ierr = PetscMalloc1(m*n,&Aback);CHKERRQ(ierr); 212 ierr = PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 213 } 214 215 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 216 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 217 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 218 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 219 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 220 LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info); 221 ierr = PetscFPTrapPop();CHKERRQ(ierr); 222 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error"); 223 R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */ 224 225 /* Extract an explicit representation of Q */ 226 Q = Ainv; 227 ierr = PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));CHKERRQ(ierr); 228 K = N; /* full rank */ 229 LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info); 230 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error"); 231 232 /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */ 233 Alpha = 1.0; 234 ldb = lda; 235 BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb); 236 /* Ainv is Q, overwritten with inverse */ 237 238 if (debug) { /* Check that pseudo-inverse worked */ 239 PetscScalar Beta = 0.0; 240 PetscInt ldc; 241 K = N; 242 ldc = N; 243 BLASgemm_("ConjugateTranspose","Normal",&N,&K,&M,&Alpha,Ainv,&lda,Aback,&ldb,&Beta,work,&ldc); 244 ierr = PetscScalarView(n*n,work,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 245 ierr = PetscFree(Aback);CHKERRQ(ierr); 246 } 247 PetscFunctionReturn(0); 248 } 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "PseudoInverseGetWorkRequired" 252 static PetscErrorCode PseudoInverseGetWorkRequired(PetscInt maxFaces, PetscInt *work) 253 { 254 PetscInt m,n,nrhs,minwork; 255 256 PetscFunctionBegin; 257 m = maxFaces; 258 n = 3; /* spatial dimension */ 259 nrhs = maxFaces; 260 minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */ 261 *work = 5*minwork; /* We can afford to be extra generous */ 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNCT__ 266 #define __FUNCT__ "PseudoInverseSVD" 267 /* Overwrites A. Can handle degenerate problems and m<n. */ 268 static PetscErrorCode PseudoInverseSVD(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work) 269 { 270 PetscBool debug = PETSC_FALSE; 271 PetscErrorCode ierr; 272 PetscInt i,j,maxmn; 273 PetscBLASInt M,N,nrhs,lda,ldb,irank,ldwork,info; 274 PetscScalar rcond,*tmpwork,*Brhs,*Aback; 275 276 PetscFunctionBegin; 277 if (debug) { 278 ierr = PetscMalloc1(m*n,&Aback);CHKERRQ(ierr); 279 ierr = PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 280 } 281 282 /* initialize to identity */ 283 tmpwork = Ainv; 284 Brhs = work; 285 maxmn = PetscMax(m,n); 286 for (j=0; j<maxmn; j++) { 287 for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j); 288 } 289 290 ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr); 291 ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr); 292 nrhs = M; 293 ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr); 294 ierr = PetscBLASIntCast(maxmn,&ldb);CHKERRQ(ierr); 295 ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr); 296 rcond = -1; 297 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 298 LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb,tau,&rcond,&irank,tmpwork,&ldwork,&info); 299 ierr = PetscFPTrapPop();CHKERRQ(ierr); 300 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGELSS error"); 301 /* The following check should be turned into a diagnostic as soon as someone wants to do this intentionally */ 302 if (irank < PetscMin(M,N)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Rank deficient least squares fit, indicates an isolated cell with two colinear points"); 303 304 /* Brhs shaped (M,nrhs) column-major coldim=mstride was overwritten by Ainv shaped (N,nrhs) column-major coldim=maxmn. 305 * Here we transpose to (N,nrhs) row-major rowdim=mstride. */ 306 for (i=0; i<n; i++) { 307 for (j=0; j<nrhs; j++) Ainv[i*mstride+j] = Brhs[i + j*maxmn]; 308 } 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "BuildLeastSquares" 314 /* Build least squares gradient reconstruction operators */ 315 static PetscErrorCode BuildLeastSquares(DM dm,PetscInt cEndInterior,DM dmFace,PetscScalar *fgeom,DM dmCell,PetscScalar *cgeom) 316 { 317 DMLabel ghostLabel, bdLabel; 318 PetscScalar *B,*Binv,*work,*tau,**gref; 319 PetscInt dim, c,cStart,cEnd,maxNumFaces,worksize; 320 PetscBool useSVD = PETSC_TRUE; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 325 ierr = DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 326 ierr = DMPlexGetMaxSizes(dm,&maxNumFaces,NULL);CHKERRQ(ierr); 327 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 328 /* TODO: Get this from the BC */ 329 ierr = DMPlexGetLabel(dm, "Face Sets", &bdLabel);CHKERRQ(ierr); 330 ierr = PseudoInverseGetWorkRequired(maxNumFaces,&worksize);CHKERRQ(ierr); 331 ierr = PetscMalloc5(maxNumFaces*dim,&B,worksize,&Binv,worksize,&work,maxNumFaces,&tau,maxNumFaces,&gref);CHKERRQ(ierr); 332 for (c=cStart; c<cEndInterior; c++) { 333 const PetscInt *faces; 334 PetscInt numFaces,usedFaces,f,i,j; 335 const CellGeom *cg; 336 PetscInt ghost, boundary; 337 ierr = DMPlexGetConeSize(dm,c,&numFaces);CHKERRQ(ierr); 338 if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction",c,numFaces); 339 ierr = DMPlexGetCone(dm,c,&faces);CHKERRQ(ierr); 340 ierr = DMPlexPointLocalRead(dmCell,c,cgeom,&cg);CHKERRQ(ierr); 341 for (f=0,usedFaces=0; f<numFaces; f++) { 342 const PetscInt *fcells; 343 PetscInt ncell,side; 344 FaceGeom *fg; 345 const CellGeom *cg1; 346 347 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 348 ierr = DMLabelGetValue(bdLabel, faces[f], &boundary);CHKERRQ(ierr); 349 if ((ghost >= 0) || (boundary >= 0)) continue; 350 ierr = DMPlexGetSupport(dm,faces[f],&fcells);CHKERRQ(ierr); 351 side = (c != fcells[0]); /* c is on left=0 or right=1 of face */ 352 ncell = fcells[!side]; /* the neighbor */ 353 ierr = DMPlexPointLocalRef(dmFace,faces[f],fgeom,&fg);CHKERRQ(ierr); 354 ierr = DMPlexPointLocalRead(dmCell,ncell,cgeom,&cg1);CHKERRQ(ierr); 355 for (j=0; j<dim; j++) B[j*numFaces+usedFaces] = cg1->centroid[j] - cg->centroid[j]; 356 gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */ 357 } 358 if (!usedFaces) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Mesh contains isolated cell (no neighbors). Is it intentional?"); 359 /* Overwrites B with garbage, returns Binv in row-major format */ 360 if (useSVD) { 361 ierr = PseudoInverseSVD(usedFaces,numFaces,dim,B,Binv,tau,worksize,work);CHKERRQ(ierr); 362 } else { 363 ierr = PseudoInverse(usedFaces,numFaces,dim,B,Binv,tau,worksize,work);CHKERRQ(ierr); 364 } 365 for (f=0,i=0; f<numFaces; f++) { 366 ierr = DMLabelGetValue(ghostLabel, faces[f], &ghost);CHKERRQ(ierr); 367 ierr = DMLabelGetValue(bdLabel, faces[f], &boundary);CHKERRQ(ierr); 368 if ((ghost >= 0) || (boundary >= 0)) continue; 369 for (j=0; j<dim; j++) gref[i][j] = Binv[j*numFaces+i]; 370 i++; 371 } 372 373 if (0) { 374 PetscReal grad[2] = {0,0}; 375 for (f=0; f<numFaces; f++) { 376 const PetscInt *fcells; 377 const CellGeom *cg1; 378 const FaceGeom *fg; 379 ierr = DMPlexGetSupport(dm,faces[f],&fcells);CHKERRQ(ierr); 380 ierr = DMPlexPointLocalRead(dmFace,faces[f],fgeom,&fg);CHKERRQ(ierr); 381 for (i=0; i<2; i++) { 382 if (fcells[i] == c) continue; 383 ierr = DMPlexPointLocalRead(dmCell,fcells[i],cgeom,&cg1);CHKERRQ(ierr); 384 PetscScalar du = cg1->centroid[0] + 3*cg1->centroid[1] - (cg->centroid[0] + 3*cg->centroid[1]); 385 grad[0] += fg->grad[!i][0] * du; 386 grad[1] += fg->grad[!i][1] * du; 387 } 388 } 389 printf("cell[%d] grad (%g,%g)\n",c,grad[0],grad[1]); 390 } 391 } 392 ierr = PetscFree5(B,Binv,work,tau,gref);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 #undef __FUNCT__ 397 #define __FUNCT__ "DMPlexTSSetupGradient" 398 static PetscErrorCode DMPlexTSSetupGradient(DM dm, PetscFV fvm, DMTS_Plex *dmplexts) 399 { 400 DM dmFace, dmCell; 401 PetscScalar *fgeom, *cgeom; 402 PetscSection sectionGrad; 403 PetscInt dim, pdim, cStart, cEnd, cEndInterior, c; 404 PetscErrorCode ierr; 405 406 PetscFunctionBegin; 407 if (dmplexts->setupGrad) PetscFunctionReturn(0); 408 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 409 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 410 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 411 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 412 /* Construct the interpolant corresponding to each face from the leat-square solution over the cell neighborhood */ 413 ierr = VecGetDM(dmplexts->facegeom, &dmFace);CHKERRQ(ierr); 414 ierr = VecGetDM(dmplexts->cellgeom, &dmCell);CHKERRQ(ierr); 415 ierr = VecGetArray(dmplexts->facegeom, &fgeom);CHKERRQ(ierr); 416 ierr = VecGetArray(dmplexts->cellgeom, &cgeom);CHKERRQ(ierr); 417 ierr = BuildLeastSquares(dm, cEndInterior, dmFace, fgeom, dmCell, cgeom);CHKERRQ(ierr); 418 ierr = VecRestoreArray(dmplexts->facegeom, &fgeom);CHKERRQ(ierr); 419 ierr = VecRestoreArray(dmplexts->cellgeom, &cgeom);CHKERRQ(ierr); 420 /* Create storage for gradients */ 421 ierr = DMClone(dm, &dmplexts->dmGrad);CHKERRQ(ierr); 422 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ionGrad);CHKERRQ(ierr); 423 ierr = PetscSectionSetChart(sectionGrad, cStart, cEnd);CHKERRQ(ierr); 424 for (c = cStart; c < cEnd; ++c) {ierr = PetscSectionSetDof(sectionGrad, c, pdim*dim);CHKERRQ(ierr);} 425 ierr = PetscSectionSetUp(sectionGrad);CHKERRQ(ierr); 426 ierr = DMSetDefaultSection(dmplexts->dmGrad, sectionGrad);CHKERRQ(ierr); 427 ierr = PetscSectionDestroy(§ionGrad);CHKERRQ(ierr); 428 dmplexts->setupGrad = PETSC_TRUE; 429 PetscFunctionReturn(0); 430 } 431 432 #undef __FUNCT__ 433 #define __FUNCT__ "DMPlexInsertBoundaryValuesFVM_Static" 434 static PetscErrorCode DMPlexInsertBoundaryValuesFVM_Static(DM dm, PetscFV fvm, PetscReal time, Vec locX, Vec Grad, DMTS_Plex *dmplexts) 435 { 436 Vec faceGeometry = dmplexts->facegeom; 437 Vec cellGeometry = dmplexts->cellgeom; 438 DM dmFace, dmCell; 439 DMLabel label; 440 const PetscScalar *facegeom, *cellgeom, *grad; 441 PetscScalar *x, *fx; 442 PetscInt numBd, b, dim, pdim, fStart, fEnd; 443 PetscErrorCode ierr; 444 445 PetscFunctionBegin; 446 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 447 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 448 /* TODO Get from BC data */ 449 ierr = DMPlexGetLabel(dm, "Face Sets", &label);CHKERRQ(ierr); 450 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 451 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 452 if (Grad) { 453 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 454 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 455 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 456 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 457 } 458 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 459 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 460 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 461 for (b = 0; b < numBd; ++b) { 462 PetscErrorCode (*func)(PetscReal,const PetscScalar*,const PetscScalar*,const PetscScalar*,PetscScalar*,void*); 463 const PetscInt *ids; 464 PetscInt numids, i; 465 void *ctx; 466 467 ierr = DMPlexGetBoundary(dm, b, NULL, NULL, NULL, (void (**)()) &func, &numids, &ids, &ctx);CHKERRQ(ierr); 468 for (i = 0; i < numids; ++i) { 469 IS faceIS; 470 const PetscInt *faces; 471 PetscInt numFaces, f; 472 473 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 474 if (!faceIS) continue; /* No points with that id on this process */ 475 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 476 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 477 for (f = 0; f < numFaces; ++f) { 478 const PetscInt face = faces[f], *cells; 479 const FaceGeom *fg; 480 481 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 482 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 483 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 484 if (Grad) { 485 const CellGeom *cg; 486 const PetscScalar *cx, *cgrad; 487 PetscScalar *xG, dx[3]; 488 PetscInt d; 489 490 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 491 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 492 ierr = DMPlexPointLocalRead(dmplexts->dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 493 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 494 WaxpyD(dim, -1, cg->centroid, fg->centroid, dx); 495 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DotD(dim, &cgrad[d*dim], dx); 496 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 497 } else { 498 const PetscScalar *xI; 499 PetscScalar *xG; 500 501 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 502 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 503 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 504 } 505 } 506 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 507 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 508 } 509 } 510 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 511 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 512 if (Grad) { 513 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 514 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 515 } 516 PetscFunctionReturn(0); 517 } 518 519 #undef __FUNCT__ 520 #define __FUNCT__ "TSComputeRHSFunction_DMPlex" 521 static PetscErrorCode TSComputeRHSFunction_DMPlex(TS ts, PetscReal time, Vec X, Vec F, void *ctx) 522 { 523 DM dm; 524 DMTS_Plex *dmplexts = (DMTS_Plex *) ctx; 525 void (*riemann)(const PetscReal[], const PetscReal[], const PetscScalar[], const PetscScalar[], PetscScalar[], void *) = dmplexts->riemann; 526 PetscFV fvm; 527 PetscLimiter lim; 528 Vec faceGeometry = dmplexts->facegeom; 529 Vec cellGeometry = dmplexts->cellgeom; 530 Vec Grad = NULL, locGrad, locX; 531 DM dmFace, dmCell; 532 DMLabel ghostLabel, bdLabel; 533 PetscCellGeometry fgeom, cgeom; 534 const PetscScalar *facegeom, *cellgeom, *x, *lgrad; 535 PetscScalar *grad, *f, *uL, *uR, *fluxL, *fluxR; 536 PetscReal *centroid, *normal, *vol, *cellPhi; 537 PetscBool computeGradients; 538 PetscInt Nf, dim, pdim, fStart, fEnd, numFaces = 0, face, iface, cell, cStart, cEnd, cEndInterior; 539 PetscErrorCode ierr; 540 541 PetscFunctionBegin; 542 PetscValidHeaderSpecific(ts,TS_CLASSID,1); 543 PetscValidHeaderSpecific(X,VEC_CLASSID,3); 544 PetscValidHeaderSpecific(F,VEC_CLASSID,5); 545 ierr = TSGetDM(ts, &dm);CHKERRQ(ierr); 546 ierr = DMGetLocalVector(dm, &locX);CHKERRQ(ierr); 547 ierr = VecZeroEntries(locX);CHKERRQ(ierr); 548 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 549 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);CHKERRQ(ierr); 550 ierr = VecZeroEntries(F);CHKERRQ(ierr); 551 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 552 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 553 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 554 ierr = PetscFVGetLimiter(fvm, &lim);CHKERRQ(ierr); 555 ierr = PetscFVGetNumComponents(fvm, &pdim);CHKERRQ(ierr); 556 ierr = PetscFVGetComputeGradients(fvm, &computeGradients);CHKERRQ(ierr); 557 if (computeGradients) { 558 ierr = DMGetGlobalVector(dmplexts->dmGrad, &Grad);CHKERRQ(ierr); 559 ierr = VecZeroEntries(Grad);CHKERRQ(ierr); 560 ierr = VecGetArray(Grad, &grad);CHKERRQ(ierr); 561 } 562 ierr = DMPlexGetLabel(dm, "ghost", &ghostLabel);CHKERRQ(ierr); 563 /* TODO: Get this from the BC */ 564 ierr = DMPlexGetLabel(dm, "Face Sets", &bdLabel);CHKERRQ(ierr); 565 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 566 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 567 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 568 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 569 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 570 ierr = VecGetArrayRead(locX, &x);CHKERRQ(ierr); 571 /* Count faces and reconstruct gradients */ 572 for (face = fStart; face < fEnd; ++face) { 573 const PetscInt *cells; 574 const FaceGeom *fg; 575 const PetscScalar *cx[2]; 576 PetscScalar *cgrad[2]; 577 PetscInt ghost, boundary, c, pd, d; 578 579 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 580 if (ghost >= 0) continue; 581 ++numFaces; 582 if (!computeGradients) continue; 583 ierr = DMLabelGetValue(bdLabel, face, &boundary);CHKERRQ(ierr); 584 if (boundary >= 0) continue; 585 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 586 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 587 for (c = 0; c < 2; ++c) { 588 ierr = DMPlexPointLocalRead(dm, cells[c], x, &cx[c]);CHKERRQ(ierr); 589 ierr = DMPlexPointGlobalRef(dmplexts->dmGrad, cells[c], grad, &cgrad[c]);CHKERRQ(ierr); 590 } 591 for (pd = 0; pd < pdim; ++pd) { 592 PetscScalar delta = cx[1][pd] - cx[0][pd]; 593 594 for (d = 0; d < dim; ++d) { 595 if (cgrad[0]) cgrad[0][pd*dim+d] += fg->grad[0][d] * delta; 596 if (cgrad[1]) cgrad[1][pd*dim+d] -= fg->grad[1][d] * delta; 597 } 598 } 599 } 600 /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */ 601 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 602 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 603 ierr = DMGetWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 604 for (cell = computeGradients && lim ? cStart : cEnd; cell < cEndInterior; ++cell) { 605 const PetscInt *faces; 606 const PetscScalar *cx; 607 const CellGeom *cg; 608 PetscScalar *cgrad; 609 PetscInt coneSize, f, pd, d; 610 611 ierr = DMPlexGetConeSize(dm, cell, &coneSize);CHKERRQ(ierr); 612 ierr = DMPlexGetCone(dm, cell, &faces);CHKERRQ(ierr); 613 ierr = DMPlexPointLocalRead(dm, cell, x, &cx);CHKERRQ(ierr); 614 ierr = DMPlexPointLocalRead(dmCell, cell, cellgeom, &cg);CHKERRQ(ierr); 615 ierr = DMPlexPointGlobalRef(dmplexts->dmGrad, cell, grad, &cgrad);CHKERRQ(ierr); 616 if (!cgrad) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Supposedly ghost cell %d, but this should be impossible", cell); 617 /* Limiter will be minimum value over all neighbors */ 618 for (d = 0; d < pdim; ++d) cellPhi[d] = PETSC_MAX_REAL; 619 for (f = 0; f < coneSize; ++f) { 620 const PetscScalar *ncx; 621 const CellGeom *ncg; 622 const PetscInt *fcells; 623 PetscInt face = faces[f], ncell; 624 PetscScalar v[3]; 625 PetscInt ghost, boundary; 626 627 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 628 ierr = DMLabelGetValue(bdLabel, face, &boundary);CHKERRQ(ierr); 629 if ((ghost >= 0) || (boundary >= 0)) continue; 630 ierr = DMPlexGetSupport(dm, face, &fcells);CHKERRQ(ierr); 631 ncell = cell == fcells[0] ? fcells[1] : fcells[0]; 632 ierr = DMPlexPointLocalRead(dm, ncell, x, &ncx);CHKERRQ(ierr); 633 ierr = DMPlexPointLocalRead(dmCell, ncell, cellgeom, &ncg);CHKERRQ(ierr); 634 WaxpyD(dim, -1, cg->centroid, ncg->centroid, v); 635 for (d = 0; d < pdim; ++d) { 636 /* We use the symmetric slope limited form of Berger, Aftosmis, and Murman 2005 */ 637 PetscScalar phi, flim = 0.5 * (ncx[d] - cx[d]) / DotD(dim, &cgrad[d*dim], v); 638 639 ierr = PetscLimiterLimit(lim, flim, &phi);CHKERRQ(ierr); 640 cellPhi[d] = PetscMin(cellPhi[d], phi); 641 } 642 } 643 /* Apply limiter to gradient */ 644 for (pd = 0; pd < pdim; ++pd) 645 /* Scalar limiter applied to each component separately */ 646 for (d = 0; d < dim; ++d) cgrad[pd*dim+d] *= cellPhi[pd]; 647 } 648 ierr = DMRestoreWorkArray(dm, pdim, PETSC_REAL, &cellPhi);CHKERRQ(ierr); 649 ierr = DMPlexInsertBoundaryValuesFVM_Static(dm, fvm, time, locX, Grad, dmplexts);CHKERRQ(ierr); 650 if (computeGradients) { 651 ierr = VecRestoreArray(Grad, &grad);CHKERRQ(ierr); 652 ierr = DMGetLocalVector(dmplexts->dmGrad, &locGrad);CHKERRQ(ierr); 653 ierr = DMGlobalToLocalBegin(dmplexts->dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 654 ierr = DMGlobalToLocalEnd(dmplexts->dmGrad, Grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 655 ierr = DMRestoreGlobalVector(dmplexts->dmGrad, &Grad);CHKERRQ(ierr); 656 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 657 } 658 ierr = PetscMalloc7(numFaces*dim,¢roid,numFaces*dim,&normal,numFaces*2,&vol,numFaces*pdim,&uL,numFaces*pdim,&uR,numFaces*pdim,&fluxL,numFaces*pdim,&fluxR);CHKERRQ(ierr); 659 /* Read out values */ 660 for (face = fStart, iface = 0; face < fEnd; ++face) { 661 const PetscInt *cells; 662 const FaceGeom *fg; 663 const CellGeom *cgL, *cgR; 664 const PetscScalar *xL, *xR, *gL, *gR; 665 PetscInt ghost, d; 666 667 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 668 if (ghost >= 0) continue; 669 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 670 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 671 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cgL);CHKERRQ(ierr); 672 ierr = DMPlexPointLocalRead(dmCell, cells[1], cellgeom, &cgR);CHKERRQ(ierr); 673 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xL);CHKERRQ(ierr); 674 ierr = DMPlexPointLocalRead(dm, cells[1], x, &xR);CHKERRQ(ierr); 675 if (computeGradients) { 676 PetscScalar dxL[3], dxR[3]; 677 678 ierr = DMPlexPointLocalRead(dmplexts->dmGrad, cells[0], lgrad, &gL);CHKERRQ(ierr); 679 ierr = DMPlexPointLocalRead(dmplexts->dmGrad, cells[1], lgrad, &gR);CHKERRQ(ierr); 680 WaxpyD(dim, -1, cgL->centroid, fg->centroid, dxL); 681 WaxpyD(dim, -1, cgR->centroid, fg->centroid, dxR); 682 for (d = 0; d < pdim; ++d) { 683 uL[iface*pdim+d] = xL[d] + DotD(dim, &gL[d*dim], dxL); 684 uR[iface*pdim+d] = xR[d] + DotD(dim, &gR[d*dim], dxR); 685 } 686 } else { 687 for (d = 0; d < pdim; ++d) { 688 uL[iface*pdim+d] = xL[d]; 689 uR[iface*pdim+d] = xR[d]; 690 } 691 } 692 for (d = 0; d < dim; ++d) { 693 centroid[iface*dim+d] = fg->centroid[d]; 694 normal[iface*dim+d] = fg->normal[d]; 695 } 696 vol[iface*2+0] = cgL->volume; 697 vol[iface*2+1] = cgR->volume; 698 ++iface; 699 } 700 if (computeGradients) { 701 ierr = VecRestoreArrayRead(locGrad,&lgrad);CHKERRQ(ierr); 702 ierr = DMRestoreLocalVector(dmplexts->dmGrad, &locGrad);CHKERRQ(ierr); 703 } 704 ierr = VecRestoreArrayRead(locX, &x);CHKERRQ(ierr); 705 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 706 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 707 fgeom.v0 = centroid; 708 fgeom.n = normal; 709 cgeom.vol = vol; 710 /* Riemann solve */ 711 ierr = PetscFVIntegrateRHSFunction(fvm, numFaces, Nf, &fvm, 0, fgeom, cgeom, uL, uR, riemann, fluxL, fluxR, dmplexts->rhsfunctionlocalctx);CHKERRQ(ierr); 712 /* Insert fluxes */ 713 ierr = VecGetArray(F, &f);CHKERRQ(ierr); 714 for (face = fStart, iface = 0; face < fEnd; ++face) { 715 const PetscInt *cells; 716 PetscScalar *fL, *fR; 717 PetscInt ghost, d; 718 719 ierr = DMLabelGetValue(ghostLabel, face, &ghost);CHKERRQ(ierr); 720 if (ghost >= 0) continue; 721 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 722 ierr = DMPlexPointGlobalRef(dm, cells[0], f, &fL);CHKERRQ(ierr); 723 ierr = DMPlexPointGlobalRef(dm, cells[1], f, &fR);CHKERRQ(ierr); 724 for (d = 0; d < pdim; ++d) { 725 if (fL) fL[d] -= fluxL[iface*pdim+d]; 726 if (fR) fR[d] += fluxR[iface*pdim+d]; 727 } 728 ++iface; 729 } 730 ierr = VecRestoreArray(F, &f);CHKERRQ(ierr); 731 ierr = PetscFree7(centroid,normal,vol,uL,uR,fluxL,fluxR);CHKERRQ(ierr); 732 ierr = DMRestoreLocalVector(dm, &locX);CHKERRQ(ierr); 733 PetscFunctionReturn(0); 734 } 735 736 #undef __FUNCT__ 737 #define __FUNCT__ "DMPlexTSSetRHSFunctionLocal" 738 /*@C 739 DMPlexTSSetRHSFunctionLocal - set a local residual evaluation function 740 741 Logically Collective 742 743 Input Arguments: 744 + dm - DM to associate callback with 745 . riemann - Riemann solver 746 - ctx - optional context for Riemann solve 747 748 Calling sequence for riemann: 749 750 $ riemann(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx) 751 752 + x - The coordinates at a point on the interface 753 . n - The normal vector to the interface 754 . uL - The state vector to the left of the interface 755 . uR - The state vector to the right of the interface 756 . flux - output array of flux through the interface 757 - ctx - optional user context 758 759 Level: beginner 760 761 .seealso: DMTSSetRHSFunctionLocal() 762 @*/ 763 PetscErrorCode DMPlexTSSetRHSFunctionLocal(DM dm, void (*riemann)(const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx), void *ctx) 764 { 765 DMTS dmts; 766 DMTS_Plex *dmplexts; 767 PetscFV fvm; 768 PetscInt Nf; 769 PetscBool computeGradients; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 774 ierr = DMGetDMTSWrite(dm, &dmts);CHKERRQ(ierr); 775 ierr = DMPlexTSGetContext(dm, dmts, &dmplexts);CHKERRQ(ierr); 776 dmplexts->riemann = riemann; 777 dmplexts->rhsfunctionlocalctx = ctx; 778 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 779 ierr = DMGetField(dm, 0, (PetscObject *) &fvm);CHKERRQ(ierr); 780 ierr = DMPlexTSSetupGeometry(dm, fvm, dmplexts);CHKERRQ(ierr); 781 ierr = PetscFVGetComputeGradients(fvm, &computeGradients);CHKERRQ(ierr); 782 if (computeGradients) {ierr = DMPlexTSSetupGradient(dm, fvm, dmplexts);CHKERRQ(ierr);} 783 ierr = DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMPlex, dmplexts);CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786