Actual source code: ex30.c

petsc-master 2016-12-08
Report Typos and Errors
  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";



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

The problem is modeled by the partial differential equation system

\begin{eqnarray}
-\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\
\nabla \cdot v & = & 0 \\
dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\
\end{eqnarray}

\begin{eqnarray}
\eta(T,Eps\_dot) & = & \hbox{constant } \hbox{if ivisc} ==0 \\
& = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\
& = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\
& = & \hbox{mantle viscosity (difn and disl) } \hbox{if ivisc} ==3
\end{eqnarray}

which is uniformly discretized on a staggered mesh:
-------$w_{ij}$------
$u_{i-1j}$ $P,T_{ij}$ $u_{ij}$
------$w_{ij-1}$-----

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

 59: #define VISC_CONST   0
 60: #define VISC_DIFN    1
 61: #define VISC_DISL    2
 62: #define VISC_FULL    3
 63: #define CELL_CENTER  0
 64: #define CELL_CORNER  1
 65: #define BC_ANALYTIC  0
 66: #define BC_NOSTRESS  1
 67: #define BC_EXPERMNT  2
 68: #define ADVECT_FV    0
 69: #define ADVECT_FROMM 1
 70: #define PLATE_SLAB   0
 71: #define PLATE_LID    1
 72: #define EPS_ZERO     0.00000001

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

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

 82: typedef struct { /* physical and miscelaneous parameters */
 83:   PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT;
 84:   PetscReal slab_dip, slab_age, slab_velocity, kappa, z_scale;
 85:   PetscReal c, d, sb, cb, skt, visc_cutoff, lid_age, eta0, continuation;
 86:   PetscReal L, V, lid_depth, fault_depth;
 87:   ViscParam diffusion, dislocation;
 88:   PetscInt  ivisc, adv_scheme, ibound, output_ivisc;
 89:   PetscBool quiet, param_test, output_to_file, pv_analytic;
 90:   PetscBool interrupted, stop_solve, toggle_kspmon, kspmon;
 91:   char      filename[PETSC_MAX_PATH_LEN];
 92: } Parameter;

 94: typedef struct { /* grid parameters */
 95:   DMBoundaryType   bx,by;
 96:   DMDAStencilType  stencil;
 97:   PetscInt         corner,ni,nj,jlid,jfault,inose;
 98:   PetscInt         dof,stencil_width,mglevels;
 99:   PetscReal        dx,dz;
100: } GridInfo;

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

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

111: /* Main routines */
112: extern PetscErrorCode SetParams(Parameter*, GridInfo*);
113: extern PetscErrorCode ReportParams(Parameter*, GridInfo*);
114: extern PetscErrorCode Initialize(DM);
115: extern PetscErrorCode UpdateSolution(SNES,AppCtx*, PetscInt*);
116: extern PetscErrorCode DoOutput(SNES,PetscInt);

118: /* Post-processing & misc */
119: extern PetscErrorCode ViscosityField(DM,Vec,Vec);
120: extern PetscErrorCode StressField(DM);
121: extern PetscErrorCode SNESConverged_Interactive(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason*, void*);
122: extern PetscErrorCode InteractiveHandler(int, void*);
123: extern PetscBool  OptionsHasName(const char pre[],const char name[]);

125: /*-----------------------------------------------------------------------*/
128: int main(int argc,char **argv)
129: /*-----------------------------------------------------------------------*/
130: {
131:   SNES           snes;
132:   AppCtx         *user;               /* user-defined work context */
133:   Parameter      param;
134:   GridInfo       grid;
135:   PetscInt       nits;
137:   MPI_Comm       comm;
138:   DM             da;

140:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
141:   PetscOptionsSetValue(NULL,"-file","ex30_output");
142:   PetscOptionsSetValue(NULL,"-snes_monitor_short",NULL);
143:   PetscOptionsSetValue(NULL,"-snes_max_it","20");
144:   PetscOptionsSetValue(NULL,"-ksp_max_it","1500");
145:   PetscOptionsSetValue(NULL,"-ksp_gmres_restart","300");
146:   PetscOptionsInsert(NULL,&argc,&argv,NULL);

148:   comm = PETSC_COMM_WORLD;

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Set up the problem parameters.
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   SetParams(&param,&grid);
154:   ReportParams(&param,&grid);

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158:   SNESCreate(comm,&snes);
159:   DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
160:   DMSetFromOptions(da);
161:   DMSetUp(da);
162:   SNESSetDM(snes,da);
163:   DMDASetFieldName(da,0,"x-velocity");
164:   DMDASetFieldName(da,1,"y-velocity");
165:   DMDASetFieldName(da,2,"pressure");
166:   DMDASetFieldName(da,3,"temperature");


169:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170:      Create user context, set problem data, create vector data structures.
171:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172:   PetscNew(&user);
173:   user->param = &param;
174:   user->grid  = &grid;
175:   DMSetApplicationContext(da,user);
176:   DMCreateGlobalVector(da,&(user->Xguess));


179:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180:      Set up the SNES solver with callback functions.
181:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,(void*)user);
183:   SNESSetFromOptions(snes);


186:   SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,NULL);
187:   PetscPushSignalHandler(InteractiveHandler,(void*)user);

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:      Initialize and solve the nonlinear system
191:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192:   Initialize(da);
193:   UpdateSolution(snes,user,&nits);

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Output variables.
197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198:   DoOutput(snes,nits);

200:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201:      Free work space.
202:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203:   VecDestroy(&user->Xguess);
204:   VecDestroy(&user->x);
205:   PetscFree(user);
206:   SNESDestroy(&snes);
207:   DMDestroy(&da);
208:   PetscPopSignalHandler();
209:   PetscFinalize();
210:   return ierr;
211: }

213: /*=====================================================================
214:   PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
215:   =====================================================================*/

217: /*---------------------------------------------------------------------*/
220: /*  manages solve: adaptive continuation method  */
221: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
222: {
223:   KSP                 ksp;
224:   PC                  pc;
225:   SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
226:   Parameter           *param   = user->param;
227:   PetscReal           cont_incr=0.3;
228:   PetscInt            its;
229:   PetscErrorCode      ierr;
230:   PetscBool           q = PETSC_FALSE;
231:   DM                  dm;

234:   SNESGetDM(snes,&dm);
235:   DMCreateGlobalVector(dm,&user->x);
236:   SNESGetKSP(snes,&ksp);
237:   KSPGetPC(ksp, &pc);
238:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);

