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