Actual source code: ex52.c

petsc-master 2020-07-10
Report Typos and Errors
  1: static char help[] = "Simple Advection-diffusion equation solved using FVM in DMPLEX\n";

  3: /*
  4:    Solves the simple advection equation given by

  6:    q_t + u (q_x) + v (q_y) - D (q_xx + q_yy) = 0 using FVM and First Order Upwind discretization.

  8:    with a user defined initial condition.

 10:    with dirichlet/neumann conditions on the four boundaries of the domain.

 12:    User can define the mesh parameters either in the command line or inside
 13:    the ProcessOptions() routine.

 15:    Contributed by: Mukkund Sunjii, Domenico Lahaye
 16: */

 18:  #include <petscdmplex.h>
 19:  #include <petscts.h>
 20:  #include <petscblaslapack.h>

 22: #if defined(PETSC_HAVE_CGNS)
 23: #undef I
 24: #include <cgnslib.h>
 25: #endif
 26: /*
 27:    User-defined routines
 28: */
 29: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec);
 30: extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *);
 31: extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *);

 33: /* Defining the usr defined context */
 34: typedef struct {
 35:     DM da;
 36:     PetscBool interpolate;                  /* Generate intermediate mesh elements */
 37:     char filename[PETSC_MAX_PATH_LEN]; /* Mesh filename */
 38:     PetscInt dim;
 39:     PetscScalar diffusion;
 40:     PetscReal u, v;
 41:     PetscScalar delta_x, delta_y;
 42:     PetscInt cells[2];
 43: } AppCtx;

 45: /* Options for the scenario */
 46: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {

 50:     options->interpolate = PETSC_TRUE;
 51:     options->filename[0] = '\0';
 52:     options->dim = 2;
 53:     options->u = 2.5;
 54:     options->v = 0.0;
 55:     options->cells[0] = 20;
 56:     options->cells[1] = 20;
 57:     options->diffusion = 0.0;

 59:     PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
 60:     PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "advection_DMPLEX.c",options->interpolate, &options->interpolate, NULL);
 61:     PetscOptionsString("-filename", "The mesh file", "advection_DMPLEX.c", options->filename, options->filename, sizeof(options->filename), NULL);
 62:     PetscOptionsRangeInt("-dim", "Problem dimension used for the non-file mesh.", "ex7.c", options->dim, &options->dim, NULL,1,3);
 63:     PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL);
 64:     PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL);
 65:     PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL);
 66:     PetscOptionsEnd();
 67:     return(0);
 68: }

 70: /*
 71:   User can provide the file containing the mesh.
 72:   Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
 73: */
 74: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm) {
 75:     size_t len;
 78:     PetscStrlen(user->filename, &len);
 79:     if (!len) {
 80:         DMLabel label;
 81:         DMPlexCreateBoxMesh(comm, user->dim, PETSC_FALSE, user->cells, NULL, NULL, NULL, user->interpolate, dm);
 82:         /* Mark boundary and set BC */
 83:         DMCreateLabel(*dm, "boundary");
 84:         DMGetLabel(*dm, "boundary", &label);
 85:         DMPlexMarkBoundaryFaces(*dm, 1, label);
 86:         DMPlexLabelComplete(*dm, label);
 87:     } else {
 88:         DMPlexCreateFromFile(comm, user->filename, user->interpolate, dm);
 89:     }
 90:     PetscObjectSetName((PetscObject) * dm, "Mesh");
 91:     DMSetFromOptions(*dm);
 92:     DMViewFromOptions(*dm, NULL, "-dm_view");
 93:     return(0);
 94: }

 96:     /* This routine is responsible for defining the local solution vector x
 97:     with a given initial solution.
 98:     The initial solution can be modified accordingly inside the loops.
 99:     No need for a local vector because there is exchange of information
100:     across the processors. Unlike for FormFunction which depends on the neighbours */
101: PetscErrorCode FormInitialSolution(DM da, Vec U) {
103:     PetscScalar *u;
104:     PetscInt cell, cStart, cEnd;
105:     PetscReal cellvol, centroid[3], normal[3];

108:     /* Get pointers to vector data */
109:     VecGetArray(U, &u);
110:     /* Get local grid boundaries */
111:     DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
112:     /* Assigning the values at the cell centers based on x and y directions */
113:     for (cell = cStart; cell < cEnd; cell++) {
114:         DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal);
115:         if (centroid[0] > 0.9 && centroid[0] < 0.95) {
116:             if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0;
117:         }
118:         else u[cell] = 0;
119:     }
120:     VecRestoreArray(U, &u);
121:     return(0);
122: }

124: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx) {
126:     PetscReal norm;
127:     MPI_Comm comm;
129:     if (step < 0) return(0); /* step of -1 indicates an interpolated solution */
130:     VecNorm(v, NORM_2, &norm);
131:     PetscObjectGetComm((PetscObject) ts, &comm);
132:     PetscPrintf(comm, "timestep %D time %g norm %g\n", step, (double) ptime, (double) norm);
133:     return(0);
134: }