240:   *nits=0;

242:   /* Isoviscous solve */
243:   if (param->ivisc == VISC_CONST && !param->stop_solve) {
244:     param->ivisc = VISC_CONST;

246:     SNESSolve(snes,0,user->x);
247:     SNESGetConvergedReason(snes,&reason);
248:     SNESGetIterationNumber(snes,&its);
249:     *nits += its;
250:     VecCopy(user->x,user->Xguess);
251:     if (param->stop_solve) goto done;
252:   }

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

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

262:       /* solve the non-linear system */
263:       VecCopy(user->Xguess,user->x);
264:       SNESSolve(snes,0,user->x);
265:       SNESGetConvergedReason(snes,&reason);
266:       SNESGetIterationNumber(snes,&its);
267:       *nits += its;
268:       if (!q) PetscPrintf(PETSC_COMM_WORLD," SNES iterations: %D, Cumulative: %D\n", its, *nits);
269:       if (param->stop_solve) goto done;

271:       if (reason<0) {
272:         /* NOT converged */
273:         cont_incr = -fabs(cont_incr)/2.0;
274:         if (fabs(cont_incr)<0.01) goto done;

276:       } else {
277:         /* converged */
278:         VecCopy(user->x,user->Xguess);
279:         if (param->continuation >= 1.0) goto done;
280:         if (its<=3)      cont_incr = 0.30001;
281:         else if (its<=8) cont_incr = 0.15001;
282:         else             cont_incr = 0.10001;

284:         if (param->continuation+cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
285:       } /* endif reason<0 */
286:     }
287:   }
288: done:
289:   if (param->stop_solve && !q) {PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");}
290:   if (reason<0 && !q) {PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");}
291:   return(0);
292: }


295: /*=====================================================================
296:   PHYSICS FUNCTIONS (compute the discrete residual)
297:   =====================================================================*/


300: /*---------------------------------------------------------------------*/
303: PETSC_STATIC_INLINE PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
304: /*---------------------------------------------------------------------*/
305: {
306:   return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
307: }

309: /*---------------------------------------------------------------------*/
312: PETSC_STATIC_INLINE PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
313: /*---------------------------------------------------------------------*/
314: {
315:   return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
316: }

318: /*---------------------------------------------------------------------*/
321: PETSC_STATIC_INLINE PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
322: /*---------------------------------------------------------------------*/
323: {
324:   return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
325: }

327: /*---------------------------------------------------------------------*/
330: PETSC_STATIC_INLINE PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
331: /*---------------------------------------------------------------------*/
332: {
333:   return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
334: }

336: /*---------------------------------------------------------------------*/
339: /*  isoviscous analytic solution for IC */
340: PETSC_STATIC_INLINE PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
341: /*---------------------------------------------------------------------*/
342: {
343:   Parameter   *param = user->param;
344:   GridInfo    *grid  = user->grid;
345:   PetscScalar st, ct, th, c=param->c, d=param->d;
346:   PetscReal   x, z,r;

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

356: /*---------------------------------------------------------------------*/
359: /*  isoviscous analytic solution for IC */
360: PETSC_STATIC_INLINE PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
361: /*---------------------------------------------------------------------*/
362: {
363:   Parameter   *param = user->param;
364:   GridInfo    *grid  = user->grid;
365:   PetscScalar st, ct, th, c=param->c, d=param->d;
366:   PetscReal   x, z, r;

368:   x = (i - grid->jlid - 0.5)*grid->dx;  z = (j - grid->jlid)*grid->dz;
369:   r = PetscSqrtReal(x*x+z*z); st = z/r;  ct = x/r;  th = atan(z/x);
370:   return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
371: }

373: /*---------------------------------------------------------------------*/
376: /*  isoviscous analytic solution for IC */
377: PETSC_STATIC_INLINE PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
378: /*---------------------------------------------------------------------*/
379: {
380:   Parameter   *param = user->param;
381:   GridInfo    *grid  = user->grid;
382:   PetscScalar x, z, r, st, ct, c=param->c, d=param->d;

384:   x = (i - grid->jlid - 0.5)*grid->dx;  z = (j - grid->jlid - 0.5)*grid->dz;
385:   r = PetscSqrtReal(x*x+z*z);  st = z/r;  ct = x/r;
386:   return (-2.0*(c*ct-d*st)/r);
387: }

391: /*  computes the second invariant of the strain rate tensor */
392: PETSC_STATIC_INLINE PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
393: /*---------------------------------------------------------------------*/
394: {
395:   Parameter   *param = user->param;
396:   GridInfo    *grid  = user->grid;
397:   PetscInt    ilim   =grid->ni-1, jlim=grid->nj-1;
398:   PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
399:   PetscScalar eps11, eps12, eps22;

401:   if (i<j) return EPS_ZERO;
402:   if (i==ilim) i--;
403:   if (j==jlim) j--;

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

408:     uE = x[j][i].u; uW = x[j][i-1].u;
409:     wN = x[j][i].w; wS = x[j-1][i].w;
410:     wE = WInterp(x,i,j-1);
411:     if (i==j) {
412:       uN = param->cb; wW = param->sb;
413:     } else {
414:       uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1);
415:     }

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

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

423:     uN = x[j+1][i].u;  uS = x[j][i].u;
424:     wE = x[j][i+1].w;  wW = x[j][i].w;
425:     if (i==j) {
426:       wN = param->sb;
427:       uW = param->cb;
428:     } else {
429:       wN = WInterp(x,i,j);
430:       uW = UInterp(x,i-1,j);
431:     }

433:     if (j==grid->jlid) {
434:       uE = 0.0;  uW = 0.0;
435:       uS = -uN;
436:       wS = -wN;
437:     } else {
438:       uE = UInterp(x,i,j);
439:       wS = WInterp(x,i,j-1);
440:     }
441:   }

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

446:   return PetscSqrtReal(0.5*(eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22));
447: }

