Actual source code: ex30.c

      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";

     27: /*F-----------------------------------------------------------------------

     29:     This PETSc 2.2.0 example by Richard F. Katz
     30:     http://www.ldeo.columbia.edu/~katz/

     32:     The problem is modeled by the partial differential equation system

     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}

     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}

     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}$-----

     52:   ------------------------------------------------------------------------F*/

     54: #include <petscsnes.h>
     55: #include <petscdm.h>
     56: #include <petscdmda.h>

     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

     73: typedef struct { /* holds the variables to be solved for */
     74:   PetscScalar u, w, p, T;
     75: } Field;

     77: typedef struct { /* parameters needed to compute viscosity */
     78:   PetscReal A, n, Estar, Vstar;
     79: } ViscParam;

     81: typedef struct { /* physical and miscellaneous 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;

     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;

    101: typedef struct { /* application context */
    102:   Vec        x, Xguess;
    103:   Parameter *param;
    104:   GridInfo  *grid;
    105: } AppCtx;

    107: /* Callback functions (static interface) */
    108: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field **, Field **, void *);

    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);

    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 *);

    123: int main(int argc, char **argv)
    124: {
    125:   SNES      snes;
    126:   AppCtx   *user; /* user-defined work context */
    127:   Parameter param;
    128:   GridInfo  grid;
    129:   PetscInt  nits;
    130:   MPI_Comm  comm;
    131:   DM        da;

    133:   PetscFunctionBeginUser;
    134:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
    135:   PetscCall(PetscOptionsSetValue(NULL, "-file", "ex30_output"));
    136:   PetscCall(PetscOptionsSetValue(NULL, "-snes_monitor_short", NULL));
    137:   PetscCall(PetscOptionsSetValue(NULL, "-snes_max_it", "20"));
    138:   PetscCall(PetscOptionsSetValue(NULL, "-ksp_max_it", "1500"));
    139:   PetscCall(PetscOptionsSetValue(NULL, "-ksp_gmres_restart", "300"));
    140:   PetscCall(PetscOptionsInsert(NULL, &argc, &argv, NULL));

    142:   comm = PETSC_COMM_WORLD;

    144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    145:      Set up the problem parameters.
    146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    147:   PetscCall(SetParams(&param, &grid));
    148:   PetscCall(ReportParams(&param, &grid));

    150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    151:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    152:   PetscCall(SNESCreate(comm, &snes));
    153:   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));
    154:   PetscCall(DMSetFromOptions(da));
    155:   PetscCall(DMSetUp(da));
    156:   PetscCall(SNESSetDM(snes, da));
    157:   PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
    158:   PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
    159:   PetscCall(DMDASetFieldName(da, 2, "pressure"));
    160:   PetscCall(DMDASetFieldName(da, 3, "temperature"));

    162:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    163:      Create user context, set problem data, create vector data structures.
    164:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    165:   PetscCall(PetscNew(&user));
    166:   user->param = &param;
    167:   user->grid  = &grid;
    168:   PetscCall(DMSetApplicationContext(da, user));
    169:   PetscCall(DMCreateGlobalVector(da, &user->Xguess));

    171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:      Set up the SNES solver with callback functions.
    173:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    174:   PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode (*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, (void *)user));
    175:   PetscCall(SNESSetFromOptions(snes));

    177:   PetscCall(SNESSetConvergenceTest(snes, SNESConverged_Interactive, (void *)user, NULL));
    178:   PetscCall(PetscPushSignalHandler(InteractiveHandler, (void *)user));

    180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    181:      Initialize and solve the nonlinear system
    182:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    183:   PetscCall(Initialize(da));
    184:   PetscCall(UpdateSolution(snes, user, &nits));

    186:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    187:      Output variables.
    188:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    189:   PetscCall(DoOutput(snes, nits));

    191:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    192:      Free work space.
    193:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    194:   PetscCall(VecDestroy(&user->Xguess));
    195:   PetscCall(VecDestroy(&user->x));
    196:   PetscCall(PetscFree(user));
    197:   PetscCall(SNESDestroy(&snes));
    198:   PetscCall(DMDestroy(&da));
    199:   PetscCall(PetscPopSignalHandler());
    200:   PetscCall(PetscFinalize());
    201:   return 0;
    202: }

    204: /*=====================================================================
    205:   PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
    206:   =====================================================================*/

    208: /*  manages solve: adaptive continuation method  */
    209: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
    210: {
    211:   KSP                 ksp;
    212:   PC                  pc;
    213:   SNESConvergedReason reason    = SNES_CONVERGED_ITERATING;
    214:   Parameter          *param     = user->param;
    215:   PetscReal           cont_incr = 0.3;
    216:   PetscInt            its;
    217:   PetscBool           q = PETSC_FALSE;
    218:   DM                  dm;

    220:   PetscFunctionBeginUser;
    221:   PetscCall(SNESGetDM(snes, &dm));
    222:   PetscCall(DMCreateGlobalVector(dm, &user->x));
    223:   PetscCall(SNESGetKSP(snes, &ksp));
    224:   PetscCall(KSPGetPC(ksp, &pc));
    225:   PetscCall(KSPSetComputeSingularValues(ksp, PETSC_TRUE));

    227:   *nits = 0;

    229:   /* Isoviscous solve */
    230:   if (param->ivisc == VISC_CONST && !param->stop_solve) {
    231:     param->ivisc = VISC_CONST;

    233:     PetscCall(SNESSolve(snes, 0, user->x));
    234:     PetscCall(SNESGetConvergedReason(snes, &reason));
    235:     PetscCall(SNESGetIterationNumber(snes, &its));
    236:     *nits += its;
    237:     PetscCall(VecCopy(user->x, user->Xguess));
    238:     if (param->stop_solve) goto done;
    239:   }

    241:   /* Olivine diffusion creep */
    242:   if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
    243:     if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Computing Variable Viscosity Solution\n"));

    245:     /* continuation method on viscosity cutoff */
    246:     for (param->continuation = 0.0;; param->continuation += cont_incr) {
    247:       if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Continuation parameter = %g\n", (double)param->continuation));

    249:       /* solve the non-linear system */
    250:       PetscCall(VecCopy(user->Xguess, user->x));
    251:       PetscCall(SNESSolve(snes, 0, user->x));
    252:       PetscCall(SNESGetConvergedReason(snes, &reason));
    253:       PetscCall(SNESGetIterationNumber(snes, &its));
    254:       *nits += its;
    255:       if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SNES iterations: %" PetscInt_FMT ", Cumulative: %" PetscInt_FMT "\n", its, *nits));
    256:       if (param->stop_solve) goto done;

    258:       if (reason < 0) {
    259:         /* NOT converged */
    260:         cont_incr = -PetscAbsReal(cont_incr) / 2.0;
    261:         if (PetscAbsReal(cont_incr) < 0.01) goto done;

    263:       } else {
    264:         /* converged */
    265:         PetscCall(VecCopy(user->x, user->Xguess));
    266:         if (param->continuation >= 1.0) goto done;
    267:         if (its <= 3) cont_incr = 0.30001;
    268:         else if (its <= 8) cont_incr = 0.15001;
    269:         else cont_incr = 0.10001;

    271:         if (param->continuation + cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
    272:       } /* endif reason<0 */
    273:     }
    274:   }
    275: done:
    276:   if (param->stop_solve && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: stopping solve.\n"));
    277:   if (reason < 0 && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "FAILED TO CONVERGE: stopping solve.\n"));
    278:   PetscFunctionReturn(PETSC_SUCCESS);
    279: }

    281: /*=====================================================================
    282:   PHYSICS FUNCTIONS (compute the discrete residual)
    283:   =====================================================================*/

    285: static inline PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
    286: {
    287:   return 0.25 * (x[j][i].u + x[j + 1][i].u + x[j][i + 1].u + x[j + 1][i + 1].u);
    288: }

    290: static inline PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
    291: {
    292:   return 0.25 * (x[j][i].w + x[j + 1][i].w + x[j][i + 1].w + x[j + 1][i + 1].w);
    293: }

    295: static inline PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
    296: {
    297:   return 0.25 * (x[j][i].p + x[j + 1][i].p + x[j][i + 1].p + x[j + 1][i + 1].p);
    298: }

    300: static inline PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
    301: {
    302:   return 0.25 * (x[j][i].T + x[j + 1][i].T + x[j][i + 1].T + x[j + 1][i + 1].T);
    303: }

    305: /*  isoviscous analytic solution for IC */
    306: static inline PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
    307: {
    308:   Parameter  *param = user->param;
    309:   GridInfo   *grid  = user->grid;
    310:   PetscScalar st, ct, th, c = param->c, d = param->d;
    311:   PetscReal   x, z, r;

    313:   x  = (i - grid->jlid) * grid->dx;
    314:   z  = (j - grid->jlid - 0.5) * grid->dz;
    315:   r  = PetscSqrtReal(x * x + z * z);
    316:   st = z / r;
    317:   ct = x / r;
    318:   th = PetscAtanReal(z / x);
    319:   return ct * (c * th * st + d * (st + th * ct)) + st * (c * (st - th * ct) + d * th * st);
    320: }

    322: /*  isoviscous analytic solution for IC */
    323: static inline PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)

    325: {
    326:   Parameter  *param = user->param;
    327:   GridInfo   *grid  = user->grid;
    328:   PetscScalar st, ct, th, c = param->c, d = param->d;
    329:   PetscReal   x, z, r;

    331:   x  = (i - grid->jlid - 0.5) * grid->dx;
    332:   z  = (j - grid->jlid) * grid->dz;
    333:   r  = PetscSqrtReal(x * x + z * z);
    334:   st = z / r;
    335:   ct = x / r;
    336:   th = PetscAtanReal(z / x);
    337:   return st * (c * th * st + d * (st + th * ct)) - ct * (c * (st - th * ct) + d * th * st);
    338: }

    340: /*  isoviscous analytic solution for IC */
    341: static inline PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
    342: {
    343:   Parameter  *param = user->param;
    344:   GridInfo   *grid  = user->grid;
    345:   PetscScalar x, z, r, st, ct, c = param->c, d = param->d;

    347:   x  = (i - grid->jlid - 0.5) * grid->dx;
    348:   z  = (j - grid->jlid - 0.5) * grid->dz;
    349:   r  = PetscSqrtReal(x * x + z * z);
    350:   st = z / r;
    351:   ct = x / r;
    352:   return -2.0 * (c * ct - d * st) / r;
    353: }

    355: /*  computes the second invariant of the strain rate tensor */
    356: static inline PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    357: {
    358:   Parameter  *param = user->param;
    359:   GridInfo   *grid  = user->grid;
    360:   PetscInt    ilim = grid->ni - 1, jlim = grid->nj - 1;
    361:   PetscScalar uN, uS, uE, uW, wN, wS, wE, wW;
    362:   PetscScalar eps11, eps12, eps22;

    364:   if (i < j) return EPS_ZERO;
    365:   if (i == ilim) i--;
    366:   if (j == jlim) j--;

    368:   if (ipos == CELL_CENTER) { /* on cell center */
    369:     if (j <= grid->jlid) return EPS_ZERO;

    371:     uE = x[j][i].u;
    372:     uW = x[j][i - 1].u;
    373:     wN = x[j][i].w;
    374:     wS = x[j - 1][i].w;
    375:     wE = WInterp(x, i, j - 1);
    376:     if (i == j) {
    377:       uN = param->cb;
    378:       wW = param->sb;
    379:     } else {
    380:       uN = UInterp(x, i - 1, j);
    381:       wW = WInterp(x, i - 1, j - 1);
    382:     }

    384:     if (j == grid->jlid + 1) uS = 0.0;
    385:     else uS = UInterp(x, i - 1, j - 1);

    387:   } else { /* on CELL_CORNER */
    388:     if (j < grid->jlid) return EPS_ZERO;

    390:     uN = x[j + 1][i].u;
    391:     uS = x[j][i].u;
    392:     wE = x[j][i + 1].w;
    393:     wW = x[j][i].w;
    394:     if (i == j) {
    395:       wN = param->sb;
    396:       uW = param->cb;
    397:     } else {
    398:       wN = WInterp(x, i, j);
    399:       uW = UInterp(x, i - 1, j);
    400:     }

    402:     if (j == grid->jlid) {
    403:       uE = 0.0;
    404:       uW = 0.0;
    405:       uS = -uN;
    406:       wS = -wN;
    407:     } else {
    408:       uE = UInterp(x, i, j);
    409:       wS = WInterp(x, i, j - 1);
    410:     }
    411:   }

    413:   eps11 = (uE - uW) / grid->dx;
    414:   eps22 = (wN - wS) / grid->dz;
    415:   eps12 = 0.5 * ((uN - uS) / grid->dz + (wE - wW) / grid->dx);

    417:   return PetscSqrtReal(0.5 * (eps11 * eps11 + 2.0 * eps12 * eps12 + eps22 * eps22));
    418: }

    420: /*  computes the shear viscosity */
    421: static inline PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
    422: {
    423:   PetscReal   result = 0.0;
    424:   ViscParam   difn = param->diffusion, disl = param->dislocation;
    425:   PetscInt    iVisc     = param->ivisc;
    426:   PetscScalar eps_scale = param->V / (param->L * 1000.0);
    427:   PetscScalar strain_power, v1, v2, P;
    428:   PetscScalar rho_g = 32340.0, R = 8.3144;

    430:   P = rho_g * (z * param->L * 1000.0); /* Pa */

    432:   if (iVisc == VISC_CONST) {
    433:     /* constant viscosity */
    434:     return 1.0;
    435:   } else if (iVisc == VISC_DIFN) {
    436:     /* diffusion creep rheology */
    437:     result = PetscRealPart(difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0);
    438:   } else if (iVisc == VISC_DISL) {
    439:     /* dislocation creep rheology */
    440:     strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);

    442:     result = PetscRealPart(disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0);
    443:   } else if (iVisc == VISC_FULL) {
    444:     /* dislocation/diffusion creep rheology */
    445:     strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);

    447:     v1 = difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0;
    448:     v2 = disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0;

    450:     result = PetscRealPart(1.0 / (1.0 / v1 + 1.0 / v2));
    451:   }

    453:   /* max viscosity is param->eta0 */
    454:   result = PetscMin(result, 1.0);
    455:   /* min viscosity is param->visc_cutoff */
    456:   result = PetscMax(result, param->visc_cutoff);
    457:   /* continuation method */
    458:   result = PetscPowReal(result, param->continuation);
    459:   return result;
    460: }

    462: /*  computes the residual of the x-component of eqn (1) above */
    463: static inline PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    464: {
    465:   Parameter  *param = user->param;
    466:   GridInfo   *grid  = user->grid;
    467:   PetscScalar dx = grid->dx, dz = grid->dz;
    468:   PetscScalar etaN, etaS, etaE, etaW, epsN = 0.0, epsS = 0.0, epsE = 0.0, epsW = 0.0;
    469:   PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdx, residual, z_scale;
    470:   PetscScalar dudxW, dudxE, dudzN, dudzS, dwdxN, dwdxS;
    471:   PetscInt    jlim = grid->nj - 1;

    473:   z_scale = param->z_scale;

    475:   if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
    476:     TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale);
    477:     if (j == jlim) TN = TS;
    478:     else TN = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
    479:     TW = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
    480:     TE = param->potentialT * x[j][i + 1].T * PetscExpScalar((j - 0.5) * dz * z_scale);
    481:     if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
    482:       epsN = CalcSecInv(x, i, j, CELL_CORNER, user);
    483:       epsS = CalcSecInv(x, i, j - 1, CELL_CORNER, user);
    484:       epsE = CalcSecInv(x, i + 1, j, CELL_CENTER, user);
    485:       epsW = CalcSecInv(x, i, j, CELL_CENTER, user);
    486:     }
    487:   }
    488:   etaN = Viscosity(TN, epsN, dz * (j + 0.5), param);
    489:   etaS = Viscosity(TS, epsS, dz * (j - 0.5), param);
    490:   etaW = Viscosity(TW, epsW, dz * j, param);
    491:   etaE = Viscosity(TE, epsE, dz * j, param);

    493:   dPdx = (x[j][i + 1].p - x[j][i].p) / dx;
    494:   if (j == jlim) dudzN = etaN * (x[j][i].w - x[j][i + 1].w) / dx;
    495:   else dudzN = etaN * (x[j + 1][i].u - x[j][i].u) / dz;
    496:   dudzS = etaS * (x[j][i].u - x[j - 1][i].u) / dz;
    497:   dudxE = etaE * (x[j][i + 1].u - x[j][i].u) / dx;
    498:   dudxW = etaW * (x[j][i].u - x[j][i - 1].u) / dx;

    500:   residual = -dPdx /* X-MOMENTUM EQUATION*/
    501:            + (dudxE - dudxW) / dx + (dudzN - dudzS) / dz;

    503:   if (param->ivisc != VISC_CONST) {
    504:     dwdxN = etaN * (x[j][i + 1].w - x[j][i].w) / dx;
    505:     dwdxS = etaS * (x[j - 1][i + 1].w - x[j - 1][i].w) / dx;

    507:     residual += (dudxE - dudxW) / dx + (dwdxN - dwdxS) / dz;
    508:   }

    510:   return residual;
    511: }

    513: /*  computes the residual of the z-component of eqn (1) above */
    514: static inline PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)

    516: {
    517:   Parameter  *param = user->param;
    518:   GridInfo   *grid  = user->grid;
    519:   PetscScalar dx = grid->dx, dz = grid->dz;
    520:   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;
    521:   PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdz, residual, z_scale;
    522:   PetscScalar dudzE, dudzW, dwdxW, dwdxE, dwdzN, dwdzS;
    523:   PetscInt    ilim = grid->ni - 1;

    525:   /* geometric and other parameters */
    526:   z_scale = param->z_scale;

    528:   /* viscosity */
    529:   if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
    530:     TN = param->potentialT * x[j + 1][i].T * PetscExpScalar((j + 0.5) * dz * z_scale);
    531:     TS = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
    532:     TW = param->potentialT * TInterp(x, i - 1, j) * PetscExpScalar(j * dz * z_scale);
    533:     if (i == ilim) TE = TW;
    534:     else TE = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
    535:     if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
    536:       epsN = CalcSecInv(x, i, j + 1, CELL_CENTER, user);
    537:       epsS = CalcSecInv(x, i, j, CELL_CENTER, user);
    538:       epsE = CalcSecInv(x, i, j, CELL_CORNER, user);
    539:       epsW = CalcSecInv(x, i - 1, j, CELL_CORNER, user);
    540:     }
    541:   }
    542:   etaN = Viscosity(TN, epsN, dz * (j + 1.0), param);
    543:   etaS = Viscosity(TS, epsS, dz * (j + 0.0), param);
    544:   etaW = Viscosity(TW, epsW, dz * (j + 0.5), param);
    545:   etaE = Viscosity(TE, epsE, dz * (j + 0.5), param);

    547:   dPdz  = (x[j + 1][i].p - x[j][i].p) / dz;
    548:   dwdzN = etaN * (x[j + 1][i].w - x[j][i].w) / dz;
    549:   dwdzS = etaS * (x[j][i].w - x[j - 1][i].w) / dz;
    550:   if (i == ilim) dwdxE = etaE * (x[j][i].u - x[j + 1][i].u) / dz;
    551:   else dwdxE = etaE * (x[j][i + 1].w - x[j][i].w) / dx;
    552:   dwdxW = 2.0 * etaW * (x[j][i].w - x[j][i - 1].w) / dx;

    554:   /* Z-MOMENTUM */
    555:   residual = -dPdz /* constant viscosity terms */
    556:            + (dwdzN - dwdzS) / dz + (dwdxE - dwdxW) / dx;

    558:   if (param->ivisc != VISC_CONST) {
    559:     dudzE = etaE * (x[j + 1][i].u - x[j][i].u) / dz;
    560:     dudzW = etaW * (x[j + 1][i - 1].u - x[j][i - 1].u) / dz;

    562:     residual += (dwdzN - dwdzS) / dz + (dudzE - dudzW) / dx;
    563:   }

    565:   return residual;
    566: }

    568: /*  computes the residual of eqn (2) above */
    569: static inline PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    570: {
    571:   GridInfo   *grid = user->grid;
    572:   PetscScalar uE, uW, wN, wS, dudx, dwdz;

    574:   uW   = x[j][i - 1].u;
    575:   uE   = x[j][i].u;
    576:   dudx = (uE - uW) / grid->dx;
    577:   wS   = x[j - 1][i].w;
    578:   wN   = x[j][i].w;
    579:   dwdz = (wN - wS) / grid->dz;

    581:   return dudx + dwdz;
    582: }

    584: /*  computes the residual of eqn (3) above */
    585: static inline PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    586: {
    587:   Parameter  *param = user->param;
    588:   GridInfo   *grid  = user->grid;
    589:   PetscScalar dx = grid->dx, dz = grid->dz;
    590:   PetscInt    ilim = grid->ni - 1, jlim = grid->nj - 1, jlid = grid->jlid;
    591:   PetscScalar TE, TN, TS, TW, residual;
    592:   PetscScalar uE, uW, wN, wS;
    593:   PetscScalar fN, fS, fE, fW, dTdxW, dTdxE, dTdzN, dTdzS;

    595:   dTdzN = (x[j + 1][i].T - x[j][i].T) / dz;
    596:   dTdzS = (x[j][i].T - x[j - 1][i].T) / dz;
    597:   dTdxE = (x[j][i + 1].T - x[j][i].T) / dx;
    598:   dTdxW = (x[j][i].T - x[j][i - 1].T) / dx;

    600:   residual = ((dTdzN - dTdzS) / dz + /* diffusion term */
    601:               (dTdxE - dTdxW) / dx) *
    602:              dx * dz / param->peclet;

    604:   if (j <= jlid && i >= j) {
    605:     /* don't advect in the lid */
    606:     return residual;
    607:   } else if (i < j) {
    608:     /* beneath the slab sfc */
    609:     uW = uE = param->cb;
    610:     wS = wN = param->sb;
    611:   } else {
    612:     /* advect in the slab and wedge */
    613:     uW = x[j][i - 1].u;
    614:     uE = x[j][i].u;
    615:     wS = x[j - 1][i].w;
    616:     wN = x[j][i].w;
    617:   }

    619:   if (param->adv_scheme == ADVECT_FV || i == ilim - 1 || j == jlim - 1 || i == 1 || j == 1) {
    620:     /* finite volume advection */
    621:     TS = (x[j][i].T + x[j - 1][i].T) / 2.0;
    622:     TN = (x[j][i].T + x[j + 1][i].T) / 2.0;
    623:     TE = (x[j][i].T + x[j][i + 1].T) / 2.0;
    624:     TW = (x[j][i].T + x[j][i - 1].T) / 2.0;
    625:     fN = wN * TN * dx;
    626:     fS = wS * TS * dx;
    627:     fE = uE * TE * dz;
    628:     fW = uW * TW * dz;

    630:   } else {
    631:     /* Fromm advection scheme */
    632:     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 - 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;
    633:     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 - 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;
    634:     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 - 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;
    635:     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 - 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;
    636:   }

    638:   residual -= (fE - fW + fN - fS);

    640:   return residual;
    641: }

    643: /*  computes the shear stress---used on the boundaries */
    644: static inline PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    645: {
    646:   Parameter  *param = user->param;
    647:   GridInfo   *grid  = user->grid;
    648:   PetscInt    ilim = grid->ni - 1, jlim = grid->nj - 1;
    649:   PetscScalar uN, uS, wE, wW;

    651:   if (j <= grid->jlid || i < j || i == ilim || j == jlim) return EPS_ZERO;

    653:   if (ipos == CELL_CENTER) { /* on cell center */

    655:     wE = WInterp(x, i, j - 1);
    656:     if (i == j) {
    657:       wW = param->sb;
    658:       uN = param->cb;
    659:     } else {
    660:       wW = WInterp(x, i - 1, j - 1);
    661:       uN = UInterp(x, i - 1, j);
    662:     }
    663:     if (j == grid->jlid + 1) uS = 0.0;
    664:     else uS = UInterp(x, i - 1, j - 1);

    666:   } else { /* on cell corner */

    668:     uN = x[j + 1][i].u;
    669:     uS = x[j][i].u;
    670:     wW = x[j][i].w;
    671:     wE = x[j][i + 1].w;
    672:   }

    674:   return (uN - uS) / grid->dz + (wE - wW) / grid->dx;
    675: }

    677: /*  computes the normal stress---used on the boundaries */
    678: static inline PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    679: {
    680:   Parameter  *param = user->param;
    681:   GridInfo   *grid  = user->grid;
    682:   PetscScalar dx = grid->dx, dz = grid->dz;
    683:   PetscInt    ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
    684:   PetscScalar epsC = 0.0, etaC, TC, uE, uW, pC, z_scale;
    685:   if (i < j || j <= grid->jlid) return EPS_ZERO;

    687:   ivisc   = param->ivisc;
    688:   z_scale = param->z_scale;

    690:   if (ipos == CELL_CENTER) { /* on cell center */

    692:     TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
    693:     if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
    694:     etaC = Viscosity(TC, epsC, dz * j, param);

    696:     uW = x[j][i - 1].u;
    697:     uE = x[j][i].u;
    698:     pC = x[j][i].p;

    700:   } else { /* on cell corner */
    701:     if (i == ilim || j == jlim) return EPS_ZERO;

    703:     TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
    704:     if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
    705:     etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);

    707:     if (i == j) uW = param->sb;
    708:     else uW = UInterp(x, i - 1, j);
    709:     uE = UInterp(x, i, j);
    710:     pC = PInterp(x, i, j);
    711:   }

    713:   return 2.0 * etaC * (uE - uW) / dx - pC;
    714: }

    716: /*  computes the normal stress---used on the boundaries */
    717: static inline PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    718: {
    719:   Parameter  *param = user->param;
    720:   GridInfo   *grid  = user->grid;
    721:   PetscScalar dz    = grid->dz;
    722:   PetscInt    ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
    723:   PetscScalar epsC = 0.0, etaC, TC;
    724:   PetscScalar pC, wN, wS, z_scale;
    725:   if (i < j || j <= grid->jlid) return EPS_ZERO;

    727:   ivisc   = param->ivisc;
    728:   z_scale = param->z_scale;

    730:   if (ipos == CELL_CENTER) { /* on cell center */

    732:     TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
    733:     if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
    734:     etaC = Viscosity(TC, epsC, dz * j, param);
    735:     wN   = x[j][i].w;
    736:     wS   = x[j - 1][i].w;
    737:     pC   = x[j][i].p;

    739:   } else { /* on cell corner */
    740:     if ((i == ilim) || (j == jlim)) return EPS_ZERO;

    742:     TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
    743:     if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
    744:     etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);
    745:     if (i == j) wN = param->sb;
    746:     else wN = WInterp(x, i, j);
    747:     wS = WInterp(x, i, j - 1);
    748:     pC = PInterp(x, i, j);
    749:   }

    751:   return 2.0 * etaC * (wN - wS) / dz - pC;
    752: }

    754: /*=====================================================================
    755:   INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
    756:   =====================================================================*/

    758: /* initializes the problem parameters and checks for
    759:    command line changes */
    760: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
    761: {
    762:   PetscReal SEC_PER_YR                     = 3600.00 * 24.00 * 365.2500;
    763:   PetscReal alpha_g_on_cp_units_inverse_km = 4.0e-5 * 9.8;

    765:   PetscFunctionBeginUser;
    766:   /* domain geometry */
    767:   param->slab_dip    = 45.0;
    768:   param->width       = 320.0; /* km */
    769:   param->depth       = 300.0; /* km */
    770:   param->lid_depth   = 35.0;  /* km */
    771:   param->fault_depth = 35.0;  /* km */

    773:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_dip", &param->slab_dip, NULL));
    774:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", &param->width, NULL));
    775:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-depth", &param->depth, NULL));
    776:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_depth", &param->lid_depth, NULL));
    777:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-fault_depth", &param->fault_depth, NULL));

    779:   param->slab_dip = param->slab_dip * PETSC_PI / 180.0; /* radians */

    781:   /* grid information */
    782:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-jfault", &grid->jfault, NULL));
    783:   grid->ni = 82;
    784:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ni", &grid->ni, NULL));

    786:   grid->dx     = param->width / ((PetscReal)(grid->ni - 2)); /* km */
    787:   grid->dz     = grid->dx * PetscTanReal(param->slab_dip);   /* km */
    788:   grid->nj     = (PetscInt)(param->depth / grid->dz + 3.0);  /* gridpoints*/
    789:   param->depth = grid->dz * (grid->nj - 2);                  /* km */
    790:   grid->inose  = 0;                                          /* gridpoints*/
    791:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-inose", &grid->inose, NULL));
    792:   grid->bx            = DM_BOUNDARY_NONE;
    793:   grid->by            = DM_BOUNDARY_NONE;
    794:   grid->stencil       = DMDA_STENCIL_BOX;
    795:   grid->dof           = 4;
    796:   grid->stencil_width = 2;
    797:   grid->mglevels      = 1;

    799:   /* boundary conditions */
    800:   param->pv_analytic = PETSC_FALSE;
    801:   param->ibound      = BC_NOSTRESS;
    802:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ibound", &param->ibound, NULL));

    804:   /* physical constants */
    805:   param->slab_velocity = 5.0;       /* cm/yr */
    806:   param->slab_age      = 50.0;      /* Ma */
    807:   param->lid_age       = 50.0;      /* Ma */
    808:   param->kappa         = 0.7272e-6; /* m^2/sec */
    809:   param->potentialT    = 1300.0;    /* degrees C */

    811:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_velocity", &param->slab_velocity, NULL));
    812:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_age", &param->slab_age, NULL));
    813:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_age", &param->lid_age, NULL));
    814:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &param->kappa, NULL));
    815:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-potentialT", &param->potentialT, NULL));

    817:   /* viscosity */
    818:   param->ivisc        = 3;    /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
    819:   param->eta0         = 1e24; /* Pa-s */
    820:   param->visc_cutoff  = 0.0;  /* factor of eta_0 */
    821:   param->continuation = 1.0;

    823:   /* constants for diffusion creep */
    824:   param->diffusion.A     = 1.8e7; /* Pa-s */
    825:   param->diffusion.n     = 1.0;   /* dim'less */
    826:   param->diffusion.Estar = 375e3; /* J/mol */
    827:   param->diffusion.Vstar = 5e-6;  /* m^3/mol */

    829:   /* constants for param->dislocationocation creep */
    830:   param->dislocation.A     = 2.8969e4; /* Pa-s */
    831:   param->dislocation.n     = 3.5;      /* dim'less */
    832:   param->dislocation.Estar = 530e3;    /* J/mol */
    833:   param->dislocation.Vstar = 14e-6;    /* m^3/mol */

    835:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ivisc", &param->ivisc, NULL));
    836:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-visc_cutoff", &param->visc_cutoff, NULL));

    838:   param->output_ivisc = param->ivisc;

    840:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-output_ivisc", &param->output_ivisc, NULL));
    841:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-vstar", &param->dislocation.Vstar, NULL));

    843:   /* output options */
    844:   param->quiet      = PETSC_FALSE;
    845:   param->param_test = PETSC_FALSE;

    847:   PetscCall(PetscOptionsHasName(NULL, NULL, "-quiet", &param->quiet));
    848:   PetscCall(PetscOptionsHasName(NULL, NULL, "-test", &param->param_test));
    849:   PetscCall(PetscOptionsGetString(NULL, NULL, "-file", param->filename, sizeof(param->filename), &param->output_to_file));

    851:   /* advection */
    852:   param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */

    854:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-adv_scheme", &param->adv_scheme, NULL));

    856:   /* misc. flags */
    857:   param->stop_solve    = PETSC_FALSE;
    858:   param->interrupted   = PETSC_FALSE;
    859:   param->kspmon        = PETSC_FALSE;
    860:   param->toggle_kspmon = PETSC_FALSE;

    862:   /* derived parameters for slab angle */
    863:   param->sb = PetscSinReal(param->slab_dip);
    864:   param->cb = PetscCosReal(param->slab_dip);
    865:   param->c  = param->slab_dip * param->sb / (param->slab_dip * param->slab_dip - param->sb * param->sb);
    866:   param->d  = (param->slab_dip * param->cb - param->sb) / (param->slab_dip * param->slab_dip - param->sb * param->sb);

    868:   /* length, velocity and time scale for non-dimensionalization */
    869:   param->L = PetscMin(param->width, param->depth);      /* km */
    870:   param->V = param->slab_velocity / 100.0 / SEC_PER_YR; /* m/sec */

    872:   /* other unit conversions and derived parameters */
    873:   param->scaled_width = param->width / param->L;                   /* dim'less */
    874:   param->scaled_depth = param->depth / param->L;                   /* dim'less */
    875:   param->lid_depth    = param->lid_depth / param->L;               /* dim'less */
    876:   param->fault_depth  = param->fault_depth / param->L;             /* dim'less */
    877:   grid->dx            = grid->dx / param->L;                       /* dim'less */
    878:   grid->dz            = grid->dz / param->L;                       /* dim'less */
    879:   grid->jlid          = (PetscInt)(param->lid_depth / grid->dz);   /* gridcells */
    880:   grid->jfault        = (PetscInt)(param->fault_depth / grid->dz); /* gridcells */
    881:   param->lid_depth    = grid->jlid * grid->dz;                     /* dim'less */
    882:   param->fault_depth  = grid->jfault * grid->dz;                   /* dim'less */
    883:   grid->corner        = grid->jlid + 1;                            /* gridcells */
    884:   param->peclet       = param->V                                   /* m/sec */
    885:                 * param->L * 1000.0                                /* m */
    886:                 / param->kappa;                                    /* m^2/sec */
    887:   param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
    888:   param->skt     = PetscSqrtReal(param->kappa * param->slab_age * SEC_PER_YR);
    889:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-peclet", &param->peclet, NULL));
    890:   PetscFunctionReturn(PETSC_SUCCESS);
    891: }

    893: /*  prints a report of the problem parameters to stdout */
    894: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
    895: {
    896:   char date[30];

    898:   PetscFunctionBeginUser;
    899:   PetscCall(PetscGetDate(date, 30));

    901:   if (!param->quiet) {
    902:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------BEGIN ex30 PARAM REPORT-------------------\n"));
    903:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Domain: \n"));
    904:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Width = %g km,         Depth = %g km\n", (double)param->width, (double)param->depth));
    905:     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));
    906:     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)));

    908:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nGrid: \n"));
    909:     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)));
    910:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  jlid = %3" PetscInt_FMT "              jfault = %3" PetscInt_FMT " \n", grid->jlid, grid->jfault));
    911:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Pe = %g\n", (double)param->peclet));

    913:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nRheology:"));
    914:     if (param->ivisc == VISC_CONST) {
    915:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                 Isoviscous \n"));
    916:       if (param->pv_analytic) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                          Pressure and Velocity prescribed! \n"));
    917:     } else if (param->ivisc == VISC_DIFN) {
    918:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                 Diffusion Creep (T-Dependent Newtonian) \n"));
    919:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                          Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
    920:     } else if (param->ivisc == VISC_DISL) {
    921:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                 Dislocation Creep (T-Dependent Non-Newtonian) \n"));
    922:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                          Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
    923:     } else if (param->ivisc == VISC_FULL) {
    924:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                 Full Rheology \n"));
    925:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                          Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
    926:     } else {
    927:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                 Invalid! \n"));
    928:       PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
    929:     }

    931:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Boundary condition:"));
    932:     if (param->ibound == BC_ANALYTIC) {
    933:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "       Isoviscous Analytic Dirichlet \n"));
    934:     } else if (param->ibound == BC_NOSTRESS) {
    935:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "       Stress-Free (normal & shear stress)\n"));
    936:     } else if (param->ibound == BC_EXPERMNT) {
    937:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "       Experimental boundary condition \n"));
    938:     } else {
    939:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "       Invalid! \n"));
    940:       PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
    941:     }

    943:     if (param->output_to_file) {
    944: #if defined(PETSC_HAVE_MATLAB)
    945:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination:       Mat file \"%s\"\n", param->filename));
    946: #else
    947:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination:       PETSc binary file \"%s\"\n", param->filename));
    948: #endif
    949:     }
    950:     if (param->output_ivisc != param->ivisc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "                          Output viscosity: -ivisc %" PetscInt_FMT "\n", param->output_ivisc));

    952:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------END ex30 PARAM REPORT---------------------\n"));
    953:   }
    954:   if (param->param_test) PetscCall(PetscEnd());
    955:   PetscFunctionReturn(PETSC_SUCCESS);
    956: }

    958: /* ------------------------------------------------------------------- */
    959: /*  generates an initial guess using the analytic solution for isoviscous
    960:     corner flow */
    961: PetscErrorCode Initialize(DM da)
    962: /* ------------------------------------------------------------------- */
    963: {
    964:   AppCtx    *user;
    965:   Parameter *param;
    966:   GridInfo  *grid;
    967:   PetscInt   i, j, is, js, im, jm;
    968:   Field    **x;
    969:   Vec        Xguess;

    971:   PetscFunctionBeginUser;
    972:   /* Get the fine grid */
    973:   PetscCall(DMGetApplicationContext(da, &user));
    974:   Xguess = user->Xguess;
    975:   param  = user->param;
    976:   grid   = user->grid;
    977:   PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
    978:   PetscCall(DMDAVecGetArray(da, Xguess, (void **)&x));

    980:   /* Compute initial guess */
    981:   for (j = js; j < js + jm; j++) {
    982:     for (i = is; i < is + im; i++) {
    983:       if (i < j) x[j][i].u = param->cb;
    984:       else if (j <= grid->jlid) x[j][i].u = 0.0;
    985:       else x[j][i].u = HorizVelocity(i, j, user);

    987:       if (i <= j) x[j][i].w = param->sb;
    988:       else if (j <= grid->jlid) x[j][i].w = 0.0;
    989:       else x[j][i].w = VertVelocity(i, j, user);

    991:       if (i < j || j <= grid->jlid) x[j][i].p = 0.0;
    992:       else x[j][i].p = Pressure(i, j, user);

    994:       x[j][i].T = PetscMin(grid->dz * (j - 0.5), 1.0);
    995:     }
    996:   }

    998:   /* Restore x to Xguess */
    999:   PetscCall(DMDAVecRestoreArray(da, Xguess, (void **)&x));
    1000:   PetscFunctionReturn(PETSC_SUCCESS);
    1001: }

    1003: /*  controls output to a file */
    1004: PetscErrorCode DoOutput(SNES snes, PetscInt its)
    1005: {
    1006:   AppCtx     *user;
    1007:   Parameter  *param;
    1008:   GridInfo   *grid;
    1009:   PetscInt    ivt;
    1010:   PetscMPIInt rank;
    1011:   PetscViewer viewer;
    1012:   Vec         res, pars;
    1013:   MPI_Comm    comm;
    1014:   DM          da;

    1016:   PetscFunctionBeginUser;
    1017:   PetscCall(SNESGetDM(snes, &da));
    1018:   PetscCall(DMGetApplicationContext(da, &user));
    1019:   param = user->param;
    1020:   grid  = user->grid;
    1021:   ivt   = param->ivisc;

    1023:   param->ivisc = param->output_ivisc;

    1025:   /* compute final residual and final viscosity/strain rate fields */
    1026:   PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
    1027:   PetscCall(ViscosityField(da, user->x, user->Xguess));

    1029:   /* get the communicator and the rank of the processor */
    1030:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
    1031:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    1033:   if (param->output_to_file) { /* send output to binary file */
    1034:     PetscCall(VecCreate(comm, &pars));
    1035:     if (rank == 0) { /* on processor 0 */
    1036:       PetscCall(VecSetSizes(pars, 20, PETSC_DETERMINE));
    1037:       PetscCall(VecSetFromOptions(pars));
    1038:       PetscCall(VecSetValue(pars, 0, (PetscScalar)grid->ni, INSERT_VALUES));
    1039:       PetscCall(VecSetValue(pars, 1, (PetscScalar)grid->nj, INSERT_VALUES));
    1040:       PetscCall(VecSetValue(pars, 2, (PetscScalar)grid->dx, INSERT_VALUES));
    1041:       PetscCall(VecSetValue(pars, 3, (PetscScalar)grid->dz, INSERT_VALUES));
    1042:       PetscCall(VecSetValue(pars, 4, (PetscScalar)param->L, INSERT_VALUES));
    1043:       PetscCall(VecSetValue(pars, 5, (PetscScalar)param->V, INSERT_VALUES));
    1044:       /* skipped 6 intentionally */
    1045:       PetscCall(VecSetValue(pars, 7, (PetscScalar)param->slab_dip, INSERT_VALUES));
    1046:       PetscCall(VecSetValue(pars, 8, (PetscScalar)grid->jlid, INSERT_VALUES));
    1047:       PetscCall(VecSetValue(pars, 9, (PetscScalar)param->lid_depth, INSERT_VALUES));
    1048:       PetscCall(VecSetValue(pars, 10, (PetscScalar)grid->jfault, INSERT_VALUES));
    1049:       PetscCall(VecSetValue(pars, 11, (PetscScalar)param->fault_depth, INSERT_VALUES));
    1050:       PetscCall(VecSetValue(pars, 12, (PetscScalar)param->potentialT, INSERT_VALUES));
    1051:       PetscCall(VecSetValue(pars, 13, (PetscScalar)param->ivisc, INSERT_VALUES));
    1052:       PetscCall(VecSetValue(pars, 14, (PetscScalar)param->visc_cutoff, INSERT_VALUES));
    1053:       PetscCall(VecSetValue(pars, 15, (PetscScalar)param->ibound, INSERT_VALUES));
    1054:       PetscCall(VecSetValue(pars, 16, (PetscScalar)its, INSERT_VALUES));
    1055:     } else { /* on some other processor */
    1056:       PetscCall(VecSetSizes(pars, 0, PETSC_DETERMINE));
    1057:       PetscCall(VecSetFromOptions(pars));
    1058:     }
    1059:     PetscCall(VecAssemblyBegin(pars));
    1060:     PetscCall(VecAssemblyEnd(pars));

    1062:     /* create viewer */
    1063: #if defined(PETSC_HAVE_MATLAB)
    1064:     PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
    1065: #else
    1066:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
    1067: #endif

    1069:     /* send vectors to viewer */
    1070:     PetscCall(PetscObjectSetName((PetscObject)res, "res"));
    1071:     PetscCall(VecView(res, viewer));
    1072:     PetscCall(PetscObjectSetName((PetscObject)user->x, "out"));
    1073:     PetscCall(VecView(user->x, viewer));
    1074:     PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "aux"));
    1075:     PetscCall(VecView(user->Xguess, viewer));
    1076:     PetscCall(StressField(da)); /* compute stress fields */
    1077:     PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "str"));
    1078:     PetscCall(VecView(user->Xguess, viewer));
    1079:     PetscCall(PetscObjectSetName((PetscObject)pars, "par"));
    1080:     PetscCall(VecView(pars, viewer));

    1082:     /* destroy viewer and vector */
    1083:     PetscCall(PetscViewerDestroy(&viewer));
    1084:     PetscCall(VecDestroy(&pars));
    1085:   }

    1087:   param->ivisc = ivt;
    1088:   PetscFunctionReturn(PETSC_SUCCESS);
    1089: }

    1091: /* ------------------------------------------------------------------- */
    1092: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
    1093: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
    1094: /* ------------------------------------------------------------------- */
    1095: {
    1096:   AppCtx    *user;
    1097:   Parameter *param;
    1098:   GridInfo  *grid;
    1099:   Vec        localX;
    1100:   Field    **v, **x;
    1101:   PetscReal  eps, /* dx,*/ dz, T, epsC, TC;
    1102:   PetscInt   i, j, is, js, im, jm, ilim, jlim, ivt;

    1104:   PetscFunctionBeginUser;
    1105:   PetscCall(DMGetApplicationContext(da, &user));
    1106:   param        = user->param;
    1107:   grid         = user->grid;
    1108:   ivt          = param->ivisc;
    1109:   param->ivisc = param->output_ivisc;

    1111:   PetscCall(DMGetLocalVector(da, &localX));
    1112:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    1113:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
    1114:   PetscCall(DMDAVecGetArray(da, localX, (void **)&x));
    1115:   PetscCall(DMDAVecGetArray(da, V, (void **)&v));

    1117:   /* Parameters */
    1118:   /* dx = grid->dx; */ dz = grid->dz;

    1120:   ilim = grid->ni - 1;
    1121:   jlim = grid->nj - 1;

    1123:   /* Compute real temperature, strain rate and viscosity */
    1124:   PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
    1125:   for (j = js; j < js + jm; j++) {
    1126:     for (i = is; i < is + im; i++) {
    1127:       T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale));
    1128:       if (i < ilim && j < jlim) {
    1129:         TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale));
    1130:       } else {
    1131:         TC = T;
    1132:       }
    1133:       eps  = PetscRealPart(CalcSecInv(x, i, j, CELL_CENTER, user));
    1134:       epsC = PetscRealPart(CalcSecInv(x, i, j, CELL_CORNER, user));

    1136:       v[j][i].u = eps;
    1137:       v[j][i].w = epsC;
    1138:       v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param);
    1139:       v[j][i].T = Viscosity(TC, epsC, dz * j, param);
    1140:     }
    1141:   }
    1142:   PetscCall(DMDAVecRestoreArray(da, V, (void **)&v));
    1143:   PetscCall(DMDAVecRestoreArray(da, localX, (void **)&x));
    1144:   PetscCall(DMRestoreLocalVector(da, &localX));

    1146:   param->ivisc = ivt;
    1147:   PetscFunctionReturn(PETSC_SUCCESS);
    1148: }

    1150: /* ------------------------------------------------------------------- */
    1151: /* post-processing: compute stress everywhere */
    1152: PetscErrorCode StressField(DM da)
    1153: /* ------------------------------------------------------------------- */
    1154: {
    1155:   AppCtx  *user;
    1156:   PetscInt i, j, is, js, im, jm;
    1157:   Vec      locVec;
    1158:   Field  **x, **y;

    1160:   PetscFunctionBeginUser;
    1161:   PetscCall(DMGetApplicationContext(da, &user));

    1163:   /* Get the fine grid of Xguess and X */
    1164:   PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
    1165:   PetscCall(DMDAVecGetArray(da, user->Xguess, (void **)&x));

    1167:   PetscCall(DMGetLocalVector(da, &locVec));
    1168:   PetscCall(DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec));
    1169:   PetscCall(DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec));
    1170:   PetscCall(DMDAVecGetArray(da, locVec, (void **)&y));

    1172:   /* Compute stress on the corner points */
    1173:   for (j = js; j < js + jm; j++) {
    1174:     for (i = is; i < is + im; i++) {
    1175:       x[j][i].u = ShearStress(y, i, j, CELL_CENTER, user);
    1176:       x[j][i].w = ShearStress(y, i, j, CELL_CORNER, user);
    1177:       x[j][i].p = XNormalStress(y, i, j, CELL_CENTER, user);
    1178:       x[j][i].T = ZNormalStress(y, i, j, CELL_CENTER, user);
    1179:     }
    1180:   }

    1182:   /* Restore the fine grid of Xguess and X */
    1183:   PetscCall(DMDAVecRestoreArray(da, user->Xguess, (void **)&x));
    1184:   PetscCall(DMDAVecRestoreArray(da, locVec, (void **)&y));
    1185:   PetscCall(DMRestoreLocalVector(da, &locVec));
    1186:   PetscFunctionReturn(PETSC_SUCCESS);
    1187: }

    1189: /*=====================================================================
    1190:   UTILITY FUNCTIONS
    1191:   =====================================================================*/

    1193: /* returns the velocity of the subducting slab and handles fault nodes for BC */
    1194: static inline PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
    1195: {
    1196:   Parameter *param = user->param;
    1197:   GridInfo  *grid  = user->grid;

    1199:   if (c == 'U' || c == 'u') {
    1200:     if (i < j - 1) return param->cb;
    1201:     else if (j <= grid->jfault) return 0.0;
    1202:     else return param->cb;

    1204:   } else {
    1205:     if (i < j) return param->sb;
    1206:     else if (j <= grid->jfault) return 0.0;
    1207:     else return param->sb;
    1208:   }
    1209: }

    1211: /*  solution to diffusive half-space cooling model for BC */
    1212: static inline PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
    1213: {
    1214:   Parameter  *param = user->param;
    1215:   PetscScalar z;
    1216:   if (plate == PLATE_LID) z = (j - 0.5) * user->grid->dz;
    1217:   else z = (j - 0.5) * user->grid->dz * param->cb; /* PLATE_SLAB */
    1218: #if defined(PETSC_HAVE_ERF)
    1219:   return (PetscReal)(erf((double)PetscRealPart(z * param->L / 2.0 / param->skt)));
    1220: #else
    1221:   (*PetscErrorPrintf)("erf() not available on this machine\n");
    1222:   MPI_Abort(PETSC_COMM_SELF, 1);
    1223: #endif
    1224: }

    1226: /*=====================================================================
    1227:   INTERACTIVE SIGNAL HANDLING
    1228:   =====================================================================*/

    1230: /* ------------------------------------------------------------------- */
    1231: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it, PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
    1232: /* ------------------------------------------------------------------- */
    1233: {
    1234:   AppCtx    *user  = (AppCtx *)ctx;
    1235:   Parameter *param = user->param;
    1236:   KSP        ksp;

    1238:   PetscFunctionBeginUser;
    1239:   if (param->interrupted) {
    1240:     param->interrupted = PETSC_FALSE;
    1241:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: exiting SNES solve. \n"));
    1242:     *reason = SNES_CONVERGED_FNORM_ABS;
    1243:     PetscFunctionReturn(PETSC_SUCCESS);
    1244:   } else if (param->toggle_kspmon) {
    1245:     param->toggle_kspmon = PETSC_FALSE;

    1247:     PetscCall(SNESGetKSP(snes, &ksp));

    1249:     if (param->kspmon) {
    1250:       PetscCall(KSPMonitorCancel(ksp));

    1252:       param->kspmon = PETSC_FALSE;
    1253:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: deactivating ksp singular value monitor. \n"));
    1254:     } else {
    1255:       PetscViewerAndFormat *vf;
    1256:       PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
    1257:       PetscCall(KSPMonitorSet(ksp, (PetscErrorCode (*)(KSP, PetscInt, PetscReal, void *))KSPMonitorSingularValue, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));

    1259:       param->kspmon = PETSC_TRUE;
    1260:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: activating ksp singular value monitor. \n"));
    1261:     }
    1262:   }
    1263:   PetscCall(SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason, ctx));
    1264:   PetscFunctionReturn(PETSC_SUCCESS);
    1265: }

    1267: /* ------------------------------------------------------------------- */
    1268: #include <signal.h>
    1269: PetscErrorCode InteractiveHandler(int signum, void *ctx)
    1270: /* ------------------------------------------------------------------- */
    1271: {
    1272:   AppCtx    *user  = (AppCtx *)ctx;
    1273:   Parameter *param = user->param;

    1275:   if (signum == SIGILL) {
    1276:     param->toggle_kspmon = PETSC_TRUE;
    1277: #if !defined(PETSC_MISSING_SIGCONT)
    1278:   } else if (signum == SIGCONT) {
    1279:     param->interrupted = PETSC_TRUE;
    1280: #endif
    1281: #if !defined(PETSC_MISSING_SIGURG)
    1282:   } else if (signum == SIGURG) {
    1283:     param->stop_solve = PETSC_TRUE;
    1284: #endif
    1285:   }
    1286:   return PETSC_SUCCESS;
    1287: }

    1289: /*  main call-back function that computes the processor-local piece of the residual */
    1290: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr)
    1291: {
    1292:   AppCtx     *user  = (AppCtx *)ptr;
    1293:   Parameter  *param = user->param;
    1294:   GridInfo   *grid  = user->grid;
    1295:   PetscScalar mag_w, mag_u;
    1296:   PetscInt    i, j, mx, mz, ilim, jlim;
    1297:   PetscInt    is, ie, js, je, ibound; /* ,ivisc */

    1299:   PetscFunctionBeginUser;
    1300:   /* Define global and local grid parameters */
    1301:   mx   = info->mx;
    1302:   mz   = info->my;
    1303:   ilim = mx - 1;
    1304:   jlim = mz - 1;
    1305:   is   = info->xs;
    1306:   ie   = info->xs + info->xm;
    1307:   js   = info->ys;
    1308:   je   = info->ys + info->ym;

    1310:   /* Define geometric and numeric parameters */
    1311:   /* ivisc = param->ivisc; */ ibound = param->ibound;

    1313:   for (j = js; j < je; j++) {
    1314:     for (i = is; i < ie; i++) {
    1315:       /************* X-MOMENTUM/VELOCITY *************/
    1316:       if (i < j) f[j][i].u = x[j][i].u - SlabVel('U', i, j, user);
    1317:       else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
    1318:         /* in the lithospheric lid */
    1319:         f[j][i].u = x[j][i].u - 0.0;
    1320:       } else if (i == ilim) {
    1321:         /* on the right side boundary */
    1322:         if (ibound == BC_ANALYTIC) {
    1323:           f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
    1324:         } else {
    1325:           f[j][i].u = XNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
    1326:         }

    1328:       } else if (j == jlim) {
    1329:         /* on the bottom boundary */
    1330:         if (ibound == BC_ANALYTIC) {
    1331:           f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
    1332:         } else if (ibound == BC_NOSTRESS) {
    1333:           f[j][i].u = XMomentumResidual(x, i, j, user);
    1334:         } else {
    1335:           /* experimental boundary condition */
    1336:         }

    1338:       } else {
    1339:         /* in the mantle wedge */
    1340:         f[j][i].u = XMomentumResidual(x, i, j, user);
    1341:       }

    1343:       /************* Z-MOMENTUM/VELOCITY *************/
    1344:       if (i <= j) {
    1345:         f[j][i].w = x[j][i].w - SlabVel('W', i, j, user);

    1347:       } else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
    1348:         /* in the lithospheric lid */
    1349:         f[j][i].w = x[j][i].w - 0.0;

    1351:       } else if (j == jlim) {
    1352:         /* on the bottom boundary */
    1353:         if (ibound == BC_ANALYTIC) {
    1354:           f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
    1355:         } else {
    1356:           f[j][i].w = ZNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
    1357:         }

    1359:       } else if (i == ilim) {
    1360:         /* on the right side boundary */
    1361:         if (ibound == BC_ANALYTIC) {
    1362:           f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
    1363:         } else if (ibound == BC_NOSTRESS) {
    1364:           f[j][i].w = ZMomentumResidual(x, i, j, user);
    1365:         } else {
    1366:           /* experimental boundary condition */
    1367:         }

    1369:       } else {
    1370:         /* in the mantle wedge */
    1371:         f[j][i].w = ZMomentumResidual(x, i, j, user);
    1372:       }

    1374:       /************* CONTINUITY/PRESSURE *************/
    1375:       if (i < j || j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
    1376:         /* in the lid or slab */
    1377:         f[j][i].p = x[j][i].p;

    1379:       } else if ((i == ilim || j == jlim) && ibound == BC_ANALYTIC) {
    1380:         /* on an analytic boundary */
    1381:         f[j][i].p = x[j][i].p - Pressure(i, j, user);

    1383:       } else {
    1384:         /* in the mantle wedge */
    1385:         f[j][i].p = ContinuityResidual(x, i, j, user);
    1386:       }

    1388:       /************* TEMPERATURE *************/
    1389:       if (j == 0) {
    1390:         /* on the surface */
    1391:         f[j][i].T = x[j][i].T + x[j + 1][i].T + PetscMax(PetscRealPart(x[j][i].T), 0.0);

    1393:       } else if (i == 0) {
    1394:         /* slab inflow boundary */
    1395:         f[j][i].T = x[j][i].T - PlateModel(j, PLATE_SLAB, user);

    1397:       } else if (i == ilim) {
    1398:         /* right side boundary */
    1399:         mag_u     = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j][i - 1].u) / param->cb, 1.0), 0.0), 5);
    1400:         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);

    1402:       } else if (j == jlim) {
    1403:         /* bottom boundary */
    1404:         mag_w     = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j - 1][i].w) / param->sb, 1.0), 0.0), 5);
    1405:         f[j][i].T = x[j][i].T - mag_w * x[j - 1][i - 1].T - (1.0 - mag_w);

    1407:       } else {
    1408:         /* in the mantle wedge */
    1409:         f[j][i].T = EnergyResidual(x, i, j, user);
    1410:       }
    1411:     }
    1412:   }
    1413:   PetscFunctionReturn(PETSC_SUCCESS);
    1414: }

    1416: /*TEST

    1418:    build:
    1419:       requires: !complex erf

    1421:    test:
    1422:       args: -ni 18 -fp_trap 0
    1423:       filter: grep -v Destination
    1424:       requires: !single

    1426: TEST*/