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