Actual source code: ex45.c
petsc-3.14.6 2021-03-30
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 = -1 * dim
15: Exact 2D solution:
17: u = 2t + x^2 + y^2
19: 2 - (2 + 2) + 2 = 0
21: Exact 3D solution:
23: u = 3t + x^2 + y^2 + z^2
25: 3 - (2 + 2 + 2) + 3 = 0
26: */
28: typedef struct {
29: PetscInt dim;
30: PetscBool simplex;
31: PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
32: } AppCtx;
34: static PetscErrorCode analytic_temp(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
35: {
36: PetscInt d;
38: *u = dim*time;
39: for (d = 0; d < dim; ++d) *u += x[d]*x[d];
40: return 0;
41: }
43: static void f0_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
44: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
45: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
46: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
47: {
48: f0[0] = u_t[0] + (PetscScalar) dim;
49: }
51: static void f1_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
52: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
53: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
54: PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
55: {
56: PetscInt d;
57: for (d = 0; d < dim; ++d) {
58: f1[d] = u_x[d];
59: }
60: }
62: static void g3_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
63: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
64: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
65: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
66: {
67: PetscInt d;
68: for (d = 0; d < dim; ++d) {
69: g3[d*dim+d] = 1.0;
70: }
71: }
73: static void g0_temp(PetscInt dim, PetscInt Nf, PetscInt NfAux,
74: const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
75: const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
76: PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])
77: {
78: g0[0] = u_tShift*1.0;
79: }
81: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
82: {
86: options->dim = 2;
87: options->simplex = PETSC_TRUE;
89: PetscOptionsBegin(comm, "", "Heat Equation Options", "DMPLEX");
90: PetscOptionsInt("-dim", "The topological mesh dimension", "ex45.c", options->dim, &options->dim, NULL);
91: PetscOptionsBool("-simplex", "Simplicial (true) or tensor (false) mesh", "ex45.c", options->simplex, &options->simplex, NULL);
92: PetscOptionsEnd();
93: return(0);
94: }
96: static PetscErrorCode CreateBCLabel(DM dm, const char name[])
97: {
98: DM plex;
99: DMLabel label;
103: DMCreateLabel(dm, name);
104: DMGetLabel(dm, name, &label);
105: DMConvert(dm, DMPLEX, &plex);
106: DMPlexMarkBoundaryFaces(plex, 1, label);
107: DMDestroy(&plex);
108: return(0);
109: }
111: static PetscErrorCode CreateMesh(MPI_Comm comm, DM *dm, AppCtx *ctx)
112: {
113: DM pdm = NULL;
114: const PetscInt dim = ctx->dim;
115: PetscBool hasLabel;
119: DMPlexCreateBoxMesh(comm, dim, ctx->simplex, NULL, NULL, NULL, NULL, PETSC_TRUE, dm);
120: PetscObjectSetName((PetscObject) *dm, "Mesh");
121: /* If no boundary marker exists, mark the whole boundary */
122: DMHasLabel(*dm, "marker", &hasLabel);
123: if (!hasLabel) {CreateBCLabel(*dm, "marker");}
124: /* Distribute mesh over processes */
125: DMPlexDistribute(*dm, 0, NULL, &pdm);
126: if (pdm) {
127: DMDestroy(dm);
128: *dm = pdm;
129: }
130: DMSetFromOptions(*dm);
131: DMViewFromOptions(*dm, NULL, "-dm_view");
132: return(0);
133: }
135: static PetscErrorCode SetupProblem(DM dm, AppCtx *ctx)
136: {
137: PetscDS prob;
138: const PetscInt id = 1;
142: DMGetDS(dm, &prob);
143: PetscDSSetResidual(prob, 0, f0_temp, f1_temp);
144: PetscDSSetJacobian(prob, 0, 0, g0_temp, NULL, NULL, g3_temp);
145: ctx->exactFuncs[0] = analytic_temp;
146: DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", "marker", 0, 0, NULL, (void (*)(void)) ctx->exactFuncs[0], NULL, 1, &id, ctx);
147: return(0);
148: }
150: static PetscErrorCode SetupDiscretization(DM dm, AppCtx* ctx)
151: {
152: DM cdm = dm;
153: const PetscInt dim = ctx->dim;
154: PetscFE fe;
158: /* Create finite element */
159: PetscFECreateDefault(PetscObjectComm((PetscObject) dm), dim, 1, ctx->simplex, "temp_", -1, &fe);
160: PetscObjectSetName((PetscObject) fe, "temperature");
161: /* Set discretization and boundary conditions for each mesh */
162: DMSetField(dm, 0, NULL, (PetscObject) fe);
163: DMCreateDS(dm);
164: SetupProblem(dm, ctx);
165: while (cdm) {
166: PetscBool hasLabel;
168: DMHasLabel(cdm, "marker", &hasLabel);
169: if (!hasLabel) {CreateBCLabel(cdm, "marker");}
170: DMCopyDisc(dm, cdm);
171: DMGetCoarseDM(cdm, &cdm);
172: }
173: PetscFEDestroy(&fe);
174: return(0);
175: }
177: int main(int argc, char **argv)
178: {
179: AppCtx ctx;
180: DM dm;
181: TS ts;
182: Vec u, r;
183: PetscReal t = 0.0;
184: PetscReal L2error = 0.0;
187: PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr;
188: ProcessOptions(PETSC_COMM_WORLD, &ctx);
189: CreateMesh(PETSC_COMM_WORLD, &dm, &ctx);
190: DMSetApplicationContext(dm, &ctx);
191: PetscMalloc1(1, &ctx.exactFuncs);
192: SetupDiscretization(dm, &ctx);
194: DMCreateGlobalVector(dm, &u);
195: PetscObjectSetName((PetscObject) u, "temperature");
196: VecDuplicate(u, &r);
198: TSCreate(PETSC_COMM_WORLD, &ts);
199: TSSetDM(ts, dm);
200: DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, &ctx);
201: DMTSSetIFunctionLocal(dm, DMPlexTSComputeIFunctionFEM, &ctx);
202: DMTSSetIJacobianLocal(dm, DMPlexTSComputeIJacobianFEM, &ctx);
203: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
204: TSSetFromOptions(ts);
206: DMProjectFunction(dm, t, ctx.exactFuncs, NULL, INSERT_ALL_VALUES, u);
207: TSSolve(ts, u);
209: TSGetTime(ts, &t);
210: DMComputeL2Diff(dm, t, ctx.exactFuncs, NULL, u, &L2error);
211: if (L2error < 1.0e-11) {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");}
212: else {PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)L2error);}
213: VecViewFromOptions(u, NULL, "-sol_vec_view");
215: VecDestroy(&u);
216: VecDestroy(&r);
217: TSDestroy(&ts);
218: DMDestroy(&dm);
219: PetscFree(ctx.exactFuncs);
220: PetscFinalize();
221: return ierr;
222: }
224: /*TEST
226: # Full solves
227: test:
228: suffix: 2d_p1_r1
229: requires: triangle
230: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
231: args: -dm_refine 1 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
232: test:
233: suffix: 2d_p1_r3
234: requires: triangle
235: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
236: args: -dm_refine 3 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
237: test:
238: suffix: 2d_p2_r1
239: requires: triangle
240: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
241: args: -dm_refine 1 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
242: test:
243: suffix: 2d_p2_r3
244: requires: triangle
245: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
246: args: -dm_refine 3 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
247: test:
248: suffix: 2d_q1_r1
249: requires: !single
250: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
251: args: -simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
252: test:
253: suffix: 2d_q1_r3
254: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
255: args: -simplex 0 -dm_refine 3 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
256: test:
257: suffix: 2d_q2_r1
258: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
259: args: -simplex 0 -dm_refine 1 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
260: test:
261: suffix: 2d_q2_r3
262: requires: !single
263: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
264: args: -simplex 0 -dm_refine 3 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
265: test:
266: suffix: 3d_p1_r1
267: requires: ctetgen
268: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
269: args: -dim 3 -dm_refine 1 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
270: test:
271: suffix: 3d_p1_r2
272: requires: ctetgen
273: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
274: args: -dim 3 -dm_refine 2 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
275: test:
276: suffix: 3d_p2_r1
277: requires: ctetgen
278: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
279: args: -dim 3 -dm_refine 1 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
280: test:
281: suffix: 3d_p2_r2
282: requires: ctetgen
283: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
284: args: -dim 3 -dm_refine 2 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
285: test:
286: suffix: 3d_q1_r1
287: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
288: args: -dim 3 -simplex 0 -dm_refine 1 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
289: test:
290: suffix: 3d_q1_r2
291: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
292: args: -dim 3 -simplex 0 -dm_refine 2 -temp_petscspace_degree 1 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
293: test:
294: suffix: 3d_q2_r1
295: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
296: args: -dim 3 -simplex 0 -dm_refine 1 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
297: test:
298: suffix: 3d_q2_r2
299: filter: sed -e "s~ATOL~RTOL~g" -e "s~ABS~RELATIVE~g"
300: args: -dim 3 -simplex 0 -dm_refine 2 -temp_petscspace_degree 2 -ts_type beuler -ts_max_steps 10 -ts_dt 0.1 -pc_type lu -ksp_monitor_short -ksp_converged_reason -snes_monitor_short -snes_converged_reason -ts_monitor
302: TEST*/