1 static const char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\ 2 The flow is driven by the subducting slab.\n\ 3 ---------------------------------ex30 help---------------------------------\n\ 4 -OPTION <DEFAULT> = (UNITS) DESCRIPTION.\n\n\ 5 -width <320> = (km) width of domain.\n\ 6 -depth <300> = (km) depth of domain.\n\ 7 -slab_dip <45> = (degrees) dip angle of the slab (determines the grid aspect ratio).\n\ 8 -lid_depth <35> = (km) depth of the static conductive lid.\n\ 9 -fault_depth <35> = (km) depth of slab-wedge mechanical coupling\n\ 10 (fault dept >= lid depth).\n\ 11 \n\ 12 -ni <82> = grid cells in x-direction. (nj adjusts to accommodate\n\ 13 the slab dip & depth). DO NOT USE -da_grid_x option!!!\n\ 14 -ivisc <3> = rheology option.\n\ 15 0 --- constant viscosity.\n\ 16 1 --- olivine diffusion creep rheology (T&P-dependent, newtonian).\n\ 17 2 --- olivine dislocation creep rheology (T&P-dependent, non-newtonian).\n\ 18 3 --- Full mantle rheology, combination of 1 & 2.\n\ 19 \n\ 20 -slab_velocity <5> = (cm/year) convergence rate of slab into subduction zone.\n\ 21 -slab_age <50> = (million yrs) age of slab for thermal profile boundary condition.\n\ 22 -lid_age <50> = (million yrs) age of lid for thermal profile boundary condition.\n\ 23 \n\ 24 FOR OTHER PARAMETER OPTIONS AND THEIR DEFAULT VALUES, see SetParams() in ex30.c.\n\ 25 ---------------------------------ex30 help---------------------------------\n"; 26 27 28 /*F----------------------------------------------------------------------- 29 30 This PETSc 2.2.0 example by Richard F. Katz 31 http://www.ldeo.columbia.edu/~katz/ 32 33 The problem is modeled by the partial differential equation system 34 35 \begin{eqnarray} 36 -\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\ 37 \nabla \cdot v & = & 0 \\ 38 dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\ 39 \end{eqnarray} 40 41 \begin{eqnarray} 42 \eta(T,Eps\_dot) & = & \hbox{constant } \hbox{if ivisc} ==0 \\ 43 & = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\ 44 & = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\ 45 & = & \hbox{mantle viscosity (difn and disl) } \hbox{if ivisc} ==3 46 \end{eqnarray} 47 48 which is uniformly discretized on a staggered mesh: 49 -------$w_{ij}$------ 50 $u_{i-1j}$ $P,T_{ij}$ $u_{ij}$ 51 ------$w_{ij-1}$----- 52 53 ------------------------------------------------------------------------F*/ 54 55 #include <petscsnes.h> 56 #include <petscdm.h> 57 #include <petscdmda.h> 58 59 #define VISC_CONST 0 60 #define VISC_DIFN 1 61 #define VISC_DISL 2 62 #define VISC_FULL 3 63 #define CELL_CENTER 0 64 #define CELL_CORNER 1 65 #define BC_ANALYTIC 0 66 #define BC_NOSTRESS 1 67 #define BC_EXPERMNT 2 68 #define ADVECT_FV 0 69 #define ADVECT_FROMM 1 70 #define PLATE_SLAB 0 71 #define PLATE_LID 1 72 #define EPS_ZERO 0.00000001 73 74 typedef struct { /* holds the variables to be solved for */ 75 PetscScalar u,w,p,T; 76 } Field; 77 78 typedef struct { /* parameters needed to compute viscosity */ 79 PetscReal A,n,Estar,Vstar; 80 } ViscParam; 81 82 typedef struct { /* physical and miscelaneous parameters */ 83 PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT; 84 PetscReal slab_dip, slab_age, slab_velocity, kappa, z_scale; 85 PetscReal c, d, sb, cb, skt, visc_cutoff, lid_age, eta0, continuation; 86 PetscReal L, V, lid_depth, fault_depth; 87 ViscParam diffusion, dislocation; 88 PetscInt ivisc, adv_scheme, ibound, output_ivisc; 89 PetscBool quiet, param_test, output_to_file, pv_analytic; 90 PetscBool interrupted, stop_solve, toggle_kspmon, kspmon; 91 char filename[PETSC_MAX_PATH_LEN]; 92 } Parameter; 93 94 typedef struct { /* grid parameters */ 95 DMBoundaryType bx,by; 96 DMDAStencilType stencil; 97 PetscInt corner,ni,nj,jlid,jfault,inose; 98 PetscInt dof,stencil_width,mglevels; 99 PetscReal dx,dz; 100 } GridInfo; 101 102 typedef struct { /* application context */ 103 Vec x,Xguess; 104 Parameter *param; 105 GridInfo *grid; 106 } AppCtx; 107 108 /* Callback functions (static interface) */ 109 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*); 110 111 /* Main routines */ 112 extern PetscErrorCode SetParams(Parameter*, GridInfo*); 113 extern PetscErrorCode ReportParams(Parameter*, GridInfo*); 114 extern PetscErrorCode Initialize(DM); 115 extern PetscErrorCode UpdateSolution(SNES,AppCtx*, PetscInt*); 116 extern PetscErrorCode DoOutput(SNES,PetscInt); 117 118 /* Post-processing & misc */ 119 extern PetscErrorCode ViscosityField(DM,Vec,Vec); 120 extern PetscErrorCode StressField(DM); 121 extern PetscErrorCode SNESConverged_Interactive(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason*, void*); 122 extern PetscErrorCode InteractiveHandler(int, void*); 123 124 /*-----------------------------------------------------------------------*/ 125 int main(int argc,char **argv) 126 /*-----------------------------------------------------------------------*/ 127 { 128 SNES snes; 129 AppCtx *user; /* user-defined work context */ 130 Parameter param; 131 GridInfo grid; 132 PetscInt nits; 133 PetscErrorCode ierr; 134 MPI_Comm comm; 135 DM da; 136 137 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 138 PetscOptionsSetValue(NULL,"-file","ex30_output"); 139 PetscOptionsSetValue(NULL,"-snes_monitor_short",NULL); 140 PetscOptionsSetValue(NULL,"-snes_max_it","20"); 141 PetscOptionsSetValue(NULL,"-ksp_max_it","1500"); 142 PetscOptionsSetValue(NULL,"-ksp_gmres_restart","300"); 143 PetscOptionsInsert(NULL,&argc,&argv,NULL); 144 145 comm = PETSC_COMM_WORLD; 146 147 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 148 Set up the problem parameters. 149 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 150 ierr = SetParams(¶m,&grid);CHKERRQ(ierr); 151 ierr = ReportParams(¶m,&grid);CHKERRQ(ierr); 152 153 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 155 ierr = SNESCreate(comm,&snes);CHKERRQ(ierr); 156 ierr = DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);CHKERRQ(ierr); 157 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 158 ierr = DMSetUp(da);CHKERRQ(ierr); 159 ierr = SNESSetDM(snes,da);CHKERRQ(ierr); 160 ierr = DMDASetFieldName(da,0,"x-velocity");CHKERRQ(ierr); 161 ierr = DMDASetFieldName(da,1,"y-velocity");CHKERRQ(ierr); 162 ierr = DMDASetFieldName(da,2,"pressure");CHKERRQ(ierr); 163 ierr = DMDASetFieldName(da,3,"temperature");CHKERRQ(ierr); 164 165 166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 167 Create user context, set problem data, create vector data structures. 168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 169 ierr = PetscNew(&user);CHKERRQ(ierr); 170 user->param = ¶m; 171 user->grid = &grid; 172 ierr = DMSetApplicationContext(da,user);CHKERRQ(ierr); 173 ierr = DMCreateGlobalVector(da,&(user->Xguess));CHKERRQ(ierr); 174 175 176 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177 Set up the SNES solver with callback functions. 178 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179 ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,(void*)user);CHKERRQ(ierr); 180 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 181 182 183 ierr = SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,NULL);CHKERRQ(ierr); 184 ierr = PetscPushSignalHandler(InteractiveHandler,(void*)user);CHKERRQ(ierr); 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Initialize and solve the nonlinear system 188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 189 ierr = Initialize(da);CHKERRQ(ierr); 190 ierr = UpdateSolution(snes,user,&nits);CHKERRQ(ierr); 191 192 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 193 Output variables. 194 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 195 ierr = DoOutput(snes,nits);CHKERRQ(ierr); 196 197 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 198 Free work space. 199 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 200 ierr = VecDestroy(&user->Xguess);CHKERRQ(ierr); 201 ierr = VecDestroy(&user->x);CHKERRQ(ierr); 202 ierr = PetscFree(user);CHKERRQ(ierr); 203 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 204 ierr = DMDestroy(&da);CHKERRQ(ierr); 205 ierr = PetscPopSignalHandler();CHKERRQ(ierr); 206 ierr = PetscFinalize(); 207 return ierr; 208 } 209 210 /*===================================================================== 211 PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve) 212 =====================================================================*/ 213 214 /*---------------------------------------------------------------------*/ 215 /* manages solve: adaptive continuation method */ 216 PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits) 217 { 218 KSP ksp; 219 PC pc; 220 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 221 Parameter *param = user->param; 222 PetscReal cont_incr=0.3; 223 PetscInt its; 224 PetscErrorCode ierr; 225 PetscBool q = PETSC_FALSE; 226 DM dm; 227 228 PetscFunctionBeginUser; 229 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 230 ierr = DMCreateGlobalVector(dm,&user->x);CHKERRQ(ierr); 231 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 232 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 233 ierr = KSPSetComputeSingularValues(ksp, PETSC_TRUE);CHKERRQ(ierr); 234 235 *nits=0; 236 237 /* Isoviscous solve */ 238 if (param->ivisc == VISC_CONST && !param->stop_solve) { 239 param->ivisc = VISC_CONST; 240 241 ierr = SNESSolve(snes,0,user->x);CHKERRQ(ierr); 242 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 243 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 244 *nits += its; 245 ierr = VecCopy(user->x,user->Xguess);CHKERRQ(ierr); 246 if (param->stop_solve) goto done; 247 } 248 249 /* Olivine diffusion creep */ 250 if (param->ivisc >= VISC_DIFN && !param->stop_solve) { 251 if (!q) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Computing Variable Viscosity Solution\n");CHKERRQ(ierr);} 252 253 /* continuation method on viscosity cutoff */ 254 for (param->continuation=0.0;; param->continuation+=cont_incr) { 255 if (!q) {ierr = PetscPrintf(PETSC_COMM_WORLD," Continuation parameter = %g\n", (double)param->continuation);CHKERRQ(ierr);} 256 257 /* solve the non-linear system */ 258 ierr = VecCopy(user->Xguess,user->x);CHKERRQ(ierr); 259 ierr = SNESSolve(snes,0,user->x);CHKERRQ(ierr); 260 ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr); 261 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 262 *nits += its; 263 if (!q) {ierr = PetscPrintf(PETSC_COMM_WORLD," SNES iterations: %D, Cumulative: %D\n", its, *nits);CHKERRQ(ierr);} 264 if (param->stop_solve) goto done; 265 266 if (reason<0) { 267 /* NOT converged */ 268 cont_incr = -PetscAbsReal(cont_incr)/2.0; 269 if (PetscAbsReal(cont_incr)<0.01) goto done; 270 271 } else { 272 /* converged */ 273 ierr = VecCopy(user->x,user->Xguess);CHKERRQ(ierr); 274 if (param->continuation >= 1.0) goto done; 275 if (its<=3) cont_incr = 0.30001; 276 else if (its<=8) cont_incr = 0.15001; 277 else cont_incr = 0.10001; 278 279 if (param->continuation+cont_incr > 1.0) cont_incr = 1.0 - param->continuation; 280 } /* endif reason<0 */ 281 } 282 } 283 done: 284 if (param->stop_solve && !q) {ierr = PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");CHKERRQ(ierr);} 285 if (reason<0 && !q) {ierr = PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");CHKERRQ(ierr);} 286 PetscFunctionReturn(0); 287 } 288 289 290 /*===================================================================== 291 PHYSICS FUNCTIONS (compute the discrete residual) 292 =====================================================================*/ 293 294 295 /*---------------------------------------------------------------------*/ 296 PETSC_STATIC_INLINE PetscScalar UInterp(Field **x, PetscInt i, PetscInt j) 297 /*---------------------------------------------------------------------*/ 298 { 299 return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u); 300 } 301 302 /*---------------------------------------------------------------------*/ 303 PETSC_STATIC_INLINE PetscScalar WInterp(Field **x, PetscInt i, PetscInt j) 304 /*---------------------------------------------------------------------*/ 305 { 306 return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w); 307 } 308 309 /*---------------------------------------------------------------------*/ 310 PETSC_STATIC_INLINE PetscScalar PInterp(Field **x, PetscInt i, PetscInt j) 311 /*---------------------------------------------------------------------*/ 312 { 313 return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p); 314 } 315 316 /*---------------------------------------------------------------------*/ 317 PETSC_STATIC_INLINE PetscScalar TInterp(Field **x, PetscInt i, PetscInt j) 318 /*---------------------------------------------------------------------*/ 319 { 320 return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T); 321 } 322 323 /*---------------------------------------------------------------------*/ 324 /* isoviscous analytic solution for IC */ 325 PETSC_STATIC_INLINE PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user) 326 /*---------------------------------------------------------------------*/ 327 { 328 Parameter *param = user->param; 329 GridInfo *grid = user->grid; 330 PetscScalar st, ct, th, c=param->c, d=param->d; 331 PetscReal x, z,r; 332 333 x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz; 334 r = PetscSqrtReal(x*x+z*z); 335 st = z/r; 336 ct = x/r; 337 th = PetscAtanReal(z/x); 338 return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st); 339 } 340 341 /*---------------------------------------------------------------------*/ 342 /* isoviscous analytic solution for IC */ 343 PETSC_STATIC_INLINE PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user) 344 /*---------------------------------------------------------------------*/ 345 { 346 Parameter *param = user->param; 347 GridInfo *grid = user->grid; 348 PetscScalar st, ct, th, c=param->c, d=param->d; 349 PetscReal x, z, r; 350 351 x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz; 352 r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r; th = PetscAtanReal(z/x); 353 return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st); 354 } 355 356 /*---------------------------------------------------------------------*/ 357 /* isoviscous analytic solution for IC */ 358 PETSC_STATIC_INLINE PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user) 359 /*---------------------------------------------------------------------*/ 360 { 361 Parameter *param = user->param; 362 GridInfo *grid = user->grid; 363 PetscScalar x, z, r, st, ct, c=param->c, d=param->d; 364 365 x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz; 366 r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r; 367 return (-2.0*(c*ct-d*st)/r); 368 } 369 370 /* computes the second invariant of the strain rate tensor */ 371 PETSC_STATIC_INLINE PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user) 372 /*---------------------------------------------------------------------*/ 373 { 374 Parameter *param = user->param; 375 GridInfo *grid = user->grid; 376 PetscInt ilim =grid->ni-1, jlim=grid->nj-1; 377 PetscScalar uN,uS,uE,uW,wN,wS,wE,wW; 378 PetscScalar eps11, eps12, eps22; 379 380 if (i<j) return EPS_ZERO; 381 if (i==ilim) i--; 382 if (j==jlim) j--; 383 384 if (ipos==CELL_CENTER) { /* on cell center */ 385 if (j<=grid->jlid) return EPS_ZERO; 386 387 uE = x[j][i].u; uW = x[j][i-1].u; 388 wN = x[j][i].w; wS = x[j-1][i].w; 389 wE = WInterp(x,i,j-1); 390 if (i==j) { 391 uN = param->cb; wW = param->sb; 392 } else { 393 uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1); 394 } 395 396 if (j==grid->jlid+1) uS = 0.0; 397 else uS = UInterp(x,i-1,j-1); 398 399 } else { /* on CELL_CORNER */ 400 if (j<grid->jlid) return EPS_ZERO; 401 402 uN = x[j+1][i].u; uS = x[j][i].u; 403 wE = x[j][i+1].w; wW = x[j][i].w; 404 if (i==j) { 405 wN = param->sb; 406 uW = param->cb; 407 } else { 408 wN = WInterp(x,i,j); 409 uW = UInterp(x,i-1,j); 410 } 411 412 if (j==grid->jlid) { 413 uE = 0.0; uW = 0.0; 414 uS = -uN; 415 wS = -wN; 416 } else { 417 uE = UInterp(x,i,j); 418 wS = WInterp(x,i,j-1); 419 } 420 } 421 422 eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz; 423 eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx); 424 425 return PetscSqrtReal(0.5*(eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22)); 426 } 427 428 /*---------------------------------------------------------------------*/ 429 /* computes the shear viscosity */ 430 PETSC_STATIC_INLINE PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param) 431 /*---------------------------------------------------------------------*/ 432 { 433 PetscReal result =0.0; 434 ViscParam difn =param->diffusion, disl=param->dislocation; 435 PetscInt iVisc =param->ivisc; 436 PetscScalar eps_scale=param->V/(param->L*1000.0); 437 PetscScalar strain_power, v1, v2, P; 438 PetscScalar rho_g = 32340.0, R=8.3144; 439 440 P = rho_g*(z*param->L*1000.0); /* Pa */ 441 442 if (iVisc==VISC_CONST) { 443 /* constant viscosity */ 444 return 1.0; 445 } else if (iVisc==VISC_DIFN) { 446 /* diffusion creep rheology */ 447 result = PetscRealPart((difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0)); 448 } else if (iVisc==VISC_DISL) { 449 /* dislocation creep rheology */ 450 strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n); 451 452 result = PetscRealPart(disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0); 453 } else if (iVisc==VISC_FULL) { 454 /* dislocation/diffusion creep rheology */ 455 strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n); 456 457 v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0; 458 v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0; 459 460 result = PetscRealPart(1.0/(1.0/v1 + 1.0/v2)); 461 } 462 463 /* max viscosity is param->eta0 */ 464 result = PetscMin(result, 1.0); 465 /* min viscosity is param->visc_cutoff */ 466 result = PetscMax(result, param->visc_cutoff); 467 /* continuation method */ 468 result = PetscPowReal(result,param->continuation); 469 return result; 470 } 471 472 /*---------------------------------------------------------------------*/ 473 /* computes the residual of the x-component of eqn (1) above */ 474 PETSC_STATIC_INLINE PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) 475 /*---------------------------------------------------------------------*/ 476 { 477 Parameter *param=user->param; 478 GridInfo *grid =user->grid; 479 PetscScalar dx = grid->dx, dz=grid->dz; 480 PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0; 481 PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale; 482 PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS; 483 PetscInt jlim = grid->nj-1; 484 485 z_scale = param->z_scale; 486 487 if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */ 488 TS = param->potentialT * TInterp(x,i,j-1) * PetscExpScalar((j-1.0)*dz*z_scale); 489 if (j==jlim) TN = TS; 490 else TN = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale); 491 TW = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale); 492 TE = param->potentialT * x[j][i+1].T * PetscExpScalar((j-0.5)*dz*z_scale); 493 if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */ 494 epsN = CalcSecInv(x,i,j, CELL_CORNER,user); 495 epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user); 496 epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user); 497 epsW = CalcSecInv(x,i,j, CELL_CENTER,user); 498 } 499 } 500 etaN = Viscosity(TN,epsN,dz*(j+0.5),param); 501 etaS = Viscosity(TS,epsS,dz*(j-0.5),param); 502 etaW = Viscosity(TW,epsW,dz*j,param); 503 etaE = Viscosity(TE,epsE,dz*j,param); 504 505 dPdx = (x[j][i+1].p - x[j][i].p)/dx; 506 if (j==jlim) dudzN = etaN * (x[j][i].w - x[j][i+1].w)/dx; 507 else dudzN = etaN * (x[j+1][i].u - x[j][i].u) /dz; 508 dudzS = etaS * (x[j][i].u - x[j-1][i].u)/dz; 509 dudxE = etaE * (x[j][i+1].u - x[j][i].u) /dx; 510 dudxW = etaW * (x[j][i].u - x[j][i-1].u)/dx; 511 512 residual = -dPdx /* X-MOMENTUM EQUATION*/ 513 +(dudxE - dudxW)/dx 514 +(dudzN - dudzS)/dz; 515 516 if (param->ivisc!=VISC_CONST) { 517 dwdxN = etaN * (x[j][i+1].w - x[j][i].w) /dx; 518 dwdxS = etaS * (x[j-1][i+1].w - x[j-1][i].w)/dx; 519 520 residual += (dudxE - dudxW)/dx + (dwdxN - dwdxS)/dz; 521 } 522 523 return residual; 524 } 525 526 /*---------------------------------------------------------------------*/ 527 /* computes the residual of the z-component of eqn (1) above */ 528 PETSC_STATIC_INLINE PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) 529 /*---------------------------------------------------------------------*/ 530 { 531 Parameter *param=user->param; 532 GridInfo *grid =user->grid; 533 PetscScalar dx = grid->dx, dz=grid->dz; 534 PetscScalar etaN =0.0,etaS=0.0,etaE=0.0,etaW=0.0,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0; 535 PetscScalar TE =0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale; 536 PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS; 537 PetscInt ilim = grid->ni-1; 538 539 /* geometric and other parameters */ 540 z_scale = param->z_scale; 541 542 /* viscosity */ 543 if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */ 544 TN = param->potentialT * x[j+1][i].T * PetscExpScalar((j+0.5)*dz*z_scale); 545 TS = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale); 546 TW = param->potentialT * TInterp(x,i-1,j) * PetscExpScalar(j*dz*z_scale); 547 if (i==ilim) TE = TW; 548 else TE = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale); 549 if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */ 550 epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user); 551 epsS = CalcSecInv(x,i,j, CELL_CENTER,user); 552 epsE = CalcSecInv(x,i,j, CELL_CORNER,user); 553 epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user); 554 } 555 } 556 etaN = Viscosity(TN,epsN,dz*(j+1.0),param); 557 etaS = Viscosity(TS,epsS,dz*(j+0.0),param); 558 etaW = Viscosity(TW,epsW,dz*(j+0.5),param); 559 etaE = Viscosity(TE,epsE,dz*(j+0.5),param); 560 561 dPdz = (x[j+1][i].p - x[j][i].p)/dz; 562 dwdzN = etaN * (x[j+1][i].w - x[j][i].w)/dz; 563 dwdzS = etaS * (x[j][i].w - x[j-1][i].w)/dz; 564 if (i==ilim) dwdxE = etaE * (x[j][i].u - x[j+1][i].u)/dz; 565 else dwdxE = etaE * (x[j][i+1].w - x[j][i].w) /dx; 566 dwdxW = 2.0*etaW * (x[j][i].w - x[j][i-1].w)/dx; 567 568 /* Z-MOMENTUM */ 569 residual = -dPdz /* constant viscosity terms */ 570 +(dwdzN - dwdzS)/dz 571 +(dwdxE - dwdxW)/dx; 572 573 if (param->ivisc!=VISC_CONST) { 574 dudzE = etaE * (x[j+1][i].u - x[j][i].u)/dz; 575 dudzW = etaW * (x[j+1][i-1].u - x[j][i-1].u)/dz; 576 577 residual += (dwdzN - dwdzS)/dz + (dudzE - dudzW)/dx; 578 } 579 580 return residual; 581 } 582 583 /*---------------------------------------------------------------------*/ 584 /* computes the residual of eqn (2) above */ 585 PETSC_STATIC_INLINE PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) 586 /*---------------------------------------------------------------------*/ 587 { 588 GridInfo *grid =user->grid; 589 PetscScalar uE,uW,wN,wS,dudx,dwdz; 590 591 uW = x[j][i-1].u; uE = x[j][i].u; dudx = (uE - uW)/grid->dx; 592 wS = x[j-1][i].w; wN = x[j][i].w; dwdz = (wN - wS)/grid->dz; 593 594 return dudx + dwdz; 595 } 596 597 /*---------------------------------------------------------------------*/ 598 /* computes the residual of eqn (3) above */ 599 PETSC_STATIC_INLINE PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user) 600 /*---------------------------------------------------------------------*/ 601 { 602 Parameter *param=user->param; 603 GridInfo *grid =user->grid; 604 PetscScalar dx = grid->dx, dz=grid->dz; 605 PetscInt ilim =grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid; 606 PetscScalar TE, TN, TS, TW, residual; 607 PetscScalar uE,uW,wN,wS; 608 PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS; 609 610 dTdzN = (x[j+1][i].T - x[j][i].T) /dz; 611 dTdzS = (x[j][i].T - x[j-1][i].T)/dz; 612 dTdxE = (x[j][i+1].T - x[j][i].T) /dx; 613 dTdxW = (x[j][i].T - x[j][i-1].T)/dx; 614 615 residual = ((dTdzN - dTdzS)/dz + /* diffusion term */ 616 (dTdxE - dTdxW)/dx)*dx*dz/param->peclet; 617 618 if (j<=jlid && i>=j) { 619 /* don't advect in the lid */ 620 return residual; 621 } else if (i<j) { 622 /* beneath the slab sfc */ 623 uW = uE = param->cb; 624 wS = wN = param->sb; 625 } else { 626 /* advect in the slab and wedge */ 627 uW = x[j][i-1].u; uE = x[j][i].u; 628 wS = x[j-1][i].w; wN = x[j][i].w; 629 } 630 631 if (param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1) { 632 /* finite volume advection */ 633 TS = (x[j][i].T + x[j-1][i].T)/2.0; 634 TN = (x[j][i].T + x[j+1][i].T)/2.0; 635 TE = (x[j][i].T + x[j][i+1].T)/2.0; 636 TW = (x[j][i].T + x[j][i-1].T)/2.0; 637 fN = wN*TN*dx; fS = wS*TS*dx; 638 fE = uE*TE*dz; fW = uW*TW*dz; 639 640 } else { 641 /* Fromm advection scheme */ 642 fE = (uE *(-x[j][i+2].T + 5.0*(x[j][i+1].T+x[j][i].T)-x[j][i-1].T)/8.0 643 - PetscAbsScalar(uE)*(-x[j][i+2].T + 3.0*(x[j][i+1].T-x[j][i].T)+x[j][i-1].T)/8.0)*dz; 644 fW = (uW *(-x[j][i+1].T + 5.0*(x[j][i].T+x[j][i-1].T)-x[j][i-2].T)/8.0 645 - PetscAbsScalar(uW)*(-x[j][i+1].T + 3.0*(x[j][i].T-x[j][i-1].T)+x[j][i-2].T)/8.0)*dz; 646 fN = (wN *(-x[j+2][i].T + 5.0*(x[j+1][i].T+x[j][i].T)-x[j-1][i].T)/8.0 647 - PetscAbsScalar(wN)*(-x[j+2][i].T + 3.0*(x[j+1][i].T-x[j][i].T)+x[j-1][i].T)/8.0)*dx; 648 fS = (wS *(-x[j+1][i].T + 5.0*(x[j][i].T+x[j-1][i].T)-x[j-2][i].T)/8.0 649 - PetscAbsScalar(wS)*(-x[j+1][i].T + 3.0*(x[j][i].T-x[j-1][i].T)+x[j-2][i].T)/8.0)*dx; 650 } 651 652 residual -= (fE - fW + fN - fS); 653 654 return residual; 655 } 656 657 /*---------------------------------------------------------------------*/ 658 /* computes the shear stress---used on the boundaries */ 659 PETSC_STATIC_INLINE PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user) 660 /*---------------------------------------------------------------------*/ 661 { 662 Parameter *param=user->param; 663 GridInfo *grid =user->grid; 664 PetscInt ilim =grid->ni-1, jlim=grid->nj-1; 665 PetscScalar uN, uS, wE, wW; 666 667 if (j<=grid->jlid || i<j || i==ilim || j==jlim) return EPS_ZERO; 668 669 if (ipos==CELL_CENTER) { /* on cell center */ 670 671 wE = WInterp(x,i,j-1); 672 if (i==j) { 673 wW = param->sb; 674 uN = param->cb; 675 } else { 676 wW = WInterp(x,i-1,j-1); 677 uN = UInterp(x,i-1,j); 678 } 679 if (j==grid->jlid+1) uS = 0.0; 680 else uS = UInterp(x,i-1,j-1); 681 682 } else { /* on cell corner */ 683 684 uN = x[j+1][i].u; uS = x[j][i].u; 685 wW = x[j][i].w; wE = x[j][i+1].w; 686 687 } 688 689 return (uN-uS)/grid->dz + (wE-wW)/grid->dx; 690 } 691 692 /*---------------------------------------------------------------------*/ 693 /* computes the normal stress---used on the boundaries */ 694 PETSC_STATIC_INLINE PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user) 695 /*---------------------------------------------------------------------*/ 696 { 697 Parameter *param=user->param; 698 GridInfo *grid =user->grid; 699 PetscScalar dx = grid->dx, dz=grid->dz; 700 PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc; 701 PetscScalar epsC =0.0, etaC, TC, uE, uW, pC, z_scale; 702 if (i<j || j<=grid->jlid) return EPS_ZERO; 703 704 ivisc=param->ivisc; z_scale = param->z_scale; 705 706 if (ipos==CELL_CENTER) { /* on cell center */ 707 708 TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale); 709 if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user); 710 etaC = Viscosity(TC,epsC,dz*j,param); 711 712 uW = x[j][i-1].u; uE = x[j][i].u; 713 pC = x[j][i].p; 714 715 } else { /* on cell corner */ 716 if (i==ilim || j==jlim) return EPS_ZERO; 717 718 TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale); 719 if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user); 720 etaC = Viscosity(TC,epsC,dz*(j+0.5),param); 721 722 if (i==j) uW = param->sb; 723 else uW = UInterp(x,i-1,j); 724 uE = UInterp(x,i,j); pC = PInterp(x,i,j); 725 } 726 727 return 2.0*etaC*(uE-uW)/dx - pC; 728 } 729 730 /*---------------------------------------------------------------------*/ 731 /* computes the normal stress---used on the boundaries */ 732 PETSC_STATIC_INLINE PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user) 733 /*---------------------------------------------------------------------*/ 734 { 735 Parameter *param=user->param; 736 GridInfo *grid =user->grid; 737 PetscScalar dz =grid->dz; 738 PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc; 739 PetscScalar epsC =0.0, etaC, TC; 740 PetscScalar pC, wN, wS, z_scale; 741 if (i<j || j<=grid->jlid) return EPS_ZERO; 742 743 ivisc=param->ivisc; z_scale = param->z_scale; 744 745 if (ipos==CELL_CENTER) { /* on cell center */ 746 747 TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale); 748 if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user); 749 etaC = Viscosity(TC,epsC,dz*j,param); 750 wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p; 751 752 } else { /* on cell corner */ 753 if ((i==ilim) || (j==jlim)) return EPS_ZERO; 754 755 TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale); 756 if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user); 757 etaC = Viscosity(TC,epsC,dz*(j+0.5),param); 758 if (i==j) wN = param->sb; 759 else wN = WInterp(x,i,j); 760 wS = WInterp(x,i,j-1); pC = PInterp(x,i,j); 761 } 762 763 return 2.0*etaC*(wN-wS)/dz - pC; 764 } 765 766 /*---------------------------------------------------------------------*/ 767 768 /*===================================================================== 769 INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS 770 =====================================================================*/ 771 772 /*---------------------------------------------------------------------*/ 773 /* initializes the problem parameters and checks for 774 command line changes */ 775 PetscErrorCode SetParams(Parameter *param, GridInfo *grid) 776 /*---------------------------------------------------------------------*/ 777 { 778 PetscErrorCode ierr, ierr_out=0; 779 PetscReal SEC_PER_YR = 3600.00*24.00*365.2500; 780 PetscReal alpha_g_on_cp_units_inverse_km=4.0e-5*9.8; 781 782 /* domain geometry */ 783 param->slab_dip = 45.0; 784 param->width = 320.0; /* km */ 785 param->depth = 300.0; /* km */ 786 param->lid_depth = 35.0; /* km */ 787 param->fault_depth = 35.0; /* km */ 788 789 ierr = PetscOptionsGetReal(NULL,NULL,"-slab_dip",&(param->slab_dip),NULL);CHKERRQ(ierr); 790 ierr = PetscOptionsGetReal(NULL,NULL,"-width",&(param->width),NULL);CHKERRQ(ierr); 791 ierr = PetscOptionsGetReal(NULL,NULL,"-depth",&(param->depth),NULL);CHKERRQ(ierr); 792 ierr = PetscOptionsGetReal(NULL,NULL,"-lid_depth",&(param->lid_depth),NULL);CHKERRQ(ierr); 793 ierr = PetscOptionsGetReal(NULL,NULL,"-fault_depth",&(param->fault_depth),NULL);CHKERRQ(ierr); 794 795 param->slab_dip = param->slab_dip*PETSC_PI/180.0; /* radians */ 796 797 /* grid information */ 798 ierr = PetscOptionsGetInt(NULL,NULL, "-jfault",&(grid->jfault),NULL);CHKERRQ(ierr); 799 grid->ni = 82; 800 ierr = PetscOptionsGetInt(NULL,NULL, "-ni",&(grid->ni),NULL);CHKERRQ(ierr); 801 802 grid->dx = param->width/((PetscReal)(grid->ni-2)); /* km */ 803 grid->dz = grid->dx*PetscTanReal(param->slab_dip); /* km */ 804 grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/ 805 param->depth = grid->dz*(grid->nj-2); /* km */ 806 grid->inose = 0; /* gridpoints*/ 807 ierr = PetscOptionsGetInt(NULL,NULL,"-inose",&(grid->inose),NULL);CHKERRQ(ierr); 808 grid->bx = DM_BOUNDARY_NONE; 809 grid->by = DM_BOUNDARY_NONE; 810 grid->stencil = DMDA_STENCIL_BOX; 811 grid->dof = 4; 812 grid->stencil_width = 2; 813 grid->mglevels = 1; 814 815 /* boundary conditions */ 816 param->pv_analytic = PETSC_FALSE; 817 param->ibound = BC_NOSTRESS; 818 ierr = PetscOptionsGetInt(NULL,NULL,"-ibound",&(param->ibound),NULL);CHKERRQ(ierr); 819 820 /* physical constants */ 821 param->slab_velocity = 5.0; /* cm/yr */ 822 param->slab_age = 50.0; /* Ma */ 823 param->lid_age = 50.0; /* Ma */ 824 param->kappa = 0.7272e-6; /* m^2/sec */ 825 param->potentialT = 1300.0; /* degrees C */ 826 827 ierr = PetscOptionsGetReal(NULL,NULL,"-slab_velocity",&(param->slab_velocity),NULL);CHKERRQ(ierr); 828 ierr = PetscOptionsGetReal(NULL,NULL,"-slab_age",&(param->slab_age),NULL);CHKERRQ(ierr); 829 ierr = PetscOptionsGetReal(NULL,NULL,"-lid_age",&(param->lid_age),NULL);CHKERRQ(ierr); 830 ierr = PetscOptionsGetReal(NULL,NULL,"-kappa",&(param->kappa),NULL);CHKERRQ(ierr); 831 ierr = PetscOptionsGetReal(NULL,NULL,"-potentialT",&(param->potentialT),NULL);CHKERRQ(ierr); 832 833 /* viscosity */ 834 param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */ 835 param->eta0 = 1e24; /* Pa-s */ 836 param->visc_cutoff = 0.0; /* factor of eta_0 */ 837 param->continuation = 1.0; 838 839 /* constants for diffusion creep */ 840 param->diffusion.A = 1.8e7; /* Pa-s */ 841 param->diffusion.n = 1.0; /* dim'less */ 842 param->diffusion.Estar = 375e3; /* J/mol */ 843 param->diffusion.Vstar = 5e-6; /* m^3/mol */ 844 845 /* constants for param->dislocationocation creep */ 846 param->dislocation.A = 2.8969e4; /* Pa-s */ 847 param->dislocation.n = 3.5; /* dim'less */ 848 param->dislocation.Estar = 530e3; /* J/mol */ 849 param->dislocation.Vstar = 14e-6; /* m^3/mol */ 850 851 ierr = PetscOptionsGetInt(NULL,NULL, "-ivisc",&(param->ivisc),NULL);CHKERRQ(ierr); 852 ierr = PetscOptionsGetReal(NULL,NULL,"-visc_cutoff",&(param->visc_cutoff),NULL);CHKERRQ(ierr); 853 854 param->output_ivisc = param->ivisc; 855 856 ierr = PetscOptionsGetInt(NULL,NULL,"-output_ivisc",&(param->output_ivisc),NULL);CHKERRQ(ierr); 857 ierr = PetscOptionsGetReal(NULL,NULL,"-vstar",&(param->dislocation.Vstar),NULL);CHKERRQ(ierr); 858 859 /* output options */ 860 param->quiet = PETSC_FALSE; 861 param->param_test = PETSC_FALSE; 862 863 ierr = PetscOptionsHasName(NULL,NULL,"-quiet",&(param->quiet));CHKERRQ(ierr); 864 ierr = PetscOptionsHasName(NULL,NULL,"-test",&(param->param_test));CHKERRQ(ierr); 865 ierr = PetscOptionsGetString(NULL,NULL,"-file",param->filename,sizeof(param->filename),&(param->output_to_file));CHKERRQ(ierr); 866 867 /* advection */ 868 param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */ 869 870 ierr = PetscOptionsGetInt(NULL,NULL,"-adv_scheme",&(param->adv_scheme),NULL);CHKERRQ(ierr); 871 872 /* misc. flags */ 873 param->stop_solve = PETSC_FALSE; 874 param->interrupted = PETSC_FALSE; 875 param->kspmon = PETSC_FALSE; 876 param->toggle_kspmon = PETSC_FALSE; 877 878 /* derived parameters for slab angle */ 879 param->sb = PetscSinReal(param->slab_dip); 880 param->cb = PetscCosReal(param->slab_dip); 881 param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb); 882 param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb); 883 884 /* length, velocity and time scale for non-dimensionalization */ 885 param->L = PetscMin(param->width,param->depth); /* km */ 886 param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */ 887 888 /* other unit conversions and derived parameters */ 889 param->scaled_width = param->width/param->L; /* dim'less */ 890 param->scaled_depth = param->depth/param->L; /* dim'less */ 891 param->lid_depth = param->lid_depth/param->L; /* dim'less */ 892 param->fault_depth = param->fault_depth/param->L; /* dim'less */ 893 grid->dx = grid->dx/param->L; /* dim'less */ 894 grid->dz = grid->dz/param->L; /* dim'less */ 895 grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */ 896 grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */ 897 param->lid_depth = grid->jlid*grid->dz; /* dim'less */ 898 param->fault_depth = grid->jfault*grid->dz; /* dim'less */ 899 grid->corner = grid->jlid+1; /* gridcells */ 900 param->peclet = param->V /* m/sec */ 901 * param->L*1000.0 /* m */ 902 / param->kappa; /* m^2/sec */ 903 param->z_scale = param->L * alpha_g_on_cp_units_inverse_km; 904 param->skt = PetscSqrtReal(param->kappa*param->slab_age*SEC_PER_YR); 905 ierr = PetscOptionsGetReal(NULL,NULL,"-peclet",&(param->peclet),NULL);CHKERRQ(ierr); 906 907 return ierr_out; 908 } 909 910 /*---------------------------------------------------------------------*/ 911 /* prints a report of the problem parameters to stdout */ 912 PetscErrorCode ReportParams(Parameter *param, GridInfo *grid) 913 /*---------------------------------------------------------------------*/ 914 { 915 PetscErrorCode ierr, ierr_out=0; 916 char date[30]; 917 918 ierr = PetscGetDate(date,30);CHKERRQ(ierr); 919 920 if (!(param->quiet)) { 921 ierr = PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");CHKERRQ(ierr); 922 ierr = PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");CHKERRQ(ierr); 923 ierr = PetscPrintf(PETSC_COMM_WORLD," Width = %g km, Depth = %g km\n",(double)param->width,(double)param->depth);CHKERRQ(ierr); 924 ierr = PetscPrintf(PETSC_COMM_WORLD," Slab dip = %g degrees, Slab velocity = %g cm/yr\n",(double)(param->slab_dip*180.0/PETSC_PI),(double)param->slab_velocity);CHKERRQ(ierr); 925 ierr = PetscPrintf(PETSC_COMM_WORLD," Lid depth = %5.2f km, Fault depth = %5.2f km\n",(double)(param->lid_depth*param->L),(double)(param->fault_depth*param->L));CHKERRQ(ierr); 926 927 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");CHKERRQ(ierr); 928 ierr = PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %D, %D [dx,dz] = %g, %g km\n",grid->ni,grid->nj,(double)(grid->dx*param->L),(double)(grid->dz*param->L));CHKERRQ(ierr); 929 ierr = PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);CHKERRQ(ierr); 930 ierr = PetscPrintf(PETSC_COMM_WORLD," Pe = %g\n",(double)param->peclet);CHKERRQ(ierr); 931 932 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");CHKERRQ(ierr); 933 if (param->ivisc==VISC_CONST) { 934 ierr = PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");CHKERRQ(ierr); 935 if (param->pv_analytic) { 936 ierr = PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");CHKERRQ(ierr); 937 } 938 } else if (param->ivisc==VISC_DIFN) { 939 ierr = PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");CHKERRQ(ierr); 940 ierr = PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));CHKERRQ(ierr); 941 } else if (param->ivisc==VISC_DISL) { 942 ierr = PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");CHKERRQ(ierr); 943 ierr = PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));CHKERRQ(ierr); 944 } else if (param->ivisc==VISC_FULL) { 945 ierr = PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");CHKERRQ(ierr); 946 ierr = PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));CHKERRQ(ierr); 947 } else { 948 ierr = PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");CHKERRQ(ierr); 949 ierr_out = 1; 950 } 951 952 ierr = PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");CHKERRQ(ierr); 953 if (param->ibound==BC_ANALYTIC) { 954 ierr = PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");CHKERRQ(ierr); 955 } else if (param->ibound==BC_NOSTRESS) { 956 ierr = PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");CHKERRQ(ierr); 957 } else if (param->ibound==BC_EXPERMNT) { 958 ierr = PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");CHKERRQ(ierr); 959 } else { 960 ierr = PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");CHKERRQ(ierr); 961 ierr_out = 1; 962 } 963 964 if (param->output_to_file) 965 #if defined(PETSC_HAVE_MATLAB_ENGINE) 966 ierr = PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);CHKERRQ(ierr); 967 #else 968 ierr = PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);CHKERRQ(ierr); 969 #endif 970 if (param->output_ivisc != param->ivisc) { 971 ierr = PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);CHKERRQ(ierr); 972 } 973 974 ierr = PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");CHKERRQ(ierr); 975 } 976 if (param->param_test) PetscEnd(); 977 return ierr_out; 978 } 979 980 /* ------------------------------------------------------------------- */ 981 /* generates an inital guess using the analytic solution for isoviscous 982 corner flow */ 983 PetscErrorCode Initialize(DM da) 984 /* ------------------------------------------------------------------- */ 985 { 986 AppCtx *user; 987 Parameter *param; 988 GridInfo *grid; 989 PetscInt i,j,is,js,im,jm; 990 PetscErrorCode ierr; 991 Field **x; 992 Vec Xguess; 993 994 /* Get the fine grid */ 995 ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 996 Xguess = user->Xguess; 997 param = user->param; 998 grid = user->grid; 999 ierr = DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);CHKERRQ(ierr); 1000 ierr = DMDAVecGetArray(da,Xguess,(void**)&x);CHKERRQ(ierr); 1001 1002 /* Compute initial guess */ 1003 for (j=js; j<js+jm; j++) { 1004 for (i=is; i<is+im; i++) { 1005 if (i<j) x[j][i].u = param->cb; 1006 else if (j<=grid->jlid) x[j][i].u = 0.0; 1007 else x[j][i].u = HorizVelocity(i,j,user); 1008 1009 if (i<=j) x[j][i].w = param->sb; 1010 else if (j<=grid->jlid) x[j][i].w = 0.0; 1011 else x[j][i].w = VertVelocity(i,j,user); 1012 1013 if (i<j || j<=grid->jlid) x[j][i].p = 0.0; 1014 else x[j][i].p = Pressure(i,j,user); 1015 1016 x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0); 1017 } 1018 } 1019 1020 /* Restore x to Xguess */ 1021 ierr = DMDAVecRestoreArray(da,Xguess,(void**)&x);CHKERRQ(ierr); 1022 1023 return 0; 1024 } 1025 1026 /*---------------------------------------------------------------------*/ 1027 /* controls output to a file */ 1028 PetscErrorCode DoOutput(SNES snes, PetscInt its) 1029 /*---------------------------------------------------------------------*/ 1030 { 1031 AppCtx *user; 1032 Parameter *param; 1033 GridInfo *grid; 1034 PetscInt ivt; 1035 PetscErrorCode ierr; 1036 PetscMPIInt rank; 1037 PetscViewer viewer; 1038 Vec res, pars; 1039 MPI_Comm comm; 1040 DM da; 1041 1042 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 1043 ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 1044 param = user->param; 1045 grid = user->grid; 1046 ivt = param->ivisc; 1047 1048 param->ivisc = param->output_ivisc; 1049 1050 /* compute final residual and final viscosity/strain rate fields */ 1051 ierr = SNESGetFunction(snes, &res, NULL, NULL);CHKERRQ(ierr); 1052 ierr = ViscosityField(da, user->x, user->Xguess);CHKERRQ(ierr); 1053 1054 /* get the communicator and the rank of the processor */ 1055 ierr = PetscObjectGetComm((PetscObject)snes, &comm);CHKERRQ(ierr); 1056 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 1057 1058 if (param->output_to_file) { /* send output to binary file */ 1059 ierr = VecCreate(comm, &pars);CHKERRQ(ierr); 1060 if (!rank) { /* on processor 0 */ 1061 ierr = VecSetSizes(pars, 20, PETSC_DETERMINE);CHKERRQ(ierr); 1062 ierr = VecSetFromOptions(pars);CHKERRQ(ierr); 1063 ierr = VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);CHKERRQ(ierr); 1064 ierr = VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);CHKERRQ(ierr); 1065 ierr = VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);CHKERRQ(ierr); 1066 ierr = VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);CHKERRQ(ierr); 1067 ierr = VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);CHKERRQ(ierr); 1068 ierr = VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);CHKERRQ(ierr); 1069 /* skipped 6 intentionally */ 1070 ierr = VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);CHKERRQ(ierr); 1071 ierr = VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);CHKERRQ(ierr); 1072 ierr = VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);CHKERRQ(ierr); 1073 ierr = VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);CHKERRQ(ierr); 1074 ierr = VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);CHKERRQ(ierr); 1075 ierr = VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);CHKERRQ(ierr); 1076 ierr = VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);CHKERRQ(ierr); 1077 ierr = VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);CHKERRQ(ierr); 1078 ierr = VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);CHKERRQ(ierr); 1079 ierr = VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);CHKERRQ(ierr); 1080 } else { /* on some other processor */ 1081 ierr = VecSetSizes(pars, 0, PETSC_DETERMINE);CHKERRQ(ierr); 1082 ierr = VecSetFromOptions(pars);CHKERRQ(ierr); 1083 } 1084 ierr = VecAssemblyBegin(pars);CHKERRQ(ierr); ierr = VecAssemblyEnd(pars);CHKERRQ(ierr); 1085 1086 /* create viewer */ 1087 #if defined(PETSC_HAVE_MATLAB_ENGINE) 1088 ierr = PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 1089 #else 1090 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 1091 #endif 1092 1093 /* send vectors to viewer */ 1094 ierr = PetscObjectSetName((PetscObject)res,"res");CHKERRQ(ierr); 1095 ierr = VecView(res,viewer);CHKERRQ(ierr); 1096 ierr = PetscObjectSetName((PetscObject)user->x,"out");CHKERRQ(ierr); 1097 ierr = VecView(user->x, viewer);CHKERRQ(ierr); 1098 ierr = PetscObjectSetName((PetscObject)(user->Xguess),"aux");CHKERRQ(ierr); 1099 ierr = VecView(user->Xguess, viewer);CHKERRQ(ierr); 1100 ierr = StressField(da);CHKERRQ(ierr); /* compute stress fields */ 1101 ierr = PetscObjectSetName((PetscObject)(user->Xguess),"str");CHKERRQ(ierr); 1102 ierr = VecView(user->Xguess, viewer);CHKERRQ(ierr); 1103 ierr = PetscObjectSetName((PetscObject)pars,"par");CHKERRQ(ierr); 1104 ierr = VecView(pars, viewer);CHKERRQ(ierr); 1105 1106 /* destroy viewer and vector */ 1107 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 1108 ierr = VecDestroy(&pars);CHKERRQ(ierr); 1109 } 1110 1111 param->ivisc = ivt; 1112 return 0; 1113 } 1114 1115 /* ------------------------------------------------------------------- */ 1116 /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */ 1117 PetscErrorCode ViscosityField(DM da, Vec X, Vec V) 1118 /* ------------------------------------------------------------------- */ 1119 { 1120 AppCtx *user; 1121 Parameter *param; 1122 GridInfo *grid; 1123 Vec localX; 1124 Field **v, **x; 1125 PetscReal eps, /* dx,*/ dz, T, epsC, TC; 1126 PetscInt i,j,is,js,im,jm,ilim,jlim,ivt; 1127 PetscErrorCode ierr; 1128 1129 PetscFunctionBeginUser; 1130 ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 1131 param = user->param; 1132 grid = user->grid; 1133 ivt = param->ivisc; 1134 param->ivisc = param->output_ivisc; 1135 1136 ierr = DMGetLocalVector(da, &localX);CHKERRQ(ierr); 1137 ierr = DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1138 ierr = DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1139 ierr = DMDAVecGetArray(da,localX,(void**)&x);CHKERRQ(ierr); 1140 ierr = DMDAVecGetArray(da,V,(void**)&v);CHKERRQ(ierr); 1141 1142 /* Parameters */ 1143 /* dx = grid->dx; */ dz = grid->dz; 1144 1145 ilim = grid->ni-1; jlim = grid->nj-1; 1146 1147 /* Compute real temperature, strain rate and viscosity */ 1148 ierr = DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);CHKERRQ(ierr); 1149 for (j=js; j<js+jm; j++) { 1150 for (i=is; i<is+im; i++) { 1151 T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*param->z_scale)); 1152 if (i<ilim && j<jlim) { 1153 TC = PetscRealPart(param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*param->z_scale)); 1154 } else { 1155 TC = T; 1156 } 1157 eps = PetscRealPart((CalcSecInv(x,i,j,CELL_CENTER,user))); 1158 epsC = PetscRealPart(CalcSecInv(x,i,j,CELL_CORNER,user)); 1159 1160 v[j][i].u = eps; 1161 v[j][i].w = epsC; 1162 v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param); 1163 v[j][i].T = Viscosity(TC,epsC,dz*j,param); 1164 } 1165 } 1166 ierr = DMDAVecRestoreArray(da,V,(void**)&v);CHKERRQ(ierr); 1167 ierr = DMDAVecRestoreArray(da,localX,(void**)&x);CHKERRQ(ierr); 1168 ierr = DMRestoreLocalVector(da, &localX);CHKERRQ(ierr); 1169 1170 param->ivisc = ivt; 1171 PetscFunctionReturn(0); 1172 } 1173 1174 /* ------------------------------------------------------------------- */ 1175 /* post-processing: compute stress everywhere */ 1176 PetscErrorCode StressField(DM da) 1177 /* ------------------------------------------------------------------- */ 1178 { 1179 AppCtx *user; 1180 PetscInt i,j,is,js,im,jm; 1181 PetscErrorCode ierr; 1182 Vec locVec; 1183 Field **x, **y; 1184 1185 ierr = DMGetApplicationContext(da,&user);CHKERRQ(ierr); 1186 1187 /* Get the fine grid of Xguess and X */ 1188 ierr = DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);CHKERRQ(ierr); 1189 ierr = DMDAVecGetArray(da,user->Xguess,(void**)&x);CHKERRQ(ierr); 1190 1191 ierr = DMGetLocalVector(da, &locVec);CHKERRQ(ierr); 1192 ierr = DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);CHKERRQ(ierr); 1193 ierr = DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);CHKERRQ(ierr); 1194 ierr = DMDAVecGetArray(da,locVec,(void**)&y);CHKERRQ(ierr); 1195 1196 /* Compute stress on the corner points */ 1197 for (j=js; j<js+jm; j++) { 1198 for (i=is; i<is+im; i++) { 1199 x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user); 1200 x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user); 1201 x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user); 1202 x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user); 1203 } 1204 } 1205 1206 /* Restore the fine grid of Xguess and X */ 1207 ierr = DMDAVecRestoreArray(da,user->Xguess,(void**)&x);CHKERRQ(ierr); 1208 ierr = DMDAVecRestoreArray(da,locVec,(void**)&y);CHKERRQ(ierr); 1209 ierr = DMRestoreLocalVector(da, &locVec);CHKERRQ(ierr); 1210 return 0; 1211 } 1212 1213 /*===================================================================== 1214 UTILITY FUNCTIONS 1215 =====================================================================*/ 1216 1217 /*---------------------------------------------------------------------*/ 1218 /* returns the velocity of the subducting slab and handles fault nodes 1219 for BC */ 1220 PETSC_STATIC_INLINE PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user) 1221 /*---------------------------------------------------------------------*/ 1222 { 1223 Parameter *param = user->param; 1224 GridInfo *grid = user->grid; 1225 1226 if (c=='U' || c=='u') { 1227 if (i<j-1) return param->cb; 1228 else if (j<=grid->jfault) return 0.0; 1229 else return param->cb; 1230 1231 } else { 1232 if (i<j) return param->sb; 1233 else if (j<=grid->jfault) return 0.0; 1234 else return param->sb; 1235 } 1236 } 1237 1238 /*---------------------------------------------------------------------*/ 1239 /* solution to diffusive half-space cooling model for BC */ 1240 PETSC_STATIC_INLINE PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user) 1241 /*---------------------------------------------------------------------*/ 1242 { 1243 Parameter *param = user->param; 1244 PetscScalar z; 1245 if (plate==PLATE_LID) z = (j-0.5)*user->grid->dz; 1246 else z = (j-0.5)*user->grid->dz*param->cb; /* PLATE_SLAB */ 1247 #if defined(PETSC_HAVE_ERF) 1248 return (PetscReal)(erf((double)PetscRealPart(z*param->L/2.0/param->skt))); 1249 #else 1250 (*PetscErrorPrintf)("erf() not available on this machine\n"); 1251 MPI_Abort(PETSC_COMM_SELF,1); 1252 #endif 1253 } 1254 1255 /*===================================================================== 1256 INTERACTIVE SIGNAL HANDLING 1257 =====================================================================*/ 1258 1259 /* ------------------------------------------------------------------- */ 1260 PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx) 1261 /* ------------------------------------------------------------------- */ 1262 { 1263 AppCtx *user = (AppCtx*) ctx; 1264 Parameter *param = user->param; 1265 KSP ksp; 1266 PetscErrorCode ierr; 1267 1268 PetscFunctionBeginUser; 1269 if (param->interrupted) { 1270 param->interrupted = PETSC_FALSE; 1271 ierr = PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");CHKERRQ(ierr); 1272 *reason = SNES_CONVERGED_FNORM_ABS; 1273 PetscFunctionReturn(0); 1274 } else if (param->toggle_kspmon) { 1275 param->toggle_kspmon = PETSC_FALSE; 1276 1277 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr); 1278 1279 if (param->kspmon) { 1280 ierr = KSPMonitorCancel(ksp);CHKERRQ(ierr); 1281 1282 param->kspmon = PETSC_FALSE; 1283 ierr = PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");CHKERRQ(ierr); 1284 } else { 1285 PetscViewerAndFormat *vf; 1286 ierr = PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);CHKERRQ(ierr); 1287 ierr = KSPMonitorSet(ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSingularValue,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);CHKERRQ(ierr); 1288 1289 param->kspmon = PETSC_TRUE; 1290 ierr = PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");CHKERRQ(ierr); 1291 } 1292 } 1293 PetscFunctionReturn(SNESConvergedDefault(snes,it,xnorm,snorm,fnorm,reason,ctx)); 1294 } 1295 1296 /* ------------------------------------------------------------------- */ 1297 #include <signal.h> 1298 PetscErrorCode InteractiveHandler(int signum, void *ctx) 1299 /* ------------------------------------------------------------------- */ 1300 { 1301 AppCtx *user = (AppCtx*) ctx; 1302 Parameter *param = user->param; 1303 1304 if (signum == SIGILL) { 1305 param->toggle_kspmon = PETSC_TRUE; 1306 #if !defined(PETSC_MISSING_SIGCONT) 1307 } else if (signum == SIGCONT) { 1308 param->interrupted = PETSC_TRUE; 1309 #endif 1310 #if !defined(PETSC_MISSING_SIGURG) 1311 } else if (signum == SIGURG) { 1312 param->stop_solve = PETSC_TRUE; 1313 #endif 1314 } 1315 return 0; 1316 } 1317 1318 /*---------------------------------------------------------------------*/ 1319 /* main call-back function that computes the processor-local piece 1320 of the residual */ 1321 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr) 1322 /*---------------------------------------------------------------------*/ 1323 { 1324 AppCtx *user = (AppCtx*)ptr; 1325 Parameter *param = user->param; 1326 GridInfo *grid = user->grid; 1327 PetscScalar mag_w, mag_u; 1328 PetscInt i,j,mx,mz,ilim,jlim; 1329 PetscInt is,ie,js,je,ibound; /* ,ivisc */ 1330 1331 PetscFunctionBeginUser; 1332 /* Define global and local grid parameters */ 1333 mx = info->mx; mz = info->my; 1334 ilim = mx-1; jlim = mz-1; 1335 is = info->xs; ie = info->xs+info->xm; 1336 js = info->ys; je = info->ys+info->ym; 1337 1338 /* Define geometric and numeric parameters */ 1339 /* ivisc = param->ivisc; */ ibound = param->ibound; 1340 1341 for (j=js; j<je; j++) { 1342 for (i=is; i<ie; i++) { 1343 1344 /************* X-MOMENTUM/VELOCITY *************/ 1345 if (i<j) f[j][i].u = x[j][i].u - SlabVel('U',i,j,user); 1346 else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) { 1347 /* in the lithospheric lid */ 1348 f[j][i].u = x[j][i].u - 0.0; 1349 } else if (i==ilim) { 1350 /* on the right side boundary */ 1351 if (ibound==BC_ANALYTIC) { 1352 f[j][i].u = x[j][i].u - HorizVelocity(i,j,user); 1353 } else { 1354 f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO; 1355 } 1356 1357 } else if (j==jlim) { 1358 /* on the bottom boundary */ 1359 if (ibound==BC_ANALYTIC) { 1360 f[j][i].u = x[j][i].u - HorizVelocity(i,j,user); 1361 } else if (ibound==BC_NOSTRESS) { 1362 f[j][i].u = XMomentumResidual(x,i,j,user); 1363 } else { 1364 /* experimental boundary condition */ 1365 } 1366 1367 } else { 1368 /* in the mantle wedge */ 1369 f[j][i].u = XMomentumResidual(x,i,j,user); 1370 } 1371 1372 /************* Z-MOMENTUM/VELOCITY *************/ 1373 if (i<=j) { 1374 f[j][i].w = x[j][i].w - SlabVel('W',i,j,user); 1375 1376 } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) { 1377 /* in the lithospheric lid */ 1378 f[j][i].w = x[j][i].w - 0.0; 1379 1380 } else if (j==jlim) { 1381 /* on the bottom boundary */ 1382 if (ibound==BC_ANALYTIC) { 1383 f[j][i].w = x[j][i].w - VertVelocity(i,j,user); 1384 } else { 1385 f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO; 1386 } 1387 1388 } else if (i==ilim) { 1389 /* on the right side boundary */ 1390 if (ibound==BC_ANALYTIC) { 1391 f[j][i].w = x[j][i].w - VertVelocity(i,j,user); 1392 } else if (ibound==BC_NOSTRESS) { 1393 f[j][i].w = ZMomentumResidual(x,i,j,user); 1394 } else { 1395 /* experimental boundary condition */ 1396 } 1397 1398 } else { 1399 /* in the mantle wedge */ 1400 f[j][i].w = ZMomentumResidual(x,i,j,user); 1401 } 1402 1403 /************* CONTINUITY/PRESSURE *************/ 1404 if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) { 1405 /* in the lid or slab */ 1406 f[j][i].p = x[j][i].p; 1407 1408 } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) { 1409 /* on an analytic boundary */ 1410 f[j][i].p = x[j][i].p - Pressure(i,j,user); 1411 1412 } else { 1413 /* in the mantle wedge */ 1414 f[j][i].p = ContinuityResidual(x,i,j,user); 1415 } 1416 1417 /************* TEMPERATURE *************/ 1418 if (j==0) { 1419 /* on the surface */ 1420 f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0); 1421 1422 } else if (i==0) { 1423 /* slab inflow boundary */ 1424 f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user); 1425 1426 } else if (i==ilim) { 1427 /* right side boundary */ 1428 mag_u = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j][i-1].u)/param->cb,1.0),0.0)), 5); 1429 f[j][i].T = x[j][i].T - mag_u*x[j-1][i-1].T - (1.0-mag_u)*PlateModel(j,PLATE_LID,user); 1430 1431 } else if (j==jlim) { 1432 /* bottom boundary */ 1433 mag_w = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j-1][i].w)/param->sb,1.0),0.0)), 5); 1434 f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w); 1435 1436 } else { 1437 /* in the mantle wedge */ 1438 f[j][i].T = EnergyResidual(x,i,j,user); 1439 } 1440 } 1441 } 1442 PetscFunctionReturn(0); 1443 } 1444 1445 1446 /*TEST 1447 1448 build: 1449 requires: !complex erf 1450 1451 test: 1452 args: -ni 18 1453 filter: grep -v Destination 1454 requires: !single 1455 1456 TEST*/ 1457