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 /*T 46 Concepts: TS^solving a system of nonlinear equations (parallel multicomponent example); 47 Concepts: DMDA^using distributed arrays; 48 Concepts: TS^multicomponent 49 Concepts: TS^differential-algebraic equation 50 Processors: n 51 T*/ 52 /* ------------------------------------------------------------------------ 53 54 We thank David E. Keyes for contributing the driven cavity discretization 55 within this example code. 56 57 This problem is modeled by the partial differential equation system 58 59 - Lap(U) - Grad_y(Omega) = 0 60 - Lap(V) + Grad_x(Omega) = 0 61 Omega_t - Lap(Omega) + Div([U*Omega,V*Omega]) - GR*Grad_x(T) = 0 62 T_t - Lap(T) + PR*Div([U*T,V*T]) = 0 63 64 in the unit square, which is uniformly discretized in each of x and 65 y in this simple encoding. 66 67 No-slip, rigid-wall Dirichlet conditions are used for [U,V]. 68 Dirichlet conditions are used for Omega, based on the definition of 69 vorticity: Omega = - Grad_y(U) + Grad_x(V), where along each 70 constant coordinate boundary, the tangential derivative is zero. 71 Dirichlet conditions are used for T on the left and right walls, 72 and insulation homogeneous Neumann conditions are used for T on 73 the top and bottom walls. 74 75 A finite difference approximation with the usual 5-point stencil 76 is used to discretize the boundary value problem to obtain a 77 nonlinear system of equations. Upwinding is used for the divergence 78 (convective) terms and central for the gradient (source) terms. 79 80 The Jacobian can be either 81 * formed via finite differencing using coloring (the default), or 82 * applied matrix-free via the option -snes_mf 83 (for larger grid problems this variant may not converge 84 without a preconditioner due to ill-conditioning). 85 86 ------------------------------------------------------------------------- */ 87 88 /* 89 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 90 Include "petscts.h" so that we can use TS solvers. Note that this 91 file automatically includes: 92 petscsys.h - base PETSc routines petscvec.h - vectors 93 petscmat.h - matrices 94 petscis.h - index sets petscksp.h - Krylov subspace methods 95 petscviewer.h - viewers petscpc.h - preconditioners 96 petscksp.h - linear solvers petscsnes.h - nonlinear solvers 97 */ 98 #include <petscts.h> 99 #include <petscdm.h> 100 #include <petscdmda.h> 101 102 /* 103 User-defined routines and data structures 104 */ 105 typedef struct { 106 PetscScalar u,v,omega,temp; 107 } Field; 108 109 PetscErrorCode FormIFunctionLocal(DMDALocalInfo*,PetscReal,Field**,Field**,Field**,void*); 110 111 typedef struct { 112 PetscReal lidvelocity,prandtl,grashof; /* physical parameters */ 113 PetscBool parabolic; /* allow a transient term corresponding roughly to artificial compressibility */ 114 PetscReal cfl_initial; /* CFL for first time step */ 115 } AppCtx; 116 117 PetscErrorCode FormInitialSolution(TS,Vec,AppCtx*); 118 119 int main(int argc,char **argv) 120 { 121 AppCtx user; /* user-defined work context */ 122 PetscInt mx,my,steps; 123 PetscErrorCode ierr; 124 TS ts; 125 DM da; 126 Vec X; 127 PetscReal ftime; 128 TSConvergedReason reason; 129 130 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 131 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 132 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da);CHKERRQ(ierr); 133 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 134 ierr = DMSetUp(da);CHKERRQ(ierr); 135 ierr = TSSetDM(ts,(DM)da);CHKERRQ(ierr); 136 137 ierr = DMDAGetInfo(da,0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 138 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 139 /* 140 Problem parameters (velocity of lid, prandtl, and grashof numbers) 141 */ 142 user.lidvelocity = 1.0/(mx*my); 143 user.prandtl = 1.0; 144 user.grashof = 1.0; 145 user.parabolic = PETSC_FALSE; 146 user.cfl_initial = 50.; 147 148 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Driven cavity/natural convection options","");CHKERRQ(ierr); 149 ierr = PetscOptionsReal("-lidvelocity","Lid velocity, related to Reynolds number","",user.lidvelocity,&user.lidvelocity,NULL);CHKERRQ(ierr); 150 ierr = PetscOptionsReal("-prandtl","Ratio of viscous to thermal diffusivity","",user.prandtl,&user.prandtl,NULL);CHKERRQ(ierr); 151 ierr = PetscOptionsReal("-grashof","Ratio of bouyant to viscous forces","",user.grashof,&user.grashof,NULL);CHKERRQ(ierr); 152 ierr = PetscOptionsBool("-parabolic","Relax incompressibility to make the system parabolic instead of differential-algebraic","",user.parabolic,&user.parabolic,NULL);CHKERRQ(ierr); 153 ierr = PetscOptionsReal("-cfl_initial","Advective CFL for the first time step","",user.cfl_initial,&user.cfl_initial,NULL);CHKERRQ(ierr); 154 ierr = PetscOptionsEnd();CHKERRQ(ierr); 155 156 ierr = DMDASetFieldName(da,0,"x-velocity");CHKERRQ(ierr); 157 ierr = DMDASetFieldName(da,1,"y-velocity");CHKERRQ(ierr); 158 ierr = DMDASetFieldName(da,2,"Omega");CHKERRQ(ierr); 159 ierr = DMDASetFieldName(da,3,"temperature");CHKERRQ(ierr); 160 161 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 162 Create user context, set problem data, create vector data structures. 163 Also, compute the initial guess. 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167 Create time integration context 168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 169 ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr); 170 ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)FormIFunctionLocal,&user);CHKERRQ(ierr); 171 ierr = TSSetMaxSteps(ts,10000);CHKERRQ(ierr); 172 ierr = TSSetMaxTime(ts,1e12);CHKERRQ(ierr); 173 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 174 ierr = TSSetTimeStep(ts,user.cfl_initial/(user.lidvelocity*mx));CHKERRQ(ierr); 175 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 176 177 ierr = 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);CHKERRQ(ierr); 178 179 180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 181 Solve the nonlinear system 182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 183 184 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr); 185 ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr); 186 187 ierr = TSSolve(ts,X);CHKERRQ(ierr); 188 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 189 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 190 ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr); 191 192 ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);CHKERRQ(ierr); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Free work space. All PETSc objects should be destroyed when they 196 are no longer needed. 197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 198 ierr = VecDestroy(&X);CHKERRQ(ierr); 199 ierr = DMDestroy(&da);CHKERRQ(ierr); 200 ierr = TSDestroy(&ts);CHKERRQ(ierr); 201 202 ierr = PetscFinalize(); 203 return ierr; 204 } 205 206 /* ------------------------------------------------------------------- */ 207 208 209 /* 210 FormInitialSolution - Forms initial approximation. 211 212 Input Parameters: 213 user - user-defined application context 214 X - vector 215 216 Output Parameter: 217 X - vector 218 */ 219 PetscErrorCode FormInitialSolution(TS ts,Vec X,AppCtx *user) 220 { 221 DM da; 222 PetscInt i,j,mx,xs,ys,xm,ym; 223 PetscErrorCode ierr; 224 PetscReal grashof,dx; 225 Field **x; 226 227 grashof = user->grashof; 228 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 229 ierr = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 230 dx = 1.0/(mx-1); 231 232 /* 233 Get local grid boundaries (for 2-dimensional DMDA): 234 xs, ys - starting grid indices (no ghost points) 235 xm, ym - widths of local grid (no ghost points) 236 */ 237 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 238 239 /* 240 Get a pointer to vector data. 241 - For default PETSc vectors, VecGetArray() returns a pointer to 242 the data array. Otherwise, the routine is implementation dependent. 243 - You MUST call VecRestoreArray() when you no longer need access to 244 the array. 245 */ 246 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 247 248 /* 249 Compute initial guess over the locally owned part of the grid 250 Initial condition is motionless fluid and equilibrium temperature 251 */ 252 for (j=ys; j<ys+ym; j++) { 253 for (i=xs; i<xs+xm; i++) { 254 x[j][i].u = 0.0; 255 x[j][i].v = 0.0; 256 x[j][i].omega = 0.0; 257 x[j][i].temp = (grashof>0)*i*dx; 258 } 259 } 260 261 /* 262 Restore vector 263 */ 264 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 265 return 0; 266 } 267 268 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscReal ptime,Field **x,Field **xdot,Field **f,void *ptr) 269 { 270 AppCtx *user = (AppCtx*)ptr; 271 PetscErrorCode ierr; 272 PetscInt xints,xinte,yints,yinte,i,j; 273 PetscReal hx,hy,dhx,dhy,hxdhy,hydhx; 274 PetscReal grashof,prandtl,lid; 275 PetscScalar u,udot,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym; 276 277 PetscFunctionBeginUser; 278 grashof = user->grashof; 279 prandtl = user->prandtl; 280 lid = user->lidvelocity; 281 282 /* 283 Define mesh intervals ratios for uniform grid. 284 285 Note: FD formulae below are normalized by multiplying through by 286 local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions. 287 288 289 */ 290 dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1); 291 hx = 1.0/dhx; hy = 1.0/dhy; 292 hxdhy = hx*dhy; hydhx = hy*dhx; 293 294 xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym; 295 296 /* Test whether we are on the bottom edge of the global array */ 297 if (yints == 0) { 298 j = 0; 299 yints = yints + 1; 300 /* bottom edge */ 301 for (i=info->xs; i<info->xs+info->xm; i++) { 302 f[j][i].u = x[j][i].u; 303 f[j][i].v = x[j][i].v; 304 f[j][i].omega = x[j][i].omega + (x[j+1][i].u - x[j][i].u)*dhy; 305 f[j][i].temp = x[j][i].temp-x[j+1][i].temp; 306 } 307 } 308 309 /* Test whether we are on the top edge of the global array */ 310 if (yinte == info->my) { 311 j = info->my - 1; 312 yinte = yinte - 1; 313 /* top edge */ 314 for (i=info->xs; i<info->xs+info->xm; i++) { 315 f[j][i].u = x[j][i].u - lid; 316 f[j][i].v = x[j][i].v; 317 f[j][i].omega = x[j][i].omega + (x[j][i].u - x[j-1][i].u)*dhy; 318 f[j][i].temp = x[j][i].temp-x[j-1][i].temp; 319 } 320 } 321 322 /* Test whether we are on the left edge of the global array */ 323 if (xints == 0) { 324 i = 0; 325 xints = xints + 1; 326 /* left edge */ 327 for (j=info->ys; j<info->ys+info->ym; j++) { 328 f[j][i].u = x[j][i].u; 329 f[j][i].v = x[j][i].v; 330 f[j][i].omega = x[j][i].omega - (x[j][i+1].v - x[j][i].v)*dhx; 331 f[j][i].temp = x[j][i].temp; 332 } 333 } 334 335 /* Test whether we are on the right edge of the global array */ 336 if (xinte == info->mx) { 337 i = info->mx - 1; 338 xinte = xinte - 1; 339 /* right edge */ 340 for (j=info->ys; j<info->ys+info->ym; j++) { 341 f[j][i].u = x[j][i].u; 342 f[j][i].v = x[j][i].v; 343 f[j][i].omega = x[j][i].omega - (x[j][i].v - x[j][i-1].v)*dhx; 344 f[j][i].temp = x[j][i].temp - (PetscReal)(grashof>0); 345 } 346 } 347 348 /* Compute over the interior points */ 349 for (j=yints; j<yinte; j++) { 350 for (i=xints; i<xinte; i++) { 351 352 /* 353 convective coefficients for upwinding 354 */ 355 vx = x[j][i].u; avx = PetscAbsScalar(vx); 356 vxp = .5*(vx+avx); vxm = .5*(vx-avx); 357 vy = x[j][i].v; avy = PetscAbsScalar(vy); 358 vyp = .5*(vy+avy); vym = .5*(vy-avy); 359 360 /* U velocity */ 361 u = x[j][i].u; 362 udot = user->parabolic ? xdot[j][i].u : 0.; 363 uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u)*hydhx; 364 uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u)*hxdhy; 365 f[j][i].u = udot + uxx + uyy - .5*(x[j+1][i].omega-x[j-1][i].omega)*hx; 366 367 /* V velocity */ 368 u = x[j][i].v; 369 udot = user->parabolic ? xdot[j][i].v : 0.; 370 uxx = (2.0*u - x[j][i-1].v - x[j][i+1].v)*hydhx; 371 uyy = (2.0*u - x[j-1][i].v - x[j+1][i].v)*hxdhy; 372 f[j][i].v = udot + uxx + uyy + .5*(x[j][i+1].omega-x[j][i-1].omega)*hy; 373 374 /* Omega */ 375 u = x[j][i].omega; 376 uxx = (2.0*u - x[j][i-1].omega - x[j][i+1].omega)*hydhx; 377 uyy = (2.0*u - x[j-1][i].omega - x[j+1][i].omega)*hxdhy; 378 f[j][i].omega = (xdot[j][i].omega + uxx + uyy 379 + (vxp*(u - x[j][i-1].omega) 380 + vxm*(x[j][i+1].omega - u)) * hy 381 + (vyp*(u - x[j-1][i].omega) 382 + vym*(x[j+1][i].omega - u)) * hx 383 - .5 * grashof * (x[j][i+1].temp - x[j][i-1].temp) * hy); 384 385 /* Temperature */ 386 u = x[j][i].temp; 387 uxx = (2.0*u - x[j][i-1].temp - x[j][i+1].temp)*hydhx; 388 uyy = (2.0*u - x[j-1][i].temp - x[j+1][i].temp)*hxdhy; 389 f[j][i].temp = (xdot[j][i].temp + uxx + uyy 390 + prandtl * ((vxp*(u - x[j][i-1].temp) 391 + vxm*(x[j][i+1].temp - u)) * hy 392 + (vyp*(u - x[j-1][i].temp) 393 + vym*(x[j+1][i].temp - u)) * hx)); 394 } 395 } 396 397 /* 398 Flop count (multiply-adds are counted as 2 operations) 399 */ 400 ierr = PetscLogFlops(84.0*info->ym*info->xm);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 /*TEST 405 406 test: 407 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' 408 requires: !complex !single 409 410 test: 411 suffix: 2 412 nsize: 4 413 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' 414 requires: !complex !single 415 416 test: 417 suffix: 3 418 nsize: 4 419 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 420 requires: !complex !single 421 422 test: 423 suffix: 4 424 nsize: 2 425 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 426 requires: !complex !single 427 428 test: 429 suffix: asm 430 nsize: 4 431 args: -da_refine 1 -lidvelocity 100 -grashof 1e3 -ts_max_steps 10 -ts_rtol 1e-3 -ts_atol 1e-3 432 requires: !complex !single 433 434 TEST*/ 435