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