449: /*---------------------------------------------------------------------*/
452: /*  computes the shear viscosity */
453: PETSC_STATIC_INLINE PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
454: /*---------------------------------------------------------------------*/
455: {
456:   PetscReal   result   =0.0;
457:   ViscParam   difn     =param->diffusion, disl=param->dislocation;
458:   PetscInt    iVisc    =param->ivisc;
459:   PetscScalar eps_scale=param->V/(param->L*1000.0);
460:   PetscScalar strain_power, v1, v2, P;
461:   PetscScalar rho_g = 32340.0, R=8.3144;

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

465:   if (iVisc==VISC_CONST) {
466:     /* constant viscosity */
467:     return 1.0;
468:   } else if (iVisc==VISC_DIFN) {
469:     /* diffusion creep rheology */
470:     result = PetscRealPart((difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0));
471:   } else if (iVisc==VISC_DISL) {
472:     /* dislocation creep rheology */
473:     strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);

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

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

483:     result = PetscRealPart(1.0/(1.0/v1 + 1.0/v2));
484:   }

486:   /* max viscosity is param->eta0 */
487:   result = PetscMin(result, 1.0);
488:   /* min viscosity is param->visc_cutoff */
489:   result = PetscMax(result, param->visc_cutoff);
490:   /* continuation method */
491:   result = PetscPowReal(result,param->continuation);
492:   return result;
493: }

495: /*---------------------------------------------------------------------*/
498: /*  computes the residual of the x-component of eqn (1) above */
499: PETSC_STATIC_INLINE PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
500: /*---------------------------------------------------------------------*/
501: {
502:   Parameter   *param=user->param;
503:   GridInfo    *grid =user->grid;
504:   PetscScalar dx    = grid->dx, dz=grid->dz;
505:   PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
506:   PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
507:   PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
508:   PetscInt    jlim = grid->nj-1;

510:   z_scale = param->z_scale;

512:   if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
513:     TS = param->potentialT * TInterp(x,i,j-1) * PetscExpScalar((j-1.0)*dz*z_scale);
514:     if (j==jlim) TN = TS;
515:     else         TN = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
516:     TW = param->potentialT * x[j][i].T        * PetscExpScalar((j-0.5)*dz*z_scale);
517:     TE = param->potentialT * x[j][i+1].T      * PetscExpScalar((j-0.5)*dz*z_scale);
518:     if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
519:       epsN = CalcSecInv(x,i,j,  CELL_CORNER,user);
520:       epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
521:       epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
522:       epsW = CalcSecInv(x,i,j,  CELL_CENTER,user);
523:     }
524:   }
525:   etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
526:   etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
527:   etaW = Viscosity(TW,epsW,dz*j,param);
528:   etaE = Viscosity(TE,epsE,dz*j,param);

530:   dPdx = (x[j][i+1].p - x[j][i].p)/dx;
531:   if (j==jlim) dudzN = etaN * (x[j][i].w   - x[j][i+1].w)/dx;
532:   else         dudzN = etaN * (x[j+1][i].u - x[j][i].u)  /dz;
533:   dudzS = etaS * (x[j][i].u    - x[j-1][i].u)/dz;
534:   dudxE = etaE * (x[j][i+1].u  - x[j][i].u)  /dx;
535:   dudxW = etaW * (x[j][i].u    - x[j][i-1].u)/dx;

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

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

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

548:   return residual;
549: }

551: /*---------------------------------------------------------------------*/
554: /*  computes the residual of the z-component of eqn (1) above */
555: PETSC_STATIC_INLINE PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
556: /*---------------------------------------------------------------------*/
557: {
558:   Parameter   *param=user->param;
559:   GridInfo    *grid =user->grid;
560:   PetscScalar dx    = grid->dx, dz=grid->dz;
561:   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;
562:   PetscScalar TE    =0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
563:   PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
564:   PetscInt    ilim = grid->ni-1;

566:   /* geometric and other parameters */
567:   z_scale = param->z_scale;

569:   /* viscosity */
570:   if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
571:     TN = param->potentialT * x[j+1][i].T      * PetscExpScalar((j+0.5)*dz*z_scale);
572:     TS = param->potentialT * x[j][i].T        * PetscExpScalar((j-0.5)*dz*z_scale);
573:     TW = param->potentialT * TInterp(x,i-1,j) * PetscExpScalar(j*dz*z_scale);
574:     if (i==ilim) TE = TW;
575:     else         TE = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
576:     if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
577:       epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
578:       epsS = CalcSecInv(x,i,j,  CELL_CENTER,user);
579:       epsE = CalcSecInv(x,i,j,  CELL_CORNER,user);
580:       epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
581:     }
582:   }
583:   etaN = Viscosity(TN,epsN,dz*(j+1.0),param);
584:   etaS = Viscosity(TS,epsS,dz*(j+0.0),param);
585:   etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
586:   etaE = Viscosity(TE,epsE,dz*(j+0.5),param);

588:   dPdz  = (x[j+1][i].p - x[j][i].p)/dz;
589:   dwdzN = etaN * (x[j+1][i].w - x[j][i].w)/dz;
590:   dwdzS = etaS * (x[j][i].w - x[j-1][i].w)/dz;
591:   if (i==ilim) dwdxE = etaE * (x[j][i].u   - x[j+1][i].u)/dz;
592:   else         dwdxE = etaE * (x[j][i+1].w - x[j][i].w)  /dx;
593:   dwdxW = 2.0*etaW * (x[j][i].w - x[j][i-1].w)/dx;

595:   /* Z-MOMENTUM */
596:   residual = -dPdz                 /* constant viscosity terms */
597:              +(dwdzN - dwdzS)/dz
598:              +(dwdxE - dwdxW)/dx;

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

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

607:   return residual;
608: }

610: /*---------------------------------------------------------------------*/
613: /*  computes the residual of eqn (2) above */
614: PETSC_STATIC_INLINE PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
615: /*---------------------------------------------------------------------*/
616: {
617:   GridInfo    *grid =user->grid;
618:   PetscScalar uE,uW,wN,wS,dudx,dwdz;

620:   uW = x[j][i-1].u; uE = x[j][i].u; dudx = (uE - uW)/grid->dx;
621:   wS = x[j-1][i].w; wN = x[j][i].w; dwdz = (wN - wS)/grid->dz;

623:   return dudx + dwdz;
624: }

