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