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