626: /*---------------------------------------------------------------------*/
629: /*  computes the residual of eqn (3) above */
630: PETSC_STATIC_INLINE PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
631: /*---------------------------------------------------------------------*/
632: {
633:   Parameter   *param=user->param;
634:   GridInfo    *grid =user->grid;
635:   PetscScalar dx    = grid->dx, dz=grid->dz;
636:   PetscInt    ilim  =grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
637:   PetscScalar TE, TN, TS, TW, residual;
638:   PetscScalar uE,uW,wN,wS;
639:   PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;

641:   dTdzN = (x[j+1][i].T - x[j][i].T)  /dz;
642:   dTdzS = (x[j][i].T   - x[j-1][i].T)/dz;
643:   dTdxE = (x[j][i+1].T - x[j][i].T)  /dx;
644:   dTdxW = (x[j][i].T   - x[j][i-1].T)/dx;

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

649:   if (j<=jlid && i>=j) {
650:     /* don't advect in the lid */
651:     return residual;
652:   } else if (i<j) {
653:     /* beneath the slab sfc */
654:     uW = uE = param->cb;
655:     wS = wN = param->sb;
656:   } else {
657:     /* advect in the slab and wedge */
658:     uW = x[j][i-1].u; uE = x[j][i].u;
659:     wS = x[j-1][i].w; wN = x[j][i].w;
660:   }

662:   if (param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1) {
663:     /* finite volume advection */
664:     TS = (x[j][i].T + x[j-1][i].T)/2.0;
665:     TN = (x[j][i].T + x[j+1][i].T)/2.0;
666:     TE = (x[j][i].T + x[j][i+1].T)/2.0;
667:     TW = (x[j][i].T + x[j][i-1].T)/2.0;
668:     fN = wN*TN*dx; fS = wS*TS*dx;
669:     fE = uE*TE*dz; fW = uW*TW*dz;

671:   } else {
672:     /* Fromm advection scheme */
673:     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
674:               - 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;
675:     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
676:               - 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;
677:     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
678:               - 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;
679:     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
680:               - 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;
681:   }

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

685:   return residual;
686: }

688: /*---------------------------------------------------------------------*/
691: /*  computes the shear stress---used on the boundaries */
692: PETSC_STATIC_INLINE PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
693: /*---------------------------------------------------------------------*/
694: {
695:   Parameter   *param=user->param;
696:   GridInfo    *grid =user->grid;
697:   PetscInt    ilim  =grid->ni-1, jlim=grid->nj-1;
698:   PetscScalar uN, uS, wE, wW;

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

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

704:     wE = WInterp(x,i,j-1);
705:     if (i==j) {
706:       wW = param->sb;
707:       uN = param->cb;
708:     } else {
709:       wW = WInterp(x,i-1,j-1);
710:       uN = UInterp(x,i-1,j);
711:     }
712:     if (j==grid->jlid+1) uS = 0.0;
713:     else                 uS = UInterp(x,i-1,j-1);

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

717:     uN = x[j+1][i].u;         uS = x[j][i].u;
718:     wW = x[j][i].w;           wE = x[j][i+1].w;

720:   }

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

725: /*---------------------------------------------------------------------*/
728: /*  computes the normal stress---used on the boundaries */
729: PETSC_STATIC_INLINE PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
730: /*---------------------------------------------------------------------*/
731: {
732:   Parameter   *param=user->param;
733:   GridInfo    *grid =user->grid;
734:   PetscScalar dx    = grid->dx, dz=grid->dz;
735:   PetscInt    ilim  =grid->ni-1, jlim=grid->nj-1, ivisc;
736:   PetscScalar epsC  =0.0, etaC, TC, uE, uW, pC, z_scale;
737:   if (i<j || j<=grid->jlid) return EPS_ZERO;

739:   ivisc=param->ivisc;  z_scale = param->z_scale;

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

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

747:     uW = x[j][i-1].u;   uE = x[j][i].u;
748:     pC = x[j][i].p;

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

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

757:     if (i==j) uW = param->sb;
758:     else      uW = UInterp(x,i-1,j);
759:     uE = UInterp(x,i,j); pC = PInterp(x,i,j);
760:   }

762:   return 2.0*etaC*(uE-uW)/dx - pC;
763: }

765: /*---------------------------------------------------------------------*/
768: /*  computes the normal stress---used on the boundaries */
769: PETSC_STATIC_INLINE PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
770: /*---------------------------------------------------------------------*/
771: {
772:   Parameter   *param=user->param;
773:   GridInfo    *grid =user->grid;
774:   PetscScalar dz    =grid->dz;
775:   PetscInt    ilim  =grid->ni-1, jlim=grid->nj-1, ivisc;
776:   PetscScalar epsC  =0.0, etaC, TC;
777:   PetscScalar pC, wN, wS, z_scale;
778:   if (i<j || j<=grid->jlid) return EPS_ZERO;

780:   ivisc=param->ivisc;  z_scale = param->z_scale;

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

784:     TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
785:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
786:     etaC = Viscosity(TC,epsC,dz*j,param);
787:     wN   = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;

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

792:     TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
793:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
794:     etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
795:     if (i==j) wN = param->sb;
796:     else      wN = WInterp(x,i,j);
797:     wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
798:   }

800:   return 2.0*etaC*(wN-wS)/dz - pC;
801: }

803: /*---------------------------------------------------------------------*/

805: /*=====================================================================
806:   INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
807:   =====================================================================*/

