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 PetscFunctionBeginUser; 280 PetscCall(PetscInitialize(&argc, &argv, (char *) 0, help)); 281 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 282 /* Create distributed array (DMPLEX) to manage parallel grid and vectors */ 283 PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user)); 284 PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &da)); 285 PetscCall(DMGetDimension(da, &dim)); 286 287 /* Specifying the fields and dof for the formula through PETSc Section 288 Create a scalar field u with 1 component on cells, faces and edges. 289 Alternatively, the field information could be added through a PETSCFV object 290 using DMAddField(...).*/ 291 numComp[0] = 1; 292 293 for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0; 294 295 numDof[0 * (dim + 1)] = 1; 296 numDof[0 * (dim + 1) + dim - 1] = 1; 297 numDof[0 * (dim + 1) + dim] = 1; 298 299 /* Setup boundary conditions */ 300 numBC = 1; 301 /* Prescribe a Dirichlet condition on u on the boundary 302 Label "marker" is made by the mesh creation routine */ 303 bcField[0] = 0; 304 PetscCall(DMGetStratumIS(da, "marker", 1, &bcPointIS[0])); 305 306 /* Create a PetscSection with this data layout */ 307 PetscCall(DMSetNumFields(da, numFields)); 308 PetscCall(DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion)); 309 310 /* Name the Field variables */ 311 PetscCall(PetscSectionSetFieldName(section, 0, "u")); 312 313 /* Tell the DM to use this section (with the specified fields and dof) */ 314 PetscCall(DMSetLocalSection(da, section)); 315 316 /* Extract global vectors from DMDA; then duplicate for remaining 317 vectors that are the same types */ 318 319 /* Create a Vec with this layout and view it */ 320 PetscCall(DMGetGlobalVector(da, &x)); 321 PetscCall(VecDuplicate(x, &r)); 322 323 /* Create timestepping solver context */ 324 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 325 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 326 PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &user)); 327 328 PetscCall(TSSetMaxTime(ts, 1.0)); 329 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 330 PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL)); 331 PetscCall(TSSetDM(ts, da)); 332 333 /* Customize nonlinear solver */ 334 PetscCall(TSSetType(ts, TSEULER)); 335 PetscCall(TSGetSNES(ts, &snes)); 336 PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf)); 337 PetscCall(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy)); 338 339 /* Set initial conditions */ 340 PetscCall(FormInitialSolution(da, x)); 341 PetscCall(TSSetTimeStep(ts, .0001)); 342 PetscCall(TSSetSolution(ts, x)); 343 /* Set runtime options */ 344 PetscCall(TSSetFromOptions(ts)); 345 /* Solve nonlinear system */ 346 PetscCall(TSSolve(ts, x)); 347 348 /* Clean up routine */ 349 PetscCall(DMRestoreGlobalVector(da, &x)); 350 PetscCall(ISDestroy(&bcPointIS[0])); 351 PetscCall(PetscSectionDestroy(§ion)); 352 PetscCall(VecDestroy(&r)); 353 PetscCall(TSDestroy(&ts)); 354 PetscCall(DMDestroy(&da)); 355 PetscCall(PetscFinalize()); 356 return 0; 357 } 358 359 /*TEST 360 361 test: 362 suffix: 0 363 args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -dm_plex_boundary_label boundary -ts_max_steps 5 -ts_type rk 364 requires: !single !complex triangle ctetgen 365 366 TEST*/ 367