Actual source code: ex52.c
petsc-3.14.6 2021-03-30
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)
47: {
51: options->interpolate = PETSC_TRUE;
52: options->filename[0] = '\0';
53: options->dim = 2;
54: options->u = 2.5;
55: options->v = 0.0;
56: options->cells[0] = 20;
57: options->cells[1] = 20;
58: options->diffusion = 0.0;
60: PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
61: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "advection_DMPLEX.c",options->interpolate, &options->interpolate, NULL);
62: PetscOptionsString("-filename", "The mesh file", "advection_DMPLEX.c", options->filename, options->filename, sizeof(options->filename), NULL);
63: PetscOptionsRangeInt("-dim", "Problem dimension used for the non-file mesh.", "ex7.c", options->dim, &options->dim, NULL,1,3);
64: PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL);
65: PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL);
66: PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL);
67: PetscOptionsEnd();
68: return(0);
69: }
71: /*
72: User can provide the file containing the mesh.
73: Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
74: */
75: static PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
76: {
77: size_t len;
80: PetscStrlen(user->filename, &len);
81: if (!len) {
82: DMLabel label;
83: DMPlexCreateBoxMesh(comm, user->dim, PETSC_FALSE, user->cells, NULL, NULL, NULL, user->interpolate, dm);
84: /* Mark boundary and set BC */
85: DMCreateLabel(*dm, "boundary");
86: DMGetLabel(*dm, "boundary", &label);
87: DMPlexMarkBoundaryFaces(*dm, 1, label);
88: DMPlexLabelComplete(*dm, label);
89: } else {
90: DMPlexCreateFromFile(comm, user->filename, user->interpolate, dm);
91: }
92: PetscObjectSetName((PetscObject) * dm, "Mesh");
93: DMSetFromOptions(*dm);
94: DMViewFromOptions(*dm, NULL, "-dm_view");
95: return(0);
96: }
98: /* This routine is responsible for defining the local solution vector x
99: with a given initial solution.
100: The initial solution can be modified accordingly inside the loops.
101: No need for a local vector because there is exchange of information
102: across the processors. Unlike for FormFunction which depends on the neighbours */
103: PetscErrorCode FormInitialSolution(DM da, Vec U)
104: {
106: PetscScalar *u;
107: PetscInt cell, cStart, cEnd;
108: PetscReal cellvol, centroid[3], normal[3];
111: /* Get pointers to vector data */
112: VecGetArray(U, &u);
113: /* Get local grid boundaries */
114: DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
115: /* Assigning the values at the cell centers based on x and y directions */
116: for (cell = cStart; cell < cEnd; cell++) {
117: DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal);
118: if (centroid[0] > 0.9 && centroid[0] < 0.95) {
119: if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0;
120: }
121: else u[cell] = 0;
122: }
123: VecRestoreArray(U, &u);
124: return(0);
125: }
127: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
128: {
130: PetscReal norm;
131: MPI_Comm comm;
134: if (step < 0) return(0); /* step of -1 indicates an interpolated solution */
135: VecNorm(v, NORM_2, &norm);
136: PetscObjectGetComm((PetscObject) ts, &comm);
137: PetscPrintf(comm, "timestep %D time %g norm %g\n", step, (double) ptime, (double) norm);
138: return(0);
139: }
141: /*
142: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
143: Input Parameters:
144: snes - the SNES context
145: its - iteration number
146: fnorm - 2-norm function value (may be estimated)
147: ctx - optional user-defined context for private data for the
148: monitor routine, as set by SNESMonitorSet()
149: */
150: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
151: {
155: SNESMonitorDefaultShort(snes, its, fnorm, vf);
156: return(0);
157: }
159: /*
160: FormFunction - Evaluates nonlinear function, F(x).
162: Input Parameters:
163: . ts - the TS context
164: . X - input vector
165: . ctx - optional user-defined context, as set by SNESSetFunction()
167: Output Parameter:
168: . F - function vector
169: */
170: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ctx)
171: {
172: AppCtx *user = (AppCtx *) ctx;
173: DM da = (DM) user->da;
175: PetscScalar *x, *f;
176: Vec localX;
177: PetscInt fStart, fEnd, nF;
178: PetscInt cell, cStart, cEnd, nC;
179: DM dmFace; /* DMPLEX for face geometry */
180: PetscFV fvm; /* specify type of FVM discretization */
181: Vec cellGeom, faceGeom; /* vector of structs related to cell/face geometry*/
182: const PetscScalar *fgeom; /* values stored in the vector facegeom */
183: PetscFVFaceGeom *fgA; /* struct with face geometry information */
184: const PetscInt *cellcone, *cellsupport;
185: PetscScalar flux_east, flux_west, flux_north, flux_south, flux_centre;
186: PetscScalar centroid_x[2], centroid_y[2], boundary = 0.0;
187: PetscScalar boundary_left = 0.0;
188: PetscReal u_plus, u_minus, v_plus, v_minus, zero = 0.0;
189: PetscScalar delta_x, delta_y;
191: /* Get the local vector from the DM object. */
193: DMGetLocalVector(da, &localX);
195: /* Scatter ghost points to local vector,using the 2-step process
196: DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
197: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
198: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
199: /* Get pointers to vector data. */
200: VecGetArray(localX, &x);
201: VecGetArray(F, &f);
203: /* Obtaining local cell and face ownership */
204: DMPlexGetHeightStratum(da, 0, &cStart, &cEnd);
205: DMPlexGetHeightStratum(da, 1, &fStart, &fEnd);
207: /* Creating the PetscFV object to obtain face and cell geometry.
208: Later to be used to compute face centroid to find cell widths. */
210: PetscFVCreate(PETSC_COMM_WORLD, &fvm);
211: PetscFVSetType(fvm, PETSCFVUPWIND);
212: /*....Retrieve precomputed cell geometry....*/
213: DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL);
214: VecGetDM(faceGeom, &dmFace);
215: VecGetArrayRead(faceGeom, &fgeom);
217: /* Spanning through all the cells and an inner loop through the faces. Find the
218: face neighbors and pick the upwinded cell value for flux. */
220: u_plus = PetscMax(user->u, zero);
221: u_minus = PetscMin(user->u, zero);
222: v_plus = PetscMax(user->v, zero);
223: v_minus = PetscMin(user->v, zero);
225: for (cell = cStart; cell < cEnd; cell++) {
226: /* Obtaining the faces of the cell */
227: DMPlexGetConeSize(da, cell, &nF);
228: DMPlexGetCone(da, cell, &cellcone);
230: /* south */
231: DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA);
232: centroid_y[0] = fgA->centroid[1];
233: /* North */
234: DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA);
235: centroid_y[1] = fgA->centroid[1];
236: /* West */
237: DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA);
238: centroid_x[0] = fgA->centroid[0];
239: /* East */
240: DMPlexPointLocalRead(dmFace, cellcone[1], fgeom, &fgA);
241: centroid_x[1] = fgA->centroid[0];
243: /* Computing the cell widths in the x and y direction */
244: delta_x = centroid_x[1] - centroid_x[0];
245: delta_y = centroid_y[1] - centroid_y[0];
247: /* Getting the neighbors of each face
248: Going through the faces by the order (cellcone) */
250: /* cellcone[0] - south */
251: DMPlexGetSupportSize(da, cellcone[0], &nC);
252: DMPlexGetSupport(da, cellcone[0], &cellsupport);
253: if (nC == 2) flux_south = (x[cellsupport[0]] * (-v_plus - user->diffusion * delta_x)) / delta_y;
254: else flux_south = (boundary * (-v_plus - user->diffusion * delta_x)) / delta_y;
256: /* cellcone[1] - east */
257: DMPlexGetSupportSize(da, cellcone[1], &nC);
258: DMPlexGetSupport(da, cellcone[1], &cellsupport);
259: if (nC == 2) flux_east = (x[cellsupport[1]] * (u_minus - user->diffusion * delta_y)) / delta_x;
260: else flux_east = (boundary * (u_minus - user->diffusion * delta_y)) / delta_x;
262: /* cellcone[2] - north */
263: DMPlexGetSupportSize(da, cellcone[2], &nC);
264: DMPlexGetSupport(da, cellcone[2], &cellsupport);
265: if (nC == 2) flux_north = (x[cellsupport[1]] * (v_minus - user->diffusion * delta_x)) / delta_y;
266: else flux_north = (boundary * (v_minus - user->diffusion * delta_x)) / delta_y;
268: /* cellcone[3] - west */
269: DMPlexGetSupportSize(da, cellcone[3], &nC);
270: DMPlexGetSupport(da, cellcone[3], &cellsupport);
271: if (nC == 2) flux_west = (x[cellsupport[0]] * (-u_plus - user->diffusion * delta_y)) / delta_x;
272: else flux_west = (boundary_left * (-u_plus - user->diffusion * delta_y)) / delta_x;
274: /* Contribution by the cell to the fluxes */
275: flux_centre = x[cell] * ((u_plus - u_minus + 2 * user->diffusion * delta_y) / delta_x +
276: (v_plus - v_minus + 2 * user->diffusion * delta_x) / delta_y);
278: /* Calculating the net flux for each cell
279: and computing the RHS time derivative f[.] */
280: f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south);
281: }
282: PetscFVDestroy(&fvm);
283: VecRestoreArray(localX, &x);
284: VecRestoreArray(F, &f);
285: DMRestoreLocalVector(da, &localX);
286: return(0);
287: }
289: int main(int argc, char **argv)
290: {
291: TS ts; /* time integrator */
292: SNES snes;
293: Vec x, r; /* solution, residual vectors */
294: PetscErrorCode ierr;
295: DM da;
296: PetscMPIInt rank;
297: PetscViewerAndFormat *vf;
298: AppCtx user; /* mesh context */
299: PetscInt numFields = 1, numBC, i;
300: PetscInt numComp[1];
301: PetscInt numDof[12];
302: PetscInt bcField[1];
303: PetscSection section;
304: IS bcPointIS[1];
306: /* Initialize program */
307: PetscInitialize(&argc, &argv, (char *) 0, help);if (ierr) return ierr;
308: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
309: /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
310: ProcessOptions(PETSC_COMM_WORLD, &user);
311: CreateMesh(PETSC_COMM_WORLD, &user, &da);
313: /* Specifying the fields and dof for the formula through PETSc Section
314: Create a scalar field u with 1 component on cells, faces and edges.
315: Alternatively, the field information could be added through a PETSCFV object
316: using DMAddField(...).*/
317: numComp[0] = 1;
319: for (i = 0; i < numFields * (user.dim + 1); ++i) numDof[i] = 0;
321: numDof[0 * (user.dim + 1)] = 1;
322: numDof[0 * (user.dim + 1) + user.dim - 1] = 1;
323: numDof[0 * (user.dim + 1) + user.dim] = 1;
325: /* Setup boundary conditions */
326: numBC = 1;
327: /* Prescribe a Dirichlet condition on u on the boundary
328: Label "marker" is made by the mesh creation routine */
329: bcField[0] = 0;
330: DMGetStratumIS(da, "marker", 1, &bcPointIS[0]);
332: /* Create a PetscSection with this data layout */
333: DMSetNumFields(da, numFields);
334: DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion);
336: /* Name the Field variables */
337: PetscSectionSetFieldName(section, 0, "u");
339: /* Tell the DM to use this section (with the specified fields and dof) */
340: DMSetLocalSection(da, section);
341: user.da = da;
343: /* Extract global vectors from DMDA; then duplicate for remaining
344: vectors that are the same types */
346: /* Create a Vec with this layout and view it */
347: DMGetGlobalVector(da, &x);
348: VecDuplicate(x, &r);
350: /* Create timestepping solver context */
351: TSCreate(PETSC_COMM_WORLD, &ts);
352: TSSetProblemType(ts, TS_NONLINEAR);
353: TSSetRHSFunction(ts, NULL, FormFunction, &user);
355: TSSetMaxTime(ts, 1.0);
356: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
357: TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL);
358: TSSetDM(ts, da);
360: /* Customize nonlinear solver */
361: TSSetType(ts, TSEULER);
362: TSGetSNES(ts, &snes);
363: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf);
364: SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy);
366: /* Set initial conditions */
367: FormInitialSolution(da, x);
368: TSSetTimeStep(ts, .0001);
369: TSSetSolution(ts, x);
370: /* Set runtime options */
371: TSSetFromOptions(ts);
372: /* Solve nonlinear system */
373: TSSolve(ts, x);
375: /* Clean up routine */
376: DMRestoreGlobalVector(da, &x);
377: ISDestroy(&bcPointIS[0]);
378: PetscSectionDestroy(§ion);
379: VecDestroy(&r);
380: TSDestroy(&ts);
381: DMDestroy(&da);
382: PetscFinalize();
383: return ierr;
384: }
386: /*TEST
388: test:
389: suffix: 0
390: args: -ts_max_steps 5 -ts_type rk
391: requires: !single !complex triangle ctetgen
393: TEST*/