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, void *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, void *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, void *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, void *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, void *ctx)
117: {
118: PetscInt d;
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: }
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: }
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: }
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
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;
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: }
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: }
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: }
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: }
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: }
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: }
200: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
201: {
202: PetscInt sol;
204: PetscFunctionBeginUser;
205: options->solType = SOL_QUADRATIC_LINEAR;
206: options->expTS = PETSC_FALSE;
207: options->lumped = PETSC_TRUE;
208: options->remesh_every = 0;
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: }
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: }
248: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
249: {
250: PetscDS ds;
251: DMLabel label;
252: const PetscInt id = 1;
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: }
291: static PetscErrorCode SetupDiscretization(DM dm, AppCtx *ctx)
292: {
293: DM plex, cdm = dm;
294: PetscFE fe;
295: PetscBool simplex;
296: PetscInt dim;
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;
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: }
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;
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));
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));
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"));
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));
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));
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"));
383: DMLabel label = NULL;
384: PetscCall(DMGetLabel(dm, "marker", &label));
385: PetscCall(DMAdaptMetric(dm, metric, label, NULL, newdm));
386: PetscCall(VecDestroy(&metric));
388: PetscCall(DMViewFromOptions(*newdm, NULL, "-remesh_dmout_view"));
390: AppCtx *ctx;
391: PetscCall(DMGetApplicationContext(dm, &ctx));
392: PetscCall(DMSetApplicationContext(*newdm, ctx));
393: PetscCall(SetupDiscretization(*newdm, ctx));
395: // TODO
396: ((DM_Plex *)(*newdm)->data)->useHashLocation = ((DM_Plex *)dm->data)->useHashLocation;
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode SetInitialConditions(TS ts, Vec u)
401: {
402: DM dm;
403: PetscReal t;
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: }
413: static PetscErrorCode TransferSetUp(TS ts, PetscInt step, PetscReal time, Vec U, PetscBool *remesh, void *tctx)
414: {
415: AppCtx *ctx = (AppCtx *)tctx;
417: PetscFunctionBeginUser;
418: *remesh = PETSC_FALSE;
419: if (ctx->remesh_every > 0 && step && step % ctx->remesh_every == 0) {
420: DM dm;
422: *remesh = PETSC_TRUE;
423: PetscCall(TSGetDM(ts, &dm));
424: PetscCall(Remesh(dm, U, &ctx->remesh_dm));
425: }
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
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;
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: }
448: int main(int argc, char **argv)
449: {
450: DM dm;
451: TS ts;
452: Vec u;
453: AppCtx ctx;
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));
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));
479: PetscCall(DMCreateGlobalVector(dm, &u));
480: PetscCall(DMTSCheckFromOptions(ts, u));
481: PetscCall(SetInitialConditions(ts, u));
482: PetscCall(PetscObjectSetName((PetscObject)u, "temperature"));
484: PetscCall(TSSetSolution(ts, u));
485: PetscCall(VecViewFromOptions(u, NULL, "-u0_view"));
486: PetscCall(VecDestroy(&u));
487: PetscCall(TSSolve(ts, NULL));
489: PetscCall(TSGetSolution(ts, &u));
490: PetscCall(VecViewFromOptions(u, NULL, "-uf_view"));
491: PetscCall(DMTSCheckFromOptions(ts, u));
492: if (ctx.expTS) PetscCall(DMTSDestroyRHSMassMatrix(dm));
494: PetscCall(TSDestroy(&ts));
495: PetscCall(DMDestroy(&dm));
496: PetscCall(PetscFinalize());
497: return 0;
498: }
500: /*TEST
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
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
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
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
683: TEST*/