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