Actual source code: ex52.c

  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:   PetscScalar diffusion;
 36:   PetscReal   u, v;
 37:   PetscScalar delta_x, delta_y;
 38: } AppCtx;

 40: /* Options for the scenario */
 41: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 42: {
 44:   options->u         = 2.5;
 45:   options->v         = 0.0;
 46:   options->diffusion = 0.0;
 47:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
 48:   PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL);
 49:   PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL);
 50:   PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL);
 51:   PetscOptionsEnd();
 52:   return 0;
 53: }

 55: /*
 56:   User can provide the file containing the mesh.
 57:   Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
 58: */
 59: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 60: {
 62:   DMCreate(comm, dm);
 63:   DMSetType(*dm, DMPLEX);
 64:   DMSetFromOptions(*dm);
 65:   DMViewFromOptions(*dm, NULL, "-dm_view");
 66:   {
 67:     DMLabel label;
 68:     DMGetLabel(*dm, "boundary", &label);
 69:     DMPlexLabelComplete(*dm, label);
 70:   }
 71:   return 0;
 72: }

 74: /* This routine is responsible for defining the local solution vector x
 75:     with a given initial solution.
 76:     The initial solution can be modified accordingly inside the loops.
 77:     No need for a local vector because there is exchange of information
 78:     across the processors. Unlike for FormFunction which depends on the neighbours */
 79: PetscErrorCode FormInitialSolution(DM da, Vec U)
 80: {
 81:   PetscScalar *u;
 82:   PetscInt     cell, cStart, cEnd;
 83:   PetscReal    cellvol, centroid[3], normal[3];

 86:   /* Get pointers to vector data */
 87:   VecGetArray(U, &u);
 88:   /* Get local grid boundaries */
 89:   DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
 90:   /* Assigning the values at the cell centers based on x and y directions */
 91:   DMGetCoordinatesLocalSetUp(da);
 92:   for (cell = cStart; cell < cEnd; cell++) {
 93:     DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal);
 94:     if (centroid[0] > 0.9 && centroid[0] < 0.95) {
 95:       if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0;
 96:     } else u[cell] = 0;
 97:   }
 98:   VecRestoreArray(U, &u);
 99:   return 0;
100: }

102: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
103: {
104:   PetscReal norm;
105:   MPI_Comm  comm;

108:   if (step < 0) return 0; /* step of -1 indicates an interpolated solution */
109:   VecNorm(v, NORM_2, &norm);
110:   PetscObjectGetComm((PetscObject)ts, &comm);
111:   PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm);
112:   return 0;
113: }

115: /*
116:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
117:    Input Parameters:
118:      snes - the SNES context
119:      its - iteration number
120:      fnorm - 2-norm function value (may be estimated)
121:      ctx - optional user-defined context for private data for the
122:          monitor routine, as set by SNESMonitorSet()
123: */
124: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
125: {
127:   SNESMonitorDefaultShort(snes, its, fnorm, vf);
128:   return 0;
129: }

