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: {
43: PetscFunctionBeginUser;
44: options->u = 2.5;
45: options->v = 0.0;
46: options->diffusion = 0.0;
47: PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
48: PetscCall(PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL));
49: PetscCall(PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL));
50: PetscCall(PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL));
51: PetscOptionsEnd();
52: PetscFunctionReturn(PETSC_SUCCESS);
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: {
61: PetscFunctionBeginUser;
62: PetscCall(DMCreate(comm, dm));
63: PetscCall(DMSetType(*dm, DMPLEX));
64: PetscCall(DMSetFromOptions(*dm));
65: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
66: {
67: DMLabel label;
68: PetscCall(DMGetLabel(*dm, "boundary", &label));
69: PetscCall(DMPlexLabelComplete(*dm, label));
70: }
71: PetscFunctionReturn(PETSC_SUCCESS);
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];
85: PetscFunctionBeginUser;
86: /* Get pointers to vector data */
87: PetscCall(VecGetArray(U, &u));
88: /* Get local grid boundaries */
89: PetscCall(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd));
90: /* Assigning the values at the cell centers based on x and y directions */
91: PetscCall(DMGetCoordinatesLocalSetUp(da));
92: for (cell = cStart; cell < cEnd; cell++) {
93: PetscCall(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: PetscCall(VecRestoreArray(U, &u));
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
103: {
104: PetscReal norm;
105: MPI_Comm comm;
107: PetscFunctionBeginUser;
108: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */
109: PetscCall(VecNorm(v, NORM_2, &norm));
110: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
111: PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
112: PetscFunctionReturn(PETSC_SUCCESS);
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: {
126: PetscFunctionBeginUser;
127: PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf));
128: PetscFunctionReturn(PETSC_SUCCESS);
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. */
163: PetscFunctionBeginUser;
164: PetscCall(TSGetDM(ts, &da));
165: PetscCall(DMGetLocalVector(da, &localX));
167: /* Scatter ghost points to local vector,using the 2-step process
168: DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */
169: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
170: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
171: /* Get pointers to vector data. */
172: PetscCall(VecGetArray(localX, &x));
173: PetscCall(VecGetArray(F, &f));
175: /* Obtaining local cell and face ownership */
176: PetscCall(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd));
177: PetscCall(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: PetscCall(PetscFVCreate(PETSC_COMM_WORLD, &fvm));
183: PetscCall(PetscFVSetType(fvm, PETSCFVUPWIND));
184: /*....Retrieve precomputed cell geometry....*/
185: PetscCall(DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL));
186: PetscCall(VecGetDM(faceGeom, &dmFace));
187: PetscCall(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: PetscCall(DMPlexGetConeSize(da, cell, &nF));
200: PetscCall(DMPlexGetCone(da, cell, &cellcone));
202: /* south */
203: PetscCall(DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA));
204: centroid_y[0] = fgA->centroid[1];
205: /* North */
206: PetscCall(DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA));
207: centroid_y[1] = fgA->centroid[1];
208: /* West */
209: PetscCall(DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA));
210: centroid_x[0] = fgA->centroid[0];
211: /* East */
212: PetscCall(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: PetscCall(DMPlexGetSupportSize(da, cellcone[0], &nC));
224: PetscCall(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: PetscCall(DMPlexGetSupportSize(da, cellcone[1], &nC));
230: PetscCall(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: PetscCall(DMPlexGetSupportSize(da, cellcone[2], &nC));
236: PetscCall(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: PetscCall(DMPlexGetSupportSize(da, cellcone[3], &nC));
242: PetscCall(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: PetscCall(PetscFVDestroy(&fvm));
254: PetscCall(VecRestoreArray(localX, &x));
255: PetscCall(VecRestoreArray(F, &f));
256: PetscCall(DMRestoreLocalVector(da, &localX));
257: PetscFunctionReturn(PETSC_SUCCESS);
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 */
277: PetscFunctionBeginUser;
278: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
279: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
280: /* Create distributed array (DMPLEX) to manage parallel grid and vectors */
281: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
282: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &da));
283: PetscCall(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: PetscCall(DMGetStratumIS(da, "marker", 1, &bcPointIS[0]));
304: /* Create a PetscSection with this data layout */
305: PetscCall(DMSetNumFields(da, numFields));
306: PetscCall(DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion));
308: /* Name the Field variables */
309: PetscCall(PetscSectionSetFieldName(section, 0, "u"));
311: /* Tell the DM to use this section (with the specified fields and dof) */
312: PetscCall(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: PetscCall(DMGetGlobalVector(da, &x));
319: PetscCall(VecDuplicate(x, &r));
321: /* Create timestepping solver context */
322: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
323: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
324: PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user));
326: PetscCall(TSSetMaxTime(ts, 1.0));
327: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
328: PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL));
329: PetscCall(TSSetDM(ts, da));
331: /* Customize nonlinear solver */
332: PetscCall(TSSetType(ts, TSEULER));
333: PetscCall(TSGetSNES(ts, &snes));
334: PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
335: PetscCall(SNESMonitorSet(snes, (PetscErrorCode(*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy));
337: /* Set initial conditions */
338: PetscCall(FormInitialSolution(da, x));
339: PetscCall(TSSetTimeStep(ts, .0001));
340: PetscCall(TSSetSolution(ts, x));
341: /* Set runtime options */
342: PetscCall(TSSetFromOptions(ts));
343: /* Solve nonlinear system */
344: PetscCall(TSSolve(ts, x));
346: /* Clean up routine */
347: PetscCall(DMRestoreGlobalVector(da, &x));
348: PetscCall(ISDestroy(&bcPointIS[0]));
349: PetscCall(PetscSectionDestroy(§ion));
350: PetscCall(VecDestroy(&r));
351: PetscCall(TSDestroy(&ts));
352: PetscCall(DMDestroy(&da));
353: PetscCall(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*/