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