131: /*
132:    FormFunction - Evaluates nonlinear function, F(x).

134:    Input Parameters:
135: .  ts - the TS context
136: .  X - input vector
137: .  ctx - optional user-defined context, as set by SNESSetFunction()

139:    Output Parameter:
140: .  F - function vector
141:  */
142: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ctx)
143: {
144:   AppCtx            *user = (AppCtx *)ctx;
145:   DM                 da;
146:   PetscScalar       *x, *f;
147:   Vec                localX;
148:   PetscInt           fStart, fEnd, nF;
149:   PetscInt           cell, cStart, cEnd, nC;
150:   DM                 dmFace;             /* DMPLEX for face geometry */
151:   PetscFV            fvm;                /* specify type of FVM discretization */
152:   Vec                cellGeom, faceGeom; /* vector of structs related to cell/face geometry*/
153:   const PetscScalar *fgeom;              /* values stored in the vector facegeom */
154:   PetscFVFaceGeom   *fgA;                /* struct with face geometry information */
155:   const PetscInt    *cellcone, *cellsupport;
156:   PetscScalar        flux_east, flux_west, flux_north, flux_south, flux_centre;
157:   PetscScalar        centroid_x[2], centroid_y[2], boundary = 0.0;
158:   PetscScalar        boundary_left = 0.0;
159:   PetscReal          u_plus, u_minus, v_plus, v_minus, zero = 0.0;
160:   PetscScalar        delta_x, delta_y;

162:   /* Get the local vector from the DM object. */
164:   TSGetDM(ts, &da);
165:   DMGetLocalVector(da, &localX);

167:   /* Scatter ghost points to local vector,using the 2-step process
168:        DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
169:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
170:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
171:   /* Get pointers to vector data. */
172:   VecGetArray(localX, &x);
173:   VecGetArray(F, &f);

175:   /* Obtaining local cell and face ownership */
176:   DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
177:   DMPlexGetHeightStratum(da, 1, &fStart, &fEnd);

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

182:   PetscFVCreate(PETSC_COMM_WORLD, &fvm);
183:   PetscFVSetType(fvm, PETSCFVUPWIND);
184:   /*....Retrieve precomputed cell geometry....*/
185:   DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL);
186:   VecGetDM(faceGeom, &dmFace);
187:   VecGetArrayRead(faceGeom, &fgeom);

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

192:   u_plus  = PetscMax(user->u, zero);
193:   u_minus = PetscMin(user->u, zero);
194:   v_plus  = PetscMax(user->v, zero);
195:   v_minus = PetscMin(user->v, zero);

197:   for (cell = cStart; cell < cEnd; cell++) {
198:     /* Obtaining the faces of the cell */
199:     DMPlexGetConeSize(da, cell, &nF);
200:     DMPlexGetCone(da, cell, &cellcone);

202:     /* south */
203:     DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA);
204:     centroid_y[0] = fgA->centroid[1];
205:     /* North */
206:     DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA);
207:     centroid_y[1] = fgA->centroid[1];
208:     /* West */
209:     DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA);
210:     centroid_x[0] = fgA->centroid[0];
211:     /* East */
212:     DMPlexPointLocalRead(dmFace, cellcone[1], fgeom, &fgA);
213:     centroid_x[1] = fgA->centroid[0];

215:     /* Computing the cell widths in the x and y direction */
216:     delta_x = centroid_x[1] - centroid_x[0];
217:     delta_y = centroid_y[1] - centroid_y[0];

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

222:     /* cellcone[0] - south */
223:     DMPlexGetSupportSize(da, cellcone[0], &nC);
224:     DMPlexGetSupport(da, cellcone[0], &cellsupport);
225:     if (nC == 2) flux_south = (x[cellsupport[0]] * (-v_plus - user->diffusion * delta_x)) / delta_y;
226:     else flux_south = (boundary * (-v_plus - user->diffusion * delta_x)) / delta_y;

228:     /* cellcone[1] - east */
229:     DMPlexGetSupportSize(da, cellcone[1], &nC);
230:     DMPlexGetSupport(da, cellcone[1], &cellsupport);
231:     if (nC == 2) flux_east = (x[cellsupport[1]] * (u_minus - user->diffusion * delta_y)) / delta_x;
232:     else flux_east = (boundary * (u_minus - user->diffusion * delta_y)) / delta_x;

234:     /* cellcone[2] - north */
235:     DMPlexGetSupportSize(da, cellcone[2], &nC);
236:     DMPlexGetSupport(da, cellcone[2], &cellsupport);
237:     if (nC == 2) flux_north = (x[cellsupport[1]] * (v_minus - user->diffusion * delta_x)) / delta_y;
238:     else flux_north = (boundary * (v_minus - user->diffusion * delta_x)) / delta_y;

