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