809: /*---------------------------------------------------------------------*/
812: /* initializes the problem parameters and checks for
813:    command line changes */
814: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
815: /*---------------------------------------------------------------------*/
816: {
817:   PetscErrorCode ierr, ierr_out=0;
818:   PetscReal      SEC_PER_YR                    = 3600.00*24.00*365.2500;
819:   PetscReal      alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;

821:   /* domain geometry */
822:   param->slab_dip    = 45.0;
823:   param->width       = 320.0;                                              /* km */
824:   param->depth       = 300.0;                                              /* km */
825:   param->lid_depth   = 35.0;                                               /* km */
826:   param->fault_depth = 35.0;                                               /* km */

828:   PetscOptionsGetReal(NULL,NULL,"-slab_dip",&(param->slab_dip),NULL);
829:   PetscOptionsGetReal(NULL,NULL,"-width",&(param->width),NULL);
830:   PetscOptionsGetReal(NULL,NULL,"-depth",&(param->depth),NULL);
831:   PetscOptionsGetReal(NULL,NULL,"-lid_depth",&(param->lid_depth),NULL);
832:   PetscOptionsGetReal(NULL,NULL,"-fault_depth",&(param->fault_depth),NULL);

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

836:   /* grid information */
837:   PetscOptionsGetInt(NULL,NULL, "-jfault",&(grid->jfault),NULL);
838:   grid->ni = 82;
839:   PetscOptionsGetInt(NULL,NULL, "-ni",&(grid->ni),NULL);

841:   grid->dx            = param->width/((PetscReal)(grid->ni-2));               /* km */
842:   grid->dz            = grid->dx*tan(param->slab_dip);                     /* km */
843:   grid->nj            = (PetscInt)(param->depth/grid->dz + 3.0);         /* gridpoints*/
844:   param->depth        = grid->dz*(grid->nj-2);                             /* km */
845:   grid->inose         = 0;                                          /* gridpoints*/
846:   PetscOptionsGetInt(NULL,NULL,"-inose",&(grid->inose),NULL);
847:   grid->bx            = DM_BOUNDARY_NONE;
848:   grid->by            = DM_BOUNDARY_NONE;
849:   grid->stencil       = DMDA_STENCIL_BOX;
850:   grid->dof           = 4;
851:   grid->stencil_width = 2;
852:   grid->mglevels      = 1;

854:   /* boundary conditions */
855:   param->pv_analytic = PETSC_FALSE;
856:   param->ibound      = BC_NOSTRESS;
857:   PetscOptionsGetInt(NULL,NULL,"-ibound",&(param->ibound),NULL);

859:   /* physical constants */
860:   param->slab_velocity = 5.0;               /* cm/yr */
861:   param->slab_age      = 50.0;              /* Ma */
862:   param->lid_age       = 50.0;              /* Ma */
863:   param->kappa         = 0.7272e-6;         /* m^2/sec */
864:   param->potentialT    = 1300.0;            /* degrees C */

866:   PetscOptionsGetReal(NULL,NULL,"-slab_velocity",&(param->slab_velocity),NULL);
867:   PetscOptionsGetReal(NULL,NULL,"-slab_age",&(param->slab_age),NULL);
868:   PetscOptionsGetReal(NULL,NULL,"-lid_age",&(param->lid_age),NULL);
869:   PetscOptionsGetReal(NULL,NULL,"-kappa",&(param->kappa),NULL);
870:   PetscOptionsGetReal(NULL,NULL,"-potentialT",&(param->potentialT),NULL);

872:   /* viscosity */
873:   param->ivisc        = 3;                  /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
874:   param->eta0         = 1e24;               /* Pa-s */
875:   param->visc_cutoff  = 0.0;                /* factor of eta_0 */
876:   param->continuation = 1.0;

878:   /* constants for diffusion creep */
879:   param->diffusion.A     = 1.8e7;             /* Pa-s */
880:   param->diffusion.n     = 1.0;               /* dim'less */
881:   param->diffusion.Estar = 375e3;             /* J/mol */
882:   param->diffusion.Vstar = 5e-6;              /* m^3/mol */

884:   /* constants for param->dislocationocation creep */
885:   param->dislocation.A     = 2.8969e4;        /* Pa-s */
886:   param->dislocation.n     = 3.5;             /* dim'less */
887:   param->dislocation.Estar = 530e3;           /* J/mol */
888:   param->dislocation.Vstar = 14e-6;           /* m^3/mol */

890:   PetscOptionsGetInt(NULL,NULL, "-ivisc",&(param->ivisc),NULL);
891:   PetscOptionsGetReal(NULL,NULL,"-visc_cutoff",&(param->visc_cutoff),NULL);

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

895:   PetscOptionsGetInt(NULL,NULL,"-output_ivisc",&(param->output_ivisc),NULL);
896:   PetscOptionsGetReal(NULL,NULL,"-vstar",&(param->dislocation.Vstar),NULL);

898:   /* output options */
899:   param->quiet      = PETSC_FALSE;
900:   param->param_test = PETSC_FALSE;

902:   PetscOptionsHasName(NULL,NULL,"-quiet",&(param->quiet));
903:   PetscOptionsHasName(NULL,NULL,"-test",&(param->param_test));
904:   PetscOptionsGetString(NULL,NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));

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

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

911:   /* misc. flags */
912:   param->stop_solve    = PETSC_FALSE;
913:   param->interrupted   = PETSC_FALSE;
914:   param->kspmon        = PETSC_FALSE;
915:   param->toggle_kspmon = PETSC_FALSE;

917:   /* derived parameters for slab angle */
918:   param->sb = PetscSinReal(param->slab_dip);
919:   param->cb = PetscCosReal(param->slab_dip);
920:   param->c  =  param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
921:   param->d  = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);

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

927:   /* other unit conversions and derived parameters */
928:   param->scaled_width = param->width/param->L;                  /* dim'less */
929:   param->scaled_depth = param->depth/param->L;                  /* dim'less */
930:   param->lid_depth    = param->lid_depth/param->L;              /* dim'less */
931:   param->fault_depth  = param->fault_depth/param->L;            /* dim'less */
932:   grid->dx            = grid->dx/param->L;                      /* dim'less */
933:   grid->dz            = grid->dz/param->L;                      /* dim'less */
934:   grid->jlid          = (PetscInt)(param->lid_depth/grid->dz);       /* gridcells */
935:   grid->jfault        = (PetscInt)(param->fault_depth/grid->dz);     /* gridcells */
936:   param->lid_depth    = grid->jlid*grid->dz;                    /* dim'less */
937:   param->fault_depth  = grid->jfault*grid->dz;                  /* dim'less */
938:   grid->corner        = grid->jlid+1;                           /* gridcells */
939:   param->peclet       = param->V                                /* m/sec */
940:                         * param->L*1000.0                       /* m */
941:                         / param->kappa;                         /* m^2/sec */
942:   param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
943:   param->skt     = PetscSqrtReal(param->kappa*param->slab_age*SEC_PER_YR);
944:   PetscOptionsGetReal(NULL,NULL,"-peclet",&(param->peclet),NULL);

