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