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