946:   return ierr_out;
947: }

949: /*---------------------------------------------------------------------*/
952: /*  prints a report of the problem parameters to stdout */
953: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
954: /*---------------------------------------------------------------------*/
955: {
956:   PetscErrorCode ierr, ierr_out=0;
957:   char           date[30];

959:   PetscGetDate(date,30);

961:   if (!(param->quiet)) {
962:     PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
963:     /* PetscPrintf(PETSC_COMM_WORLD,"                   %s\n",&(date[0]));*/

965:     PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
966:     PetscPrintf(PETSC_COMM_WORLD,"  Width = %g km,         Depth = %g km\n",(double)param->width,(double)param->depth);
967:     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);
968:     PetscPrintf(PETSC_COMM_WORLD,"  Lid depth = %5.2f km,   Fault depth = %5.2f km\n",param->lid_depth*param->L,param->fault_depth*param->L);

970:     PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
971:     PetscPrintf(PETSC_COMM_WORLD,"  [ni,nj] = %D, %D       [dx,dz] = %g, %g km\n",grid->ni,grid->nj,(double)grid->dx*param->L,(double)(grid->dz*param->L));
972:     PetscPrintf(PETSC_COMM_WORLD,"  jlid = %3D              jfault = %3D \n",grid->jlid,grid->jfault);
973:     PetscPrintf(PETSC_COMM_WORLD,"  Pe = %g\n",(double)param->peclet);

975:     PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
976:     if (param->ivisc==VISC_CONST) {
977:       PetscPrintf(PETSC_COMM_WORLD,"                 Isoviscous \n");
978:       if (param->pv_analytic) {
979:         PetscPrintf(PETSC_COMM_WORLD,"                          Pressure and Velocity prescribed! \n");
980:       }
981:     } else if (param->ivisc==VISC_DIFN) {
982:       PetscPrintf(PETSC_COMM_WORLD,"                 Diffusion Creep (T-Dependent Newtonian) \n");
983:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
984:     } else if (param->ivisc==VISC_DISL) {
985:       PetscPrintf(PETSC_COMM_WORLD,"                 Dislocation Creep (T-Dependent Non-Newtonian) \n");
986:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
987:     } else if (param->ivisc==VISC_FULL) {
988:       PetscPrintf(PETSC_COMM_WORLD,"                 Full Rheology \n");
989:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
990:     } else {
991:       PetscPrintf(PETSC_COMM_WORLD,"                 Invalid! \n");
992:       ierr_out = 1;
993:     }

995:     PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
996:     if (param->ibound==BC_ANALYTIC) {
997:       PetscPrintf(PETSC_COMM_WORLD,"       Isoviscous Analytic Dirichlet \n");
998:     } else if (param->ibound==BC_NOSTRESS) {
999:       PetscPrintf(PETSC_COMM_WORLD,"       Stress-Free (normal & shear stress)\n");
1000:     } else if (param->ibound==BC_EXPERMNT) {
1001:       PetscPrintf(PETSC_COMM_WORLD,"       Experimental boundary condition \n");
1002:     } else {
1003:       PetscPrintf(PETSC_COMM_WORLD,"       Invalid! \n");
1004:       ierr_out = 1;
1005:     }

1007:     if (param->output_to_file)
1008: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1009:       PetscPrintf(PETSC_COMM_WORLD,"Output Destination:       Mat file \"%s\"\n",param->filename);
1010: #else
1011:       PetscPrintf(PETSC_COMM_WORLD,"Output Destination:       PETSc binary file \"%s\"\n",param->filename);
1012: #endif
1013:     if (param->output_ivisc != param->ivisc) {
1014:       PetscPrintf(PETSC_COMM_WORLD,"                          Output viscosity: -ivisc %D\n",param->output_ivisc);
1015:     }

1017:     PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
1018:   }
1019:   if (param->param_test) PetscEnd();
1020:   return ierr_out;
1021: }

1023: /* ------------------------------------------------------------------- */
1026: /*  generates an inital guess using the analytic solution for isoviscous
1027:     corner flow */
1028: PetscErrorCode Initialize(DM da)
1029: /* ------------------------------------------------------------------- */
1030: {
1031:   AppCtx         *user;
1032:   Parameter      *param;
1033:   GridInfo       *grid;
1034:   PetscInt       i,j,is,js,im,jm;
1036:   Field          **x;
1037:   Vec            Xguess;

1039:   /* Get the fine grid */
1040:   DMGetApplicationContext(da,&user);
1041:   Xguess = user->Xguess;
1042:   param  = user->param;
1043:   grid   = user->grid;
1044:   DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1045:   DMDAVecGetArray(da,Xguess,(void**)&x);

1047:   /* Compute initial guess */
1048:   for (j=js; j<js+jm; j++) {
1049:     for (i=is; i<is+im; i++) {
1050:       if (i<j)                x[j][i].u = param->cb;
1051:       else if (j<=grid->jlid) x[j][i].u = 0.0;
1052:       else                    x[j][i].u = HorizVelocity(i,j,user);

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

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

1061:       x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1062:     }
1063:   }

1065:   /* Restore x to Xguess */
1066:   DMDAVecRestoreArray(da,Xguess,(void**)&x);

1068:   return 0;
1069: }

