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