136: /*
137:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
138:    Input Parameters:
139:      snes - the SNES context
140:      its - iteration number
141:      fnorm - 2-norm function value (may be estimated)
142:      ctx - optional user-defined context for private data for the
143:          monitor routine, as set by SNESMonitorSet()
144: */
145: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf) {
148:     SNESMonitorDefaultShort(snes, its, fnorm, vf);
149:     return(0);
150: }

152: /*
153:    FormFunction - Evaluates nonlinear function, F(x).

155:    Input Parameters:
156: .  ts - the TS context
157: .  X - input vector
158: .  ctx - optional user-defined context, as set by SNESSetFunction()

160:    Output Parameter:
161: .  F - function vector
162:  */
163: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ctx) {
164:     AppCtx *user = (AppCtx *) ctx;
165:     DM da = (DM) user->da;
167:     PetscScalar *x, *f;
168:     Vec localX;
169:     PetscInt fStart, fEnd, nF;
170:     PetscInt cell, cStart, cEnd, nC;
171:     DM dmFace;      /* DMPLEX for face geometry */
172:     PetscFV fvm;                /* specify type of FVM discretization */
173:     Vec cellGeom, faceGeom; /* vector of structs related to cell/face geometry*/
174:     const PetscScalar *fgeom;             /* values stored in the vector facegeom */
175:     PetscFVFaceGeom *fgA;               /* struct with face geometry information */
176:     const PetscInt *cellcone, *cellsupport;
177:     PetscScalar flux_east, flux_west, flux_north, flux_south, flux_centre;
178:     PetscScalar centroid_x[2], centroid_y[2], boundary = 0.0;
179:     PetscScalar boundary_left = 0.0;
180:     PetscReal u_plus, u_minus, v_plus, v_minus, zero = 0.0;
181:     PetscScalar delta_x, delta_y;

183:     /* Get the local vector from the DM object. */
185:     DMGetLocalVector(da, &localX);

187:     /* Scatter ghost points to local vector,using the 2-step process
188:        DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
189:     DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
190:     DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
191:     /* Get pointers to vector data. */
192:     VecGetArray(localX, &x);
193:     VecGetArray(F, &f);

195:     /* Obtaining local cell and face ownership */
196:     DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
197:     DMPlexGetHeightStratum(da, 1, &fStart, &fEnd);

199:     /* Creating the PetscFV object to obtain face and cell geometry.
200:     Later to be used to compute face centroid to find cell widths. */

202:     PetscFVCreate(PETSC_COMM_WORLD, &fvm);
203:     PetscFVSetType(fvm, PETSCFVUPWIND);
204:     /*....Retrieve precomputed cell geometry....*/
205:     DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL);
206:     VecGetDM(faceGeom, &dmFace);
207:     VecGetArrayRead(faceGeom, &fgeom);

209:     /* Spanning through all the cells and an inner loop through the faces. Find the
210:     face neighbors and pick the upwinded cell value for flux. */

212:     u_plus = PetscMax(user->u, zero);
213:     u_minus = PetscMin(user->u, zero);
214:     v_plus = PetscMax(user->v, zero);
215:     v_minus = PetscMin(user->v, zero);

217:     for (cell = cStart; cell < cEnd; cell++) {
218:         /* Obtaining the faces of the cell */
219:         DMPlexGetConeSize(da, cell, &nF);
220:         DMPlexGetCone(da, cell, &cellcone);

222:         /* south */
223:         DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA);
224:         centroid_y[0] = fgA->centroid[1];
225:         /* North */
226:         DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA);
227:         centroid_y[1] = fgA->centroid[1];
228:         /* West */
229:         DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA);
230:         centroid_x[0] = fgA->centroid[0];
231:         /* East */
232:         DMPlexPointLocalRead(dmFace, cellcone[1], fgeom, &fgA);
233:         centroid_x[1] = fgA->centroid[0];

235:         /* Computing the cell widths in the x and y direction */
236:         delta_x = centroid_x[1] - centroid_x[0];
237:         delta_y = centroid_y[1] - centroid_y[0];

239:         /* Getting the neighbors of each face
240:            Going through the faces by the order (cellcone) */

242:         /* cellcone[0] - south */
243:         DMPlexGetSupportSize(da, cellcone[0], &nC);
244:         DMPlexGetSupport(da, cellcone[0], &cellsupport);
245:         if (nC == 2) flux_south = (x[cellsupport[0]] * (-v_plus - user->diffusion * delta_x)) / delta_y;
246:         else flux_south = (boundary * (-v_plus - user->diffusion * delta_x)) / delta_y;

248:         /* cellcone[1] - east */
249:         DMPlexGetSupportSize(da, cellcone[1], &nC);
250:         DMPlexGetSupport(da, cellcone[1], &cellsupport);
251:         if (nC == 2) flux_east = (x[cellsupport[1]] * (u_minus - user->diffusion * delta_y)) / delta_x;
252:         else flux_east = (boundary * (u_minus - user->diffusion * delta_y)) / delta_x;

