1 static char help[] = "Simple Advection-diffusion equation solved using FVM in DMPLEX\n";
2
3 /*
4 Solves the simple advection equation given by
5
6 q_t + u (q_x) + v (q_y) - D (q_xx + q_yy) = 0 using FVM and First Order Upwind discretization.
7
8 with a user defined initial condition.
9
10 with dirichlet/neumann conditions on the four boundaries of the domain.
11
12 User can define the mesh parameters either in the command line or inside
13 the ProcessOptions() routine.
14
15 Contributed by: Mukkund Sunjii, Domenico Lahaye
16 */
17
18 #include <petscdmplex.h>
19 #include <petscts.h>
20 #include <petscblaslapack.h>
21
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 *);
32
33 /* Defining the usr defined context */
34 typedef struct {
35 PetscScalar diffusion;
36 PetscReal u, v;
37 PetscScalar delta_x, delta_y;
38 } AppCtx;
39
40 /* Options for the scenario */
ProcessOptions(MPI_Comm comm,AppCtx * options)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 }
54
55 /*
56 User can provide the file containing the mesh.
57 Or can generate the mesh using DMPlexCreateBoxMesh with the specified options.
58 */
CreateMesh(MPI_Comm comm,AppCtx * user,DM * dm)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 }
73
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 */
FormInitialSolution(DM da,Vec U)79 PetscErrorCode FormInitialSolution(DM da, Vec U)
80 {
81 PetscScalar *u;
82 PetscInt cell, cStart, cEnd;
83 PetscReal cellvol, centroid[3], normal[3];
84
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 }
101
MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscCtx ctx)102 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscCtx ctx)
103 {
104 PetscReal norm;
105 MPI_Comm comm;
106
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 }
114
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 */
MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat * vf)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 }
130
131 /*
132 FormFunction - Evaluates nonlinear function, F(x).
133
134 Input Parameters:
135 . ts - the TS context
136 . X - input vector
137 . ctx - optional user-defined context, as set by SNESSetFunction()
138
139 Output Parameter:
140 . F - function vector
141 */
FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,PetscCtx ctx)142 PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, PetscCtx 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;
161
162 /* Get the local vector from the DM object. */
163 PetscFunctionBeginUser;
164 PetscCall(TSGetDM(ts, &da));
165 PetscCall(DMGetLocalVector(da, &localX));
166
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));
174
175 /* Obtaining local cell and face ownership */
176 PetscCall(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd));
177 PetscCall(DMPlexGetHeightStratum(da, 1, &fStart, &fEnd));
178
179 /* Creating the PetscFV object to obtain face and cell geometry.
180 Later to be used to compute face centroid to find cell widths. */
181
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));
188
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. */
191
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);
196
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));
201
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];
214
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];
218
219 /* Getting the neighbors of each face
220 Going through the faces by the order (cellcone) */
221
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;
227
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;
233
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;
239
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;
245
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);
248
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 }
259
main(int argc,char ** argv)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];
275
276 /* Initialize program */
277 PetscFunctionBeginUser;
278 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
284
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;
290
291 for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0;
292
293 numDof[0 * (dim + 1)] = 1;
294 numDof[0 * (dim + 1) + dim - 1] = 1;
295 numDof[0 * (dim + 1) + dim] = 1;
296
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]));
303
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));
307
308 /* Name the Field variables */
309 PetscCall(PetscSectionSetFieldName(section, 0, "u"));
310
311 /* Tell the DM to use this section (with the specified fields and dof) */
312 PetscCall(DMSetLocalSection(da, section));
313
314 /* Extract global vectors from DMDA; then duplicate for remaining
315 vectors that are the same types */
316
317 /* Create a Vec with this layout and view it */
318 PetscCall(DMGetGlobalVector(da, &x));
319 PetscCall(VecDuplicate(x, &r));
320
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));
325
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));
330
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, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
336
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));
345
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 }
356
357 /*TEST
358
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
363
364 TEST*/
365