1 static char help[] = "Hybrid Finite Element-Finite Volume Example.\n"; 2 /*F 3 Here we are advecting a passive tracer in a harmonic velocity field, defined by 4 a forcing function $f$: 5 \begin{align} 6 -\Delta \mathbf{u} + f &= 0 \\ 7 \frac{\partial\phi}{\partial t} + \nabla\cdot \phi \mathbf{u} &= 0 8 \end{align} 9 F*/ 10 11 #include <petscdmplex.h> 12 #include <petscds.h> 13 #include <petscts.h> 14 15 #include <petsc/private/dmpleximpl.h> /* For DotD */ 16 17 typedef enum {VEL_ZERO, VEL_CONSTANT, VEL_HARMONIC, VEL_SHEAR} VelocityDistribution; 18 19 typedef enum {ZERO, CONSTANT, GAUSSIAN, TILTED, DELTA} PorosityDistribution; 20 21 static PetscErrorCode constant_u_2d(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 22 23 /* 24 FunctionalFunc - Calculates the value of a functional of the solution at a point 25 26 Input Parameters: 27 + dm - The DM 28 . time - The TS time 29 . x - The coordinates of the evaluation point 30 . u - The field values at point x 31 - ctx - A user context, or NULL 32 33 Output Parameter: 34 . f - The value of the functional at point x 35 36 */ 37 typedef PetscErrorCode (*FunctionalFunc)(DM, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *); 38 39 typedef struct _n_Functional *Functional; 40 struct _n_Functional { 41 char *name; 42 FunctionalFunc func; 43 void *ctx; 44 PetscInt offset; 45 Functional next; 46 }; 47 48 typedef struct { 49 /* Problem definition */ 50 PetscBool useFV; /* Use a finite volume scheme for advection */ 51 PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 52 VelocityDistribution velocityDist; 53 PorosityDistribution porosityDist; 54 PetscReal inflowState; 55 PetscReal source[3]; 56 /* Monitoring */ 57 PetscInt numMonitorFuncs, maxMonitorFunc; 58 Functional *monitorFuncs; 59 PetscInt errorFunctional; 60 Functional functionalRegistry; 61 } AppCtx; 62 63 static AppCtx *globalUser; 64 65 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 66 { 67 const char *velocityDist[4] = {"zero", "constant", "harmonic", "shear"}; 68 const char *porosityDist[5] = {"zero", "constant", "gaussian", "tilted", "delta"}; 69 PetscInt vd, pd, d; 70 PetscBool flg; 71 PetscErrorCode ierr; 72 73 PetscFunctionBeginUser; 74 options->useFV = PETSC_FALSE; 75 options->velocityDist = VEL_HARMONIC; 76 options->porosityDist = ZERO; 77 options->inflowState = -2.0; 78 options->numMonitorFuncs = 0; 79 options->source[0] = 0.5; 80 options->source[1] = 0.5; 81 options->source[2] = 0.5; 82 83 ierr = PetscOptionsBegin(comm, "", "Magma Dynamics Options", "DMPLEX");PetscCall(ierr); 84 PetscCall(PetscOptionsBool("-use_fv", "Use the finite volume method for advection", "ex18.c", options->useFV, &options->useFV, NULL)); 85 vd = options->velocityDist; 86 PetscCall(PetscOptionsEList("-velocity_dist","Velocity distribution type","ex18.c",velocityDist,4,velocityDist[options->velocityDist],&vd,NULL)); 87 options->velocityDist = (VelocityDistribution) vd; 88 pd = options->porosityDist; 89 PetscCall(PetscOptionsEList("-porosity_dist","Initial porosity distribution type","ex18.c",porosityDist,5,porosityDist[options->porosityDist],&pd,NULL)); 90 options->porosityDist = (PorosityDistribution) pd; 91 PetscCall(PetscOptionsReal("-inflow_state", "The inflow state", "ex18.c", options->inflowState, &options->inflowState, NULL)); 92 d = 2; 93 PetscCall(PetscOptionsRealArray("-source_loc", "The source location", "ex18.c", options->source, &d, &flg)); 94 PetscCheck(!flg || d == 2,PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must give dim coordinates for the source location, not %d", d); 95 ierr = PetscOptionsEnd();PetscCall(ierr); 96 97 PetscFunctionReturn(0); 98 } 99 100 static PetscErrorCode ProcessMonitorOptions(MPI_Comm comm, AppCtx *options) 101 { 102 Functional func; 103 char *names[256]; 104 PetscInt f; 105 PetscErrorCode ierr; 106 107 PetscFunctionBeginUser; 108 ierr = PetscOptionsBegin(comm, "", "Simulation Monitor Options", "DMPLEX");PetscCall(ierr); 109 options->numMonitorFuncs = PETSC_STATIC_ARRAY_LENGTH(names); 110 PetscCall(PetscOptionsStringArray("-monitor", "List of functionals to monitor", "", names, &options->numMonitorFuncs, NULL)); 111 PetscCall(PetscMalloc1(options->numMonitorFuncs, &options->monitorFuncs)); 112 for (f = 0; f < options->numMonitorFuncs; ++f) { 113 for (func = options->functionalRegistry; func; func = func->next) { 114 PetscBool match; 115 116 PetscCall(PetscStrcasecmp(names[f], func->name, &match)); 117 if (match) break; 118 } 119 PetscCheck(func,comm, PETSC_ERR_USER, "No known functional '%s'", names[f]); 120 options->monitorFuncs[f] = func; 121 /* Jed inserts a de-duplication of functionals here */ 122 PetscCall(PetscFree(names[f])); 123 } 124 /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */ 125 options->maxMonitorFunc = -1; 126 for (func = options->functionalRegistry; func; func = func->next) { 127 for (f = 0; f < options->numMonitorFuncs; ++f) { 128 Functional call = options->monitorFuncs[f]; 129 130 if (func->func == call->func && func->ctx == call->ctx) options->maxMonitorFunc = PetscMax(options->maxMonitorFunc, func->offset); 131 } 132 } 133 ierr = PetscOptionsEnd();PetscCall(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 static PetscErrorCode FunctionalRegister(Functional *functionalRegistry, const char name[], PetscInt *offset, FunctionalFunc func, void *ctx) 138 { 139 Functional *ptr, f; 140 PetscInt lastoffset = -1; 141 142 PetscFunctionBeginUser; 143 for (ptr = functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset; 144 PetscCall(PetscNew(&f)); 145 PetscCall(PetscStrallocpy(name, &f->name)); 146 f->offset = lastoffset + 1; 147 f->func = func; 148 f->ctx = ctx; 149 f->next = NULL; 150 *ptr = f; 151 *offset = f->offset; 152 PetscFunctionReturn(0); 153 } 154 155 static PetscErrorCode FunctionalDestroy(Functional *link) 156 { 157 Functional next, l; 158 159 PetscFunctionBeginUser; 160 if (!link) PetscFunctionReturn(0); 161 l = *link; 162 *link = NULL; 163 for (; l; l=next) { 164 next = l->next; 165 PetscCall(PetscFree(l->name)); 166 PetscCall(PetscFree(l)); 167 } 168 PetscFunctionReturn(0); 169 } 170 171 static void f0_zero_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 172 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 173 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 174 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 175 { 176 PetscInt comp; 177 for (comp = 0; comp < dim; ++comp) f0[comp] = u[comp]; 178 } 179 180 static void f0_constant_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 181 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 182 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 183 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 184 { 185 PetscScalar wind[3] = {0.0, 0.0, 0.0}; 186 PetscInt comp; 187 188 constant_u_2d(dim, t, x, Nf, wind, NULL); 189 for (comp = 0; comp < dim && comp < 3; ++comp) f0[comp] = u[comp] - wind[comp]; 190 } 191 192 static void f1_constant_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 193 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 194 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 195 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 196 { 197 PetscInt comp; 198 for (comp = 0; comp < dim*dim; ++comp) f1[comp] = 0.0; 199 } 200 201 static void g0_constant_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 202 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 203 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 204 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 205 { 206 PetscInt d; 207 for (d = 0; d < dim; ++d) g0[d*dim+d] = 1.0; 208 } 209 210 static void g0_constant_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 211 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 212 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 213 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 214 { 215 g0[0] = 1.0; 216 } 217 218 static void f0_lap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 219 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 220 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 221 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 222 { 223 PetscInt comp; 224 for (comp = 0; comp < dim; ++comp) f0[comp] = 4.0; 225 } 226 227 static void f1_lap_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 228 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 229 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 230 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 231 { 232 PetscInt comp, d; 233 for (comp = 0; comp < dim; ++comp) { 234 for (d = 0; d < dim; ++d) { 235 f1[comp*dim+d] = u_x[comp*dim+d]; 236 } 237 } 238 } 239 240 static void f0_lap_periodic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 241 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 242 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 243 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 244 { 245 f0[0] = -PetscSinReal(2.0*PETSC_PI*x[0]); 246 f0[1] = 2.0*PETSC_PI*x[1]*PetscCosReal(2.0*PETSC_PI*x[0]); 247 } 248 249 static void f0_lap_doubly_periodic_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, 250 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 251 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 252 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 253 { 254 f0[0] = -2.0*PetscSinReal(2.0*PETSC_PI*x[0])*PetscCosReal(2.0*PETSC_PI*x[1]); 255 f0[1] = 2.0*PetscSinReal(2.0*PETSC_PI*x[1])*PetscCosReal(2.0*PETSC_PI*x[0]); 256 } 257 258 void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 259 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 260 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 261 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 262 { 263 const PetscInt Ncomp = dim; 264 PetscInt compI, d; 265 266 for (compI = 0; compI < Ncomp; ++compI) { 267 for (d = 0; d < dim; ++d) { 268 g3[((compI*Ncomp+compI)*dim+d)*dim+d] = 1.0; 269 } 270 } 271 } 272 273 /* \frac{\partial\phi}{\partial t} + \nabla\phi \cdot \mathbf{u} + \phi \nabla \cdot \mathbf{u} = 0 */ 274 static void f0_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, 275 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 276 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 277 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 278 { 279 PetscInt d; 280 f0[0] = u_t[dim]; 281 for (d = 0; d < dim; ++d) f0[0] += u[dim]*u_x[d*dim+d] + u_x[dim*dim+d]*u[d]; 282 } 283 284 static void f1_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, 285 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 286 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 287 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 288 { 289 PetscInt d; 290 for (d = 0; d < dim; ++d) f1[0] = 0.0; 291 } 292 293 void g0_adv_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 294 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 295 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 296 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 297 { 298 PetscInt d; 299 g0[0] = u_tShift; 300 for (d = 0; d < dim; ++d) g0[0] += u_x[d*dim+d]; 301 } 302 303 void g1_adv_pp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 304 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 305 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 306 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 307 { 308 PetscInt d; 309 for (d = 0; d < dim; ++d) g1[d] = u[d]; 310 } 311 312 void g0_adv_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 313 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 314 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 315 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 316 { 317 PetscInt d; 318 for (d = 0; d < dim; ++d) g0[0] += u_x[dim*dim+d]; 319 } 320 321 void g1_adv_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, 322 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 323 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 324 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]) 325 { 326 PetscInt d; 327 for (d = 0; d < dim; ++d) g1[d*dim+d] = u[dim]; 328 } 329 330 static void riemann_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx) 331 { 332 PetscReal wind[3] = {0.0, 1.0, 0.0}; 333 PetscReal wn = DMPlex_DotRealD_Internal(PetscMin(dim,3), wind, n); 334 335 flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn; 336 } 337 338 static void riemann_coupled_advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx) 339 { 340 PetscReal wn = DMPlex_DotD_Internal(dim, uL, n); 341 342 #if 1 343 flux[0] = (wn > 0 ? uL[dim] : uR[dim]) * wn; 344 #else 345 /* if (fabs(uL[0] - wind[0]) > 1.0e-7 || fabs(uL[1] - wind[1]) > 1.0e-7) PetscPrintf(PETSC_COMM_SELF, "wind (%g, %g) uL (%g, %g) uR (%g, %g)\n", wind[0], wind[1], uL[0], uL[1], uR[0], uR[1]); */ 346 /* Smear it out */ 347 flux[0] = 0.5*((uL[dim] + uR[dim]) + (uL[dim] - uR[dim])*tanh(1.0e5*wn)) * wn; 348 #endif 349 } 350 351 static PetscErrorCode zero_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 352 { 353 u[0] = 0.0; 354 u[1] = 0.0; 355 return 0; 356 } 357 358 static PetscErrorCode constant_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 359 { 360 u[0] = 0.0; 361 u[1] = 1.0; 362 return 0; 363 } 364 365 /* Coordinates of the point which was at x at t = 0 */ 366 static PetscErrorCode constant_x_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 367 { 368 const PetscReal t = *((PetscReal *) ctx); 369 u[0] = x[0]; 370 u[1] = x[1] + t; 371 #if 0 372 PetscCall(DMLocalizeCoordinate(globalUser->dm, u, PETSC_FALSE, u)); 373 #else 374 u[1] = u[1] - (int) PetscRealPart(u[1]); 375 #endif 376 return 0; 377 } 378 379 /* 380 In 2D we use the exact solution: 381 382 u = x^2 + y^2 383 v = 2 x^2 - 2xy 384 phi = h(x + y + (u + v) t) 385 f_x = f_y = 4 386 387 so that 388 389 -\Delta u + f = <-4, -4> + <4, 4> = 0 390 {\partial\phi}{\partial t} - \nabla\cdot \phi u = 0 391 h_t(x + y + (u + v) t) - u . grad phi - phi div u 392 = u h' + v h' - u h_x - v h_y 393 = 0 394 395 We will conserve phi since 396 397 \nabla \cdot u = 2x - 2x = 0 398 399 Also try h((x + ut)^2 + (y + vt)^2), so that 400 401 h_t((x + ut)^2 + (y + vt)^2) - u . grad phi - phi div u 402 = 2 h' (u (x + ut) + v (y + vt)) - u h_x - v h_y 403 = 2 h' (u (x + ut) + v (y + vt)) - u h' 2 (x + u t) - v h' 2 (y + vt) 404 = 2 h' (u (x + ut) + v (y + vt) - u (x + u t) - v (y + vt)) 405 = 0 406 407 */ 408 static PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 409 { 410 u[0] = x[0]*x[0] + x[1]*x[1]; 411 u[1] = 2.0*x[0]*x[0] - 2.0*x[0]*x[1]; 412 return 0; 413 } 414 415 /* 416 In 2D we use the exact, periodic solution: 417 418 u = sin(2 pi x)/4 pi^2 419 v = -y cos(2 pi x)/2 pi 420 phi = h(x + y + (u + v) t) 421 f_x = -sin(2 pi x) 422 f_y = 2 pi y cos(2 pi x) 423 424 so that 425 426 -\Delta u + f = <sin(2pi x), -2pi y cos(2pi x)> + <-sin(2pi x), 2pi y cos(2pi x)> = 0 427 428 We will conserve phi since 429 430 \nabla \cdot u = cos(2pi x)/2pi - cos(2pi x)/2pi = 0 431 */ 432 static PetscErrorCode periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 433 { 434 u[0] = PetscSinReal(2.0*PETSC_PI*x[0])/PetscSqr(2.0*PETSC_PI); 435 u[1] = -x[1]*PetscCosReal(2.0*PETSC_PI*x[0])/(2.0*PETSC_PI); 436 return 0; 437 } 438 439 /* 440 In 2D we use the exact, doubly periodic solution: 441 442 u = sin(2 pi x) cos(2 pi y)/4 pi^2 443 v = -sin(2 pi y) cos(2 pi x)/4 pi^2 444 phi = h(x + y + (u + v) t) 445 f_x = -2sin(2 pi x) cos(2 pi y) 446 f_y = 2sin(2 pi y) cos(2 pi x) 447 448 so that 449 450 -\Delta u + f = <2 sin(2pi x) cos(2pi y), -2 sin(2pi y) cos(2pi x)> + <-2 sin(2pi x) cos(2pi y), 2 sin(2pi y) cos(2pi x)> = 0 451 452 We will conserve phi since 453 454 \nabla \cdot u = cos(2pi x) cos(2pi y)/2pi - cos(2pi y) cos(2pi x)/2pi = 0 455 */ 456 static PetscErrorCode doubly_periodic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 457 { 458 u[0] = PetscSinReal(2.0*PETSC_PI*x[0])*PetscCosReal(2.0*PETSC_PI*x[1])/PetscSqr(2.0*PETSC_PI); 459 u[1] = -PetscSinReal(2.0*PETSC_PI*x[1])*PetscCosReal(2.0*PETSC_PI*x[0])/PetscSqr(2.0*PETSC_PI); 460 return 0; 461 } 462 463 static PetscErrorCode shear_bc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 464 { 465 u[0] = x[1] - 0.5; 466 u[1] = 0.0; 467 return 0; 468 } 469 470 static PetscErrorCode initialVelocity(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 471 { 472 PetscInt d; 473 for (d = 0; d < dim; ++d) u[d] = 0.0; 474 return 0; 475 } 476 477 static PetscErrorCode zero_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 478 { 479 u[0] = 0.0; 480 return 0; 481 } 482 483 static PetscErrorCode constant_phi(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 484 { 485 u[0] = 1.0; 486 return 0; 487 } 488 489 static PetscErrorCode delta_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 490 { 491 PetscReal x0[2]; 492 PetscScalar xn[2]; 493 494 x0[0] = globalUser->source[0]; 495 x0[1] = globalUser->source[1]; 496 constant_x_2d(dim, time, x0, Nf, xn, ctx); 497 { 498 const PetscReal xi = x[0] - PetscRealPart(xn[0]); 499 const PetscReal eta = x[1] - PetscRealPart(xn[1]); 500 const PetscReal r2 = xi*xi + eta*eta; 501 502 u[0] = r2 < 1.0e-7 ? 1.0 : 0.0; 503 } 504 return 0; 505 } 506 507 /* 508 Gaussian blob, initially centered on (0.5, 0.5) 509 510 xi = x(t) - x0, eta = y(t) - y0 511 512 where x(t), y(t) are the integral curves of v(t), 513 514 dx/dt . grad f = v . f 515 516 Check: constant v(t) = {v0, w0}, x(t) = {x0 + v0 t, y0 + w0 t} 517 518 v0 f_x + w0 f_y = v . f 519 */ 520 static PetscErrorCode gaussian_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 521 { 522 const PetscReal x0[2] = {0.5, 0.5}; 523 const PetscReal sigma = 1.0/6.0; 524 PetscScalar xn[2]; 525 526 constant_x_2d(dim, time, x0, Nf, xn, ctx); 527 { 528 /* const PetscReal xi = x[0] + (sin(2.0*PETSC_PI*x[0])/(4.0*PETSC_PI*PETSC_PI))*t - x0[0]; */ 529 /* const PetscReal eta = x[1] + (-x[1]*cos(2.0*PETSC_PI*x[0])/(2.0*PETSC_PI))*t - x0[1]; */ 530 const PetscReal xi = x[0] - PetscRealPart(xn[0]); 531 const PetscReal eta = x[1] - PetscRealPart(xn[1]); 532 const PetscReal r2 = xi*xi + eta*eta; 533 534 u[0] = PetscExpReal(-r2/(2.0*sigma*sigma))/(sigma*PetscSqrtReal(2.0*PETSC_PI)); 535 } 536 return 0; 537 } 538 539 static PetscErrorCode tilted_phi_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 540 { 541 PetscReal x0[3]; 542 const PetscReal wind[3] = {0.0, 1.0, 0.0}; 543 const PetscReal t = *((PetscReal *) ctx); 544 545 DMPlex_WaxpyD_Internal(2, -t, wind, x, x0); 546 if (x0[1] > 0) u[0] = 1.0*x[0] + 3.0*x[1]; 547 else u[0] = -2.0; /* Inflow state */ 548 return 0; 549 } 550 551 static PetscErrorCode tilted_phi_coupled_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) 552 { 553 PetscReal ur[3]; 554 PetscReal x0[3]; 555 const PetscReal t = *((PetscReal *) ctx); 556 557 ur[0] = PetscRealPart(u[0]); ur[1] = PetscRealPart(u[1]); ur[2] = PetscRealPart(u[2]); 558 DMPlex_WaxpyD_Internal(2, -t, ur, x, x0); 559 if (x0[1] > 0) u[0] = 1.0*x[0] + 3.0*x[1]; 560 else u[0] = -2.0; /* Inflow state */ 561 return 0; 562 } 563 564 static PetscErrorCode advect_inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 565 { 566 AppCtx *user = (AppCtx *) ctx; 567 568 PetscFunctionBeginUser; 569 xG[0] = user->inflowState; 570 PetscFunctionReturn(0); 571 } 572 573 static PetscErrorCode advect_outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx) 574 { 575 PetscFunctionBeginUser; 576 //xG[0] = xI[dim]; 577 xG[0] = xI[2]; 578 PetscFunctionReturn(0); 579 } 580 581 static PetscErrorCode ExactSolution(DM dm, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx) 582 { 583 AppCtx *user = (AppCtx *) ctx; 584 PetscInt dim; 585 586 PetscFunctionBeginUser; 587 PetscCall(DMGetDimension(dm, &dim)); 588 switch (user->porosityDist) { 589 case TILTED: 590 if (user->velocityDist == VEL_ZERO) tilted_phi_2d(dim, time, x, 2, u, (void *) &time); 591 else tilted_phi_coupled_2d(dim, time, x, 2, u, (void *) &time); 592 break; 593 case GAUSSIAN: 594 gaussian_phi_2d(dim, time, x, 2, u, (void *) &time); 595 break; 596 case DELTA: 597 delta_phi_2d(dim, time, x, 2, u, (void *) &time); 598 break; 599 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type"); 600 } 601 PetscFunctionReturn(0); 602 } 603 604 static PetscErrorCode Functional_Error(DM dm, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx) 605 { 606 AppCtx *user = (AppCtx *) ctx; 607 PetscScalar yexact[3]={0,0,0}; 608 609 PetscFunctionBeginUser; 610 PetscCall(ExactSolution(dm, time, x, yexact, ctx)); 611 f[user->errorFunctional] = PetscAbsScalar(y[0] - yexact[0]); 612 PetscFunctionReturn(0); 613 } 614 615 static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) 616 { 617 PetscFunctionBeginUser; 618 PetscCall(DMCreate(comm, dm)); 619 PetscCall(DMSetType(*dm, DMPLEX)); 620 #if 0 621 PetscBool periodic = user->bd[0] == DM_BOUNDARY_PERIODIC || user->bd[0] == DM_BOUNDARY_TWIST || user->bd[1] == DM_BOUNDARY_PERIODIC || user->bd[1] == DM_BOUNDARY_TWIST ? PETSC_TRUE : PETSC_FALSE; 622 const PetscReal L[3] = {1.0, 1.0, 1.0}; 623 PetscReal maxCell[3]; 624 PetscInt d; 625 626 if (periodic) {for (d = 0; d < 3; ++d) maxCell[d] = 1.1*(L[d]/cells[d]); PetscCall(DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, user->bd));} 627 #endif 628 PetscCall(DMSetFromOptions(*dm)); 629 PetscCall(DMViewFromOptions(*dm, NULL, "-orig_dm_view")); 630 PetscFunctionReturn(0); 631 } 632 633 static PetscErrorCode SetupBC(DM dm, AppCtx *user) 634 { 635 PetscErrorCode (*exactFuncs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx); 636 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 637 PetscDS prob; 638 DMLabel label; 639 PetscBool check; 640 PetscInt dim, n = 3; 641 const char *prefix; 642 643 PetscFunctionBeginUser; 644 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 645 PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *) bdt, &n, NULL)); 646 PetscCall(DMGetDimension(dm, &dim)); 647 /* Set initial guesses and exact solutions */ 648 switch (dim) { 649 case 2: 650 user->initialGuess[0] = initialVelocity; 651 switch(user->porosityDist) { 652 case ZERO: user->initialGuess[1] = zero_phi;break; 653 case CONSTANT: user->initialGuess[1] = constant_phi;break; 654 case GAUSSIAN: user->initialGuess[1] = gaussian_phi_2d;break; 655 case DELTA: user->initialGuess[1] = delta_phi_2d;break; 656 case TILTED: 657 if (user->velocityDist == VEL_ZERO) user->initialGuess[1] = tilted_phi_2d; 658 else user->initialGuess[1] = tilted_phi_coupled_2d; 659 break; 660 } 661 break; 662 default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %D not supported", dim); 663 } 664 exactFuncs[0] = user->initialGuess[0]; 665 exactFuncs[1] = user->initialGuess[1]; 666 switch (dim) { 667 case 2: 668 switch (user->velocityDist) { 669 case VEL_ZERO: 670 exactFuncs[0] = zero_u_2d; break; 671 case VEL_CONSTANT: 672 exactFuncs[0] = constant_u_2d; break; 673 case VEL_HARMONIC: 674 switch (bdt[0]) { 675 case DM_BOUNDARY_PERIODIC: 676 switch (bdt[1]) { 677 case DM_BOUNDARY_PERIODIC: 678 exactFuncs[0] = doubly_periodic_u_2d; break; 679 default: 680 exactFuncs[0] = periodic_u_2d; break; 681 } 682 break; 683 default: 684 exactFuncs[0] = quadratic_u_2d; break; 685 } 686 break; 687 case VEL_SHEAR: 688 exactFuncs[0] = shear_bc; break; 689 default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", dim); 690 } 691 break; 692 default: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dimension %D not supported", dim); 693 } 694 { 695 PetscBool isImplicit = PETSC_FALSE; 696 697 PetscCall(PetscOptionsHasName(NULL,"", "-use_implicit", &isImplicit)); 698 if (user->velocityDist == VEL_CONSTANT && !isImplicit) user->initialGuess[0] = exactFuncs[0]; 699 } 700 PetscCall(PetscOptionsHasName(NULL,NULL, "-dmts_check", &check)); 701 if (check) { 702 user->initialGuess[0] = exactFuncs[0]; 703 user->initialGuess[1] = exactFuncs[1]; 704 } 705 /* Set BC */ 706 PetscCall(DMGetDS(dm, &prob)); 707 PetscCall(DMGetLabel(dm, "marker", &label)); 708 PetscCall(PetscDSSetExactSolution(prob, 0, exactFuncs[0], user)); 709 PetscCall(PetscDSSetExactSolution(prob, 1, exactFuncs[1], user)); 710 if (label) { 711 const PetscInt id = 1; 712 713 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) exactFuncs[0], NULL, user, NULL)); 714 } 715 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 716 if (label && user->useFV) { 717 const PetscInt inflowids[] = {100,200,300}, outflowids[] = {101}; 718 719 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 1, 0, NULL, (void (*)(void)) advect_inflow, NULL, user, NULL)); 720 PetscCall(DMAddBoundary(dm, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 1, 0, NULL, (void (*)(void)) advect_outflow, NULL, user, NULL)); 721 } 722 PetscFunctionReturn(0); 723 } 724 725 static PetscErrorCode SetupProblem(DM dm, AppCtx *user) 726 { 727 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE}; 728 PetscDS prob; 729 PetscInt n = 3; 730 const char *prefix; 731 732 PetscFunctionBeginUser; 733 PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix)); 734 PetscCall(PetscOptionsGetEnumArray(NULL, prefix, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *) bdt, &n, NULL)); 735 PetscCall(DMGetDS(dm, &prob)); 736 switch (user->velocityDist) { 737 case VEL_ZERO: 738 PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_constant_u)); 739 break; 740 case VEL_CONSTANT: 741 PetscCall(PetscDSSetResidual(prob, 0, f0_constant_u, f1_constant_u)); 742 PetscCall(PetscDSSetJacobian(prob, 0, 0, g0_constant_uu, NULL, NULL, NULL)); 743 PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_constant_pp, NULL, NULL, NULL)); 744 break; 745 case VEL_HARMONIC: 746 switch (bdt[0]) { 747 case DM_BOUNDARY_PERIODIC: 748 switch (bdt[1]) { 749 case DM_BOUNDARY_PERIODIC: 750 PetscCall(PetscDSSetResidual(prob, 0, f0_lap_doubly_periodic_u, f1_lap_u)); 751 break; 752 default: 753 PetscCall(PetscDSSetResidual(prob, 0, f0_lap_periodic_u, f1_lap_u)); 754 break; 755 } 756 break; 757 default: 758 PetscCall(PetscDSSetResidual(prob, 0, f0_lap_u, f1_lap_u)); 759 break; 760 } 761 PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu)); 762 break; 763 case VEL_SHEAR: 764 PetscCall(PetscDSSetResidual(prob, 0, f0_zero_u, f1_lap_u)); 765 PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu)); 766 break; 767 } 768 PetscCall(PetscDSSetResidual(prob, 1, f0_advection, f1_advection)); 769 PetscCall(PetscDSSetJacobian(prob, 1, 1, g0_adv_pp, g1_adv_pp, NULL, NULL)); 770 PetscCall(PetscDSSetJacobian(prob, 1, 0, g0_adv_pu, g1_adv_pu, NULL, NULL)); 771 if (user->velocityDist == VEL_ZERO) PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_advection)); 772 else PetscCall(PetscDSSetRiemannSolver(prob, 1, riemann_coupled_advection)); 773 774 PetscCall(FunctionalRegister(&user->functionalRegistry, "Error", &user->errorFunctional, Functional_Error, user)); 775 PetscFunctionReturn(0); 776 } 777 778 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *user) 779 { 780 DM cdm = dm; 781 PetscQuadrature q; 782 PetscFE fe[2]; 783 PetscFV fv; 784 MPI_Comm comm; 785 PetscInt dim; 786 787 PetscFunctionBeginUser; 788 /* Create finite element */ 789 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 790 PetscCall(DMGetDimension(dm, &dim)); 791 PetscCall(PetscFECreateDefault(comm, dim, dim, PETSC_FALSE, "velocity_", PETSC_DEFAULT, &fe[0])); 792 PetscCall(PetscObjectSetName((PetscObject) fe[0], "velocity")); 793 PetscCall(PetscFECreateDefault(comm, dim, 1, PETSC_FALSE, "porosity_", PETSC_DEFAULT, &fe[1])); 794 PetscCall(PetscFECopyQuadrature(fe[0], fe[1])); 795 PetscCall(PetscObjectSetName((PetscObject) fe[1], "porosity")); 796 797 PetscCall(PetscFVCreate(PetscObjectComm((PetscObject) dm), &fv)); 798 PetscCall(PetscObjectSetName((PetscObject) fv, "porosity")); 799 PetscCall(PetscFVSetFromOptions(fv)); 800 PetscCall(PetscFVSetNumComponents(fv, 1)); 801 PetscCall(PetscFVSetSpatialDimension(fv, dim)); 802 PetscCall(PetscFEGetQuadrature(fe[0], &q)); 803 PetscCall(PetscFVSetQuadrature(fv, q)); 804 805 PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe[0])); 806 if (user->useFV) PetscCall(DMSetField(dm, 1, NULL, (PetscObject) fv)); 807 else PetscCall(DMSetField(dm, 1, NULL, (PetscObject) fe[1])); 808 PetscCall(DMCreateDS(dm)); 809 PetscCall(SetupProblem(dm, user)); 810 811 /* Set discretization and boundary conditions for each mesh */ 812 while (cdm) { 813 PetscCall(DMCopyDisc(dm, cdm)); 814 PetscCall(DMGetCoarseDM(cdm, &cdm)); 815 /* Coordinates were never localized for coarse meshes */ 816 if (cdm) PetscCall(DMLocalizeCoordinates(cdm)); 817 } 818 PetscCall(PetscFEDestroy(&fe[0])); 819 PetscCall(PetscFEDestroy(&fe[1])); 820 PetscCall(PetscFVDestroy(&fv)); 821 PetscFunctionReturn(0); 822 } 823 824 static PetscErrorCode CreateDM(MPI_Comm comm, AppCtx *user, DM *dm) 825 { 826 PetscFunctionBeginUser; 827 PetscCall(CreateMesh(comm, user, dm)); 828 /* Handle refinement, etc. */ 829 PetscCall(DMSetFromOptions(*dm)); 830 /* Construct ghost cells */ 831 if (user->useFV) { 832 DM gdm; 833 834 PetscCall(DMPlexConstructGhostCells(*dm, NULL, NULL, &gdm)); 835 PetscCall(DMDestroy(dm)); 836 *dm = gdm; 837 } 838 /* Localize coordinates */ 839 PetscCall(DMLocalizeCoordinates(*dm)); 840 PetscCall(PetscObjectSetName((PetscObject)(*dm),"Mesh")); 841 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 842 /* Setup problem */ 843 PetscCall(SetupDiscretization(*dm, user)); 844 /* Setup BC */ 845 PetscCall(SetupBC(*dm, user)); 846 PetscFunctionReturn(0); 847 } 848 849 static PetscErrorCode SetInitialConditionFVM(DM dm, Vec X, PetscInt field, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx) 850 { 851 PetscDS prob; 852 DM dmCell; 853 Vec cellgeom; 854 const PetscScalar *cgeom; 855 PetscScalar *x; 856 PetscInt dim, Nf, cStart, cEnd, c; 857 858 PetscFunctionBeginUser; 859 PetscCall(DMGetDS(dm, &prob)); 860 PetscCall(DMGetDimension(dm, &dim)); 861 PetscCall(PetscDSGetNumFields(prob, &Nf)); 862 PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL)); 863 PetscCall(VecGetDM(cellgeom, &dmCell)); 864 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 865 PetscCall(VecGetArrayRead(cellgeom, &cgeom)); 866 PetscCall(VecGetArray(X, &x)); 867 for (c = cStart; c < cEnd; ++c) { 868 PetscFVCellGeom *cg; 869 PetscScalar *xc; 870 871 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 872 PetscCall(DMPlexPointGlobalFieldRef(dm, c, field, x, &xc)); 873 if (xc) PetscCall((*func)(dim, 0.0, cg->centroid, Nf, xc, ctx)); 874 } 875 PetscCall(VecRestoreArrayRead(cellgeom, &cgeom)); 876 PetscCall(VecRestoreArray(X, &x)); 877 PetscFunctionReturn(0); 878 } 879 880 static PetscErrorCode MonitorFunctionals(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx) 881 { 882 AppCtx *user = (AppCtx *) ctx; 883 char *ftable = NULL; 884 DM dm; 885 PetscSection s; 886 Vec cellgeom; 887 const PetscScalar *x; 888 PetscScalar *a; 889 PetscReal *xnorms; 890 PetscInt pStart, pEnd, p, Nf, f; 891 892 PetscFunctionBeginUser; 893 PetscCall(VecViewFromOptions(X, (PetscObject) ts, "-view_solution")); 894 PetscCall(VecGetDM(X, &dm)); 895 PetscCall(DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL)); 896 PetscCall(DMGetLocalSection(dm, &s)); 897 PetscCall(PetscSectionGetNumFields(s, &Nf)); 898 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd)); 899 PetscCall(PetscCalloc1(Nf*2, &xnorms)); 900 PetscCall(VecGetArrayRead(X, &x)); 901 for (p = pStart; p < pEnd; ++p) { 902 for (f = 0; f < Nf; ++f) { 903 PetscInt dof, cdof, d; 904 905 PetscCall(PetscSectionGetFieldDof(s, p, f, &dof)); 906 PetscCall(PetscSectionGetFieldConstraintDof(s, p, f, &cdof)); 907 PetscCall(DMPlexPointGlobalFieldRead(dm, p, f, x, &a)); 908 /* TODO Use constrained indices here */ 909 for (d = 0; d < dof-cdof; ++d) xnorms[f*2+0] = PetscMax(xnorms[f*2+0], PetscAbsScalar(a[d])); 910 for (d = 0; d < dof-cdof; ++d) xnorms[f*2+1] += PetscAbsScalar(a[d]); 911 } 912 } 913 PetscCall(VecRestoreArrayRead(X, &x)); 914 if (stepnum >= 0) { /* No summary for final time */ 915 DM dmCell, *fdm; 916 Vec *fv; 917 const PetscScalar *cgeom; 918 PetscScalar **fx; 919 PetscReal *fmin, *fmax, *fint, *ftmp, t; 920 PetscInt cStart, cEnd, c, fcount, f, num; 921 922 size_t ftableused,ftablealloc; 923 924 /* Functionals have indices after registering, this is an upper bound */ 925 fcount = user->numMonitorFuncs; 926 PetscCall(PetscMalloc4(fcount,&fmin,fcount,&fmax,fcount,&fint,fcount,&ftmp)); 927 PetscCall(PetscMalloc3(fcount,&fdm,fcount,&fv,fcount,&fx)); 928 for (f = 0; f < fcount; ++f) { 929 PetscSection fs; 930 const char *name = user->monitorFuncs[f]->name; 931 932 fmin[f] = PETSC_MAX_REAL; 933 fmax[f] = PETSC_MIN_REAL; 934 fint[f] = 0; 935 /* Make monitor vecs */ 936 PetscCall(DMClone(dm, &fdm[f])); 937 PetscCall(DMGetOutputSequenceNumber(dm, &num, &t)); 938 PetscCall(DMSetOutputSequenceNumber(fdm[f], num, t)); 939 PetscCall(PetscSectionClone(s, &fs)); 940 PetscCall(PetscSectionSetFieldName(fs, 0, NULL)); 941 PetscCall(PetscSectionSetFieldName(fs, 1, name)); 942 PetscCall(DMSetLocalSection(fdm[f], fs)); 943 PetscCall(PetscSectionDestroy(&fs)); 944 PetscCall(DMGetGlobalVector(fdm[f], &fv[f])); 945 PetscCall(PetscObjectSetName((PetscObject) fv[f], name)); 946 PetscCall(VecGetArray(fv[f], &fx[f])); 947 } 948 PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd)); 949 PetscCall(VecGetDM(cellgeom, &dmCell)); 950 PetscCall(VecGetArrayRead(cellgeom, &cgeom)); 951 PetscCall(VecGetArrayRead(X, &x)); 952 for (c = cStart; c < cEnd; ++c) { 953 PetscFVCellGeom *cg; 954 PetscScalar *cx; 955 956 PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg)); 957 PetscCall(DMPlexPointGlobalFieldRead(dm, c, 1, x, &cx)); 958 if (!cx) continue; /* not a global cell */ 959 for (f = 0; f < user->numMonitorFuncs; ++f) { 960 Functional func = user->monitorFuncs[f]; 961 PetscScalar *fxc; 962 963 PetscCall(DMPlexPointGlobalFieldRef(dm, c, 1, fx[f], &fxc)); 964 /* I need to make it easier to get interpolated values here */ 965 PetscCall((*func->func)(dm, time, cg->centroid, cx, ftmp, func->ctx)); 966 fxc[0] = ftmp[user->monitorFuncs[f]->offset]; 967 } 968 for (f = 0; f < fcount; ++f) { 969 fmin[f] = PetscMin(fmin[f], ftmp[f]); 970 fmax[f] = PetscMax(fmax[f], ftmp[f]); 971 fint[f] += cg->volume * ftmp[f]; 972 } 973 } 974 PetscCall(VecRestoreArrayRead(cellgeom, &cgeom)); 975 PetscCall(VecRestoreArrayRead(X, &x)); 976 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts))); 977 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts))); 978 PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, fint, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts))); 979 /* Output functional data */ 980 ftablealloc = fcount * 100; 981 ftableused = 0; 982 PetscCall(PetscCalloc1(ftablealloc, &ftable)); 983 for (f = 0; f < user->numMonitorFuncs; ++f) { 984 Functional func = user->monitorFuncs[f]; 985 PetscInt id = func->offset; 986 char newline[] = "\n"; 987 char buffer[256], *p, *prefix; 988 size_t countused, len; 989 990 /* Create string with functional outputs */ 991 if (f % 3) { 992 PetscCall(PetscArraycpy(buffer, " ", 2)); 993 p = buffer + 2; 994 } else if (f) { 995 PetscCall(PetscArraycpy(buffer, newline, sizeof(newline)-1)); 996 p = buffer + sizeof(newline) - 1; 997 } else { 998 p = buffer; 999 } 1000 PetscCall(PetscSNPrintfCount(p, sizeof buffer-(p-buffer), "%12s [%12.6g,%12.6g] int %12.6g", &countused, func->name, (double) fmin[id], (double) fmax[id], (double) fint[id])); 1001 countused += p - buffer; 1002 /* reallocate */ 1003 if (countused > ftablealloc-ftableused-1) { 1004 char *ftablenew; 1005 1006 ftablealloc = 2*ftablealloc + countused; 1007 PetscCall(PetscMalloc1(ftablealloc, &ftablenew)); 1008 PetscCall(PetscArraycpy(ftablenew, ftable, ftableused)); 1009 PetscCall(PetscFree(ftable)); 1010 ftable = ftablenew; 1011 } 1012 PetscCall(PetscArraycpy(ftable+ftableused, buffer, countused)); 1013 ftableused += countused; 1014 ftable[ftableused] = 0; 1015 /* Output vecs */ 1016 PetscCall(VecRestoreArray(fv[f], &fx[f])); 1017 PetscCall(PetscStrlen(func->name, &len)); 1018 PetscCall(PetscMalloc1(len+2,&prefix)); 1019 PetscCall(PetscStrcpy(prefix, func->name)); 1020 PetscCall(PetscStrcat(prefix, "_")); 1021 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fv[f], prefix)); 1022 PetscCall(VecViewFromOptions(fv[f], NULL, "-vec_view")); 1023 PetscCall(PetscFree(prefix)); 1024 PetscCall(DMRestoreGlobalVector(fdm[f], &fv[f])); 1025 PetscCall(DMDestroy(&fdm[f])); 1026 } 1027 PetscCall(PetscFree4(fmin, fmax, fint, ftmp)); 1028 PetscCall(PetscFree3(fdm, fv, fx)); 1029 PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), "% 3D time %8.4g |x| (", stepnum, (double) time)); 1030 for (f = 0; f < Nf; ++f) { 1031 if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), ", ")); 1032 PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), "%8.4g", (double) xnorms[f*2+0])); 1033 } 1034 PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), ") |x|_1 (")); 1035 for (f = 0; f < Nf; ++f) { 1036 if (f > 0) PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), ", ")); 1037 PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), "%8.4g", (double) xnorms[f*2+1])); 1038 } 1039 PetscCall(PetscPrintf(PetscObjectComm((PetscObject) ts), ") %s\n", ftable ? ftable : "")); 1040 PetscCall(PetscFree(ftable)); 1041 } 1042 PetscCall(PetscFree(xnorms)); 1043 PetscFunctionReturn(0); 1044 } 1045 1046 int main(int argc, char **argv) 1047 { 1048 MPI_Comm comm; 1049 TS ts; 1050 DM dm; 1051 Vec u; 1052 AppCtx user; 1053 PetscReal t0, t = 0.0; 1054 void *ctxs[2]; 1055 1056 ctxs[0] = &t; 1057 ctxs[1] = &t; 1058 PetscCall(PetscInitialize(&argc, &argv, (char*) 0, help)); 1059 comm = PETSC_COMM_WORLD; 1060 user.functionalRegistry = NULL; 1061 globalUser = &user; 1062 PetscCall(ProcessOptions(comm, &user)); 1063 PetscCall(TSCreate(comm, &ts)); 1064 PetscCall(TSSetType(ts, TSBEULER)); 1065 PetscCall(CreateDM(comm, &user, &dm)); 1066 PetscCall(TSSetDM(ts, dm)); 1067 PetscCall(ProcessMonitorOptions(comm, &user)); 1068 1069 PetscCall(DMCreateGlobalVector(dm, &u)); 1070 PetscCall(PetscObjectSetName((PetscObject) u, "solution")); 1071 if (user.useFV) { 1072 PetscBool isImplicit = PETSC_FALSE; 1073 1074 PetscCall(PetscOptionsHasName(NULL,"", "-use_implicit", &isImplicit)); 1075 if (isImplicit) { 1076 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user)); 1077 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user)); 1078 } 1079 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user)); 1080 PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &user)); 1081 } else { 1082 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &user)); 1083 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &user)); 1084 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &user)); 1085 } 1086 if (user.useFV) PetscCall(TSMonitorSet(ts, MonitorFunctionals, &user, NULL)); 1087 PetscCall(TSSetMaxSteps(ts, 1)); 1088 PetscCall(TSSetMaxTime(ts, 2.0)); 1089 PetscCall(TSSetTimeStep(ts,0.01)); 1090 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1091 PetscCall(TSSetFromOptions(ts)); 1092 1093 PetscCall(DMProjectFunction(dm, 0.0, user.initialGuess, ctxs, INSERT_VALUES, u)); 1094 if (user.useFV) PetscCall(SetInitialConditionFVM(dm, u, 1, user.initialGuess[1], ctxs[1])); 1095 PetscCall(VecViewFromOptions(u, NULL, "-init_vec_view")); 1096 PetscCall(TSGetTime(ts, &t)); 1097 t0 = t; 1098 PetscCall(DMTSCheckFromOptions(ts, u)); 1099 PetscCall(TSSolve(ts, u)); 1100 PetscCall(TSGetTime(ts, &t)); 1101 if (t > t0) PetscCall(DMTSCheckFromOptions(ts, u)); 1102 PetscCall(VecViewFromOptions(u, NULL, "-sol_vec_view")); 1103 { 1104 PetscReal ftime; 1105 PetscInt nsteps; 1106 TSConvergedReason reason; 1107 1108 PetscCall(TSGetSolveTime(ts, &ftime)); 1109 PetscCall(TSGetStepNumber(ts, &nsteps)); 1110 PetscCall(TSGetConvergedReason(ts, &reason)); 1111 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %D steps\n", TSConvergedReasons[reason], (double) ftime, nsteps)); 1112 } 1113 1114 PetscCall(VecDestroy(&u)); 1115 PetscCall(DMDestroy(&dm)); 1116 PetscCall(TSDestroy(&ts)); 1117 PetscCall(PetscFree(user.monitorFuncs)); 1118 PetscCall(FunctionalDestroy(&user.functionalRegistry)); 1119 PetscCall(PetscFinalize()); 1120 return 0; 1121 } 1122 1123 /*TEST 1124 1125 testset: 1126 args: -dm_plex_simplex 0 -dm_plex_box_faces 3,3,3 1127 1128 # 2D harmonic velocity, no porosity 1129 test: 1130 suffix: p1p1 1131 requires: !complex !single 1132 args: -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_factor_shift_type nonzero -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1133 test: 1134 suffix: p1p1_xper 1135 requires: !complex !single 1136 args: -dm_refine 1 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1137 test: 1138 suffix: p1p1_xper_ref 1139 requires: !complex !single 1140 args: -dm_refine 2 -dm_plex_box_bd periodic,none -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1141 test: 1142 suffix: p1p1_xyper 1143 requires: !complex !single 1144 args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1145 test: 1146 suffix: p1p1_xyper_ref 1147 requires: !complex !single 1148 args: -dm_refine 2 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 1 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1149 test: 1150 suffix: p2p1 1151 requires: !complex !single 1152 args: -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1153 test: 1154 suffix: p2p1_xyper 1155 requires: !complex !single 1156 args: -dm_refine 1 -dm_plex_box_bd periodic,periodic -velocity_petscspace_degree 2 -porosity_petscspace_degree 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -pc_factor_shift_type nonzero -ksp_rtol 1.0e-8 -ts_monitor -snes_error_if_not_converged -ksp_error_if_not_converged -dmts_check 1157 1158 test: 1159 suffix: adv_1 1160 requires: !complex !single 1161 args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -ts_view -dm_view 1162 1163 test: 1164 suffix: adv_2 1165 requires: !complex 1166 TODO: broken memory corruption 1167 args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3,4 -bc_outflow 1,2 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason 1168 1169 test: 1170 suffix: adv_3 1171 requires: !complex 1172 TODO: broken memory corruption 1173 args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason 1174 1175 test: 1176 suffix: adv_3_ex 1177 requires: !complex 1178 args: -dm_plex_box_bd periodic,none -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.1 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -ksp_max_it 100 -ts_view -dm_view 1179 1180 test: 1181 suffix: adv_4 1182 requires: !complex 1183 TODO: broken memory corruption 1184 args: -use_fv -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 3 -bc_outflow 1 -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason 1185 1186 # 2D Advection, box, delta 1187 test: 1188 suffix: adv_delta_yper_0 1189 requires: !complex 1190 TODO: broken 1191 args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error 1192 1193 test: 1194 suffix: adv_delta_yper_1 1195 requires: !complex 1196 TODO: broken 1197 args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 1 -source_loc 0.416666,0.416666 1198 1199 test: 1200 suffix: adv_delta_yper_2 1201 requires: !complex 1202 TODO: broken 1203 args: -dm_plex_box_bd none,periodic -use_fv -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type euler -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -dm_view -monitor Error -dm_refine 2 -source_loc 0.458333,0.458333 1204 1205 test: 1206 suffix: adv_delta_yper_fim_0 1207 requires: !complex 1208 TODO: broken 1209 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason 1210 1211 test: 1212 suffix: adv_delta_yper_fim_1 1213 requires: !complex 1214 TODO: broken 1215 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic 1216 1217 test: 1218 suffix: adv_delta_yper_fim_2 1219 requires: !complex 1220 TODO: broken 1221 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -mat_coloring_greedy_symmetric 0 -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -snes_linesearch_type basic 1222 1223 test: 1224 suffix: adv_delta_yper_im_0 1225 requires: !complex 1226 TODO: broken 1227 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason 1228 1229 test: 1230 suffix: adv_delta_yper_im_1 1231 requires: !complex 1232 TODO: broken 1233 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666 1234 1235 test: 1236 suffix: adv_delta_yper_im_2 1237 requires: !complex 1238 TODO: broken 1239 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 0 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333 1240 1241 test: 1242 suffix: adv_delta_yper_im_3 1243 requires: !complex 1244 TODO: broken 1245 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason 1246 1247 # I believe the nullspace is sin(pi y) 1248 test: 1249 suffix: adv_delta_yper_im_4 1250 requires: !complex 1251 TODO: broken 1252 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 40 -ts_dt 0.166666 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 1 -source_loc 0.416666,0.416666 1253 1254 test: 1255 suffix: adv_delta_yper_im_5 1256 requires: !complex 1257 TODO: broken 1258 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_mimex_version 0 -ts_max_time 5.0 -ts_max_steps 80 -ts_dt 0.083333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -dm_refine 2 -source_loc 0.458333,0.458333 1259 1260 test: 1261 suffix: adv_delta_yper_im_6 1262 requires: !complex 1263 TODO: broken 1264 args: -dm_plex_box_bd none,periodic -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 2 -bc_outflow 4 -ts_view -monitor Error -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type svd -snes_converged_reason 1265 # 2D Advection, magma benchmark 1 1266 1267 test: 1268 suffix: adv_delta_shear_im_0 1269 requires: !complex 1270 TODO: broken 1271 args: -dm_plex_box_bd periodic,none -dm_refine 2 -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist shear -porosity_dist delta -inflow_state 0.0 -ts_type mimex -ts_max_time 5.0 -ts_max_steps 20 -ts_dt 0.333333 -bc_inflow 1,3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 -pc_type lu -snes_converged_reason -source_loc 0.458333,0.708333 1272 # 2D Advection, box, gaussian 1273 1274 test: 1275 suffix: adv_gauss 1276 requires: !complex 1277 TODO: broken 1278 args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view 1279 1280 test: 1281 suffix: adv_gauss_im 1282 requires: !complex 1283 TODO: broken 1284 args: -use_fv -use_implicit -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 1285 1286 test: 1287 suffix: adv_gauss_im_1 1288 requires: !complex 1289 TODO: broken 1290 args: -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 1291 1292 test: 1293 suffix: adv_gauss_im_2 1294 requires: !complex 1295 TODO: broken 1296 args: -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type beuler -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 3 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -snes_rtol 1.0e-7 1297 1298 test: 1299 suffix: adv_gauss_corner 1300 requires: !complex 1301 TODO: broken 1302 args: -use_fv -velocity_dist constant -porosity_dist gaussian -inflow_state 0.0 -ts_type ssp -ts_max_time 2.0 -ts_max_steps 100 -ts_dt 0.01 -bc_inflow 1 -bc_outflow 2 -ts_view -dm_view 1303 1304 # 2D Advection+Harmonic 12- 1305 test: 1306 suffix: adv_harm_0 1307 requires: !complex 1308 TODO: broken memory corruption 1309 args: -velocity_petscspace_degree 2 -use_fv -velocity_dist harmonic -porosity_dist gaussian -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -bc_inflow 1,2,4 -bc_outflow 3 -use_implicit -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -ksp_max_it 100 -ts_view -dm_view -snes_converged_reason -ksp_converged_reason -snes_monitor -dmts_check 1310 1311 # Must check that FV BCs propagate to coarse meshes 1312 # Must check that FV BC ids propagate to coarse meshes 1313 # Must check that FE+FV BCs work at the same time 1314 # 2D Advection, matching wind in ex11 8-11 1315 # NOTE implicit solves are limited by accuracy of FD Jacobian 1316 test: 1317 suffix: adv_0 1318 requires: !complex !single exodusii 1319 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -velocity_dist zero -porosity_dist tilted -ts_type ssp -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view 1320 1321 test: 1322 suffix: adv_0_im 1323 requires: !complex exodusii 1324 TODO: broken memory corruption 1325 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist zero -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu 1326 1327 test: 1328 suffix: adv_0_im_2 1329 requires: !complex exodusii 1330 TODO: broken 1331 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type lu -snes_rtol 1.0e-7 1332 1333 test: 1334 suffix: adv_0_im_3 1335 requires: !complex exodusii 1336 TODO: broken 1337 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 1 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7 1338 1339 test: 1340 suffix: adv_0_im_4 1341 requires: !complex exodusii 1342 TODO: broken 1343 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -use_fv -use_implicit -velocity_petscspace_degree 2 -velocity_dist constant -porosity_dist tilted -ts_type beuler -ts_max_time 2.0 -ts_max_steps 1000 -ts_dt 0.993392 -ts_view -dm_view -snes_fd_color -snes_fd_color_use_mat -mat_coloring_type greedy -pc_type svd -snes_rtol 1.0e-7 1344 # 2D Advection, misc 1345 1346 TEST*/ 1347