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