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