240:     /* cellcone[3] - west */
241:     DMPlexGetSupportSize(da, cellcone[3], &nC);
242:     DMPlexGetSupport(da, cellcone[3], &cellsupport);
243:     if (nC == 2) flux_west = (x[cellsupport[0]] * (-u_plus - user->diffusion * delta_y)) / delta_x;
244:     else flux_west = (boundary_left * (-u_plus - user->diffusion * delta_y)) / delta_x;

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

249:     /* Calculating the net flux for each cell
250:            and computing the RHS time derivative f[.] */
251:     f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south);
252:   }
253:   PetscFVDestroy(&fvm);
254:   VecRestoreArray(localX, &x);
255:   VecRestoreArray(F, &f);
256:   DMRestoreLocalVector(da, &localX);
257:   return 0;
258: }

260: int main(int argc, char **argv)
261: {
262:   TS                    ts; /* time integrator */
263:   SNES                  snes;
264:   Vec                   x, r; /* solution, residual vectors */
265:   DM                    da;
266:   PetscMPIInt           rank;
267:   PetscViewerAndFormat *vf;
268:   AppCtx                user; /* mesh context */
269:   PetscInt              dim, numFields = 1, numBC, i;
270:   PetscInt              numComp[1];
271:   PetscInt              numDof[12];
272:   PetscInt              bcField[1];
273:   PetscSection          section;
274:   IS                    bcPointIS[1];

276:   /* Initialize program */
278:   PetscInitialize(&argc, &argv, (char *)0, help);
279:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
280:   /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
281:   ProcessOptions(PETSC_COMM_WORLD, &user);
282:   CreateMesh(PETSC_COMM_WORLD, &user, &da);
283:   DMGetDimension(da, &dim);

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

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

293:   numDof[0 * (dim + 1)]           = 1;
294:   numDof[0 * (dim + 1) + dim - 1] = 1;
295:   numDof[0 * (dim + 1) + dim]     = 1;

297:   /* Setup boundary conditions */
298:   numBC = 1;
299:   /* Prescribe a Dirichlet condition on u on the boundary
300:        Label "marker" is made by the mesh creation routine  */
301:   bcField[0] = 0;
302:   DMGetStratumIS(da, "marker", 1, &bcPointIS[0]);

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

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

311:   /* Tell the DM to use this section (with the specified fields and dof) */
312:   DMSetLocalSection(da, section);

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

317:   /* Create a Vec with this layout and view it */
318:   DMGetGlobalVector(da, &x);
319:   VecDuplicate(x, &r);

321:   /* Create timestepping solver context */
322:   TSCreate(PETSC_COMM_WORLD, &ts);
323:   TSSetProblemType(ts, TS_NONLINEAR);
324:   TSSetRHSFunction(ts, NULL, FormFunction, &user);

326:   TSSetMaxTime(ts, 1.0);
327:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
328:   TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL);
329:   TSSetDM(ts, da);

331:   /* Customize nonlinear solver */
332:   TSSetType(ts, TSEULER);
333:   TSGetSNES(ts, &snes);
334:   PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf);
335:   SNESMonitorSet(snes, (PetscErrorCode(*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);

337:   /* Set initial conditions */
338:   FormInitialSolution(da, x);
339:   TSSetTimeStep(ts, .0001);
340:   TSSetSolution(ts, x);
341:   /* Set runtime options */
342:   TSSetFromOptions(ts);
343:   /* Solve nonlinear system */
344:   TSSolve(ts, x);

346:   /* Clean up routine */
347:   DMRestoreGlobalVector(da, &x);
348:   ISDestroy(&bcPointIS[0]);
349:   PetscSectionDestroy(&section);
350:   VecDestroy(&r);
351:   TSDestroy(&ts);
352:   DMDestroy(&da);
353:   PetscFinalize();
354:   return 0;
355: }

357: /*TEST

359:     test:
360:       suffix: 0
361:       args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -dm_plex_boundary_label boundary -ts_max_steps 5 -ts_type rk
362:       requires: !single !complex triangle ctetgen

364: TEST*/