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");CHKERRQ(ierr); 51 CHKERRQ(PetscOptionsReal("-u", "The x component of the convective coefficient", "advection_DMPLEX.c", options->u, &options->u, NULL)); 52 CHKERRQ(PetscOptionsReal("-v", "The y component of the convective coefficient", "advection_DMPLEX.c", options->v, &options->v, NULL)); 53 CHKERRQ(PetscOptionsScalar("-diffus", "The diffusive coefficient", "advection_DMPLEX.c", options->diffusion, &options->diffusion, NULL)); 54 ierr = PetscOptionsEnd();CHKERRQ(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 CHKERRQ(DMCreate(comm, dm)); 66 CHKERRQ(DMSetType(*dm, DMPLEX)); 67 CHKERRQ(DMSetFromOptions(*dm)); 68 CHKERRQ(DMViewFromOptions(*dm, NULL, "-dm_view")); 69 { 70 DMLabel label; 71 CHKERRQ(DMGetLabel(*dm, "boundary", &label)); 72 CHKERRQ(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 CHKERRQ(VecGetArray(U, &u)); 91 /* Get local grid boundaries */ 92 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecNorm(v, NORM_2, &norm)); 113 CHKERRQ(PetscObjectGetComm((PetscObject) ts, &comm)); 114 CHKERRQ(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 CHKERRQ(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 CHKERRQ(TSGetDM(ts, &da)); 168 CHKERRQ(DMGetLocalVector(da, &localX)); 169 170 /* Scatter ghost points to local vector,using the 2-step process 171 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). */ 172 CHKERRQ(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 173 CHKERRQ(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 174 /* Get pointers to vector data. */ 175 CHKERRQ(VecGetArray(localX, &x)); 176 CHKERRQ(VecGetArray(F, &f)); 177 178 /* Obtaining local cell and face ownership */ 179 CHKERRQ(DMPlexGetHeightStratum(da, 0, &cStart, &cEnd)); 180 CHKERRQ(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 CHKERRQ(PetscFVCreate(PETSC_COMM_WORLD, &fvm)); 186 CHKERRQ(PetscFVSetType(fvm, PETSCFVUPWIND)); 187 /*....Retrieve precomputed cell geometry....*/ 188 CHKERRQ(DMPlexGetDataFVM(da, fvm, &cellGeom, &faceGeom, NULL)); 189 CHKERRQ(VecGetDM(faceGeom, &dmFace)); 190 CHKERRQ(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 CHKERRQ(DMPlexGetConeSize(da, cell, &nF)); 203 CHKERRQ(DMPlexGetCone(da, cell, &cellcone)); 204 205 /* south */ 206 CHKERRQ(DMPlexPointLocalRead(dmFace, cellcone[0], fgeom, &fgA)); 207 centroid_y[0] = fgA->centroid[1]; 208 /* North */ 209 CHKERRQ(DMPlexPointLocalRead(dmFace, cellcone[2], fgeom, &fgA)); 210 centroid_y[1] = fgA->centroid[1]; 211 /* West */ 212 CHKERRQ(DMPlexPointLocalRead(dmFace, cellcone[3], fgeom, &fgA)); 213 centroid_x[0] = fgA->centroid[0]; 214 /* East */ 215 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(da, cellcone[0], &nC)); 227 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(da, cellcone[1], &nC)); 233 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(da, cellcone[2], &nC)); 239 CHKERRQ(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 CHKERRQ(DMPlexGetSupportSize(da, cellcone[3], &nC)); 245 CHKERRQ(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 CHKERRQ(PetscFVDestroy(&fvm)); 258 CHKERRQ(VecRestoreArray(localX, &x)); 259 CHKERRQ(VecRestoreArray(F, &f)); 260 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc, &argv, (char *) 0, help)); 282 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 283 /* Create distributed array (DMPLEX) to manage parallel grid and vectors */ 284 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 285 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &da)); 286 CHKERRQ(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 CHKERRQ(DMGetStratumIS(da, "marker", 1, &bcPointIS[0])); 306 307 /* Create a PetscSection with this data layout */ 308 CHKERRQ(DMSetNumFields(da, numFields)); 309 CHKERRQ(DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion)); 310 311 /* Name the Field variables */ 312 CHKERRQ(PetscSectionSetFieldName(section, 0, "u")); 313 314 /* Tell the DM to use this section (with the specified fields and dof) */ 315 CHKERRQ(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 CHKERRQ(DMGetGlobalVector(da, &x)); 322 CHKERRQ(VecDuplicate(x, &r)); 323 324 /* Create timestepping solver context */ 325 CHKERRQ(TSCreate(PETSC_COMM_WORLD, &ts)); 326 CHKERRQ(TSSetProblemType(ts, TS_NONLINEAR)); 327 CHKERRQ(TSSetRHSFunction(ts, NULL, FormFunction, &user)); 328 329 CHKERRQ(TSSetMaxTime(ts, 1.0)); 330 CHKERRQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 331 CHKERRQ(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL)); 332 CHKERRQ(TSSetDM(ts, da)); 333 334 /* Customize nonlinear solver */ 335 CHKERRQ(TSSetType(ts, TSEULER)); 336 CHKERRQ(TSGetSNES(ts, &snes)); 337 CHKERRQ(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf)); 338 CHKERRQ(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy)); 339 340 /* Set initial conditions */ 341 CHKERRQ(FormInitialSolution(da, x)); 342 CHKERRQ(TSSetTimeStep(ts, .0001)); 343 CHKERRQ(TSSetSolution(ts, x)); 344 /* Set runtime options */ 345 CHKERRQ(TSSetFromOptions(ts)); 346 /* Solve nonlinear system */ 347 CHKERRQ(TSSolve(ts, x)); 348 349 /* Clean up routine */ 350 CHKERRQ(DMRestoreGlobalVector(da, &x)); 351 CHKERRQ(ISDestroy(&bcPointIS[0])); 352 CHKERRQ(PetscSectionDestroy(§ion)); 353 CHKERRQ(VecDestroy(&r)); 354 CHKERRQ(TSDestroy(&ts)); 355 CHKERRQ(DMDestroy(&da)); 356 CHKERRQ(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