1 static char help[] = "Second Order TVD Finite Volume Example.\n"; 2 /*F 3 4 We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops 5 over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values, 6 \begin{equation} 7 f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R) 8 \end{equation} 9 and then update the cell values given the cell volume. 10 \begin{eqnarray} 11 f^L_i &-=& \frac{f_i}{vol^L} \\ 12 f^R_i &+=& \frac{f_i}{vol^R} 13 \end{eqnarray} 14 15 As an example, we can consider the shallow water wave equation, 16 \begin{eqnarray} 17 h_t + \nabla\cdot \left( uh \right) &=& 0 \\ 18 (uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0 19 \end{eqnarray} 20 where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity. 21 22 A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function, 23 \begin{eqnarray} 24 f^{L,R}_h &=& uh^{L,R} \cdot \hat n \\ 25 f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\ 26 c^{L,R} &=& \sqrt{g h^{L,R}} \\ 27 s &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\ 28 f_i &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right) 29 \end{eqnarray} 30 where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux. 31 32 The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial 33 over a neighborhood of the given element. 34 35 The mesh is read in from an ExodusII file, usually generated by Cubit. 36 37 The example also shows how to handle AMR in a time-dependent TS solver. 38 F*/ 39 #include <petscdmplex.h> 40 #include <petscdmforest.h> 41 #include <petscds.h> 42 #include <petscts.h> 43 44 #include "ex11.h" 45 46 static PetscFunctionList PhysicsList, PhysicsRiemannList_SW, PhysicsRiemannList_Euler; 47 48 /* 'User' implements a discretization of a continuous model. */ 49 typedef struct _n_User *User; 50 typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *); 51 typedef PetscErrorCode (*SetUpBCFunction)(DM, PetscDS, Physics); 52 typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *); 53 typedef PetscErrorCode (*SetupFields)(Physics, PetscSection); 54 static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *); 55 static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *); 56 static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *); 57 58 typedef struct _n_FunctionalLink *FunctionalLink; 59 struct _n_FunctionalLink { 60 char *name; 61 FunctionalFunction func; 62 void *ctx; 63 PetscInt offset; 64 FunctionalLink next; 65 }; 66 67 struct _n_Model { 68 MPI_Comm comm; /* Does not do collective communication, but some error conditions can be collective */ 69 Physics physics; 70 FunctionalLink functionalRegistry; 71 PetscInt maxComputed; 72 PetscInt numMonitored; 73 FunctionalLink *functionalMonitored; 74 PetscInt numCall; 75 FunctionalLink *functionalCall; 76 SolutionFunction solution; 77 SetUpBCFunction setupbc; 78 void *solutionctx; 79 PetscReal maxspeed; /* estimate of global maximum speed (for CFL calculation) */ 80 PetscReal bounds[2 * DIM]; 81 PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *); 82 void *errorCtx; 83 PetscErrorCode (*setupCEED)(DM, Physics); 84 }; 85 86 struct _n_User { 87 PetscInt vtkInterval; /* For monitor */ 88 char outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */ 89 PetscInt monitorStepOffset; 90 Model model; 91 PetscBool vtkmon; 92 }; 93 94 #ifdef PETSC_HAVE_LIBCEED 95 // Free a plain data context that was allocated using PETSc; returning libCEED error codes 96 static int FreeContextPetsc(void *data) 97 { 98 if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed"); 99 return CEED_ERROR_SUCCESS; 100 } 101 #endif 102 103 /******************* Advect ********************/ 104 typedef enum { 105 ADVECT_SOL_TILTED, 106 ADVECT_SOL_BUMP, 107 ADVECT_SOL_BUMP_CAVITY 108 } AdvectSolType; 109 static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0}; 110 typedef enum { 111 ADVECT_SOL_BUMP_CONE, 112 ADVECT_SOL_BUMP_COS 113 } AdvectSolBumpType; 114 static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0}; 115 116 typedef struct { 117 PetscReal wind[DIM]; 118 } Physics_Advect_Tilted; 119 typedef struct { 120 PetscReal center[DIM]; 121 PetscReal radius; 122 AdvectSolBumpType type; 123 } Physics_Advect_Bump; 124 125 typedef struct { 126 PetscReal inflowState; 127 AdvectSolType soltype; 128 union 129 { 130 Physics_Advect_Tilted tilted; 131 Physics_Advect_Bump bump; 132 } sol; 133 struct { 134 PetscInt Solution; 135 PetscInt Error; 136 } functional; 137 } Physics_Advect; 138 139 static const struct FieldDescription PhysicsFields_Advect[] = { 140 {"U", 1}, 141 {NULL, 0} 142 }; 143 144 static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 145 { 146 Physics phys = (Physics)ctx; 147 Physics_Advect *advect = (Physics_Advect *)phys->data; 148 149 PetscFunctionBeginUser; 150 xG[0] = advect->inflowState; 151 PetscFunctionReturn(PETSC_SUCCESS); 152 } 153 154 static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 155 { 156 PetscFunctionBeginUser; 157 xG[0] = xI[0]; 158 PetscFunctionReturn(PETSC_SUCCESS); 159 } 160 161 static void PhysicsRiemann_Advect(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys) 162 { 163 Physics_Advect *advect = (Physics_Advect *)phys->data; 164 PetscReal wind[DIM], wn; 165 166 switch (advect->soltype) { 167 case ADVECT_SOL_TILTED: { 168 Physics_Advect_Tilted *tilted = &advect->sol.tilted; 169 wind[0] = tilted->wind[0]; 170 wind[1] = tilted->wind[1]; 171 } break; 172 case ADVECT_SOL_BUMP: 173 wind[0] = -qp[1]; 174 wind[1] = qp[0]; 175 break; 176 case ADVECT_SOL_BUMP_CAVITY: { 177 PetscInt i; 178 PetscReal comp2[3] = {0., 0., 0.}, rad2; 179 180 rad2 = 0.; 181 for (i = 0; i < dim; i++) { 182 comp2[i] = qp[i] * qp[i]; 183 rad2 += comp2[i]; 184 } 185 186 wind[0] = -qp[1]; 187 wind[1] = qp[0]; 188 if (rad2 > 1.) { 189 PetscInt maxI = 0; 190 PetscReal maxComp2 = comp2[0]; 191 192 for (i = 1; i < dim; i++) { 193 if (comp2[i] > maxComp2) { 194 maxI = i; 195 maxComp2 = comp2[i]; 196 } 197 } 198 wind[maxI] = 0.; 199 } 200 } break; 201 default: { 202 PetscInt i; 203 for (i = 0; i < DIM; ++i) wind[i] = 0.0; 204 } 205 /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */ 206 } 207 wn = Dot2Real(wind, n); 208 flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn; 209 } 210 211 static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx) 212 { 213 Physics phys = (Physics)ctx; 214 Physics_Advect *advect = (Physics_Advect *)phys->data; 215 216 PetscFunctionBeginUser; 217 switch (advect->soltype) { 218 case ADVECT_SOL_TILTED: { 219 PetscReal x0[DIM]; 220 Physics_Advect_Tilted *tilted = &advect->sol.tilted; 221 Waxpy2Real(-time, tilted->wind, x, x0); 222 if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1]; 223 else u[0] = advect->inflowState; 224 } break; 225 case ADVECT_SOL_BUMP_CAVITY: 226 case ADVECT_SOL_BUMP: { 227 Physics_Advect_Bump *bump = &advect->sol.bump; 228 PetscReal x0[DIM], v[DIM], r, cost, sint; 229 cost = PetscCosReal(time); 230 sint = PetscSinReal(time); 231 x0[0] = cost * x[0] + sint * x[1]; 232 x0[1] = -sint * x[0] + cost * x[1]; 233 Waxpy2Real(-1, bump->center, x0, v); 234 r = Norm2Real(v); 235 switch (bump->type) { 236 case ADVECT_SOL_BUMP_CONE: 237 u[0] = PetscMax(1 - r / bump->radius, 0); 238 break; 239 case ADVECT_SOL_BUMP_COS: 240 u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI); 241 break; 242 } 243 } break; 244 default: 245 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type"); 246 } 247 PetscFunctionReturn(PETSC_SUCCESS); 248 } 249 250 static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx) 251 { 252 Physics phys = (Physics)ctx; 253 Physics_Advect *advect = (Physics_Advect *)phys->data; 254 PetscScalar yexact[1] = {0.0}; 255 256 PetscFunctionBeginUser; 257 PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys)); 258 f[advect->functional.Solution] = PetscRealPart(y[0]); 259 f[advect->functional.Error] = PetscAbsScalar(y[0] - yexact[0]); 260 PetscFunctionReturn(PETSC_SUCCESS); 261 } 262 263 static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys) 264 { 265 const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101}; 266 DMLabel label; 267 268 PetscFunctionBeginUser; 269 /* Register "canned" boundary conditions and defaults for where to apply. */ 270 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 271 PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Inflow, NULL, phys, NULL)); 272 PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Outflow, NULL, phys, NULL)); 273 PetscFunctionReturn(PETSC_SUCCESS); 274 } 275 276 static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject) 277 { 278 Physics_Advect *advect; 279 280 PetscFunctionBeginUser; 281 phys->field_desc = PhysicsFields_Advect; 282 phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Advect; 283 PetscCall(PetscNew(&advect)); 284 phys->data = advect; 285 mod->setupbc = SetUpBC_Advect; 286 287 PetscOptionsHeadBegin(PetscOptionsObject, "Advect options"); 288 { 289 PetscInt two = 2, dof = 1; 290 advect->soltype = ADVECT_SOL_TILTED; 291 PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL)); 292 switch (advect->soltype) { 293 case ADVECT_SOL_TILTED: { 294 Physics_Advect_Tilted *tilted = &advect->sol.tilted; 295 two = 2; 296 tilted->wind[0] = 0.0; 297 tilted->wind[1] = 1.0; 298 PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL)); 299 advect->inflowState = -2.0; 300 PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL)); 301 phys->maxspeed = Norm2Real(tilted->wind); 302 } break; 303 case ADVECT_SOL_BUMP_CAVITY: 304 case ADVECT_SOL_BUMP: { 305 Physics_Advect_Bump *bump = &advect->sol.bump; 306 two = 2; 307 bump->center[0] = 2.; 308 bump->center[1] = 0.; 309 PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL)); 310 bump->radius = 0.9; 311 PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL)); 312 bump->type = ADVECT_SOL_BUMP_CONE; 313 PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL)); 314 phys->maxspeed = 3.; /* radius of mesh, kludge */ 315 } break; 316 } 317 } 318 PetscOptionsHeadEnd(); 319 /* Initial/transient solution with default boundary conditions */ 320 PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys)); 321 /* Register "canned" functionals */ 322 PetscCall(ModelFunctionalRegister(mod, "Solution", &advect->functional.Solution, PhysicsFunctional_Advect, phys)); 323 PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys)); 324 PetscFunctionReturn(PETSC_SUCCESS); 325 } 326 327 /******************* Shallow Water ********************/ 328 static const struct FieldDescription PhysicsFields_SW[] = { 329 {"Height", 1 }, 330 {"Momentum", DIM}, 331 {NULL, 0 } 332 }; 333 334 static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 335 { 336 PetscFunctionBeginUser; 337 xG[0] = xI[0]; 338 xG[1] = -xI[1]; 339 xG[2] = -xI[2]; 340 PetscFunctionReturn(PETSC_SUCCESS); 341 } 342 343 static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx) 344 { 345 PetscReal dx[2], r, sigma; 346 347 PetscFunctionBeginUser; 348 PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time); 349 dx[0] = x[0] - 1.5; 350 dx[1] = x[1] - 1.0; 351 r = Norm2Real(dx); 352 sigma = 0.5; 353 u[0] = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma))); 354 u[1] = 0.0; 355 u[2] = 0.0; 356 PetscFunctionReturn(PETSC_SUCCESS); 357 } 358 359 static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx) 360 { 361 Physics phys = (Physics)ctx; 362 Physics_SW *sw = (Physics_SW *)phys->data; 363 const SWNode *x = (const SWNode *)xx; 364 PetscReal u[2]; 365 PetscReal h; 366 367 PetscFunctionBeginUser; 368 h = x->h; 369 Scale2Real(1. / x->h, x->uh, u); 370 f[sw->functional.Height] = h; 371 f[sw->functional.Speed] = Norm2Real(u) + PetscSqrtReal(sw->gravity * h); 372 f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h)); 373 PetscFunctionReturn(PETSC_SUCCESS); 374 } 375 376 static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys) 377 { 378 const PetscInt wallids[] = {100, 101, 200, 300}; 379 DMLabel label; 380 381 PetscFunctionBeginUser; 382 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 383 PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_SW_Wall, NULL, phys, NULL)); 384 PetscFunctionReturn(PETSC_SUCCESS); 385 } 386 387 #ifdef PETSC_HAVE_LIBCEED 388 static PetscErrorCode CreateQFunctionContext_SW(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx) 389 { 390 Physics_SW *in = (Physics_SW *)phys->data; 391 Physics_SW *sw; 392 393 PetscFunctionBeginUser; 394 PetscCall(PetscCalloc1(1, &sw)); 395 396 sw->gravity = in->gravity; 397 398 PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx)); 399 PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*sw), sw)); 400 PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc)); 401 PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gravity", offsetof(Physics_SW, gravity), 1, "Accelaration due to gravity")); 402 PetscFunctionReturn(0); 403 } 404 #endif 405 406 static PetscErrorCode SetupCEED_SW(DM dm, Physics physics) 407 { 408 #ifdef PETSC_HAVE_LIBCEED 409 Ceed ceed; 410 CeedQFunctionContext qfCtx; 411 #endif 412 413 PetscFunctionBegin; 414 #ifdef PETSC_HAVE_LIBCEED 415 PetscCall(DMGetCeed(dm, &ceed)); 416 PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx)); 417 PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx)); 418 PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx)); 419 #endif 420 PetscFunctionReturn(PETSC_SUCCESS); 421 } 422 423 static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject) 424 { 425 Physics_SW *sw; 426 char sw_riemann[64] = "rusanov"; 427 428 PetscFunctionBeginUser; 429 phys->field_desc = PhysicsFields_SW; 430 PetscCall(PetscNew(&sw)); 431 phys->data = sw; 432 mod->setupbc = SetUpBC_SW; 433 mod->setupCEED = SetupCEED_SW; 434 435 PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov)); 436 PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL)); 437 #ifdef PETSC_HAVE_LIBCEED 438 PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED)); 439 #endif 440 441 PetscOptionsHeadBegin(PetscOptionsObject, "SW options"); 442 { 443 void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics); 444 sw->gravity = 1.0; 445 PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL)); 446 PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL)); 447 PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW)); 448 phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW; 449 } 450 PetscOptionsHeadEnd(); 451 phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */ 452 453 PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys)); 454 PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys)); 455 PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys)); 456 PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys)); 457 PetscFunctionReturn(PETSC_SUCCESS); 458 } 459 460 /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/ 461 /* An initial-value and self-similar solutions of the compressible Euler equations */ 462 /* Ravi Samtaney and D. I. Pullin */ 463 /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */ 464 static const struct FieldDescription PhysicsFields_Euler[] = { 465 {"Density", 1 }, 466 {"Momentum", DIM}, 467 {"Energy", 1 }, 468 {NULL, 0 } 469 }; 470 471 /* initial condition */ 472 int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx); 473 static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx) 474 { 475 PetscInt i; 476 Physics phys = (Physics)ctx; 477 Physics_Euler *eu = (Physics_Euler *)phys->data; 478 EulerNode *uu = (EulerNode *)u; 479 PetscReal p0, gamma, c = 0.; 480 481 PetscFunctionBeginUser; 482 PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time); 483 484 for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */ 485 /* set E and rho */ 486 gamma = eu->gamma; 487 488 if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) { 489 /******************* Euler Density Shock ********************/ 490 /* On initial-value and self-similar solutions of the compressible Euler equations */ 491 /* Ravi Samtaney and D. I. Pullin */ 492 /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */ 493 /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity, */ 494 p0 = 1.; 495 if (x[0] < 0.0 + x[1] * eu->itana) { 496 if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */ 497 PetscReal amach, rho, press, gas1, p1; 498 amach = eu->amach; 499 rho = 1.; 500 press = p0; 501 p1 = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0)); 502 gas1 = (gamma - 1.0) / (gamma + 1.0); 503 uu->r = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0); 504 uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach); 505 uu->E = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0]; 506 } else { /* left of discontinuity (0) */ 507 uu->r = 1.; /* rho = 1 */ 508 uu->E = p0 / (gamma - 1.0); 509 } 510 } else { /* right of discontinuity (2) */ 511 uu->r = eu->rhoR; 512 uu->E = p0 / (gamma - 1.0); 513 } 514 } else if (eu->type == EULER_SHOCK_TUBE) { 515 /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */ 516 if (x[0] < 0.0) { 517 uu->r = 8.; 518 uu->E = 10. / (gamma - 1.); 519 } else { 520 uu->r = 1.; 521 uu->E = 1. / (gamma - 1.); 522 } 523 } else if (eu->type == EULER_LINEAR_WAVE) { 524 initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]); 525 } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type); 526 527 /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */ 528 PetscCall(SpeedOfSound_PG(gamma, uu, &c)); 529 c = (uu->ru[0] / uu->r) + c; 530 if (c > phys->maxspeed) phys->maxspeed = c; 531 PetscFunctionReturn(PETSC_SUCCESS); 532 } 533 534 /* PetscReal* => EulerNode* conversion */ 535 static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx) 536 { 537 PetscInt i; 538 const EulerNode *xI = (const EulerNode *)a_xI; 539 EulerNode *xG = (EulerNode *)a_xG; 540 Physics phys = (Physics)ctx; 541 Physics_Euler *eu = (Physics_Euler *)phys->data; 542 543 PetscFunctionBeginUser; 544 xG->r = xI->r; /* ghost cell density - same */ 545 xG->E = xI->E; /* ghost cell energy - same */ 546 if (n[1] != 0.) { /* top and bottom */ 547 xG->ru[0] = xI->ru[0]; /* copy tang to wall */ 548 xG->ru[1] = -xI->ru[1]; /* reflect perp to t/b wall */ 549 } else { /* sides */ 550 for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */ 551 } 552 if (eu->type == EULER_LINEAR_WAVE) { /* debug */ 553 #if 0 554 PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]); 555 #endif 556 } 557 PetscFunctionReturn(PETSC_SUCCESS); 558 } 559 560 static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx) 561 { 562 Physics phys = (Physics)ctx; 563 Physics_Euler *eu = (Physics_Euler *)phys->data; 564 const EulerNode *x = (const EulerNode *)xx; 565 PetscReal p; 566 567 PetscFunctionBeginUser; 568 f[eu->monitor.Density] = x->r; 569 f[eu->monitor.Momentum] = NormDIM(x->ru); 570 f[eu->monitor.Energy] = x->E; 571 f[eu->monitor.Speed] = NormDIM(x->ru) / x->r; 572 PetscCall(Pressure_PG(eu->gamma, x, &p)); 573 f[eu->monitor.Pressure] = p; 574 PetscFunctionReturn(PETSC_SUCCESS); 575 } 576 577 static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys) 578 { 579 Physics_Euler *eu = (Physics_Euler *)phys->data; 580 DMLabel label; 581 582 PetscFunctionBeginUser; 583 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 584 if (eu->type == EULER_LINEAR_WAVE) { 585 const PetscInt wallids[] = {100, 101}; 586 PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL)); 587 } else { 588 const PetscInt wallids[] = {100, 101, 200, 300}; 589 PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL)); 590 } 591 PetscFunctionReturn(PETSC_SUCCESS); 592 } 593 594 #ifdef PETSC_HAVE_LIBCEED 595 static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx) 596 { 597 Physics_Euler *in = (Physics_Euler *)phys->data; 598 Physics_Euler *eu; 599 600 PetscFunctionBeginUser; 601 PetscCall(PetscCalloc1(1, &eu)); 602 603 eu->gamma = in->gamma; 604 605 PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx)); 606 PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu)); 607 PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc)); 608 PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio")); 609 PetscFunctionReturn(0); 610 } 611 #endif 612 613 static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics) 614 { 615 #ifdef PETSC_HAVE_LIBCEED 616 Ceed ceed; 617 CeedQFunctionContext qfCtx; 618 #endif 619 620 PetscFunctionBegin; 621 #ifdef PETSC_HAVE_LIBCEED 622 PetscCall(DMGetCeed(dm, &ceed)); 623 PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx)); 624 PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx)); 625 PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx)); 626 #endif 627 PetscFunctionReturn(PETSC_SUCCESS); 628 } 629 630 static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject) 631 { 632 Physics_Euler *eu; 633 634 PetscFunctionBeginUser; 635 phys->field_desc = PhysicsFields_Euler; 636 phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov; 637 PetscCall(PetscNew(&eu)); 638 phys->data = eu; 639 mod->setupbc = SetUpBC_Euler; 640 mod->setupCEED = SetupCEED_Euler; 641 642 PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov)); 643 #ifdef PETSC_HAVE_LIBCEED 644 PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED)); 645 #endif 646 647 PetscOptionsHeadBegin(PetscOptionsObject, "Euler options"); 648 { 649 void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics); 650 PetscReal alpha; 651 char type[64] = "linear_wave"; 652 char eu_riemann[64] = "godunov"; 653 PetscBool is; 654 eu->gamma = 1.4; 655 eu->amach = 2.02; 656 eu->rhoR = 3.0; 657 eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */ 658 PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL)); 659 PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler)); 660 phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler; 661 PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL)); 662 PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL)); 663 PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL)); 664 alpha = 60.; 665 PetscCall(PetscOptionsRangeReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL, 0.0, 90.0)); 666 eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0); 667 PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL)); 668 PetscCall(PetscStrcmp(type, "linear_wave", &is)); 669 if (is) { 670 /* Remember this should be periodic */ 671 eu->type = EULER_LINEAR_WAVE; 672 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave")); 673 } else { 674 PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type); 675 PetscCall(PetscStrcmp(type, "iv_shock", &is)); 676 if (is) { 677 eu->type = EULER_IV_SHOCK; 678 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock")); 679 } else { 680 PetscCall(PetscStrcmp(type, "ss_shock", &is)); 681 if (is) { 682 eu->type = EULER_SS_SHOCK; 683 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock")); 684 } else { 685 PetscCall(PetscStrcmp(type, "shock_tube", &is)); 686 if (is) eu->type = EULER_SHOCK_TUBE; 687 else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type); 688 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube")); 689 } 690 } 691 } 692 } 693 PetscOptionsHeadEnd(); 694 phys->maxspeed = 0.; /* will get set in solution */ 695 PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys)); 696 PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys)); 697 PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys)); 698 PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys)); 699 PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys)); 700 PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys)); 701 PetscFunctionReturn(PETSC_SUCCESS); 702 } 703 704 static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx) 705 { 706 PetscReal err = 0.; 707 PetscInt i, j; 708 709 PetscFunctionBeginUser; 710 for (i = 0; i < numComps; i++) { 711 for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j])); 712 } 713 *error = volume * err; 714 PetscFunctionReturn(PETSC_SUCCESS); 715 } 716 717 PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition) 718 { 719 PetscSF sfPoint; 720 PetscSection coordSection; 721 Vec coordinates; 722 PetscSection sectionCell; 723 PetscScalar *part; 724 PetscInt cStart, cEnd, c; 725 PetscMPIInt rank; 726 727 PetscFunctionBeginUser; 728 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 729 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 730 PetscCall(DMClone(dm, dmCell)); 731 PetscCall(DMGetPointSF(dm, &sfPoint)); 732 PetscCall(DMSetPointSF(*dmCell, sfPoint)); 733 PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection)); 734 PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates)); 735 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank)); 736 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell)); 737 PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd)); 738 PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd)); 739 for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1)); 740 PetscCall(PetscSectionSetUp(sectionCell)); 741 PetscCall(DMSetLocalSection(*dmCell, sectionCell)); 742 PetscCall(PetscSectionDestroy(§ionCell)); 743 PetscCall(DMCreateLocalVector(*dmCell, partition)); 744 PetscCall(PetscObjectSetName((PetscObject)*partition, "partition")); 745 PetscCall(VecGetArray(*partition, &part)); 746 for (c = cStart; c < cEnd; ++c) { 747 PetscScalar *p; 748 749 PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p)); 750 p[0] = rank; 751 } 752 PetscCall(VecRestoreArray(*partition, &part)); 753 PetscFunctionReturn(PETSC_SUCCESS); 754 } 755 756 PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user) 757 { 758 DM plex, dmMass, dmFace, dmCell, dmCoord; 759 PetscSection coordSection; 760 Vec coordinates, facegeom, cellgeom; 761 PetscSection sectionMass; 762 PetscScalar *m; 763 const PetscScalar *fgeom, *cgeom, *coords; 764 PetscInt vStart, vEnd, v; 765 766 PetscFunctionBeginUser; 767 PetscCall(DMConvert(dm, DMPLEX, &plex)); 768 PetscCall(DMGetCoordinateSection(dm, &coordSection)); 769 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 770 PetscCall(DMClone(dm, &dmMass)); 771 PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection)); 772 PetscCall(DMSetCoordinatesLocal(dmMass, coordinates)); 773 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionMass)); 774 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 775 PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd)); 776 for (v = vStart; v < vEnd; ++v) { 777 PetscInt numFaces; 778 779 PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces)); 780 PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces)); 781 } 782 PetscCall(PetscSectionSetUp(sectionMass)); 783 PetscCall(DMSetLocalSection(dmMass, sectionMass)); 784 PetscCall(PetscSectionDestroy(§ionMass)); 785 PetscCall(DMGetLocalVector(dmMass, massMatrix)); 786 PetscCall(VecGetArray(*massMatrix, &m)); 787 PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL)); 788 PetscCall(VecGetDM(facegeom, &dmFace)); 789 PetscCall(VecGetArrayRead(facegeom, &fgeom)); 790 PetscCall(VecGetDM(cellgeom, &dmCell)); 791 PetscCall(VecGetArrayRead(cellgeom, &cgeom)); 792 PetscCall(DMGetCoordinateDM(dm, &dmCoord)); 793 PetscCall(VecGetArrayRead(coordinates, &coords)); 794 for (v = vStart; v < vEnd; ++v) { 795 const PetscInt *faces; 796 PetscFVFaceGeom *fgA, *fgB, *cg; 797 PetscScalar *vertex; 798 PetscInt numFaces, sides[2], f, g; 799 800 PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex)); 801 PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces)); 802 PetscCall(DMPlexGetSupport(dmMass, v, &faces)); 803 for (f = 0; f < numFaces; ++f) { 804 sides[0] = faces[f]; 805 PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA)); 806 for (g = 0; g < numFaces; ++g) { 807 const PetscInt *cells = NULL; 808 PetscReal area = 0.0; 809 PetscInt numCells; 810 811 sides[1] = faces[g]; 812 PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB)); 813 PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells)); 814 PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces"); 815 PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg)); 816 area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0])); 817 area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0])); 818 m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5; 819 PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells)); 820 } 821 } 822 } 823 PetscCall(VecRestoreArrayRead(facegeom, &fgeom)); 824 PetscCall(VecRestoreArrayRead(cellgeom, &cgeom)); 825 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 826 PetscCall(VecRestoreArray(*massMatrix, &m)); 827 PetscCall(DMDestroy(&dmMass)); 828 PetscCall(DMDestroy(&plex)); 829 PetscFunctionReturn(PETSC_SUCCESS); 830 } 831 832 /* Behavior will be different for multi-physics or when using non-default boundary conditions */ 833 static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx) 834 { 835 PetscFunctionBeginUser; 836 mod->solution = func; 837 mod->solutionctx = ctx; 838 PetscFunctionReturn(PETSC_SUCCESS); 839 } 840 841 static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx) 842 { 843 FunctionalLink link, *ptr; 844 PetscInt lastoffset = -1; 845 846 PetscFunctionBeginUser; 847 for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset; 848 PetscCall(PetscNew(&link)); 849 PetscCall(PetscStrallocpy(name, &link->name)); 850 link->offset = lastoffset + 1; 851 link->func = func; 852 link->ctx = ctx; 853 link->next = NULL; 854 *ptr = link; 855 *offset = link->offset; 856 PetscFunctionReturn(PETSC_SUCCESS); 857 } 858 859 static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject) 860 { 861 PetscInt i, j; 862 FunctionalLink link; 863 char *names[256]; 864 865 PetscFunctionBeginUser; 866 mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names); 867 PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL)); 868 /* Create list of functionals that will be computed somehow */ 869 PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored)); 870 /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */ 871 PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall)); 872 mod->numCall = 0; 873 for (i = 0; i < mod->numMonitored; i++) { 874 for (link = mod->functionalRegistry; link; link = link->next) { 875 PetscBool match; 876 PetscCall(PetscStrcasecmp(names[i], link->name, &match)); 877 if (match) break; 878 } 879 PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]); 880 mod->functionalMonitored[i] = link; 881 for (j = 0; j < i; j++) { 882 if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name; 883 } 884 mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */ 885 next_name: 886 PetscCall(PetscFree(names[i])); 887 } 888 889 /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */ 890 mod->maxComputed = -1; 891 for (link = mod->functionalRegistry; link; link = link->next) { 892 for (i = 0; i < mod->numCall; i++) { 893 FunctionalLink call = mod->functionalCall[i]; 894 if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset); 895 } 896 } 897 PetscFunctionReturn(PETSC_SUCCESS); 898 } 899 900 static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link) 901 { 902 FunctionalLink l, next; 903 904 PetscFunctionBeginUser; 905 if (!link) PetscFunctionReturn(PETSC_SUCCESS); 906 l = *link; 907 *link = NULL; 908 for (; l; l = next) { 909 next = l->next; 910 PetscCall(PetscFree(l->name)); 911 PetscCall(PetscFree(l)); 912 } 913 PetscFunctionReturn(PETSC_SUCCESS); 914 } 915 916 /* put the solution callback into a functional callback */ 917 static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx) 918 { 919 Model mod; 920 921 PetscFunctionBeginUser; 922 mod = (Model)modctx; 923 PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx)); 924 PetscFunctionReturn(PETSC_SUCCESS); 925 } 926 927 PetscErrorCode SetInitialCondition(DM dm, Vec X, User user) 928 { 929 PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 930 void *ctx[1]; 931 Model mod = user->model; 932 933 PetscFunctionBeginUser; 934 func[0] = SolutionFunctional; 935 ctx[0] = (void *)mod; 936 PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X)); 937 PetscFunctionReturn(PETSC_SUCCESS); 938 } 939 940 static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer) 941 { 942 PetscFunctionBeginUser; 943 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer)); 944 PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK)); 945 PetscCall(PetscViewerFileSetName(*viewer, filename)); 946 PetscFunctionReturn(PETSC_SUCCESS); 947 } 948 949 static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx) 950 { 951 User user = (User)ctx; 952 DM dm, plex; 953 PetscViewer viewer; 954 char filename[PETSC_MAX_PATH_LEN], *ftable = NULL; 955 PetscReal xnorm; 956 957 PetscFunctionBeginUser; 958 PetscCall(PetscObjectSetName((PetscObject)X, "u")); 959 PetscCall(VecGetDM(X, &dm)); 960 PetscCall(VecNorm(X, NORM_INFINITY, &xnorm)); 961 962 if (stepnum >= 0) stepnum += user->monitorStepOffset; 963 if (stepnum >= 0) { /* No summary for final time */ 964 Model mod = user->model; 965 Vec cellgeom; 966 PetscInt c, cStart, cEnd, fcount, i; 967 size_t ftableused, ftablealloc; 968 const PetscScalar *cgeom, *x; 969 DM dmCell; 970 DMLabel vtkLabel; 971 PetscReal *fmin, *fmax, *fintegral, *ftmp; 972 973 PetscCall(DMConvert(dm, DMPLEX, &plex)); 974 PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL)); 975 fcount = mod->maxComputed + 1; 976 PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp)); 977 for (i = 0; i < fcount; i++) { 978 fmin[i] = PETSC_MAX_REAL; 979 fmax[i] = PETSC_MIN_REAL; 980 fintegral[i] = 0; 981 } 982 PetscCall(VecGetDM(cellgeom, &dmCell)); 983 PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd)); 984 PetscCall(VecGetArrayRead(cellgeom, &cgeom)); 985 PetscCall(VecGetArrayRead(X, &x)); 986 PetscCall(DMGetLabel(dm, "vtk", &vtkLabel)); 987 for (c = cStart; c < cEnd; ++c) { 988 PetscFVCellGeom *cg; 989 const PetscScalar *cx = NULL; 990 PetscInt vtkVal = 0; 991 992 /* not that these two routines as currently implemented work for any dm with a 993 * localSection/globalSection */ 994 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 995 PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx)); 996 if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal)); 997 if (!vtkVal || !cx) continue; /* ghost, or not a global cell */ 998 for (i = 0; i < mod->numCall; i++) { 999 FunctionalLink flink = mod->functionalCall[i]; 1000 PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx)); 1001 } 1002 for (i = 0; i < fcount; i++) { 1003 fmin[i] = PetscMin(fmin[i], ftmp[i]); 1004 fmax[i] = PetscMax(fmax[i], ftmp[i]); 1005 fintegral[i] += cg->volume * ftmp[i]; 1006 } 1007 } 1008 PetscCall(VecRestoreArrayRead(cellgeom, &cgeom)); 1009 PetscCall(VecRestoreArrayRead(X, &x)); 1010 PetscCall(DMDestroy(&plex)); 1011 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts))); 1012 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts))); 1013 PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts))); 1014 1015 ftablealloc = fcount * 100; 1016 ftableused = 0; 1017 PetscCall(PetscMalloc1(ftablealloc, &ftable)); 1018 for (i = 0; i < mod->numMonitored; i++) { 1019 size_t countused; 1020 char buffer[256], *p; 1021 FunctionalLink flink = mod->functionalMonitored[i]; 1022 PetscInt id = flink->offset; 1023 if (i % 3) { 1024 PetscCall(PetscArraycpy(buffer, " ", 2)); 1025 p = buffer + 2; 1026 } else if (i) { 1027 char newline[] = "\n"; 1028 PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1)); 1029 p = buffer + sizeof(newline) - 1; 1030 } else { 1031 p = buffer; 1032 } 1033 PetscCall(PetscSNPrintfCount(p, sizeof buffer - (p - buffer), "%12s [%10.7g,%10.7g] int %10.7g", &countused, flink->name, (double)fmin[id], (double)fmax[id], (double)fintegral[id])); 1034 countused--; 1035 countused += p - buffer; 1036 if (countused > ftablealloc - ftableused - 1) { /* reallocate */ 1037 char *ftablenew; 1038 ftablealloc = 2 * ftablealloc + countused; 1039 PetscCall(PetscMalloc(ftablealloc, &ftablenew)); 1040 PetscCall(PetscArraycpy(ftablenew, ftable, ftableused)); 1041 PetscCall(PetscFree(ftable)); 1042 ftable = ftablenew; 1043 } 1044 PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused)); 1045 ftableused += countused; 1046 ftable[ftableused] = 0; 1047 } 1048 PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp)); 1049 1050 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT " time %8.4g |x| %8.4g %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : "")); 1051 PetscCall(PetscFree(ftable)); 1052 } 1053 if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS); 1054 if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) { 1055 if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */ 1056 PetscCall(TSGetStepNumber(ts, &stepnum)); 1057 } 1058 PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum)); 1059 PetscCall(OutputVTK(dm, filename, &viewer)); 1060 PetscCall(VecView(X, viewer)); 1061 PetscCall(PetscViewerDestroy(&viewer)); 1062 } 1063 PetscFunctionReturn(PETSC_SUCCESS); 1064 } 1065 1066 static PetscErrorCode initializeTS(DM dm, User user, TS *ts) 1067 { 1068 PetscBool useCeed; 1069 1070 PetscFunctionBeginUser; 1071 PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts)); 1072 PetscCall(TSSetType(*ts, TSSSP)); 1073 PetscCall(TSSetDM(*ts, dm)); 1074 if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL)); 1075 PetscCall(DMPlexGetUseCeed(dm, &useCeed)); 1076 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user)); 1077 if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user)); 1078 else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user)); 1079 PetscCall(TSSetMaxTime(*ts, 2.0)); 1080 PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER)); 1081 PetscFunctionReturn(PETSC_SUCCESS); 1082 } 1083 1084 typedef struct { 1085 PetscFV fvm; 1086 VecTagger refineTag; 1087 VecTagger coarsenTag; 1088 DM adaptedDM; 1089 User user; 1090 PetscReal cfl; 1091 PetscLimiter limiter; 1092 PetscLimiter noneLimiter; 1093 } TransferCtx; 1094 1095 static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx) 1096 { 1097 TransferCtx *tctx = (TransferCtx *)ctx; 1098 PetscFV fvm = tctx->fvm; 1099 VecTagger refineTag = tctx->refineTag; 1100 VecTagger coarsenTag = tctx->coarsenTag; 1101 User user = tctx->user; 1102 DM dm, gradDM, plex, cellDM, adaptedDM = NULL; 1103 Vec cellGeom, faceGeom; 1104 PetscBool isForest, computeGradient; 1105 Vec grad, locGrad, locX, errVec; 1106 PetscInt cStart, cEnd, c, dim, nRefine, nCoarsen; 1107 PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2]; 1108 PetscScalar *errArray; 1109 const PetscScalar *pointVals; 1110 const PetscScalar *pointGrads; 1111 const PetscScalar *pointGeom; 1112 DMLabel adaptLabel = NULL; 1113 IS refineIS, coarsenIS; 1114 1115 PetscFunctionBeginUser; 1116 *resize = PETSC_FALSE; 1117 PetscCall(VecGetDM(sol, &dm)); 1118 PetscCall(DMGetDimension(dm, &dim)); 1119 PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter)); 1120 PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient)); 1121 PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE)); 1122 PetscCall(DMIsForest(dm, &isForest)); 1123 PetscCall(DMConvert(dm, DMPLEX, &plex)); 1124 PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM)); 1125 PetscCall(DMCreateLocalVector(plex, &locX)); 1126 PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL)); 1127 PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX)); 1128 PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX)); 1129 PetscCall(DMCreateGlobalVector(gradDM, &grad)); 1130 PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad)); 1131 PetscCall(DMCreateLocalVector(gradDM, &locGrad)); 1132 PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad)); 1133 PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad)); 1134 PetscCall(VecDestroy(&grad)); 1135 PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd)); 1136 PetscCall(VecGetArrayRead(locGrad, &pointGrads)); 1137 PetscCall(VecGetArrayRead(cellGeom, &pointGeom)); 1138 PetscCall(VecGetArrayRead(locX, &pointVals)); 1139 PetscCall(VecGetDM(cellGeom, &cellDM)); 1140 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel)); 1141 PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec)); 1142 PetscCall(VecSetUp(errVec)); 1143 PetscCall(VecGetArray(errVec, &errArray)); 1144 for (c = cStart; c < cEnd; c++) { 1145 PetscReal errInd = 0.; 1146 PetscScalar *pointGrad; 1147 PetscScalar *pointVal; 1148 PetscFVCellGeom *cg; 1149 1150 PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad)); 1151 PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg)); 1152 PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal)); 1153 1154 PetscCall(user->model->errorIndicator(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx)); 1155 errArray[c - cStart] = errInd; 1156 minMaxInd[0] = PetscMin(minMaxInd[0], errInd); 1157 minMaxInd[1] = PetscMax(minMaxInd[1], errInd); 1158 } 1159 PetscCall(VecRestoreArray(errVec, &errArray)); 1160 PetscCall(VecRestoreArrayRead(locX, &pointVals)); 1161 PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom)); 1162 PetscCall(VecRestoreArrayRead(locGrad, &pointGrads)); 1163 PetscCall(VecDestroy(&locGrad)); 1164 PetscCall(VecDestroy(&locX)); 1165 PetscCall(DMDestroy(&plex)); 1166 1167 PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL)); 1168 PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL)); 1169 PetscCall(ISGetSize(refineIS, &nRefine)); 1170 PetscCall(ISGetSize(coarsenIS, &nCoarsen)); 1171 if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS)); 1172 if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS)); 1173 PetscCall(ISDestroy(&coarsenIS)); 1174 PetscCall(ISDestroy(&refineIS)); 1175 PetscCall(VecDestroy(&errVec)); 1176 1177 PetscCall(PetscFVSetComputeGradients(fvm, computeGradient)); 1178 PetscCall(PetscFVSetLimiter(fvm, tctx->limiter)); 1179 minMaxInd[1] = -minMaxInd[1]; 1180 PetscCall(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm))); 1181 PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1]))); 1182 if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */ 1183 PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM)); 1184 } 1185 PetscCall(DMLabelDestroy(&adaptLabel)); 1186 if (adaptedDM) { 1187 tctx->adaptedDM = adaptedDM; 1188 PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen)); 1189 *resize = PETSC_TRUE; 1190 } else { 1191 PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep)); 1192 } 1193 PetscFunctionReturn(PETSC_SUCCESS); 1194 } 1195 1196 static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx) 1197 { 1198 TransferCtx *tctx = (TransferCtx *)ctx; 1199 DM dm; 1200 PetscReal time; 1201 1202 PetscFunctionBeginUser; 1203 PetscCall(TSGetDM(ts, &dm)); 1204 PetscCall(TSGetTime(ts, &time)); 1205 PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM"); 1206 for (PetscInt i = 0; i < nv; i++) { 1207 const char *name; 1208 1209 PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i])); 1210 PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time)); 1211 PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name)); 1212 PetscCall(PetscObjectSetName((PetscObject)vecsout[i], name)); 1213 } 1214 PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */ 1215 1216 Model mod = tctx->user->model; 1217 Physics phys = mod->physics; 1218 PetscReal minRadius; 1219 1220 PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius)); 1221 PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts))); 1222 PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed"); 1223 1224 PetscReal dt = tctx->cfl * minRadius / mod->maxspeed; 1225 PetscCall(TSSetTimeStep(ts, dt)); 1226 1227 PetscCall(TSSetDM(ts, tctx->adaptedDM)); 1228 PetscCall(DMDestroy(&tctx->adaptedDM)); 1229 PetscFunctionReturn(PETSC_SUCCESS); 1230 } 1231 1232 int main(int argc, char **argv) 1233 { 1234 MPI_Comm comm; 1235 PetscDS prob; 1236 PetscFV fvm; 1237 PetscLimiter limiter = NULL, noneLimiter = NULL; 1238 User user; 1239 Model mod; 1240 Physics phys; 1241 DM dm, plex; 1242 PetscReal ftime, cfl, dt, minRadius; 1243 PetscInt dim, nsteps; 1244 TS ts; 1245 TSConvergedReason reason; 1246 Vec X; 1247 PetscViewer viewer; 1248 PetscBool vtkCellGeom, useAMR; 1249 PetscInt adaptInterval; 1250 char physname[256] = "advect"; 1251 VecTagger refineTag = NULL, coarsenTag = NULL; 1252 TransferCtx tctx; 1253 1254 PetscFunctionBeginUser; 1255 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 1256 comm = PETSC_COMM_WORLD; 1257 1258 PetscCall(PetscNew(&user)); 1259 PetscCall(PetscNew(&user->model)); 1260 PetscCall(PetscNew(&user->model->physics)); 1261 mod = user->model; 1262 phys = mod->physics; 1263 mod->comm = comm; 1264 useAMR = PETSC_FALSE; 1265 adaptInterval = 1; 1266 1267 /* Register physical models to be available on the command line */ 1268 PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect)); 1269 PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW)); 1270 PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler)); 1271 1272 PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", ""); 1273 { 1274 cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */ 1275 PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL)); 1276 user->vtkInterval = 1; 1277 PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL)); 1278 user->vtkmon = PETSC_TRUE; 1279 PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL)); 1280 vtkCellGeom = PETSC_FALSE; 1281 PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename))); 1282 PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL)); 1283 PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL)); 1284 PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL)); 1285 PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL)); 1286 } 1287 PetscOptionsEnd(); 1288 1289 if (useAMR) { 1290 VecTaggerBox refineBox, coarsenBox; 1291 1292 refineBox.min = refineBox.max = PETSC_MAX_REAL; 1293 coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL; 1294 1295 PetscCall(VecTaggerCreate(comm, &refineTag)); 1296 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_")); 1297 PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE)); 1298 PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox)); 1299 PetscCall(VecTaggerSetFromOptions(refineTag)); 1300 PetscCall(VecTaggerSetUp(refineTag)); 1301 PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view")); 1302 1303 PetscCall(VecTaggerCreate(comm, &coarsenTag)); 1304 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_")); 1305 PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE)); 1306 PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox)); 1307 PetscCall(VecTaggerSetFromOptions(coarsenTag)); 1308 PetscCall(VecTaggerSetUp(coarsenTag)); 1309 PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view")); 1310 } 1311 1312 PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", ""); 1313 { 1314 PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *); 1315 PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL)); 1316 PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate)); 1317 PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics))); 1318 PetscCall((*physcreate)(mod, phys, PetscOptionsObject)); 1319 /* Count number of fields and dofs */ 1320 for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof; 1321 PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname); 1322 PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject)); 1323 } 1324 PetscOptionsEnd(); 1325 1326 /* Create mesh */ 1327 { 1328 PetscInt i; 1329 1330 PetscCall(DMCreate(comm, &dm)); 1331 PetscCall(DMSetType(dm, DMPLEX)); 1332 PetscCall(DMSetFromOptions(dm)); 1333 for (i = 0; i < DIM; i++) { 1334 mod->bounds[2 * i] = 0.; 1335 mod->bounds[2 * i + 1] = 1.; 1336 }; 1337 dim = DIM; 1338 { /* a null name means just do a hex box */ 1339 PetscInt cells[3] = {1, 1, 1}, n = 3; 1340 PetscBool flg2, skew = PETSC_FALSE; 1341 PetscInt nret2 = 2 * DIM; 1342 PetscOptionsBegin(comm, NULL, "Rectangular mesh options", ""); 1343 PetscCall(PetscOptionsRealArray("-grid_bounds", "bounds of the mesh in each direction (i.e., x_min,x_max,y_min,y_max", "", mod->bounds, &nret2, &flg2)); 1344 PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL)); 1345 PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL)); 1346 PetscOptionsEnd(); 1347 /* TODO Rewrite this with Mark, and remove grid_bounds at that time */ 1348 if (flg2) { 1349 PetscInt dimEmbed, i; 1350 PetscInt nCoords; 1351 PetscScalar *coords; 1352 Vec coordinates; 1353 1354 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 1355 PetscCall(DMGetCoordinateDim(dm, &dimEmbed)); 1356 PetscCall(VecGetLocalSize(coordinates, &nCoords)); 1357 PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size"); 1358 PetscCall(VecGetArray(coordinates, &coords)); 1359 for (i = 0; i < nCoords; i += dimEmbed) { 1360 PetscInt j; 1361 1362 PetscScalar *coord = &coords[i]; 1363 for (j = 0; j < dimEmbed; j++) { 1364 coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]); 1365 if (dim == 2 && cells[1] == 1 && j == 0 && skew) { 1366 if (cells[0] == 2 && i == 8) { 1367 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */ 1368 } else if (cells[0] == 3) { 1369 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.; 1370 else if (i == 4) coord[j] = mod->bounds[1] / 2.; 1371 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.; 1372 } 1373 } 1374 } 1375 } 1376 PetscCall(VecRestoreArray(coordinates, &coords)); 1377 PetscCall(DMSetCoordinatesLocal(dm, coordinates)); 1378 } 1379 } 1380 } 1381 PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view")); 1382 PetscCall(DMGetDimension(dm, &dim)); 1383 1384 /* set up BCs, functions, tags */ 1385 PetscCall(DMCreateLabel(dm, "Face Sets")); 1386 mod->errorIndicator = ErrorIndicator_Simple; 1387 1388 { 1389 DM gdm; 1390 1391 PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm)); 1392 PetscCall(DMDestroy(&dm)); 1393 dm = gdm; 1394 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 1395 } 1396 1397 PetscCall(PetscFVCreate(comm, &fvm)); 1398 PetscCall(PetscFVSetFromOptions(fvm)); 1399 PetscCall(PetscFVSetNumComponents(fvm, phys->dof)); 1400 PetscCall(PetscFVSetSpatialDimension(fvm, dim)); 1401 PetscCall(PetscObjectSetName((PetscObject)fvm, "")); 1402 { 1403 PetscInt f, dof; 1404 for (f = 0, dof = 0; f < phys->nfields; f++) { 1405 PetscInt newDof = phys->field_desc[f].dof; 1406 1407 if (newDof == 1) { 1408 PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name)); 1409 } else { 1410 PetscInt j; 1411 1412 for (j = 0; j < newDof; j++) { 1413 char compName[256] = "Unknown"; 1414 1415 PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j)); 1416 PetscCall(PetscFVSetComponentName(fvm, dof + j, compName)); 1417 } 1418 } 1419 dof += newDof; 1420 } 1421 } 1422 /* FV is now structured with one field having all physics as components */ 1423 PetscCall(DMAddField(dm, NULL, (PetscObject)fvm)); 1424 PetscCall(DMCreateDS(dm)); 1425 PetscCall(DMGetDS(dm, &prob)); 1426 PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann)); 1427 PetscCall(PetscDSSetContext(prob, 0, user->model->physics)); 1428 PetscCall((*mod->setupbc)(dm, prob, phys)); 1429 PetscCall(PetscDSSetFromOptions(prob)); 1430 { 1431 char convType[256]; 1432 PetscBool flg; 1433 1434 PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX"); 1435 PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg)); 1436 PetscOptionsEnd(); 1437 if (flg) { 1438 DM dmConv; 1439 1440 PetscCall(DMConvert(dm, convType, &dmConv)); 1441 if (dmConv) { 1442 PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view")); 1443 PetscCall(DMDestroy(&dm)); 1444 dm = dmConv; 1445 PetscCall(DMSetFromOptions(dm)); 1446 } 1447 } 1448 } 1449 #ifdef PETSC_HAVE_LIBCEED 1450 { 1451 PetscBool useCeed; 1452 PetscCall(DMPlexGetUseCeed(dm, &useCeed)); 1453 if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics)); 1454 } 1455 #endif 1456 1457 PetscCall(initializeTS(dm, user, &ts)); 1458 1459 PetscCall(DMCreateGlobalVector(dm, &X)); 1460 PetscCall(PetscObjectSetName((PetscObject)X, "solution")); 1461 PetscCall(SetInitialCondition(dm, X, user)); 1462 if (useAMR) { 1463 PetscInt adaptIter; 1464 1465 /* use no limiting when reconstructing gradients for adaptivity */ 1466 PetscCall(PetscFVGetLimiter(fvm, &limiter)); 1467 PetscCall(PetscObjectReference((PetscObject)limiter)); 1468 PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter)); 1469 PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE)); 1470 1471 /* Refinement context */ 1472 tctx.fvm = fvm; 1473 tctx.refineTag = refineTag; 1474 tctx.coarsenTag = coarsenTag; 1475 tctx.adaptedDM = NULL; 1476 tctx.user = user; 1477 tctx.noneLimiter = noneLimiter; 1478 tctx.limiter = limiter; 1479 tctx.cfl = cfl; 1480 1481 /* Do some initial refinement steps */ 1482 for (adaptIter = 0;; ++adaptIter) { 1483 PetscLogDouble bytes; 1484 PetscBool resize; 1485 1486 PetscCall(PetscMemoryGetCurrentUsage(&bytes)); 1487 PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes)); 1488 PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view")); 1489 PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view")); 1490 1491 PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx)); 1492 if (!resize) break; 1493 PetscCall(DMDestroy(&dm)); 1494 PetscCall(VecDestroy(&X)); 1495 dm = tctx.adaptedDM; 1496 tctx.adaptedDM = NULL; 1497 PetscCall(TSSetDM(ts, dm)); 1498 PetscCall(DMCreateGlobalVector(dm, &X)); 1499 PetscCall(PetscObjectSetName((PetscObject)X, "solution")); 1500 PetscCall(SetInitialCondition(dm, X, user)); 1501 } 1502 } 1503 1504 PetscCall(DMConvert(dm, DMPLEX, &plex)); 1505 if (vtkCellGeom) { 1506 DM dmCell; 1507 Vec cellgeom, partition; 1508 1509 PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL)); 1510 PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer)); 1511 PetscCall(VecView(cellgeom, viewer)); 1512 PetscCall(PetscViewerDestroy(&viewer)); 1513 PetscCall(CreatePartitionVec(dm, &dmCell, &partition)); 1514 PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer)); 1515 PetscCall(VecView(partition, viewer)); 1516 PetscCall(PetscViewerDestroy(&viewer)); 1517 PetscCall(VecDestroy(&partition)); 1518 PetscCall(DMDestroy(&dmCell)); 1519 } 1520 /* collect max maxspeed from all processes -- todo */ 1521 PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius)); 1522 PetscCall(DMDestroy(&plex)); 1523 PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts))); 1524 PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname); 1525 dt = cfl * minRadius / mod->maxspeed; 1526 PetscCall(TSSetTimeStep(ts, dt)); 1527 PetscCall(TSSetFromOptions(ts)); 1528 1529 /* When using adaptive mesh refinement 1530 specify callbacks to refine the solution 1531 and interpolate data from old to new mesh */ 1532 if (useAMR) { PetscCall(TSSetResize(ts, adaptToleranceFVMSetUp, Transfer, &tctx)); } 1533 PetscCall(TSSetSolution(ts, X)); 1534 PetscCall(VecDestroy(&X)); 1535 PetscCall(TSSolve(ts, NULL)); 1536 PetscCall(TSGetSolveTime(ts, &ftime)); 1537 PetscCall(TSGetStepNumber(ts, &nsteps)); 1538 PetscCall(TSGetConvergedReason(ts, &reason)); 1539 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps)); 1540 PetscCall(TSDestroy(&ts)); 1541 1542 PetscCall(VecTaggerDestroy(&refineTag)); 1543 PetscCall(VecTaggerDestroy(&coarsenTag)); 1544 PetscCall(PetscFunctionListDestroy(&PhysicsList)); 1545 PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW)); 1546 PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler)); 1547 PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry)); 1548 PetscCall(PetscFree(user->model->functionalMonitored)); 1549 PetscCall(PetscFree(user->model->functionalCall)); 1550 PetscCall(PetscFree(user->model->physics->data)); 1551 PetscCall(PetscFree(user->model->physics)); 1552 PetscCall(PetscFree(user->model)); 1553 PetscCall(PetscFree(user)); 1554 PetscCall(PetscLimiterDestroy(&limiter)); 1555 PetscCall(PetscLimiterDestroy(&noneLimiter)); 1556 PetscCall(PetscFVDestroy(&fvm)); 1557 PetscCall(DMDestroy(&dm)); 1558 PetscCall(PetscFinalize()); 1559 return 0; 1560 } 1561 1562 /* Subroutine to set up the initial conditions for the */ 1563 /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */ 1564 /* ----------------------------------------------------------------------- */ 1565 int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3]) 1566 { 1567 int j, k; 1568 /* Wc=matmul(lv,Ueq) 3 vars */ 1569 for (k = 0; k < 3; ++k) { 1570 wc[k] = 0.; 1571 for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j]; 1572 } 1573 return 0; 1574 } 1575 /* ----------------------------------------------------------------------- */ 1576 int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3]) 1577 { 1578 int k, j; 1579 /* V=matmul(rv,WC) 3 vars */ 1580 for (k = 0; k < 3; ++k) { 1581 v[k] = 0.; 1582 for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j]; 1583 } 1584 return 0; 1585 } 1586 /* ---------------------------------------------------------------------- */ 1587 int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma) 1588 { 1589 int j, k; 1590 PetscReal rho, csnd, p0; 1591 /* PetscScalar u; */ 1592 1593 for (k = 0; k < 3; ++k) 1594 for (j = 0; j < 3; ++j) { 1595 lv[k][j] = 0.; 1596 rv[k][j] = 0.; 1597 } 1598 rho = ueq[0]; 1599 /* u = ueq[1]; */ 1600 p0 = ueq[2]; 1601 csnd = PetscSqrtReal(gamma * p0 / rho); 1602 lv[0][1] = rho * .5; 1603 lv[0][2] = -.5 / csnd; 1604 lv[1][0] = csnd; 1605 lv[1][2] = -1. / csnd; 1606 lv[2][1] = rho * .5; 1607 lv[2][2] = .5 / csnd; 1608 rv[0][0] = -1. / csnd; 1609 rv[1][0] = 1. / rho; 1610 rv[2][0] = -csnd; 1611 rv[0][1] = 1. / csnd; 1612 rv[0][2] = 1. / csnd; 1613 rv[1][2] = 1. / rho; 1614 rv[2][2] = csnd; 1615 return 0; 1616 } 1617 1618 int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx) 1619 { 1620 PetscReal p0, u0, wcp[3], wc[3]; 1621 PetscReal lv[3][3]; 1622 PetscReal vp[3]; 1623 PetscReal rv[3][3]; 1624 PetscReal eps, ueq[3], rho0, twopi; 1625 1626 /* Function Body */ 1627 twopi = 2. * PETSC_PI; 1628 eps = 1e-4; /* perturbation */ 1629 rho0 = 1e3; /* density of water */ 1630 p0 = 101325.; /* init pressure of 1 atm (?) */ 1631 u0 = 0.; 1632 ueq[0] = rho0; 1633 ueq[1] = u0; 1634 ueq[2] = p0; 1635 /* Project initial state to characteristic variables */ 1636 eigenvectors(rv, lv, ueq, gamma); 1637 projecteqstate(wc, ueq, lv); 1638 wcp[0] = wc[0]; 1639 wcp[1] = wc[1]; 1640 wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx); 1641 projecttoprim(vp, wcp, rv); 1642 ux->r = vp[0]; /* density */ 1643 ux->ru[0] = vp[0] * vp[1]; /* x momentum */ 1644 ux->ru[1] = 0.; 1645 #if defined DIM > 2 1646 if (dim > 2) ux->ru[2] = 0.; 1647 #endif 1648 /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */ 1649 ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1]; 1650 return 0; 1651 } 1652 1653 /*TEST 1654 1655 testset: 1656 args: -dm_plex_adj_cone -dm_plex_adj_closure 0 1657 1658 test: 1659 suffix: adv_2d_tri_0 1660 requires: triangle 1661 TODO: how did this ever get in main when there is no support for this 1662 args: -ufv_vtk_interval 0 -simplex -dm_refine 3 -dm_plex_faces 1,1 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3 1663 1664 test: 1665 suffix: adv_2d_tri_1 1666 requires: triangle 1667 TODO: how did this ever get in main when there is no support for this 1668 args: -ufv_vtk_interval 0 -simplex -dm_refine 5 -dm_plex_faces 1,1 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 1669 1670 test: 1671 suffix: tut_1 1672 requires: exodusii 1673 nsize: 1 1674 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo - 1675 1676 test: 1677 suffix: tut_2 1678 requires: exodusii 1679 nsize: 1 1680 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw 1681 1682 test: 1683 suffix: tut_3 1684 requires: exodusii 1685 nsize: 4 1686 args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -monitor Error -advect_sol_type bump -petscfv_type leastsquares -petsclimiter_type sin 1687 1688 test: 1689 suffix: tut_4 1690 requires: exodusii 1691 nsize: 4 1692 args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -physics sw -monitor Height,Energy -petscfv_type leastsquares -petsclimiter_type minmod 1693 1694 testset: 1695 args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 1696 1697 # 2D Advection 0-10 1698 test: 1699 suffix: 0 1700 requires: exodusii 1701 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo 1702 1703 test: 1704 suffix: 1 1705 requires: exodusii 1706 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo 1707 1708 test: 1709 suffix: 2 1710 requires: exodusii 1711 nsize: 2 1712 args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo 1713 1714 test: 1715 suffix: 3 1716 requires: exodusii 1717 nsize: 2 1718 args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo 1719 1720 test: 1721 suffix: 4 1722 requires: exodusii 1723 nsize: 4 1724 args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple 1725 1726 test: 1727 suffix: 5 1728 requires: exodusii 1729 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1 1730 1731 test: 1732 suffix: 7 1733 requires: exodusii 1734 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1 1735 1736 test: 1737 suffix: 8 1738 requires: exodusii 1739 nsize: 2 1740 args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1 1741 1742 test: 1743 suffix: 9 1744 requires: exodusii 1745 nsize: 8 1746 args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1 1747 1748 test: 1749 suffix: 10 1750 requires: exodusii 1751 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo 1752 1753 # 2D Shallow water 1754 testset: 1755 args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0 1756 1757 test: 1758 suffix: sw_0 1759 requires: exodusii 1760 args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \ 1761 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \ 1762 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1763 -monitor height,energy 1764 1765 test: 1766 suffix: sw_ceed 1767 requires: exodusii libceed 1768 args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \ 1769 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \ 1770 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1771 -monitor height,energy 1772 1773 test: 1774 suffix: sw_ceed_small 1775 requires: exodusii libceed 1776 args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \ 1777 -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_lower 0,1 -dm_plex_box_upper 6.28,3 -dm_plex_box_faces 8,2 \ 1778 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1779 -monitor height,energy 1780 1781 test: 1782 suffix: sw_1 1783 nsize: 2 1784 args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \ 1785 -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_faces 24,12 -dm_plex_box_lower 0,1 -dm_plex_box_upper 1,3 -dm_distribute_overlap 1 \ 1786 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1787 -monitor height,energy 1788 1789 test: 1790 suffix: sw_hll 1791 args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \ 1792 -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \ 1793 -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1794 -monitor height,energy 1795 1796 # 2D Euler 1797 testset: 1798 args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \ 1799 -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy 1800 1801 test: 1802 suffix: euler_0 1803 requires: exodusii !complex 1804 args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \ 1805 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \ 1806 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 1807 1808 test: 1809 suffix: euler_ceed 1810 requires: exodusii libceed 1811 args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \ 1812 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \ 1813 -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 1814 1815 testset: 1816 args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1 1817 1818 # 2D Advection: p4est 1819 test: 1820 suffix: p4est_advec_2d 1821 requires: p4est 1822 args: -ufv_vtk_interval 0 -dm_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 2 -dm_p4est_refine_pattern hash -dm_forest_maximum_refinement 5 1823 1824 # Advection in a box 1825 test: 1826 suffix: adv_2d_quad_0 1827 args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3 1828 1829 test: 1830 suffix: adv_2d_quad_1 1831 args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 1832 timeoutfactor: 3 1833 1834 test: 1835 suffix: adv_2d_quad_p4est_0 1836 requires: p4est 1837 args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3 1838 1839 test: 1840 suffix: adv_2d_quad_p4est_1 1841 requires: p4est 1842 args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 1843 timeoutfactor: 3 1844 1845 test: # broken for quad precision 1846 suffix: adv_2d_quad_p4est_adapt_0 1847 requires: p4est !__float128 1848 args: -ufv_vtk_interval 0 -dm_refine 3 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 -ufv_use_amr -refine_vec_tagger_box 0.005,inf -coarsen_vec_tagger_box 0,1.e-5 -petscfv_type leastsquares -ts_max_time 0.01 1849 timeoutfactor: 3 1850 1851 test: 1852 suffix: adv_0 1853 requires: exodusii 1854 args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201 1855 1856 test: 1857 suffix: shock_0 1858 requires: p4est !single !complex 1859 args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \ 1860 -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \ 1861 -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \ 1862 -bc_wall 1,2,3,4 -physics euler -eu_type iv_shock -ufv_cfl 10 -eu_alpha 60. -eu_gamma 1.4 -eu_amach 2.02 -eu_rho2 3. \ 1863 -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \ 1864 -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \ 1865 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy 1866 timeoutfactor: 3 1867 1868 # Test GLVis visualization of PetscFV fields 1869 test: 1870 suffix: glvis_adv_2d_tet 1871 args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \ 1872 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \ 1873 -ts_monitor_solution glvis: -ts_max_steps 0 1874 1875 test: 1876 suffix: glvis_adv_2d_quad 1877 args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \ 1878 -dm_refine 5 -dm_plex_separate_marker \ 1879 -ts_monitor_solution glvis: -ts_max_steps 0 1880 1881 TEST*/ 1882