1071: /*---------------------------------------------------------------------*/
1074: /*  controls output to a file */
1075: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1076: /*---------------------------------------------------------------------*/
1077: {
1078:   AppCtx         *user;
1079:   Parameter      *param;
1080:   GridInfo       *grid;
1081:   PetscInt       ivt;
1083:   PetscMPIInt    rank;
1084:   PetscViewer    viewer;
1085:   Vec            res, pars;
1086:   MPI_Comm       comm;
1087:   DM             da;

1089:   SNESGetDM(snes,&da);
1090:   DMGetApplicationContext(da,&user);
1091:   param = user->param;
1092:   grid  = user->grid;
1093:   ivt   = param->ivisc;

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

1097:   /* compute final residual and final viscosity/strain rate fields */
1098:   SNESGetFunction(snes, &res, NULL, NULL);
1099:   ViscosityField(da, user->x, user->Xguess);

1101:   /* get the communicator and the rank of the processor */
1102:   PetscObjectGetComm((PetscObject)snes, &comm);
1103:   MPI_Comm_rank(comm, &rank);

1105:   if (param->output_to_file) { /* send output to binary file */
1106:     VecCreate(comm, &pars);
1107:     if (!rank) { /* on processor 0 */
1108:       VecSetSizes(pars, 20, PETSC_DETERMINE);
1109:       VecSetFromOptions(pars);
1110:       VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1111:       VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1112:       VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1113:       VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1114:       VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1115:       VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1116:       /* skipped 6 intentionally */
1117:       VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1118:       VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1119:       VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1120:       VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1121:       VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1122:       VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1123:       VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1124:       VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1125:       VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1126:       VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1127:     } else { /* on some other processor */
1128:       VecSetSizes(pars, 0, PETSC_DETERMINE);
1129:       VecSetFromOptions(pars);
1130:     }
1131:     VecAssemblyBegin(pars); VecAssemblyEnd(pars);

1133:     /* create viewer */
1134: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1135:     PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1136: #else
1137:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1138: #endif

1140:     /* send vectors to viewer */
1141:     PetscObjectSetName((PetscObject)res,"res");
1142:     VecView(res,viewer);
1143:     PetscObjectSetName((PetscObject)user->x,"out");
1144:     VecView(user->x, viewer);
1145:     PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1146:     VecView(user->Xguess, viewer);
1147:     StressField(da); /* compute stress fields */
1148:     PetscObjectSetName((PetscObject)(user->Xguess),"str");
1149:     VecView(user->Xguess, viewer);
1150:     PetscObjectSetName((PetscObject)pars,"par");
1151:     VecView(pars, viewer);

1153:     /* destroy viewer and vector */
1154:     PetscViewerDestroy(&viewer);
1155:     VecDestroy(&pars);
1156:   }

1158:   param->ivisc = ivt;
1159:   return 0;
1160: }

1162: /* ------------------------------------------------------------------- */
1165: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1166: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1167: /* ------------------------------------------------------------------- */
1168: {
1169:   AppCtx         *user;
1170:   Parameter      *param;
1171:   GridInfo       *grid;
1172:   Vec            localX;
1173:   Field          **v, **x;
1174:   PetscReal      eps, /* dx,*/ dz, T, epsC, TC;
1175:   PetscInt       i,j,is,js,im,jm,ilim,jlim,ivt;

1179:   DMGetApplicationContext(da,&user);
1180:   param        = user->param;
1181:   grid         = user->grid;
1182:   ivt          = param->ivisc;
1183:   param->ivisc = param->output_ivisc;

1185:   DMGetLocalVector(da, &localX);
1186:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1187:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1188:   DMDAVecGetArray(da,localX,(void**)&x);
1189:   DMDAVecGetArray(da,V,(void**)&v);

1191:   /* Parameters */
1192:   /* dx = grid->dx; */ dz = grid->dz;

1194:   ilim = grid->ni-1; jlim = grid->nj-1;

1196:   /* Compute real temperature, strain rate and viscosity */
1197:   DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1198:   for (j=js; j<js+jm; j++) {
1199:     for (i=is; i<is+im; i++) {
1200:       T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*param->z_scale));
1201:       if (i<ilim && j<jlim) {
1202:         TC = PetscRealPart(param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*param->z_scale));
1203:       } else {
1204:         TC = T;
1205:       }
1206:       eps  = PetscRealPart((CalcSecInv(x,i,j,CELL_CENTER,user)));
1207:       epsC = PetscRealPart(CalcSecInv(x,i,j,CELL_CORNER,user));

1209:       v[j][i].u = eps;
1210:       v[j][i].w = epsC;
1211:       v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1212:       v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1213:     }
1214:   }
1215:   DMDAVecRestoreArray(da,V,(void**)&v);
1216:   DMDAVecRestoreArray(da,localX,(void**)&x);
1217:   DMRestoreLocalVector(da, &localX);

1219:   param->ivisc = ivt;
1220:   return(0);
1221: }

1223: /* ------------------------------------------------------------------- */
1226: /* post-processing: compute stress everywhere */
1227: PetscErrorCode StressField(DM da)
1228: /* ------------------------------------------------------------------- */
1229: {
1230:   AppCtx         *user;
1231:   PetscInt       i,j,is,js,im,jm;
1233:   Vec            locVec;
1234:   Field          **x, **y;

1236:   DMGetApplicationContext(da,&user);

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

1242:   DMGetLocalVector(da, &locVec);
1243:   DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);
1244:   DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);
1245:   DMDAVecGetArray(da,locVec,(void**)&y);

1247:   /* Compute stress on the corner points */
1248:   for (j=js; j<js+jm; j++) {
1249:     for (i=is; i<is+im; i++) {
1250:       x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1251:       x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1252:       x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1253:       x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1254:     }
1255:   }

1257:   /* Restore the fine grid of Xguess and X */
1258:   DMDAVecRestoreArray(da,user->Xguess,(void**)&x);
1259:   DMDAVecRestoreArray(da,locVec,(void**)&y);
1260:   DMRestoreLocalVector(da, &locVec);
1261:   return 0;
1262: }

1264: /*=====================================================================
1265:   UTILITY FUNCTIONS
1266:   =====================================================================*/

1268: /*---------------------------------------------------------------------*/
1271: /* returns the velocity of the subducting slab and handles fault nodes
1272:    for BC */
1273: PETSC_STATIC_INLINE PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1274: /*---------------------------------------------------------------------*/
1275: {
1276:   Parameter *param = user->param;
1277:   GridInfo  *grid  = user->grid;

1279:   if (c=='U' || c=='u') {
1280:     if (i<j-1) return param->cb;
1281:     else if (j<=grid->jfault) return 0.0;
1282:     else return param->cb;

1284:   } else {
1285:     if (i<j) return param->sb;
1286:     else if (j<=grid->jfault) return 0.0;
1287:     else return param->sb;
1288:   }
1289: }

