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