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