1291: /*---------------------------------------------------------------------*/
1294: /*  solution to diffusive half-space cooling model for BC */
1295: PETSC_STATIC_INLINE PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1296: /*---------------------------------------------------------------------*/
1297: {
1298:   Parameter     *param = user->param;
1299:   PetscScalar   z;
1300:   if (plate==PLATE_LID) z = (j-0.5)*user->grid->dz;
1301:   else z = (j-0.5)*user->grid->dz*param->cb;  /* PLATE_SLAB */
1302: #if defined(PETSC_HAVE_ERF)
1303:   return (erf(PetscRealPart(z*param->L/2.0/param->skt)));
1304: #else
1305:   SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
1306: #endif
1307: }


1310: /* ------------------------------------------------------------------- */
1313: /*  utility function */
1314: PetscBool  OptionsHasName(const char pre[],const char name[])
1315: /* ------------------------------------------------------------------- */
1316: {
1317:   PetscBool      retval;
1319:   PetscOptionsHasName(NULL,pre,name,&retval);CHKERRABORT(PETSC_COMM_WORLD,ierr);
1320:   return retval;
1321: }

1323: /*=====================================================================
1324:   INTERACTIVE SIGNAL HANDLING
1325:   =====================================================================*/

1327: /* ------------------------------------------------------------------- */
1330: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1331: /* ------------------------------------------------------------------- */
1332: {
1333:   AppCtx         *user  = (AppCtx*) ctx;
1334:   Parameter      *param = user->param;
1335:   KSP            ksp;

1339:   if (param->interrupted) {
1340:     param->interrupted = PETSC_FALSE;
1341:     PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1342:     *reason = SNES_CONVERGED_FNORM_ABS;
1343:     return(0);
1344:   } else if (param->toggle_kspmon) {
1345:     param->toggle_kspmon = PETSC_FALSE;

1347:     SNESGetKSP(snes, &ksp);

1349:     if (param->kspmon) {
1350:       KSPMonitorCancel(ksp);

1352:       param->kspmon = PETSC_FALSE;
1353:       PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1354:     } else {
1355:       PetscViewerAndFormat *vf;
1356:       PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
1357:       KSPMonitorSet(ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSingularValue,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);

1359:       param->kspmon = PETSC_TRUE;
1360:       PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1361:     }
1362:   }
1363:   PetscFunctionReturn(SNESConvergedDefault(snes,it,xnorm,snorm,fnorm,reason,ctx));
1364: }

1366: /* ------------------------------------------------------------------- */
1367: #include <signal.h>
1370: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1371: /* ------------------------------------------------------------------- */
1372: {
1373:   AppCtx    *user  = (AppCtx*) ctx;
1374:   Parameter *param = user->param;

1376:   if (signum == SIGILL) {
1377:     param->toggle_kspmon = PETSC_TRUE;
1378: #if !defined(PETSC_MISSING_SIGCONT)
1379:   } else if (signum == SIGCONT) {
1380:     param->interrupted = PETSC_TRUE;
1381: #endif
1382: #if !defined(PETSC_MISSING_SIGURG)
1383:   } else if (signum == SIGURG) {
1384:     param->stop_solve = PETSC_TRUE;
1385: #endif
1386:   }
1387:   return 0;
1388: }

1390: /*---------------------------------------------------------------------*/
1393: /*  main call-back function that computes the processor-local piece
1394:     of the residual */
1395: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
1396: /*---------------------------------------------------------------------*/
1397: {
1398:   AppCtx      *user  = (AppCtx*)ptr;
1399:   Parameter   *param = user->param;
1400:   GridInfo    *grid  = user->grid;
1401:   PetscScalar mag_w, mag_u;
1402:   PetscInt    i,j,mx,mz,ilim,jlim;
1403:   PetscInt    is,ie,js,je,ibound;    /* ,ivisc */

1406:   /* Define global and local grid parameters */
1407:   mx   = info->mx;     mz   = info->my;
1408:   ilim = mx-1;         jlim = mz-1;
1409:   is   = info->xs;     ie   = info->xs+info->xm;
1410:   js   = info->ys;     je   = info->ys+info->ym;

1412:   /* Define geometric and numeric parameters */
1413:   /* ivisc = param->ivisc; */ ibound = param->ibound;

1415:   for (j=js; j<je; j++) {
1416:     for (i=is; i<ie; i++) {

1418:       /************* X-MOMENTUM/VELOCITY *************/
1419:       if (i<j) f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
1420:       else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1421:         /* in the lithospheric lid */
1422:         f[j][i].u = x[j][i].u - 0.0;
1423:       } else if (i==ilim) {
1424:         /* on the right side boundary */
1425:         if (ibound==BC_ANALYTIC) {
1426:           f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1427:         } else {
1428:           f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1429:         }

1431:       } else if (j==jlim) {
1432:         /* on the bottom boundary */
1433:         if (ibound==BC_ANALYTIC) {
1434:           f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1435:         } else if (ibound==BC_NOSTRESS) {
1436:           f[j][i].u = XMomentumResidual(x,i,j,user);
1437:         } else {
1438:           /* experimental boundary condition */
1439:         }

1441:       } else {
1442:         /* in the mantle wedge */
1443:         f[j][i].u = XMomentumResidual(x,i,j,user);
1444:       }

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

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

1454:       } else if (j==jlim) {
1455:         /* on the bottom boundary */
1456:         if (ibound==BC_ANALYTIC) {
1457:           f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1458:         } else {
1459:           f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1460:         }

1462:       } else if (i==ilim) {
1463:         /* on the right side boundary */
1464:         if (ibound==BC_ANALYTIC) {
1465:           f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1466:         } else if (ibound==BC_NOSTRESS) {
1467:           f[j][i].w = ZMomentumResidual(x,i,j,user);
1468:         } else {
1469:           /* experimental boundary condition */
1470:         }

1472:       } else {
1473:         /* in the mantle wedge */
1474:         f[j][i].w =  ZMomentumResidual(x,i,j,user);
1475:       }

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

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

1486:       } else {
1487:         /* in the mantle wedge */
1488:         f[j][i].p = ContinuityResidual(x,i,j,user);
1489:       }

1491:       /************* TEMPERATURE *************/
1492:       if (j==0) {
1493:         /* on the surface */
1494:         f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0);

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

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

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

1510:       } else {
1511:         /* in the mantle wedge */
1512:         f[j][i].T = EnergyResidual(x,i,j,user);
1513:       }
1514:     }
1515:   }
1516:   return(0);
1517: }