Actual source code: ex45.c

petsc-3.14.4 2021-02-03
Report Typos and Errors
  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*/