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 { 114 AppCtx user; /* user-defined work context */ 115 PetscInt mx,my,steps; 116 PetscErrorCode ierr; 117 TS ts; 118 DM da; 119 Vec X; 120 PetscReal ftime; 121 TSConvergedReason reason; 122 123 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 124 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 125 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)); 126 PetscCall(DMSetFromOptions(da)); 127 PetscCall(DMSetUp(da)); 128 PetscCall(TSSetDM(ts,(DM)da)); 129 130 ierr = DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 131 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);PetscCall(ierr); 132 /* 133 Problem parameters (velocity of lid, prandtl, and grashof numbers) 134 */ 135 user.lidvelocity = 1.0/(mx*my); 136 user.prandtl = 1.0; 137 user.grashof = 1.0; 138 user.parabolic = PETSC_FALSE; 139 user.cfl_initial = 50.; 140 141 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Driven cavity/natural convection options","");PetscCall(ierr); 142 PetscCall(PetscOptionsReal("-lidvelocity","Lid velocity, related to Reynolds number","",user.lidvelocity,&user.lidvelocity,NULL)); 143 PetscCall(PetscOptionsReal("-prandtl","Ratio of viscous to thermal diffusivity","",user.prandtl,&user.prandtl,NULL)); 144 PetscCall(PetscOptionsReal("-grashof","Ratio of bouyant to viscous forces","",user.grashof,&user.grashof,NULL)); 145 PetscCall(PetscOptionsBool("-parabolic","Relax incompressibility to make the system parabolic instead of differential-algebraic","",user.parabolic,&user.parabolic,NULL)); 146 PetscCall(PetscOptionsReal("-cfl_initial","Advective CFL for the first time step","",user.cfl_initial,&user.cfl_initial,NULL)); 147 ierr = PetscOptionsEnd();PetscCall(ierr); 148 149 PetscCall(DMDASetFieldName(da,0,"x-velocity")); 150 PetscCall(DMDASetFieldName(da,1,"y-velocity")); 151 PetscCall(DMDASetFieldName(da,2,"Omega")); 152 PetscCall(DMDASetFieldName(da,3,"temperature")); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Create user context, set problem data, create vector data structures. 156 Also, compute the initial guess. 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 160 Create time integration context 161 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 162 PetscCall(DMSetApplicationContext(da,&user)); 163 PetscCall(DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&user)); 164 PetscCall(TSSetMaxSteps(ts,10000)); 165 PetscCall(TSSetMaxTime(ts,1e12)); 166 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 167 PetscCall(TSSetTimeStep(ts,user.cfl_initial/(user.lidvelocity*mx))); 168 PetscCall(TSSetFromOptions(ts)); 169 170 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%Dx%D grid, lid velocity = %g, prandtl # = %g, grashof # = %g\n",mx,my,(double)user.lidvelocity,(double)user.prandtl,(double)user.grashof)); 171 172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 Solve the nonlinear system 174 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 175 176 PetscCall(DMCreateGlobalVector(da,&X)); 177 PetscCall(FormInitialSolution(ts,X,&user)); 178 179 PetscCall(TSSolve(ts,X)); 180 PetscCall(TSGetSolveTime(ts,&ftime)); 181 PetscCall(TSGetStepNumber(ts,&steps)); 182 PetscCall(TSGetConvergedReason(ts,&reason)); 183 184 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps)); 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Free work space. All PETSc objects should be destroyed when they 188 are no longer needed. 189 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 190 PetscCall(VecDestroy(&X)); 191 PetscCall(DMDestroy(&da)); 192 PetscCall(TSDestroy(&ts)); 193 194 PetscCall(PetscFinalize()); 195 return 0; 196 } 197 198 /* ------------------------------------------------------------------- */ 199 200 /* 201 FormInitialSolution - Forms initial approximation. 202 203 Input Parameters: 204 user - user-defined application context 205 X - vector 206 207 Output Parameter: 208 X - vector 209 */ 210 PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user) 211 { 212 DM da; 213 PetscInt i,j,mx,xs,ys,xm,ym; 214 PetscReal grashof,dx; 215 Field **x; 216 217 grashof = user->grashof; 218 PetscCall(TSGetDM(ts,&da)); 219 PetscCall(DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0)); 220 dx = 1.0/(mx-1); 221 222 /* 223 Get local grid boundaries (for 2-dimensional DMDA): 224 xs, ys - starting grid indices (no ghost points) 225 xm, ym - widths of local grid (no ghost points) 226 */ 227 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 228 229 /* 230 Get a pointer to vector data. 231 - For default PETSc vectors, VecGetArray() returns a pointer to 232 the data array. Otherwise, the routine is implementation dependent. 233 - You MUST call VecRestoreArray() when you no longer need access to 234 the array. 235 */ 236 PetscCall(DMDAVecGetArray(da,X,&x)); 237 238 /* 239 Compute initial guess over the locally owned part of the grid 240 Initial condition is motionless fluid and equilibrium temperature 241 */ 242 for (j=ys; j<ys+ym; j++) { 243 for (i=xs; i<xs+xm; i++) { 244 x[j][i].u = 0.0; 245 x[j][i].v = 0.0; 246 x[j][i].omega = 0.0; 247 x[j][i].temp = (grashof>0)*i*dx; 248 } 249 } 250 251 /* 252 Restore vector 253 */ 254 PetscCall(DMDAVecRestoreArray(da,X,&x)); 255 return 0; 256 } 257 258 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr) 259 { 260 AppCtx *user = (AppCtx*)ptr; 261 PetscInt xints,xinte,yints,yinte,i,j; 262 PetscReal hx,hy,dhx,dhy,hxdhy,hydhx; 263 PetscReal grashof,prandtl,lid; 264 PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym; 265 266 PetscFunctionBeginUser; 267 grashof = user->grashof; 268 prandtl = user->prandtl; 269 lid = user->lidvelocity; 270 271 /* 272 Define mesh intervals ratios for uniform grid. 273 274 Note: FD formulae below are normalized by multiplying through by 275 local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions. 276 277 */ 278 dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1); 279 hx = 1.0/dhx; hy = 1.0/dhy; 280 hxdhy = hx*dhy; hydhx = hy*dhx; 281 282 xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym; 283 284 /* Test whether we are on the bottom edge of the global array */ 285 if (yints == 0) { 286 j = 0; 287 yints = yints + 1; 288 /* bottom edge */ 289 for (i=info->xs; i<info->xs+info->xm; i++) { 290 f[j][i].u = x[j][i].u; 291 f[j][i].v = x[j][i].v; 292 f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy; 293 f[j][i].temp = x[j][i].temp-x[j+1][i].temp; 294 } 295 } 296 297 /* Test whether we are on the top edge of the global array */ 298 if (yinte == info->my) { 299 j = info->my - 1; 300 yinte = yinte - 1; 301 /* top edge */ 302 for (i=info->xs; i<info->xs+info->xm; i++) { 303 f[j][i].u = x[j][i].u - lid; 304 f[j][i].v = x[j][i].v; 305 f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy; 306 f[j][i].temp = x[j][i].temp-x[j-1][i].temp; 307 } 308 } 309 310 /* Test whether we are on the left edge of the global array */ 311 if (xints == 0) { 312 i = 0; 313 xints = xints + 1; 314 /* left edge */ 315 for (j=info->ys; j<info->ys+info->ym; j++) { 316 f[j][i].u = x[j][i].u; 317 f[j][i].v = x[j][i].v; 318 f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx; 319 f[j][i].temp = x[j][i].temp; 320 } 321 } 322 323 /* Test whether we are on the right edge of the global array */ 324 if (xinte == info->mx) { 325 i = info->mx - 1; 326 xinte = xinte - 1; 327 /* right edge */ 328 for (j=info->ys; j<info->ys+info->ym; j++) { 329 f[j][i].u = x[j][i].u; 330 f[j][i].v = x[j][i].v; 331 f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx; 332 f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0); 333 } 334 } 335 336 /* Compute over the interior points */ 337 for (j=yints; j<yinte; j++) { 338 for (i=xints; i<xinte; i++) { 339 340 /* 341 convective coefficients for upwinding 342 */ 343 vx = x[j][i].u; avx = PetscAbsScalar(vx); 344 vxp = .5*(vx+avx); vxm = .5*(vx-avx); 345 vy = x[j][i].v; avy = PetscAbsScalar(vy); 346 vyp = .5*(vy+avy); vym = .5*(vy-avy); 347 348 /* U velocity */ 349 u = x[j][i].u; 350 udot = user->parabolic ? xdot[j][i].u : 0.; 351 uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx; 352 uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy; 353 f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx; 354 355 /* V velocity */ 356 u = x[j][i].v; 357 udot = user->parabolic ? xdot[j][i].v : 0.; 358 uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx; 359 uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy; 360 f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy; 361 362 /* Omega */ 363 u = x[j][i].omega; 364 uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx; 365 uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy; 366 f[j][i].omega = (xdot[j][i].omega + uxx + uyy 367 + (vxp*(u - x[j][i-1].omega) 368 + vxm*(x[j][i+1].omega - u)) * hy 369 + (vyp*(u - x[j-1][i].omega) 370 + vym*(x[j+1][i].omega - u)) * hx 371 - .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 378 + prandtl * ((vxp*(u - x[j][i-1].temp) 379 + vxm*(x[j][i+1].temp - u)) * hy 380 + (vyp*(u - x[j-1][i].temp) 381 + vym*(x[j+1][i].temp - u)) * hx)); 382 } 383 } 384 385 /* 386 Flop count (multiply-adds are counted as 2 operations) 387 */ 388 PetscCall(PetscLogFlops(84.0*info->ym*info->xm)); 389 PetscFunctionReturn(0); 390 } 391 392 /*TEST 393 394 test: 395 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' 396 requires: !complex !single 397 398 test: 399 suffix: 2 400 nsize: 4 401 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' 402 requires: !complex !single 403 404 test: 405 suffix: 3 406 nsize: 4 407 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 408 requires: !complex !single 409 410 test: 411 suffix: 4 412 nsize: 2 413 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 414 requires: !complex !single 415 416 test: 417 suffix: asm 418 nsize: 4 419 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 420 requires: !complex !single 421 422 TEST*/ 423