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 for (cell = cStart; cell < cEnd; cell++) { 92 PetscCall(DMPlexComputeCellGeometryFVM(da, cell, &cellvol, centroid, normal)); 93 if (centroid[0] > 0.9 && centroid[0] < 0.95) { 94 if (centroid[1] > 0.9 && centroid[1] < 0.95) u[cell] = 2.0; 95 } 96 else u[cell] = 0; 97 } 98 PetscCall(VecRestoreArray(U, &u)); 99 PetscFunctionReturn(0); 100 } 101 102 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx) 103 { 104 PetscReal norm; 105 MPI_Comm comm; 106 107 PetscFunctionBeginUser; 108 if (step < 0) PetscFunctionReturn(0); /* 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 %D time %g norm %g\n", step, (double) ptime, (double) norm)); 112 PetscFunctionReturn(0); 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 */ 124 PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf) 125 { 126 PetscFunctionBeginUser; 127 PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf)); 128 PetscFunctionReturn(0); 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 */ 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; 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 + 248 (v_plus - v_minus + 2 * user->diffusion * delta_x) / delta_y); 249 250 /* Calculating the net flux for each cell 251 and computing the RHS time derivative f[.] */ 252 f[cell] = -(flux_centre + flux_east + flux_west + flux_north + flux_south); 253 } 254 PetscCall(PetscFVDestroy(&fvm)); 255 PetscCall(VecRestoreArray(localX, &x)); 256 PetscCall(VecRestoreArray(F, &f)); 257 PetscCall(DMRestoreLocalVector(da, &localX)); 258 PetscFunctionReturn(0); 259 } 260 261 int main(int argc, char **argv) 262 { 263 TS ts; /* time integrator */ 264 SNES snes; 265 Vec x, r; /* solution, residual vectors */ 266 DM da; 267 PetscMPIInt rank; 268 PetscViewerAndFormat *vf; 269 AppCtx user; /* mesh context */ 270 PetscInt dim, numFields = 1, numBC, i; 271 PetscInt numComp[1]; 272 PetscInt numDof[12]; 273 PetscInt bcField[1]; 274 PetscSection section; 275 IS bcPointIS[1]; 276 277 /* Initialize program */ 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)); 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,(PetscErrorCode (*)(void **)) 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