1 2 static char help[] = "Transient nonlinear driven cavity in 2d.\n\ 3 \n\ 4 The 2D driven cavity problem is solved in a velocity-vorticity formulation.\n\ 5 The flow can be driven with the lid or with bouyancy or both:\n\ 6 -lidvelocity <lid>, where <lid> = dimensionless velocity of lid\n\ 7 -grashof <gr>, where <gr> = dimensionless temperature gradent\n\ 8 -prandtl <pr>, where <pr> = dimensionless thermal/momentum diffusity ratio\n\ 9 -contours : draw contour plots of solution\n\n"; 10 /* 11 See src/snes/tutorials/ex19.c for the steady-state version. 12 13 There used to be a SNES example (src/snes/tutorials/ex27.c) that 14 implemented this algorithm without using TS and was used for the numerical 15 results in the paper 16 17 Todd S. Coffey and C. T. Kelley and David E. Keyes, Pseudotransient 18 Continuation and Differential-Algebraic Equations, 2003. 19 20 That example was removed because it used obsolete interfaces, but the 21 algorithms from the paper can be reproduced using this example. 22 23 Note: The paper describes the algorithm as being linearly implicit but the 24 numerical results were created using nonlinearly implicit Euler. The 25 algorithm as described (linearly implicit) is more efficient and is the 26 default when using TSPSEUDO. If you want to reproduce the numerical 27 results from the paper, you'll have to change the SNES to converge the 28 nonlinear solve (e.g., -snes_type newtonls). The DAE versus ODE variants 29 are controlled using the -parabolic option. 30 31 Comment preserved from snes/tutorials/ex27.c, since removed: 32 33 [H]owever Figure 3.1 was generated with a slightly different algorithm 34 (see targets runex27 and runex27_p) than described in the paper. In 35 particular, the described algorithm is linearly implicit, advancing to 36 the next step after one Newton step, so that the steady state residual 37 is always used, but the figure was generated by converging each step to 38 a relative tolerance of 1.e-3. On the example problem, setting 39 -snes_type ksponly has only minor impact on number of steps, but 40 significantly reduces the required number of linear solves. 41 42 See also https://lists.mcs.anl.gov/pipermail/petsc-dev/2010-March/002362.html 43 */ 44 45 /* ------------------------------------------------------------------------ 46 47 We thank David E. Keyes for contributing the driven cavity discretization 48 within this example code. 49 50 This problem is modeled by the partial differential equation system 51 52 - Lap(U) - Grad_y(Omega) = 0 53 - Lap(V) + Grad_x(Omega) = 0 54 Omega_t - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0 55 T_t - Lap(T) + PR*Div([U*T,V*T]) = 0 56 57 in the unit square, which is uniformly discretized in each of x and 58 y in this simple encoding. 59 60 No-slip, rigid-wall Dirichlet conditions are used for [U,V]. 61 Dirichlet conditions are used for Omega, based on the definition of 62 vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each 63 constant coordinate boundary, the tangential derivative is zero. 64 Dirichlet conditions are used for T on the left and right walls, 65 and insulation homogeneous Neumann conditions are used for T on 66 the top and bottom walls. 67 68 A finite difference approximation with the usual 5-point stencil 69 is used to discretize the boundary value problem to obtain a 70 nonlinear system of equations. Upwinding is used for the divergence 71 (convective) terms and central for the gradient (source) terms. 72 73 The Jacobian can be either 74 * formed via finite differencing using coloring (the default), or 75 * applied matrix-free via the option -snes_mf 76 (for larger grid problems this variant may not converge 77 without a preconditioner due to ill-conditioning). 78 79 ------------------------------------------------------------------------- */ 80 81 /* 82 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 83 Include "petscts.h" so that we can use TS solvers. Note that this 84 file automatically includes: 85 petscsys.h - base PETSc routines petscvec.h - vectors 86 petscmat.h - matrices 87 petscis.h - index sets petscksp.h - Krylov subspace methods 88 petscviewer.h - viewers petscpc.h - preconditioners 89 petscksp.h - linear solvers petscsnes.h - nonlinear solvers 90 */ 91 #include <petscts.h> 92 #include <petscdm.h> 93 #include <petscdmda.h> 94 95 /* 96 User-defined routines and data structures 97 */ 98 typedef struct { 99 PetscScalar u, v, omega, temp; 100 } Field; 101 102 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, void *); 103 104 typedef struct { 105 PetscReal lidvelocity, prandtl, grashof; /* physical parameters */ 106 PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */ 107 PetscReal cfl_initial; /* CFL for first time step */ 108 } AppCtx; 109 110 PetscErrorCode FormInitialSolution(TS, Vec, AppCtx *); 111 112 int main(int argc, char **argv) { 113 AppCtx user; /* user-defined work context */ 114 PetscInt mx, my, steps; 115 TS ts; 116 DM da; 117 Vec X; 118 PetscReal ftime; 119 TSConvergedReason reason; 120 121 PetscFunctionBeginUser; 122 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 123 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 124 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 4, 1, 0, 0, &da)); 125 PetscCall(DMSetFromOptions(da)); 126 PetscCall(DMSetUp(da)); 127 PetscCall(TSSetDM(ts, (DM)da)); 128 129 PetscCall(DMDAGetInfo(da, 0, &mx, &my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 130 /* 131 Problem parameters (velocity of lid, prandtl, and grashof numbers) 132 */ 133 user.lidvelocity = 1.0 / (mx * my); 134 user.prandtl = 1.0; 135 user.grashof = 1.0; 136 user.parabolic = PETSC_FALSE; 137 user.cfl_initial = 50.; 138 139 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Driven cavity/natural convection options", ""); 140 PetscCall(PetscOptionsReal("-lidvelocity", "Lid velocity, related to Reynolds number", "", user.lidvelocity, &user.lidvelocity, NULL)); 141 PetscCall(PetscOptionsReal("-prandtl", "Ratio of viscous to thermal diffusivity", "", user.prandtl, &user.prandtl, NULL)); 142 PetscCall(PetscOptionsReal("-grashof", "Ratio of bouyant to viscous forces", "", user.grashof, &user.grashof, NULL)); 143 PetscCall(PetscOptionsBool("-parabolic", "Relax incompressibility to make the system parabolic instead of differential-algebraic", "", user.parabolic, &user.parabolic, NULL)); 144 PetscCall(PetscOptionsReal("-cfl_initial", "Advective CFL for the first time step", "", user.cfl_initial, &user.cfl_initial, NULL)); 145 PetscOptionsEnd(); 146 147 PetscCall(DMDASetFieldName(da, 0, "x-velocity")); 148 PetscCall(DMDASetFieldName(da, 1, "y-velocity")); 149 PetscCall(DMDASetFieldName(da, 2, "Omega")); 150 PetscCall(DMDASetFieldName(da, 3, "temperature")); 151 152 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 153 Create user context, set problem data, create vector data structures. 154 Also, compute the initial guess. 155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 156 157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 Create time integration context 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 160 PetscCall(DMSetApplicationContext(da, &user)); 161 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocal)FormIFunctionLocal, &user)); 162 PetscCall(TSSetMaxSteps(ts, 10000)); 163 PetscCall(TSSetMaxTime(ts, 1e12)); 164 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 165 PetscCall(TSSetTimeStep(ts, user.cfl_initial / (user.lidvelocity * mx))); 166 PetscCall(TSSetFromOptions(ts)); 167 168 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "x%" PetscInt_FMT " grid, lid velocity = %g, prandtl # = %g, grashof # = %g\n", mx, my, (double)user.lidvelocity, (double)user.prandtl, (double)user.grashof)); 169 170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 171 Solve the nonlinear system 172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 173 174 PetscCall(DMCreateGlobalVector(da, &X)); 175 PetscCall(FormInitialSolution(ts, X, &user)); 176 177 PetscCall(TSSolve(ts, X)); 178 PetscCall(TSGetSolveTime(ts, &ftime)); 179 PetscCall(TSGetStepNumber(ts, &steps)); 180 PetscCall(TSGetConvergedReason(ts, &reason)); 181 182 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps)); 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Free work space. All PETSc objects should be destroyed when they 186 are no longer needed. 187 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 188 PetscCall(VecDestroy(&X)); 189 PetscCall(DMDestroy(&da)); 190 PetscCall(TSDestroy(&ts)); 191 192 PetscCall(PetscFinalize()); 193 return 0; 194 } 195 196 /* ------------------------------------------------------------------- */ 197 198 /* 199 FormInitialSolution - Forms initial approximation. 200 201 Input Parameters: 202 user - user-defined application context 203 X - vector 204 205 Output Parameter: 206 X - vector 207 */ 208 PetscErrorCode FormInitialSolution(TS ts, Vec X, AppCtx *user) { 209 DM da; 210 PetscInt i, j, mx, xs, ys, xm, ym; 211 PetscReal grashof, dx; 212 Field **x; 213 214 grashof = user->grashof; 215 PetscCall(TSGetDM(ts, &da)); 216 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 217 dx = 1.0 / (mx - 1); 218 219 /* 220 Get local grid boundaries (for 2-dimensional DMDA): 221 xs, ys - starting grid indices (no ghost points) 222 xm, ym - widths of local grid (no ghost points) 223 */ 224 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 225 226 /* 227 Get a pointer to vector data. 228 - For default PETSc vectors, VecGetArray() returns a pointer to 229 the data array. Otherwise, the routine is implementation dependent. 230 - You MUST call VecRestoreArray() when you no longer need access to 231 the array. 232 */ 233 PetscCall(DMDAVecGetArray(da, X, &x)); 234 235 /* 236 Compute initial guess over the locally owned part of the grid 237 Initial condition is motionless fluid and equilibrium temperature 238 */ 239 for (j = ys; j < ys + ym; j++) { 240 for (i = xs; i < xs + xm; i++) { 241 x[j][i].u = 0.0; 242 x[j][i].v = 0.0; 243 x[j][i].omega = 0.0; 244 x[j][i].temp = (grashof > 0) * i * dx; 245 } 246 } 247 248 /* 249 Restore vector 250 */ 251 PetscCall(DMDAVecRestoreArray(da, X, &x)); 252 return 0; 253 } 254 255 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal ptime, Field **x, Field **xdot, Field **f, void *ptr) { 256 AppCtx *user = (AppCtx *)ptr; 257 PetscInt xints, xinte, yints, yinte, i, j; 258 PetscReal hx, hy, dhx, dhy, hxdhy, hydhx; 259 PetscReal grashof, prandtl, lid; 260 PetscScalar u, udot, uxx, uyy, vx, vy, avx, avy, vxp, vxm, vyp, vym; 261 262 PetscFunctionBeginUser; 263 grashof = user->grashof; 264 prandtl = user->prandtl; 265 lid = user->lidvelocity; 266 267 /* 268 Define mesh intervals ratios for uniform grid. 269 270 Note: FD formulae below are normalized by multiplying through by 271 local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions. 272 273 */ 274 dhx = (PetscReal)(info->mx - 1); 275 dhy = (PetscReal)(info->my - 1); 276 hx = 1.0 / dhx; 277 hy = 1.0 / dhy; 278 hxdhy = hx * dhy; 279 hydhx = hy * dhx; 280 281 xints = info->xs; 282 xinte = info->xs + info->xm; 283 yints = info->ys; 284 yinte = info->ys + info->ym; 285 286 /* Test whether we are on the bottom edge of the global array */ 287 if (yints == 0) { 288 j = 0; 289 yints = yints + 1; 290 /* bottom edge */ 291 for (i = info->xs; i < info->xs + info->xm; i++) { 292 f[j][i].u = x[j][i].u; 293 f[j][i].v = x[j][i].v; 294 f[j][i].omega = x[j][i].omega + (x[j + 1][i].u - x[j][i].u) * dhy; 295 f[j][i].temp = x[j][i].temp - x[j + 1][i].temp; 296 } 297 } 298 299 /* Test whether we are on the top edge of the global array */ 300 if (yinte == info->my) { 301 j = info->my - 1; 302 yinte = yinte - 1; 303 /* top edge */ 304 for (i = info->xs; i < info->xs + info->xm; i++) { 305 f[j][i].u = x[j][i].u - lid; 306 f[j][i].v = x[j][i].v; 307 f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j - 1][i].u) * dhy; 308 f[j][i].temp = x[j][i].temp - x[j - 1][i].temp; 309 } 310 } 311 312 /* Test whether we are on the left edge of the global array */ 313 if (xints == 0) { 314 i = 0; 315 xints = xints + 1; 316 /* left edge */ 317 for (j = info->ys; j < info->ys + info->ym; j++) { 318 f[j][i].u = x[j][i].u; 319 f[j][i].v = x[j][i].v; 320 f[j][i].omega = x[j][i].omega - (x[j][i + 1].v - x[j][i].v) * dhx; 321 f[j][i].temp = x[j][i].temp; 322 } 323 } 324 325 /* Test whether we are on the right edge of the global array */ 326 if (xinte == info->mx) { 327 i = info->mx - 1; 328 xinte = xinte - 1; 329 /* right edge */ 330 for (j = info->ys; j < info->ys + info->ym; j++) { 331 f[j][i].u = x[j][i].u; 332 f[j][i].v = x[j][i].v; 333 f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i - 1].v) * dhx; 334 f[j][i].temp = x[j][i].temp - (PetscReal)(grashof > 0); 335 } 336 } 337 338 /* Compute over the interior points */ 339 for (j = yints; j < yinte; j++) { 340 for (i = xints; i < xinte; i++) { 341 /* 342 convective coefficients for upwinding 343 */ 344 vx = x[j][i].u; 345 avx = PetscAbsScalar(vx); 346 vxp = .5 * (vx + avx); 347 vxm = .5 * (vx - avx); 348 vy = x[j][i].v; 349 avy = PetscAbsScalar(vy); 350 vyp = .5 * (vy + avy); 351 vym = .5 * (vy - avy); 352 353 /* U velocity */ 354 u = x[j][i].u; 355 udot = user->parabolic ? xdot[j][i].u : 0.; 356 uxx = (2.0 * u - x[j][i - 1].u - x[j][i + 1].u) * hydhx; 357 uyy = (2.0 * u - x[j - 1][i].u - x[j + 1][i].u) * hxdhy; 358 f[j][i].u = udot + uxx + uyy - .5 * (x[j + 1][i].omega - x[j - 1][i].omega) * hx; 359 360 /* V velocity */ 361 u = x[j][i].v; 362 udot = user->parabolic ? xdot[j][i].v : 0.; 363 uxx = (2.0 * u - x[j][i - 1].v - x[j][i + 1].v) * hydhx; 364 uyy = (2.0 * u - x[j - 1][i].v - x[j + 1][i].v) * hxdhy; 365 f[j][i].v = udot + uxx + uyy + .5 * (x[j][i + 1].omega - x[j][i - 1].omega) * hy; 366 367 /* Omega */ 368 u = x[j][i].omega; 369 uxx = (2.0 * u - x[j][i - 1].omega - x[j][i + 1].omega) * hydhx; 370 uyy = (2.0 * u - x[j - 1][i].omega - x[j + 1][i].omega) * hxdhy; 371 f[j][i].omega = (xdot[j][i].omega + uxx + uyy + (vxp * (u - x[j][i - 1].omega) + vxm * (x[j][i + 1].omega - u)) * hy + (vyp * (u - x[j - 1][i].omega) + vym * (x[j + 1][i].omega - u)) * hx - .5 * grashof * (x[j][i + 1].temp - x[j][i - 1].temp) * hy); 372 373 /* Temperature */ 374 u = x[j][i].temp; 375 uxx = (2.0 * u - x[j][i - 1].temp - x[j][i + 1].temp) * hydhx; 376 uyy = (2.0 * u - x[j - 1][i].temp - x[j + 1][i].temp) * hxdhy; 377 f[j][i].temp = (xdot[j][i].temp + uxx + uyy + prandtl * ((vxp * (u - x[j][i - 1].temp) + vxm * (x[j][i + 1].temp - u)) * hy + (vyp * (u - x[j - 1][i].temp) + vym * (x[j + 1][i].temp - u)) * hx)); 378 } 379 } 380 381 /* 382 Flop count (multiply-adds are counted as 2 operations) 383 */ 384 PetscCall(PetscLogFlops(84.0 * info->ym * info->xm)); 385 PetscFunctionReturn(0); 386 } 387 388 /*TEST 389 390 test: 391 args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts' 392 requires: !complex !single 393 394 test: 395 suffix: 2 396 nsize: 4 397 args: -da_grid_x 20 -da_grid_y 20 -lidvelocity 100 -grashof 1e3 -ts_max_steps 100 -ts_rtol 1e-3 -ts_atol 1e-3 -ts_type rosw -ts_rosw_type ra3pw -ts_monitor -ts_monitor_solution_vtk 'foo-%03d.vts' 398 requires: !complex !single 399 400 test: 401 suffix: 3 402 nsize: 4 403 args: -da_refine 2 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 -pc_type none -ts_type beuler -ts_monitor -snes_monitor_short -snes_type aspin -da_overlap 4 404 requires: !complex !single 405 406 test: 407 suffix: 4 408 nsize: 2 409 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 410 requires: !complex !single 411 412 test: 413 suffix: asm 414 nsize: 4 415 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 416 requires: !complex !single 417 418 TEST*/ 419