254:         /* cellcone[2] - north */
255:         DMPlexGetSupportSize(da, cellcone[2], &nC);
256:         DMPlexGetSupport(da, cellcone[2], &cellsupport);
257:         if (nC == 2) flux_north = (x[cellsupport[1]] * (v_minus - user->diffusion * delta_x)) / delta_y;
258:         else flux_north = (boundary * (v_minus - user->diffusion * delta_x)) / delta_y;

260:         /* cellcone[3] - west */
261:         DMPlexGetSupportSize(da, cellcone[3], &nC);
262:         DMPlexGetSupport(da, cellcone[3], &cellsupport);
263:         if (nC == 2) flux_west = (x[cellsupport[0]] * (-u_plus - user->diffusion * delta_y)) / delta_x;
264:         else flux_west = (boundary_left * (-u_plus - user->diffusion * delta_y)) / delta_x;

266:         /* Contribution by the cell to the fluxes */
267:         flux_centre = x[cell] * ((u_plus - u_minus + 2 * user->diffusion * delta_y) / delta_x +
268:                                  (v_plus - v_minus + 2 * user->diffusion * delta_x) / delta_y);

270:         /* Calculating the net flux for each cell
271:            and computing the RHS time derivative f[.] */
272:         f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south);

274:     }

276:     PetscFVDestroy(&fvm);
277:     VecRestoreArray(localX, &x);
278:     VecRestoreArray(F, &f);
279:     DMRestoreLocalVector(da, &localX);

281:     return(0);
282: }

284: int main(int argc, char **argv) {
285:     TS ts;                         /* time integrator */
286:     SNES snes;
287:     Vec x, r;                        /* solution, residual vectors */
289:     DM da;
290:     PetscMPIInt rank;
291:     PetscViewerAndFormat *vf;
292:     AppCtx user;                             /* mesh context */
293:     PetscInt numFields = 1, numBC, i;
294:     PetscInt numComp[1];
295:     PetscInt numDof[12];
296:     PetscInt bcField[1];
297:     PetscSection section;
298:     IS bcPointIS[1];

300:     /* Initialize program */
301:     PetscInitialize(&argc, &argv, (char *) 0, help);
302:     if (ierr) return ierr;
303:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
304:     /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
305:     ProcessOptions(PETSC_COMM_WORLD, &user);
306:     CreateMesh(PETSC_COMM_WORLD, &user, &da);

308:     /* Specifying the fields and dof for the formula through PETSc Section
309:     Create a scalar field u with 1 component on cells, faces and edges.
310:     Alternatively, the field information could be added through a PETSCFV object
311:     using DMAddField(...).*/
312:     numComp[0] = 1;

314:     for (i = 0; i < numFields * (user.dim + 1); ++i) numDof[i] = 0;

316:     numDof[0 * (user.dim + 1)] = 1;
317:     numDof[0 * (user.dim + 1) + user.dim - 1] = 1;
318:     numDof[0 * (user.dim + 1) + user.dim] = 1;

320:     /* Setup boundary conditions */
321:     numBC = 1;
322:     /* Prescribe a Dirichlet condition on u on the boundary
323:        Label "marker" is made by the mesh creation routine  */
324:     bcField[0] = 0;
325:     DMGetStratumIS(da, "marker", 1, &bcPointIS[0]);

327:     /* Create a PetscSection with this data layout */
328:     DMSetNumFields(da, numFields);
329:     DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, &section);

331:     /* Name the Field variables */
332:     PetscSectionSetFieldName(section, 0, "u");

334:     /* Tell the DM to use this section (with the specified fields and dof) */
335:     DMSetLocalSection(da, section);
336:     user.da = da;

338:     /* Extract global vectors from DMDA; then duplicate for remaining
339:        vectors that are the same types */

341:     /* Create a Vec with this layout and view it */
342:     DMGetGlobalVector(da, &x);
343:     VecDuplicate(x, &r);

345:     /* Create timestepping solver context */
346:     TSCreate(PETSC_COMM_WORLD, &ts);
347:     TSSetProblemType(ts, TS_NONLINEAR);
348:     TSSetRHSFunction(ts, NULL, FormFunction, &user);

350:     TSSetMaxTime(ts, 1.0);
351:     TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
352:     TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL);
353:     TSSetDM(ts, da);

355:     /* Customize nonlinear solver */
356:     TSSetType(ts, TSEULER);
357:     TSGetSNES(ts, &snes);
358:     PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf);
359:     SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy);

361:      /* Set initial conditions */
362:     FormInitialSolution(da, x);
363:     TSSetTimeStep(ts, .0001);
364:     TSSetSolution(ts, x);
365:     /* Set runtime options */
366:     TSSetFromOptions(ts);
367:     /* Solve nonlinear system */
368:     TSSolve(ts, x);

370:     /* Clean up routine */
371:     DMRestoreGlobalVector(da, &x);
372:     ISDestroy(&bcPointIS[0]);
373:     PetscSectionDestroy(&section);
374:     VecDestroy(&r);
375:     TSDestroy(&ts);
376:     DMDestroy(&da);
377:     PetscFinalize();
378:     return ierr;
379: }

381: /*TEST

383:     test:
384:       suffix: 0
385:       args: -ts_max_steps 5 -ts_type rk
386:       requires: !single !complex triangle ctetgen

388: TEST*/