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 PetscInt remesh_every; /* Remesh every number of steps */ 31 DM remesh_dm; /* New DM after remeshing */ 32 } AppCtx; 33 34 /* 35 Exact 2D solution: 36 u = 2t + x^2 + y^2 37 u_t = 2 38 \Delta u = 2 + 2 = 4 39 f = 2 40 F(u) = 2 - (2 + 2) + 2 = 0 41 42 Exact 3D solution: 43 u = 3t + x^2 + y^2 + z^2 44 F(u) = 3 - (2 + 2 + 2) + 3 = 0 45 */ 46 static PetscErrorCode mms_quad_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 47 { 48 PetscInt d; 49 50 *u = dim * time; 51 for (d = 0; d < dim; ++d) *u += x[d] * x[d]; 52 return PETSC_SUCCESS; 53 } 54 55 static PetscErrorCode mms_quad_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 56 { 57 *u = dim; 58 return PETSC_SUCCESS; 59 } 60 61 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[]) 62 { 63 f0[0] = -(PetscScalar)dim; 64 } 65 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[]) 66 { 67 PetscScalar exp[1] = {0.}; 68 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); 69 f0[0] = u_t[0] - exp[0]; 70 } 71 72 /* 73 Exact 2D solution: 74 u = 2*cos(t) + x^2 + y^2 75 F(u) = -2*sint(t) - (2 + 2) + 2*sin(t) + 4 = 0 76 77 Exact 3D solution: 78 u = 3*cos(t) + x^2 + y^2 + z^2 79 F(u) = -3*sin(t) - (2 + 2 + 2) + 3*sin(t) + 6 = 0 80 */ 81 static PetscErrorCode mms_quad_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 82 { 83 PetscInt d; 84 85 *u = dim * PetscCosReal(time); 86 for (d = 0; d < dim; ++d) *u += x[d] * x[d]; 87 return PETSC_SUCCESS; 88 } 89 90 static PetscErrorCode mms_quad_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 91 { 92 *u = -dim * PetscSinReal(time); 93 return PETSC_SUCCESS; 94 } 95 96 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[]) 97 { 98 f0[0] = -dim * (PetscSinReal(t) + 2.0); 99 } 100 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[]) 101 { 102 PetscScalar exp[1] = {0.}; 103 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); 104 f0[0] = u_t[0] - exp[0]; 105 } 106 107 /* 108 Exact 2D solution: 109 u = 2\pi^2 t + cos(\pi x) + cos(\pi y) 110 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 111 112 Exact 3D solution: 113 u = 3\pi^2 t + cos(\pi x) + cos(\pi y) + cos(\pi z) 114 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 115 */ 116 static PetscErrorCode mms_trig_lin(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 117 { 118 PetscInt d; 119 120 *u = dim * PetscSqr(PETSC_PI) * time; 121 for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]); 122 return PETSC_SUCCESS; 123 } 124 125 static PetscErrorCode mms_trig_lin_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 126 { 127 *u = dim * PetscSqr(PETSC_PI); 128 return PETSC_SUCCESS; 129 } 130 131 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[]) 132 { 133 PetscInt d; 134 f0[0] = u_t[0]; 135 for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * (PetscCosReal(PETSC_PI * x[d]) - 1.0); 136 } 137 138 /* 139 Exact 2D solution: 140 u = pi^2 cos(t) + cos(\pi x) + cos(\pi y) 141 u_t = -pi^2 sin(t) 142 \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y)) 143 f = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y)) 144 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 145 146 Exact 3D solution: 147 u = pi^2 cos(t) + cos(\pi x) + cos(\pi y) + cos(\pi z) 148 u_t = -pi^2 sin(t) 149 \Delta u = -\pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) 150 f = pi^2 sin(t) - \pi^2 (cos(\pi x) + cos(\pi y) + cos(\pi z)) 151 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 152 */ 153 static PetscErrorCode mms_trig_trig(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 154 { 155 PetscInt d; 156 157 *u = PetscSqr(PETSC_PI) * PetscCosReal(time); 158 for (d = 0; d < dim; ++d) *u += PetscCosReal(PETSC_PI * x[d]); 159 return PETSC_SUCCESS; 160 } 161 162 static PetscErrorCode mms_trig_trig_t(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 163 { 164 *u = -PetscSqr(PETSC_PI) * PetscSinReal(time); 165 return PETSC_SUCCESS; 166 } 167 168 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[]) 169 { 170 PetscInt d; 171 f0[0] -= PetscSqr(PETSC_PI) * PetscSinReal(t); 172 for (d = 0; d < dim; ++d) f0[0] += PetscSqr(PETSC_PI) * PetscCosReal(PETSC_PI * x[d]); 173 } 174 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[]) 175 { 176 PetscScalar exp[1] = {0.}; 177 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); 178 f0[0] = u_t[0] - exp[0]; 179 } 180 181 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[]) 182 { 183 for (PetscInt d = 0; d < dim; ++d) f1[d] = -u_x[d]; 184 } 185 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[]) 186 { 187 for (PetscInt 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 for (PetscInt d = 0; d < dim; ++d) g3[d * dim + d] = 1.0; 193 } 194 195 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[]) 196 { 197 g0[0] = u_tShift * 1.0; 198 } 199 200 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 201 { 202 PetscInt sol; 203 204 PetscFunctionBeginUser; 205 options->solType = SOL_QUADRATIC_LINEAR; 206 options->expTS = PETSC_FALSE; 207 options->lumped = PETSC_TRUE; 208 options->remesh_every = 0; 209 210 PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX"); 211 sol = options->solType; 212 PetscCall(PetscOptionsEList("-sol_type", "Type of exact solution", "ex45.c", solutionTypes, NUM_SOLUTION_TYPES, solutionTypes[options->solType], &sol, NULL)); 213 options->solType = (SolutionType)sol; 214 PetscCall(PetscOptionsBool("-explicit", "Use explicit timestepping", "ex45.c", options->expTS, &options->expTS, NULL)); 215 PetscCall(PetscOptionsBool("-lumped", "Lump the mass matrix", "ex45.c", options->lumped, &options->lumped, NULL)); 216 PetscCall(PetscOptionsInt("-remesh_every", "Remesh every number of steps", "ex45.c", options->remesh_every, &options->remesh_every, NULL)); 217 PetscOptionsEnd(); 218 PetscFunctionReturn(PETSC_SUCCESS); 219 } 220 221 static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx) 222 { 223 PetscFunctionBeginUser; 224 PetscCall(DMCreate(comm, dm)); 225 PetscCall(DMSetType(*dm, DMPLEX)); 226 PetscCall(DMSetFromOptions(*dm)); 227 { 228 char convType[256]; 229 PetscBool flg; 230 PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX"); 231 PetscCall(PetscOptionsFList("-dm_plex_convert_type", "Convert DMPlex to another format", __FILE__, DMList, DMPLEX, convType, 256, &flg)); 232 PetscOptionsEnd(); 233 if (flg) { 234 DM dmConv; 235 PetscCall(DMConvert(*dm, convType, &dmConv)); 236 if (dmConv) { 237 PetscCall(DMDestroy(dm)); 238 *dm = dmConv; 239 PetscCall(DMSetFromOptions(*dm)); 240 PetscCall(DMSetUp(*dm)); 241 } 242 } 243 } 244 PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); 245 PetscFunctionReturn(PETSC_SUCCESS); 246 } 247 248 static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx) 249 { 250 PetscDS ds; 251 DMLabel label; 252 const PetscInt id = 1; 253 254 PetscFunctionBeginUser; 255 PetscCall(DMGetLabel(dm, "marker", &label)); 256 PetscCall(DMGetDS(dm, &ds)); 257 PetscCall(PetscDSSetJacobian(ds, 0, 0, g0_temp, NULL, NULL, g3_temp)); 258 switch (ctx->solType) { 259 case SOL_QUADRATIC_LINEAR: 260 PetscCall(PetscDSSetResidual(ds, 0, f0_quad_lin, f1_temp)); 261 PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_lin_exp, f1_temp_exp)); 262 PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_lin, ctx)); 263 PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_lin_t, ctx)); 264 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)); 265 break; 266 case SOL_QUADRATIC_TRIG: 267 PetscCall(PetscDSSetResidual(ds, 0, f0_quad_trig, f1_temp)); 268 PetscCall(PetscDSSetRHSResidual(ds, 0, f0_quad_trig_exp, f1_temp_exp)); 269 PetscCall(PetscDSSetExactSolution(ds, 0, mms_quad_trig, ctx)); 270 PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_quad_trig_t, ctx)); 271 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)); 272 break; 273 case SOL_TRIG_LINEAR: 274 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_lin, f1_temp)); 275 PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_lin, ctx)); 276 PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_lin_t, ctx)); 277 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)); 278 break; 279 case SOL_TRIG_TRIG: 280 PetscCall(PetscDSSetResidual(ds, 0, f0_trig_trig, f1_temp)); 281 PetscCall(PetscDSSetRHSResidual(ds, 0, f0_trig_trig_exp, f1_temp_exp)); 282 PetscCall(PetscDSSetExactSolution(ds, 0, mms_trig_trig, ctx)); 283 PetscCall(PetscDSSetExactSolutionTimeDerivative(ds, 0, mms_trig_trig_t, ctx)); 284 break; 285 default: 286 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid solution type: %s (%d)", solutionTypes[PetscMin(ctx->solType, NUM_SOLUTION_TYPES)], ctx->solType); 287 } 288 PetscFunctionReturn(PETSC_SUCCESS); 289 } 290 291 static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx) 292 { 293 DM plex, cdm = dm; 294 PetscFE fe; 295 PetscBool simplex; 296 PetscInt dim; 297 298 PetscFunctionBeginUser; 299 PetscCall(DMGetDimension(dm, &dim)); 300 PetscCall(DMConvert(dm, DMPLEX, &plex)); 301 PetscCall(DMPlexIsSimplex(plex, &simplex)); 302 PetscCall(DMDestroy(&plex)); 303 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, simplex, "temp_", -1, &fe)); 304 PetscCall(PetscObjectSetName((PetscObject)fe, "temperature")); 305 /* Set discretization and boundary conditions for each mesh */ 306 PetscCall(DMSetField(dm, 0, NULL, (PetscObject)fe)); 307 PetscCall(DMCreateDS(dm)); 308 if (ctx->expTS) { 309 PetscDS ds; 310 311 PetscCall(DMGetDS(dm, &ds)); 312 PetscCall(PetscDSSetImplicit(ds, 0, PETSC_FALSE)); 313 } 314 PetscCall(SetupProblem(dm, ctx)); 315 while (cdm) { 316 PetscCall(DMCopyDisc(dm, cdm)); 317 PetscCall(DMGetCoarseDM(cdm, &cdm)); 318 } 319 PetscCall(PetscFEDestroy(&fe)); 320 PetscFunctionReturn(PETSC_SUCCESS); 321 } 322 323 #include <petsc/private/dmpleximpl.h> 324 static PetscErrorCode Remesh(DM dm, Vec U, DM *newdm) 325 { 326 PetscFunctionBeginUser; 327 PetscCall(DMViewFromOptions(dm, NULL, "-remesh_dmin_view")); 328 *newdm = NULL; 329 330 PetscInt dim; 331 DM plex; 332 PetscBool simplex; 333 PetscCall(DMGetDimension(dm, &dim)); 334 PetscCall(DMConvert(dm, DMPLEX, &plex)); 335 PetscCall(DMPlexIsSimplex(plex, &simplex)); 336 PetscCall(DMDestroy(&plex)); 337 338 DM dmGrad; 339 PetscFE fe; 340 PetscCall(DMClone(dm, &dmGrad)); 341 PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, dim, simplex, "grad_", -1, &fe)); 342 PetscCall(DMSetField(dmGrad, 0, NULL, (PetscObject)fe)); 343 PetscCall(PetscFEDestroy(&fe)); 344 PetscCall(DMCreateDS(dmGrad)); 345 346 Vec locU, locG; 347 PetscCall(DMGetLocalVector(dm, &locU)); 348 PetscCall(DMGetLocalVector(dmGrad, &locG)); 349 PetscCall(DMGlobalToLocal(dm, U, INSERT_VALUES, locU)); 350 PetscCall(VecViewFromOptions(locU, NULL, "-remesh_input_view")); 351 PetscCall(DMPlexComputeGradientClementInterpolant(dm, locU, locG)); 352 PetscCall(VecViewFromOptions(locG, NULL, "-remesh_grad_view")); 353 354 const PetscScalar *g; 355 PetscScalar *u; 356 PetscInt n; 357 PetscCall(VecGetLocalSize(locU, &n)); 358 PetscCall(VecGetArrayWrite(locU, &u)); 359 PetscCall(VecGetArrayRead(locG, &g)); 360 for (PetscInt i = 0; i < n; i++) { 361 PetscReal norm = 0.0; 362 for (PetscInt d = 0; d < dim; d++) norm += PetscSqr(PetscRealPart(g[dim * i + d])); 363 u[i] = PetscSqrtReal(norm); 364 } 365 PetscCall(VecRestoreArrayRead(locG, &g)); 366 PetscCall(VecRestoreArrayWrite(locU, &u)); 367 368 DM dmM; 369 Vec metric; 370 PetscCall(DMClone(dm, &dmM)); 371 PetscCall(DMPlexMetricCreateIsotropic(dmM, 0, locU, &metric)); 372 PetscCall(DMDestroy(&dmM)); 373 PetscCall(DMRestoreLocalVector(dm, &locU)); 374 PetscCall(DMRestoreLocalVector(dmGrad, &locG)); 375 PetscCall(DMDestroy(&dmGrad)); 376 377 // TODO remove? 378 PetscScalar scale = 10.0; 379 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-metric_scale", &scale, NULL)); 380 PetscCall(VecScale(metric, scale)); 381 PetscCall(VecViewFromOptions(metric, NULL, "-metric_view")); 382 383 DMLabel label = NULL; 384 PetscCall(DMGetLabel(dm, "marker", &label)); 385 PetscCall(DMAdaptMetric(dm, metric, label, NULL, newdm)); 386 PetscCall(VecDestroy(&metric)); 387 388 PetscCall(DMViewFromOptions(*newdm, NULL, "-remesh_dmout_view")); 389 390 AppCtx *ctx; 391 PetscCall(DMGetApplicationContext(dm, &ctx)); 392 PetscCall(DMSetApplicationContext(*newdm, ctx)); 393 PetscCall(SetupDiscretization(*newdm, ctx)); 394 395 // TODO 396 ((DM_Plex *)(*newdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation; 397 PetscFunctionReturn(PETSC_SUCCESS); 398 } 399 400 static PetscErrorCode SetInitialConditions(TS ts, Vec u) 401 { 402 DM dm; 403 PetscReal t; 404 405 PetscFunctionBeginUser; 406 PetscCall(TSGetDM(ts, &dm)); 407 PetscCall(TSGetTime(ts, &t)); 408 PetscCall(DMComputeExactSolution(dm, t, u, NULL)); 409 PetscCall(VecSetOptionsPrefix(u, NULL)); 410 PetscFunctionReturn(PETSC_SUCCESS); 411 } 412 413 static PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec U, PetscBool *remesh, void *tctx) 414 { 415 AppCtx *ctx = (AppCtx *)tctx; 416 417 PetscFunctionBeginUser; 418 *remesh = PETSC_FALSE; 419 if (ctx->remesh_every > 0 && step && step % ctx->remesh_every == 0) { 420 DM dm; 421 422 *remesh = PETSC_TRUE; 423 PetscCall(TSGetDM(ts, &dm)); 424 PetscCall(Remesh(dm, U, &ctx->remesh_dm)); 425 } 426 PetscFunctionReturn(PETSC_SUCCESS); 427 } 428 429 static PetscErrorCode TransferVecs(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *tctx) 430 { 431 AppCtx *ctx = (AppCtx *)tctx; 432 DM dm; 433 Mat Interp; 434 435 PetscFunctionBeginUser; 436 PetscCall(TSGetDM(ts, &dm)); 437 PetscCall(DMCreateInterpolation(dm, ctx->remesh_dm, &Interp, NULL)); 438 for (PetscInt i = 0; i < nv; i++) { 439 PetscCall(DMCreateGlobalVector(ctx->remesh_dm, &vout[i])); 440 PetscCall(MatMult(Interp, vin[i], vout[i])); 441 } 442 PetscCall(MatDestroy(&Interp)); 443 PetscCall(TSSetDM(ts, ctx->remesh_dm)); 444 PetscCall(DMDestroy(&ctx->remesh_dm)); 445 PetscFunctionReturn(PETSC_SUCCESS); 446 } 447 448 int main(int argc, char **argv) 449 { 450 DM dm; 451 TS ts; 452 Vec u; 453 AppCtx ctx; 454 455 PetscFunctionBeginUser; 456 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 457 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &ctx)); 458 PetscCall(CreateMesh(PETSC_COMM_WORLD, &dm, &ctx)); 459 PetscCall(DMSetApplicationContext(dm, &ctx)); 460 PetscCall(SetupDiscretization(dm, &ctx)); 461 462 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 463 PetscCall(TSSetDM(ts, dm)); 464 PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx)); 465 if (ctx.expTS) { 466 PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFEM, &ctx)); 467 if (ctx.lumped) PetscCall(DMTSCreateRHSMassMatrixLumped(dm)); 468 else PetscCall(DMTSCreateRHSMassMatrix(dm)); 469 } else { 470 PetscCall(DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx)); 471 PetscCall(DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx)); 472 } 473 PetscCall(TSSetMaxTime(ts, 1.0)); 474 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 475 PetscCall(TSSetFromOptions(ts)); 476 PetscCall(TSSetComputeInitialCondition(ts, SetInitialConditions)); 477 PetscCall(TSSetResize(ts, PETSC_FALSE, TransferSetUp, TransferVecs, &ctx)); 478 479 PetscCall(DMCreateGlobalVector(dm, &u)); 480 PetscCall(DMTSCheckFromOptions(ts, u)); 481 PetscCall(SetInitialConditions(ts, u)); 482 PetscCall(PetscObjectSetName((PetscObject)u, "temperature")); 483 484 PetscCall(TSSetSolution(ts, u)); 485 PetscCall(VecViewFromOptions(u, NULL, "-u0_view")); 486 PetscCall(VecDestroy(&u)); 487 PetscCall(TSSolve(ts, NULL)); 488 489 PetscCall(TSGetSolution(ts, &u)); 490 PetscCall(VecViewFromOptions(u, NULL, "-uf_view")); 491 PetscCall(DMTSCheckFromOptions(ts, u)); 492 if (ctx.expTS) PetscCall(DMTSDestroyRHSMassMatrix(dm)); 493 494 PetscCall(TSDestroy(&ts)); 495 PetscCall(DMDestroy(&dm)); 496 PetscCall(PetscFinalize()); 497 return 0; 498 } 499 500 /*TEST 501 502 test: 503 suffix: 2d_p1 504 requires: triangle 505 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 506 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 507 test: 508 suffix: 2d_p1_exp 509 requires: triangle 510 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -explicit \ 511 -ts_type euler -ts_max_steps 4 -ts_dt 1e-3 -ts_monitor_error 512 test: 513 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 514 suffix: 2d_p1_sconv 515 requires: triangle 516 args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 517 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 518 test: 519 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.1] 520 suffix: 2d_p1_sconv_2 521 requires: triangle 522 args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 523 -ts_type beuler -ts_max_steps 1 -ts_dt 1e-6 -snes_error_if_not_converged -pc_type lu 524 test: 525 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 526 suffix: 2d_p1_tconv 527 requires: triangle 528 args: -sol_type quadratic_trig -temp_petscspace_degree 1 -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 test: 531 # -dm_refine 6 -convest_num_refine 3 get L_2 convergence rate: [1.0] 532 suffix: 2d_p1_tconv_2 533 requires: triangle 534 args: -sol_type trig_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 535 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 536 test: 537 # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid 538 suffix: 2d_p1_exp_tconv_2 539 requires: triangle 540 args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \ 541 -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 -lumped 0 -mass_pc_type lu 542 test: 543 # The L_2 convergence rate cannot be seen since stability of the explicit integrator requires that is be more accurate than the grid 544 suffix: 2d_p1_exp_tconv_2_lump 545 requires: triangle 546 args: -sol_type trig_trig -temp_petscspace_degree 1 -explicit -ts_convergence_estimate -convest_num_refine 1 \ 547 -ts_type euler -ts_max_steps 4 -ts_dt 1e-4 548 test: 549 suffix: 2d_p2 550 requires: triangle 551 args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 552 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 553 test: 554 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 555 suffix: 2d_p2_sconv 556 requires: triangle 557 args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 558 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 559 test: 560 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [3.1] 561 suffix: 2d_p2_sconv_2 562 requires: triangle 563 args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 564 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 565 test: 566 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 567 suffix: 2d_p2_tconv 568 requires: triangle 569 args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 570 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 571 test: 572 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 573 suffix: 2d_p2_tconv_2 574 requires: triangle 575 args: -sol_type trig_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 576 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 577 test: 578 suffix: 2d_q1 579 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 580 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 581 test: 582 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 583 suffix: 2d_q1_sconv 584 args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 585 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 586 test: 587 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 588 suffix: 2d_q1_tconv 589 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 590 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 591 test: 592 suffix: 2d_q2 593 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 594 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 595 test: 596 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 597 suffix: 2d_q2_sconv 598 args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 599 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 600 test: 601 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 602 suffix: 2d_q2_tconv 603 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 604 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 605 606 test: 607 suffix: 3d_p1 608 requires: ctetgen 609 args: -sol_type quadratic_linear -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 610 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 611 test: 612 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 613 suffix: 3d_p1_sconv 614 requires: ctetgen 615 args: -sol_type quadratic_linear -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 616 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 617 test: 618 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 619 suffix: 3d_p1_tconv 620 requires: ctetgen 621 args: -sol_type quadratic_trig -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 622 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 623 test: 624 suffix: 3d_p2 625 requires: ctetgen 626 args: -sol_type quadratic_linear -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 627 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 628 test: 629 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 630 suffix: 3d_p2_sconv 631 requires: ctetgen 632 args: -sol_type trig_linear -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 633 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 634 test: 635 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 636 suffix: 3d_p2_tconv 637 requires: ctetgen 638 args: -sol_type quadratic_trig -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 639 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 640 test: 641 suffix: 3d_q1 642 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -dmts_check .0001 \ 643 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 644 test: 645 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [1.9] 646 suffix: 3d_q1_sconv 647 args: -sol_type quadratic_linear -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 648 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00001 -snes_error_if_not_converged -pc_type lu 649 test: 650 # -dm_refine 4 -convest_num_refine 3 get L_2 convergence rate: [1.2] 651 suffix: 3d_q1_tconv 652 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 1 -ts_convergence_estimate -convest_num_refine 1 \ 653 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 654 test: 655 suffix: 3d_q2 656 args: -sol_type quadratic_linear -dm_plex_simplex 0 -dm_refine 0 -temp_petscspace_degree 2 -dmts_check .0001 \ 657 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 658 test: 659 # -dm_refine 2 -convest_num_refine 3 get L_2 convergence rate: [2.9] 660 suffix: 3d_q2_sconv 661 args: -sol_type trig_linear -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -ts_convergence_temporal 0 -convest_num_refine 1 \ 662 -ts_type beuler -ts_max_steps 1 -ts_dt 0.00000001 -snes_error_if_not_converged -pc_type lu 663 test: 664 # -dm_refine 3 -convest_num_refine 3 get L_2 convergence rate: [1.0] 665 suffix: 3d_q2_tconv 666 args: -sol_type quadratic_trig -dm_plex_simplex 0 -temp_petscspace_degree 2 -ts_convergence_estimate -convest_num_refine 1 \ 667 -ts_type beuler -ts_max_steps 4 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 668 669 test: 670 # 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 671 suffix: egads_sphere 672 requires: egads ctetgen 673 args: -sol_type quadratic_linear \ 674 -dm_plex_boundary_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sphere_example.egadslite -dm_plex_boundary_label marker \ 675 -temp_petscspace_degree 2 -dmts_check .0001 \ 676 -ts_type beuler -ts_max_steps 5 -ts_dt 0.1 -snes_error_if_not_converged -pc_type lu 677 678 test: 679 suffix: remesh 680 requires: triangle mmg 681 args: -sol_type quadratic_trig -dm_refine 2 -temp_petscspace_degree 1 -ts_type beuler -ts_dt 0.01 -snes_error_if_not_converged -pc_type lu -grad_petscspace_degree 1 -dm_adaptor mmg -dm_plex_hash_location -remesh_every 5 682 683 TEST*/ 684