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 PetscErrorCode ierr; 270 DM da; 271 PetscMPIInt rank; 272 PetscViewerAndFormat *vf; 273 AppCtx user; /* mesh context */ 274 PetscInt dim, numFields = 1, numBC, i; 275 PetscInt numComp[1]; 276 PetscInt numDof[12]; 277 PetscInt bcField[1]; 278 PetscSection section; 279 IS bcPointIS[1]; 280 281 /* Initialize program */ 282 ierr = PetscInitialize(&argc, &argv, (char *) 0, help);if (ierr) return ierr; 283 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 284 /* Create distributed array (DMPLEX) to manage parallel grid and vectors */ 285 CHKERRQ(ProcessOptions(PETSC_COMM_WORLD, &user)); 286 CHKERRQ(CreateMesh(PETSC_COMM_WORLD, &user, &da)); 287 CHKERRQ(DMGetDimension(da, &dim)); 288 289 /* Specifying the fields and dof for the formula through PETSc Section 290 Create a scalar field u with 1 component on cells, faces and edges. 291 Alternatively, the field information could be added through a PETSCFV object 292 using DMAddField(...).*/ 293 numComp[0] = 1; 294 295 for (i = 0; i < numFields * (dim + 1); ++i) numDof[i] = 0; 296 297 numDof[0 * (dim + 1)] = 1; 298 numDof[0 * (dim + 1) + dim - 1] = 1; 299 numDof[0 * (dim + 1) + dim] = 1; 300 301 /* Setup boundary conditions */ 302 numBC = 1; 303 /* Prescribe a Dirichlet condition on u on the boundary 304 Label "marker" is made by the mesh creation routine */ 305 bcField[0] = 0; 306 CHKERRQ(DMGetStratumIS(da, "marker", 1, &bcPointIS[0])); 307 308 /* Create a PetscSection with this data layout */ 309 CHKERRQ(DMSetNumFields(da, numFields)); 310 CHKERRQ(DMPlexCreateSection(da, NULL, numComp, numDof, numBC, bcField, NULL, bcPointIS, NULL, §ion)); 311 312 /* Name the Field variables */ 313 CHKERRQ(PetscSectionSetFieldName(section, 0, "u")); 314 315 /* Tell the DM to use this section (with the specified fields and dof) */ 316 CHKERRQ(DMSetLocalSection(da, section)); 317 318 /* Extract global vectors from DMDA; then duplicate for remaining 319 vectors that are the same types */ 320 321 /* Create a Vec with this layout and view it */ 322 CHKERRQ(DMGetGlobalVector(da, &x)); 323 CHKERRQ(VecDuplicate(x, &r)); 324 325 /* Create timestepping solver context */ 326 CHKERRQ(TSCreate(PETSC_COMM_WORLD, &ts)); 327 CHKERRQ(TSSetProblemType(ts, TS_NONLINEAR)); 328 CHKERRQ(TSSetRHSFunction(ts, NULL, FormFunction, &user)); 329 330 CHKERRQ(TSSetMaxTime(ts, 1.0)); 331 CHKERRQ(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 332 CHKERRQ(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL)); 333 CHKERRQ(TSSetDM(ts, da)); 334 335 /* Customize nonlinear solver */ 336 CHKERRQ(TSSetType(ts, TSEULER)); 337 CHKERRQ(TSGetSNES(ts, &snes)); 338 CHKERRQ(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf)); 339 CHKERRQ(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *)) MySNESMonitor, vf,(PetscErrorCode (*)(void **)) PetscViewerAndFormatDestroy)); 340 341 /* Set initial conditions */ 342 CHKERRQ(FormInitialSolution(da, x)); 343 CHKERRQ(TSSetTimeStep(ts, .0001)); 344 CHKERRQ(TSSetSolution(ts, x)); 345 /* Set runtime options */ 346 CHKERRQ(TSSetFromOptions(ts)); 347 /* Solve nonlinear system */ 348 CHKERRQ(TSSolve(ts, x)); 349 350 /* Clean up routine */ 351 CHKERRQ(DMRestoreGlobalVector(da, &x)); 352 CHKERRQ(ISDestroy(&bcPointIS[0])); 353 CHKERRQ(PetscSectionDestroy(§ion)); 354 CHKERRQ(VecDestroy(&r)); 355 CHKERRQ(TSDestroy(&ts)); 356 CHKERRQ(DMDestroy(&da)); 357 ierr = PetscFinalize(); 358 return ierr; 359 } 360 361 /*TEST 362 363 test: 364 suffix: 0 365 args: -dm_plex_simplex 0 -dm_plex_box_faces 20,20 -dm_plex_boundary_label boundary -ts_max_steps 5 -ts_type rk 366 requires: !single !complex triangle ctetgen 367 368 TEST*/ 369