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, (char *)0, 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));

1001:   PetscFunctionReturn(PETSC_SUCCESS);
1002: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1417: /*TEST

1419:    build:
1420:       requires: !complex erf

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

1427: TEST*/