1 static char help[] = "Heat Equation in 2d and 3d with finite elements.\n\ 2 We solve the heat equation in a rectangular\n\ 3 domain, using a parallel unstructured mesh (DMPLEX) to discretize it.\n\ 4 Contributed by: Julian Andrej <juan@tf.uni-kiel.de>\n\n\n"; 5 6 #include <petscdmplex.h> 7 #include <petscds.h> 8 #include <petscts.h> 9 10 /* 11 Heat equation: 12 13 du/dt - \Delta u + f = 0 14 */ 15 16 typedef enum {SOL_QUADRATIC_LINEAR, SOL_QUADRATIC_TRIG, SOL_TRIG_LINEAR, SOL_TRIG_TRIG, NUM_SOLUTION_TYPES} SolutionType; 17 const char *solutionTypes[NUM_SOLUTION_TYPES+1] = {"quadratic_linear", "quadratic_trig", "trig_linear", "trig_trig", "unknown"}; 18 19 typedef struct { 20 SolutionType solType; /* Type of exact solution */ 21 /* Solver setup */ 22 PetscBool expTS; /* Flag for explicit timestepping */ 23 PetscBool lumped; /* Lump the mass matrix */ 24 } AppCtx; 25 26 /* 27 Exact 2D solution: 28 u = 2t + x^2 + y^2 29 u_t = 2 30 \Delta u = 2 + 2 = 4 31 f = 2 32 F(u) = 2 - (2 + 2) + 2 = 0 33 34 Exact 3D solution: 35 u = 3t + x^2 + y^2 + z^2 36 F(u) = 3 - (2 + 2 + 2) + 3 = 0 37 */ 38 static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 39 { 40 PetscInt d; 41 42 *u = dim*time; 43 for (d = 0; d < dim; ++d) *u += x[d]*x[d]; 44 return 0; 45 } 46 47 static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 48 { 49 *u = dim; 50 return 0; 51 } 52 53 static void f0_quad_lin_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 54 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 55 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 56 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 57 { 58 f0[0] = -(PetscScalar) dim; 59 } 60 static void f0_quad_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux, 61 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 62 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 63 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 64 { 65 PetscScalar exp[1] = {0.}; 66 f0_quad_lin_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp); 67 f0[0] = u_t[0] - exp[0]; 68 } 69 70 /* 71 Exact 2D solution: 72 u = 2*cos(t) + x^2 + y^2 73 F(u) = -2*sint(t) - (2 + 2) + 2*sin(t) + 4 = 0 74 75 Exact 3D solution: 76 u = 3*cos(t) + x^2 + y^2 + z^2 77 F(u) = -3*sin(t) - (2 + 2 + 2) + 3*sin(t) + 6 = 0 78 */ 79 static PetscErrorCode mms_quad_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 80 { 81 PetscInt d; 82 83 *u = dim*PetscCosReal(time); 84 for (d = 0; d < dim; ++d) *u += x[d]*x[d]; 85 return 0; 86 } 87 88 static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 89 { 90 *u = -dim*PetscSinReal(time); 91 return 0; 92 } 93 94 static void f0_quad_trig_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 95 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 96 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 97 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 98 { 99 f0[0] = -dim*(PetscSinReal(t) + 2.0); 100 } 101 static void f0_quad_trig(PetscInt dim, PetscInt Nf, PetscInt NfAux, 102 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 103 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 104 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 105 { 106 PetscScalar exp[1] = {0.}; 107 f0_quad_trig_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp); 108 f0[0] = u_t[0] - exp[0]; 109 } 110 111 /* 112 Exact 2D solution: 113 u = 2\pi^2 t + cos(\pi x) + cos(\pi y) 114 F(u) = 2\pi^2 - \pi^2 (cos(\pi x) + cos(\pi y)) + \pi^2 (cos(\pi x) + cos(\pi y)) - 2\pi^2 = 0 115 116 Exact 3D solution: 117 u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z) 118 F(u) = 3\pi^2 - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) + \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) - 3\pi^2 = 0 119 */ 120 static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 121 { 122 PetscInt d; 123 124 *u = dim*PetscSqr(PETSC_PI)*time; 125 for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI*x[d]); 126 return 0; 127 } 128 129 static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 130 { 131 *u = dim*PetscSqr(PETSC_PI); 132 return 0; 133 } 134 135 static void f0_trig_lin(PetscInt dim, PetscInt Nf, PetscInt NfAux, 136 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 137 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 138 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 139 { 140 PetscInt d; 141 f0[0] = u_t[0]; 142 for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI)*(PetscCosReal(PETSC_PI*x[d]) - 1.0); 143 } 144 145 /* 146 Exact 2D solution: 147 u = pi^2 cos(t) + cos(\pi x) + cos(\pi y) 148 u_t = -pi^2 sin(t) 149 \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y)) 150 f = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y)) 151 F(u) = -\pi^2 sin(t) + \pi^2 (cos(\pi x) + cos(\pi y)) - \pi^2 (cos(\pi x) + cos(\pi y)) + \pi^2 sin(t) = 0 152 153 Exact 3D solution: 154 u = pi^2 cos(t) + cos(\pi x) + cos(\pi y) + cos(\pi z) 155 u_t = -pi^2 sin(t) 156 \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) 157 f = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) 158 F(u) = -\pi^2 sin(t) + \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) + \pi^2 sin(t) = 0 159 */ 160 static PetscErrorCode mms_trig_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 161 { 162 PetscInt d; 163 164 *u = PetscSqr(PETSC_PI)*PetscCosReal(time); 165 for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI*x[d]); 166 return 0; 167 } 168 169 static PetscErrorCode mms_trig_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 170 { 171 *u = -PetscSqr(PETSC_PI)*PetscSinReal(time); 172 return 0; 173 } 174 175 static void f0_trig_trig_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 176 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 177 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 178 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 179 { 180 PetscInt d; 181 f0[0] -= PetscSqr(PETSC_PI)*PetscSinReal(t); 182 for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI)*PetscCosReal(PETSC_PI*x[d]); 183 } 184 static void f0_trig_trig(PetscInt dim, PetscInt Nf, PetscInt NfAux, 185 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 186 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 187 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]) 188 { 189 PetscScalar exp[1] = {0.}; 190 f0_trig_trig_exp(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, exp); 191 f0[0] = u_t[0] - exp[0]; 192 } 193 194 static void f1_temp_exp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 195 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 196 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 197 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]) 198 { 199 PetscInt d; 200 for (d = 0; d < dim; ++d) f1[d] = -u_x[d]; 201 } 202 static void f1_temp(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 f1[]) 206 { 207 PetscInt d; 208 for (d = 0; d < dim; ++d) f1[d] = u_x[d]; 209 } 210 211 static void g3_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 212 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 213 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 214 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]) 215 { 216 PetscInt d; 217 for (d = 0; d < dim; ++d) g3[d*dim+d] = 1.0; 218 } 219 220 static void g0_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux, 221 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], 222 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], 223 PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]) 224 { 225 g0[0] = u_tShift*1.0; 226 } 227 228 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 229 { 230 PetscInt sol; 231 PetscErrorCode ierr; 232 233 PetscFunctionBeginUser; 234 options->solType = SOL_QUADRATIC_LINEAR; 235 options->expTS = PETSC_FALSE; 236 options->lumped = PETSC_TRUE; 237 238 ierr = PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");CHKERRQ(ierr); 239 CHKERRQ(PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL)); 240 options->solType = (SolutionType) sol; 241 CHKERRQ(PetscOptionsBool("-explicit", "Use explicit timestepping", "ex45.c", options->expTS, &options->expTS, NULL)); 242 CHKERRQ(PetscOptionsBool("-lumped", "Lump the mass matrix", "ex45.c", options->lumped, &options->lumped, NULL)); 243 ierr = PetscOptionsEnd();CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 248 { 249 PetscFunctionBeginUser; 250 CHKERRQ(DMCreate(comm, dm)); 251 CHKERRQ(DMSetType(*dm, DMPLEX)); 252 CHKERRQ(DMSetFromOptions(*dm)); 253 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 254 PetscFunctionReturn(0); 255 } 256 257 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 258 { 259 PetscDS ds; 260 DMLabel label; 261 const PetscInt id = 1; 262 263 PetscFunctionBeginUser; 264 CHKERRQ(DMGetLabel(dm, "marker", &label)); 265 CHKERRQ(DMGetDS(dm, &ds)); 266 CHKERRQ(PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp)); 267 switch (ctx->solType) { 268 case SOL_QUADRATIC_LINEAR: 269 CHKERRQ(PetscDSSetResidual(ds, 0, f0_quad_lin, f1_temp)); 270 CHKERRQ(PetscDSSetRHSResidual(ds, 0, f0_quad_lin_exp, f1_temp_exp)); 271 CHKERRQ(PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx)); 272 CHKERRQ(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx)); 273 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms_quad_lin, (void (*)(void)) mms_quad_lin_t, ctx, NULL)); 274 break; 275 case SOL_QUADRATIC_TRIG: 276 CHKERRQ(PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp)); 277 CHKERRQ(PetscDSSetRHSResidual(ds, 0, f0_quad_trig_exp, f1_temp_exp)); 278 CHKERRQ(PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx)); 279 CHKERRQ(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx)); 280 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms_quad_trig, (void (*)(void)) mms_quad_trig_t, ctx, NULL)); 281 break; 282 case SOL_TRIG_LINEAR: 283 CHKERRQ(PetscDSSetResidual(ds, 0, f0_trig_lin, f1_temp)); 284 CHKERRQ(PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx)); 285 CHKERRQ(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx)); 286 CHKERRQ(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)(void)) mms_trig_lin, (void (*)(void)) mms_trig_lin_t, ctx, NULL)); 287 break; 288 case SOL_TRIG_TRIG: 289 CHKERRQ(PetscDSSetResidual(ds, 0, f0_trig_trig, f1_temp)); 290 CHKERRQ(PetscDSSetRHSResidual(ds, 0, f0_trig_trig_exp, f1_temp_exp)); 291 CHKERRQ(PetscDSSetExactSolution(ds, 0, mms_trig_trig, ctx)); 292 CHKERRQ(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_trig_t, ctx)); 293 break; 294 default: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%D)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType); 295 } 296 PetscFunctionReturn(0); 297 } 298 299 static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx) 300 { 301 DM cdm = dm; 302 PetscFE fe; 303 DMPolytopeType ct; 304 PetscBool simplex; 305 PetscInt dim, cStart; 306 307 PetscFunctionBeginUser; 308 CHKERRQ(DMGetDimension(dm, &dim)); 309 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 310 CHKERRQ(DMPlexGetCellType(dm, cStart, &ct)); 311 simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE; 312 /* Create finite element */ 313 CHKERRQ(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe)); 314 CHKERRQ(PetscObjectSetName((PetscObject) fe, "temperature")); 315 /* Set discretization and boundary conditions for each mesh */ 316 CHKERRQ(DMSetField(dm, 0, NULL, (PetscObject) fe)); 317 CHKERRQ(DMCreateDS(dm)); 318 if (ctx->expTS) { 319 PetscDS ds; 320 321 CHKERRQ(DMGetDS(dm, &ds)); 322 CHKERRQ(PetscDSSetImplicit(ds, 0, PETSC_FALSE)); 323 } 324 CHKERRQ(SetupProblem(dm, ctx)); 325 while (cdm) { 326 CHKERRQ(DMCopyDisc(dm, cdm)); 327 CHKERRQ(DMGetCoarseDM(cdm, &cdm)); 328 } 329 CHKERRQ(PetscFEDestroy(&fe)); 330 PetscFunctionReturn(0); 331 } 332 333 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 334 { 335 DM dm; 336 PetscReal t; 337 338 PetscFunctionBegin; 339 CHKERRQ(TSGetDM(ts, &dm)); 340 CHKERRQ(TSGetTime(ts, &t)); 341 CHKERRQ(DMComputeExactSolution(dm, t, u, NULL)); 342 PetscFunctionReturn(0); 343 } 344 345 int main(int argc, char **argv) 346 { 347 DM dm; 348 TS ts; 349 Vec u; 350 AppCtx ctx; 351 PetscErrorCode ierr; 352 353 ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; 354 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 355 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx)); 356 CHKERRQ(DMSetApplicationContext(dm, &ctx)); 357 CHKERRQ(SetupDiscretization(dm, &ctx)); 358 359 CHKERRQ(TSCreate(PETSC_COMM_WORLD, &ts)); 360 CHKERRQ(TSSetDM(ts, dm)); 361 CHKERRQ(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 362 if (ctx.expTS) { 363 CHKERRQ(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFEM, &ctx)); 364 if (ctx.lumped) CHKERRQ(DMTSCreateRHSMassMatrixLumped(dm)); 365 else CHKERRQ(DMTSCreateRHSMassMatrix(dm)); 366 } else { 367 CHKERRQ(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 368 CHKERRQ(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 369 } 370 CHKERRQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 371 CHKERRQ(TSSetFromOptions(ts)); 372 CHKERRQ(TSSetComputeInitialCondition(ts, SetInitialConditions)); 373 374 CHKERRQ(DMCreateGlobalVector(dm, &u)); 375 CHKERRQ(DMTSCheckFromOptions(ts, u)); 376 CHKERRQ(SetInitialConditions(ts, u)); 377 CHKERRQ(PetscObjectSetName((PetscObject) u, "temperature")); 378 CHKERRQ(TSSolve(ts, u)); 379 CHKERRQ(DMTSCheckFromOptions(ts, u)); 380 if (ctx.expTS) CHKERRQ(DMTSDestroyRHSMassMatrix(dm)); 381 382 CHKERRQ(VecDestroy(&u)); 383 CHKERRQ(TSDestroy(&ts)); 384 CHKERRQ(DMDestroy(&dm)); 385 ierr = PetscFinalize(); 386 return ierr; 387 } 388 389 /*TEST 390 391 test: 392 suffix: 2d_p1 393 requires: triangle 394 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 395 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 396 test: 397 suffix: 2d_p1_exp 398 requires: triangle 399 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -explicit \ 400 -ts_type euler -ts_max_steps 4 -ts_dt 1e-3 -ts_monitor_error 401 test: 402 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 403 suffix: 2d_p1_sconv 404 requires: triangle 405 args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 406 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 407 test: 408 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.1] 409 suffix: 2d_p1_sconv_2 410 requires: triangle 411 args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 412 -ts_type beuler -ts_max_steps 1 -ts_dt 1e-6 -snes_error_if_not_converged -pc_type lu 413 test: 414 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 415 suffix: 2d_p1_tconv 416 requires: triangle 417 args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 418 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 419 test: 420 # -dm_refine 6 -convest_num_refine 3 get L_2 convergence rate: [1.0] 421 suffix: 2d_p1_tconv_2 422 requires: triangle 423 args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 424 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 425 test: 426 # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid 427 suffix: 2d_p1_exp_tconv_2 428 requires: triangle 429 args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \ 430 -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 -lumped 0 -mass_pc_type lu 431 test: 432 # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid 433 suffix: 2d_p1_exp_tconv_2_lump 434 requires: triangle 435 args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \ 436 -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 437 test: 438 suffix: 2d_p2 439 requires: triangle 440 args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 441 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 442 test: 443 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 444 suffix: 2d_p2_sconv 445 requires: triangle 446 args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 447 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 448 test: 449 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [3.1] 450 suffix: 2d_p2_sconv_2 451 requires: triangle 452 args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 453 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 454 test: 455 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 456 suffix: 2d_p2_tconv 457 requires: triangle 458 args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 459 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 460 test: 461 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 462 suffix: 2d_p2_tconv_2 463 requires: triangle 464 args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 465 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 466 test: 467 suffix: 2d_q1 468 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 469 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 470 test: 471 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 472 suffix: 2d_q1_sconv 473 args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 474 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 475 test: 476 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 477 suffix: 2d_q1_tconv 478 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 479 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 480 test: 481 suffix: 2d_q2 482 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 483 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 484 test: 485 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 486 suffix: 2d_q2_sconv 487 args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 488 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 489 test: 490 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 491 suffix: 2d_q2_tconv 492 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 493 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 494 495 test: 496 suffix: 3d_p1 497 requires: ctetgen 498 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 499 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 500 test: 501 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 502 suffix: 3d_p1_sconv 503 requires: ctetgen 504 args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 505 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 506 test: 507 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 508 suffix: 3d_p1_tconv 509 requires: ctetgen 510 args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 511 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 512 test: 513 suffix: 3d_p2 514 requires: ctetgen 515 args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 516 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 517 test: 518 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 519 suffix: 3d_p2_sconv 520 requires: ctetgen 521 args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 522 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 523 test: 524 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 525 suffix: 3d_p2_tconv 526 requires: ctetgen 527 args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 528 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 529 test: 530 suffix: 3d_q1 531 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 532 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 533 test: 534 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 535 suffix: 3d_q1_sconv 536 args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 537 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 538 test: 539 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 540 suffix: 3d_q1_tconv 541 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 542 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 543 test: 544 suffix: 3d_q2 545 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 546 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 547 test: 548 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 549 suffix: 3d_q2_sconv 550 args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 551 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 552 test: 553 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 554 suffix: 3d_q2_tconv 555 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 556 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 557 558 test: 559 # For a nice picture, -bd_dm_refine 2 -dm_refine 1 -dm_view hdf5:${PETSC_DIR}/sol.h5 -ts_monitor_solution hdf5:${PETSC_DIR}/sol.h5::append 560 suffix: egads_sphere 561 requires: egads ctetgen 562 args: -sol_type quadratic_linear \ 563 -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_plex_boundary_label marker -bd_dm_plex_scale 40 \ 564 -temp_petscspace_degree 2 -dmts_check .0001 \ 565 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 566 567 TEST*/ 568