Actual source code: ex30.c

  1: /*$Id: ex19.c,v 1.30 2001/08/07 21:31:17 bsmith Exp $*/

  3: static char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\\n\
  4: The flow is driven by the subducting slab.\n\
  5:   -ivisc <#> = rheology option.\n\
  6:       0 --- constant viscosity.\n\
  7:       1 --- olivine diffusion creep rheology (T-dependent, newtonian).\n\
  8:       2 --- weak temperature dependent rheology (1200/T, newtonian).\n\
  9:   -ibound <#> = boundary condition \n\
 10:       0 --- isoviscous analytic.\n\
 11:       1 --- stress free. \n\
 12:       2 --- stress is von neumann. \n\
 13:   -icorner <#> = i index of wedge corner point.\n\
 14:   -jcorner <#> = j index of wedge corner point.\n\
 15:   -slab_dip <#> = dip of the subducting slab in DEGREES.\n\
 16:   -back_arc <#> = distance from trench to back-arc in KM.(if unspecified then no back-arc). \n\
 17:   -u_back_arcocity <#> = full spreading rate of back arc as a factor of slab velocity. \n\
 18:   -width <#> = width of domain in KM.\n\
 19:   -depth <#> = depth of domain in KM.\n\
 20:   -lid_depth <#> = depth to the base of the lithosphere in KM.\n\
 21:   -slab_dip <#> = dip of the subducting slab in DEGREES.\n\
 22:   -slab_velocity <#> = velocity of slab in CM/YEAR.\n\
 23:   -slab_age <#> = age of slab in MILLIONS OF YEARS. \n\
 24:   -potentialT <#> = mantle potential temperature in degrees CENTIGRADE.\n\
 25:   -kappa <#> = thermal diffusivity in M^2/SEC. \n\
 26:   -peclet <#> = dimensionless Peclet number (default 111.691)\n\\n";

 28: /*T
 29:    Concepts: SNES^solving a system of nonlinear equations (parallel multicomponent example);
 30:    Concepts: DA^using distributed arrays;
 31:    Concepts: multicomponent
 32:    Processors: n
 33: T*/

 35: /* ------------------------------------------------------------------------

 37:     This problem is modeled by the partial differential equation system
 38:   
 39:          -Grad(P) + Div[eta (Grad(v) + Grad(v)^T)] = 0
 40:           Div(U,W) = 0

 42:     which is uniformly discretized on a staggered mesh:
 43:                        ------w_ij------
 44:                       /               /
 45:                   u_i-1j    P_ij    u_ij
 46:                     /               /
 47:                      ------w_ij-1----

 49:   ------------------------------------------------------------------------- */

 51: /* 
 52:    Include "petscda.h" so that we can use distributed arrays (DAs).
 53:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 54:    file automatically includes:
 55:      petsc.h       - base PETSc routines   petscvec.h - vectors
 56:      petscsys.h    - system routines       petscmat.h - matrices
 57:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 58:      petscviewer.h - viewers               petscpc.h  - preconditioners
 59:      petscksp.h   - linear solvers 
 60: */
 61:  #include petscsnes.h
 62:  #include petscda.h

 64: /* 
 65:    User-defined routines and data structures
 66: */

 68: /*
 69:     The next two structures are essentially the same. The difference is that
 70:   the first does not contain derivative information (as used by ADIC) while the
 71:   second does. The first is used only to contain the solution at the previous time-step
 72:   which is a constant in the computation of the current time step and hence passive 
 73:   (meaning not active in terms of derivatives).
 74: */

 76: typedef struct {
 77:   PetscScalar u,w,p,T;
 78: } Field;

 80: typedef struct {
 81:   int          mglevels;
 82:   PetscReal    width, depth, scaled_width, scaled_depth, peclet, potentialT;
 83:   PetscReal    slab_dip, slab_age, slab_velocity, lid_depth, kappa;
 84:   PetscReal    u_back_arc, x_back_arc;
 85:   int          icorner, jcorner, ivisc, ifromm, ibound;
 86:   PetscTruth   PreLoading, back_arc;
 87: } Parameter;

 89: typedef struct {
 90:   Vec          Xold,func;
 91:   Parameter    *param;
 92:   PetscReal    fnorm_ini, fnorm;
 93: } AppCtx;

 95: PetscReal HALFPI = 3.14159265358979323846/2.0;

 97: int SetParams(Parameter *param, int mx, int mz);
 98: extern int FormInitialGuess(SNES,Vec,void*);
 99: extern int FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
100: extern PassiveScalar HorizVelocity(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
101: extern PassiveScalar VertVelocity(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
102: extern PassiveScalar Pressure(PassiveScalar x, PassiveScalar z, PassiveScalar c, PassiveScalar d);
103: extern PassiveScalar LidVelocity(PassiveScalar x, PassiveScalar xBA, PetscTruth BA);
104: extern PetscScalar Viscosity(PetscScalar T, int iVisc);
105: extern PetscScalar UInterp(Field **x, int i, int j, PetscScalar fr);
106: extern PetscScalar WInterp(Field **x, int i, int j, PetscScalar fr);
107: extern PetscScalar PInterp(Field **x, int i, int j, PetscScalar fr);
108: extern PetscScalar TInterp(Field **x, int i, int j, PetscScalar fr);
109: extern void CalcJunk(PassiveScalar,PassiveScalar,PassiveScalar,PassiveScalar*,PassiveScalar*,PassiveScalar*,PassiveScalar*,PassiveScalar*);
110: int Initialize(DMMG*);

112: /*-----------------------------------------------------------------------*/
115: int main(int argc,char **argv)
116: /*-----------------------------------------------------------------------*/
117: {
118:   DMMG       *dmmg;               /* multilevel grid structure */
119:   AppCtx     *user;                /* user-defined work context */
120:   Parameter  param;
121:   int        mx,mz,its;
122:   int        i,ierr,tmpVisc,tmpBound;
123:   MPI_Comm   comm;
124:   SNES       snes;
125:   DA         da;
126:   Vec        res;
127:   PetscReal  SEC_PER_YR = 3600.00*24.00*356.2500;

129:   PetscInitialize(&argc,&argv,(char *)0,help);
130:   PetscOptionsSetValue("-preload","off");
131:   PetscOptionsSetValue("-mat_type","seqaij");
132:   PetscOptionsSetValue("-snes_monitor",PETSC_NULL);
133:   PetscOptionsSetValue("-ksp_truemonitor",PETSC_NULL);
134:   PetscOptionsSetValue("-pc_type","lu");
135:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);

137:   comm = PETSC_COMM_WORLD;
138:   PreLoadBegin(PETSC_TRUE,"SetUp");

140:     param.PreLoading = PreLoading;  //where is PreLoading defined?
141:     DMMGCreate(comm,1,&user,&dmmg);
142:     param.mglevels = DMMGGetLevels(dmmg);

144:     /*
145:       Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
146:       for principal unknowns (x) and governing residuals (f)
147:     */
148:     DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,-4,-4,PETSC_DECIDE,PETSC_DECIDE,4,2,0,0,&da);
149:     DMMGSetDM(dmmg,(DM)da);
150:     DADestroy(da);
151:     DAGetInfo(DMMGGetDA(dmmg),0,&mx,&mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
152:                      PETSC_IGNORE,PETSC_IGNORE);

154:     /* 
155:       Problem parameters
156:     */
157:     SetParams(&param,mx,mz);

159:     DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
160:     DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
161:     DASetFieldName(DMMGGetDA(dmmg),2,"pressure");
162:     DASetFieldName(DMMGGetDA(dmmg),3,"temperature");

164:     /*======================================================================*/
165: 
166:     PetscMalloc(param.mglevels*sizeof(AppCtx),&user);
167:     for (i=0; i<param.mglevels; i++) {
168:       VecDuplicate(dmmg[i]->x, &(user[i].Xold));
169:       VecDuplicate(dmmg[i]->x, &(user[i].func));
170:       user[i].param = &param;
171:       dmmg[i]->user = &user[i];
172:     }
173:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:        Create user context, set problem data, create vector data structures.
175:        Also, compute the initial guess.
176:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: 
178:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179:        Create nonlinear solver context
180:        Process NOT adiC(100): WInterp UInterp PInterp TInterp FormFunctionLocal Viscosity
181:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182:     /* DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,admf_FormFunctionLocal);*/
183:     DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,0,0);
184:     DMMGSetInitialGuess(dmmg,FormInitialGuess);
185: 
186:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187:        Solve the nonlinear system
188:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189:     PreLoadStage("Solve");
190:     Initialize(dmmg);

192:     if (param.ivisc>0) {
193:       PetscPrintf(PETSC_COMM_WORLD,"Doing Constant Viscosity Solve\n", its);
194:       tmpVisc = param.ivisc; param.ivisc=0; tmpBound = param.ibound; param.ibound=0;
195:       snes = DMMGGetSNES(dmmg);
196:       SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,2,PETSC_DEFAULT);
197:       DMMGSolve(dmmg);
198:       SNESSetTolerances(snes,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,50,PETSC_DEFAULT);
199:       VecCopy(DMMGGetx(dmmg),user->Xold);
200:       param.ivisc = tmpVisc; param.ibound=tmpBound;
201:       PetscPrintf(PETSC_COMM_WORLD,"Doing Variable Viscosity Solve\n", its);
202:     }
203:     if (param.ifromm==1)
204:       PetscPrintf(PETSC_COMM_WORLD,"Using Fromm advection scheme\n", its);
205:     DMMGSolve(dmmg);
206:     snes = DMMGGetSNES(dmmg);
207:     SNESGetIterationNumber(snes,&its);
208:     PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n", its);
209: 
210:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211:       Output stuff to Matlab socket.
212:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213:     SNESGetFunction(snes, &res, PETSC_NULL, PETSC_NULL);
214:     VecView(res, PETSC_VIEWER_SOCKET_WORLD);
215:     VecView(DMMGGetx(dmmg),PETSC_VIEWER_SOCKET_WORLD);
216:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&mx);
217:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&mz);
218:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.slab_dip));
219:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.scaled_width));
220:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.scaled_depth));
221:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.lid_depth));
222:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.potentialT));
223:     PetscViewerSocketPutReal(PETSC_VIEWER_SOCKET_WORLD,1,1,&(param.x_back_arc));
224:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.icorner));
225:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.jcorner));
226:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.ivisc));
227:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(param.ibound));
228:     PetscViewerSocketPutInt(PETSC_VIEWER_SOCKET_WORLD,1,&(its));

230:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231:        Free work space.  All PETSc objects should be destroyed when they
232:        are no longer needed.
233:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234:     for (i=0; i<param.mglevels; i++) {
235:       VecDestroy(user[i].Xold);
236:       VecDestroy(user[i].func);
237:     }
238:     PetscFree(user);
239:     DMMGDestroy(dmmg);
240:     PreLoadEnd();
241: 
242:   PetscFinalize();
243:   return 0;
244: }

246: /* ------------------------------------------------------------------- */
249: /* ------------------------------------------------------------------- */
250: int Initialize(DMMG *dmmg)
251: {
252:   AppCtx    *user = (AppCtx*)dmmg[dmmg[0]->nlevels-1]->user;
253:   Parameter *param = user->param;
254:   DA        da;
255:   int       i,j,mx,mz,ierr,xs,ys,xm,ym,iw,jw;
256:   int       mglevel,ic,jc;
257:   PetscReal dx, dz, c, d, beta, itb, cb, sb, sbi, skt;
258:   PetscReal xp, zp, r, st, ct, th, fr, xPrimeSlab;
259:   Field     **x;

261:   da = (DA)(dmmg[param->mglevels-1]->dm); /* getting the fine grid */

263:   beta = param->slab_dip;
264:   CalcJunk(param->kappa,param->slab_age,beta,&skt,&cb,&sb,&itb,&sbi);
265:   c = beta*sb/(beta*beta-sb*sb);
266:   d = (beta*cb-sb)/(beta*beta-sb*sb);

268:   DAGetInfo(da,0,&mx,&mz,0,0,0,0,0,0,0,0);
269:   dx  = param->scaled_width/((PetscReal)(mx-2));
270:   dz  = param->scaled_depth/((PetscReal)(mz-2));

272:   ic = param->icorner; jc = param->jcorner;
273:   xPrimeSlab = (ic-1)*dx;

275:   fr = (1.0-dz/dx*itb)/2.0;
276:   PetscPrintf(PETSC_COMM_WORLD,"interpolation fraction = %g\n", fr);
277:   if (fr<0.0)
278:     SETERRQ(1," USER ERROR: Grid shear exceeds limit! Decrease dz/dx!");

280:   /*
281:      Get local grid boundaries (for 2-dimensional DA):
282:        xs, ys   - starting grid indices (no ghost points)
283:        xm, ym   - scaled_widths of local grid (no ghost points)
284:   */
285:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

287:   /*
288:      Get a pointer to vector data.
289:        - For default PETSc vectors, VecGetArray() returns a pointer to
290:          the data array.  Otherwise, the routine is implementation dependent.
291:        - You MUST call VecRestoreArray() when you no longer need access to
292:          the array.
293:   */
294:   DAVecGetArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,(void**)&x);

296:   /*
297:      Compute initial guess (analytic soln) over the locally owned part of the grid
298:      Initial condition is isoviscous corner flow and uniform temperature in interior
299:   */
300:   for (j=ys; j<ys+ym; j++) {
301:     jw = j - jc + 1;
302:     for (i=xs; i<xs+xm; i++) {
303:       iw = i - ic + 1;

305:       if        (i<ic) { /* slab */
306:         x[j][i].p   = 0.0;
307:         x[j][i].u   = cb;
308:         x[j][i].w   = sb;
309:         xp = (i-0.5)*dx; zp = (xPrimeSlab - xp)*tan(beta);
310:         x[j][i].T   = erf(zp*param->lid_depth/2.0/skt);
311:       } else if (j<jc) { /* lid */
312:         x[j][i].p   = 0.0;
313:         x[j][i].u   = 0.0;
314:         x[j][i].w   = 0.0;
315:         zp = (j-0.5)*dz;
316:         x[j][i].T   = zp;
317:       } else { /* wedge */
318:       /* horizontal velocity */
319:       zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
320:       x[j][i].u =  HorizVelocity(xp,zp,c,d);
321:       /* vertical velocity */
322:       zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
323:       x[j][i].w =  VertVelocity(xp,zp,c,d);
324:       /* pressure */
325:       zp = (jw-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
326:       x[j][i].p =  Pressure(xp,zp,c,d);
327:       /* temperature */
328:       x[j][i].T = 1.0;

330:       }
331:     }
332:   }

334:   /* Trash initial guess */
335:   /* 
336:   for (j=ys; j<ys+ym; j++) {
337:     for (i=xs; i<xs+xm; i++) {
338:       x[j][i].u = 0.0;
339:       x[j][i].w = 0.0;
340:       x[j][i].p = 0.0;
341:       x[j][i].T = 0.0;
342:     }
343:   }
344:   */

346:   /* Restore x to Xold */
347:   DAVecRestoreArray(da,((AppCtx*)dmmg[param->mglevels-1]->user)->Xold,(void**)&x);
348: 
349:   /* Restrict Xold to coarser levels */
350:   for (mglevel=param->mglevels-1; mglevel>0; mglevel--) {
351:     MatRestrict(dmmg[mglevel]->R, ((AppCtx*)dmmg[mglevel]->user)->Xold, ((AppCtx*)dmmg[mglevel-1]->user)->Xold);
352:     VecPointwiseMult(dmmg[mglevel]->Rscale,((AppCtx*)dmmg[mglevel-1]->user)->Xold,((AppCtx*)dmmg[mglevel-1]->user)->Xold);
353:   }
354: 
355:   return 0;
356: }

358: /*---------------------------------------------------------------------*/
361: int FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
362: /*---------------------------------------------------------------------*/
363: {
364:   AppCtx         *user = (AppCtx*)ptr;
365:   Parameter      *param = user->param;
366:   PetscTruth     back_arc;
367:   int            ierr,i,j,mx,mz,sh,ish,ic,jc,iw,jw,ilim,jlim;
368:   int            xints,xinte,yints,yinte,ivisc,ifromm,ibound;
369:   PetscScalar    dudxN,dudxS,dudxW,dudxE,dudzN,dudzS,dudzE,dudzW,dudxC,dudzC;
370:   PetscScalar    dwdzE,dwdzW,dwdxW,dwdxE,dwdzN,dwdzS,dwdxN,dwdxS,dwdxC,dwdzC;
371:   PetscScalar    pN,pS,pE,pW,uE,uW,uC,uN,uS,wN,wS,wE,wW,wC;
372:   PetscScalar    uNE,uNW,uSE,uSW,wNE,wNW,wSE,wSW;
373:   PetscScalar    vE, vN, vS, vW, TE, TN, TS, TW, TC, pC;
374:   PetscScalar    fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS,TNE,TNW,TSE,TSW,dTdxN,dTdxS;
375:   PetscScalar    etaN,etaS,etaE,etaW,etaC;
376:   PassiveScalar  dx, dz, dxp, dzp, c, d, beta, itb, sbi, pe;
377:   PassiveScalar  xp, zp,  cb, sb, skt, z_scale, fr, x_back_arc, u_back_arc;
378:   PassiveScalar  eps=0.000000001, alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;



383:   /* 
384:      Define geometric and numeric parameters
385:   */
386:   back_arc = param->back_arc; x_back_arc = param->x_back_arc; u_back_arc=param->u_back_arc;
387:   pe = param->peclet;  beta = param->slab_dip;
388:   ivisc = param->ivisc; ifromm = param->ifromm;   ibound = param->ibound;
389:   z_scale = param->lid_depth * alpha_g_on_cp_units_inverse_km;
390:   ic = param->icorner;   jc = param->jcorner;
391:   CalcJunk(param->kappa,param->slab_age,beta,&skt,&cb,&sb,&itb,&sbi);
392:   c = beta*sb/(beta*beta-sb*sb);
393:   d = (beta*cb-sb)/(beta*beta-sb*sb);
394: 
395:   /* 
396:      Define global and local grid parameters
397:   */
398:   mx   = info->mx;                         mz   = info->my;
399:   ilim = mx-1;                             jlim = mz-1;
400:   dx   = param->scaled_width/((double)(mx-2));
401:   dz   = param->scaled_depth/((double)(mz-2));
402:   dxp  = dx;                                dzp = dz*sbi;
403:   xints = info->xs;                       xinte = info->xs+info->xm;
404:   yints = info->ys;                       yinte = info->ys+info->ym;
405:   fr   = (1.0-dz/dx*itb)/2.0;

407:   /* 
408:      Stokes equation, no buoyancy terms
409:      Steady state advection-diffusion equation for temperature
410:   */

412:   for (j=yints; j<yinte; j++) {
413:     jw = j-jc+1;
414:     for (i=xints; i<xinte; i++) {
415:       iw = i-ic+1;

417:       /* X-MOMENTUM/VELOCITY */
418:       if        (i<ic) { /* within the slab */
419:         f[j][i].u = x[j][i].u - cb;
420:       } else if (j<jc) { /* within the lid */
421:         f[j][i].u = x[j][i].u - u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
422:       } else if (i==ilim) { /* On the INFLOW boundary */
423:         if (ibound==0) { /* isoviscous analytic soln */
424:           zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
425:           f[j][i].u = x[j][i].u - HorizVelocity(xp,zp,c,d);
426:         } else { /* normal stress = 0 boundary condition */
427:           uE = x[j][i].u; uW = x[j][i-1].u; pC = x[j][i].p;
428:           TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
429:           etaC = Viscosity(TC,ivisc);
430:           f[j][i].u = 2.0*etaC*( uE - uW )/dxp - pC;
431:                }
432:       } else if (j==jlim) { /* On the OUTFLOW boundary */
433:         if (ibound==0) { /* isoviscous analytic soln */
434:           zp = (jw-0.5)*dz; xp = iw*dx+zp*itb;
435:           f[j][i].u = x[j][i].u - HorizVelocity(xp,zp,c,d);
436:         } else { /* shear stress = 0 boundary condition */
437:           uN = x[j][i].u;   wE = x[j-1][i+1].w;
438:           uS = x[j-1][i].u; wW = x[j-1][i].w;
439:           uW = UInterp(x,i-1,j-1,fr); uE = UInterp(x,i,j-1,fr);
440:           f[j][i].u =  sbi*( uN - uS )/dzp
441:                      - itb*( uE - uW )/dxp
442:                      +     ( wE - wW )/dxp;
443:         }
444:       } else {     /* Mantle wedge, horizontal velocity */

446:         pW = x[j][i].p; pE = x[j][i+1].p;

448:         TN = param->potentialT * TInterp(x,i,j,fr)   * exp(  j     *dz*z_scale );
449:         TS = param->potentialT * TInterp(x,i,j-1,fr) * exp( (j-1.0)*dz*z_scale );
450:         TE = param->potentialT * x[j][i+1].T         * exp( (j-0.5)*dz*z_scale );
451:         TW = param->potentialT * x[j][i].T           * exp( (j-0.5)*dz*z_scale );
452:         etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
453:         etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);
454:         /*if (j==jc) etaS = 1.0;*/

456:         /* ------ BEGIN VAR VISCOSITY USE ONLY ------- */
457:         dwdxN = etaN * ( x[j][i+1].w   - x[j][i].w   )/dxp;
458:         dwdxS = etaS * ( x[j-1][i+1].w - x[j-1][i].w )/dxp;
459:         if (i<ilim-1) { wE = WInterp(x,i+1,j-1,fr); }
460:         else { wE = ( x[j][i].w + x[j-1][i].w )/2.0; }
461:         wC = WInterp(x,i,j-1,fr);
462:         wW = WInterp(x,i-1,j-1,fr);   if (i==ic) wW = sb;
463:         dwdxE = etaE * ( wE - wC )/dxp;
464:         dwdxW = etaW * ( wC - wW )/dxp;
465:         /* ------ END VAR VISCOSITY USE ONLY ------- */

467:         /* ------ BGN ISOVISC BETA != 0 USE ONLY ------- */
468:         uNE = UInterp(x,i,j,fr);   uNW = UInterp(x,i-1,j,fr);
469:         uSE = UInterp(x,i,j-1,fr); uSW = UInterp(x,i-1,j-1,fr);
470:         if (j==jc) {
471:           xp = (iw+0.5)*dx+(j-1)*dz*itb; uSE = u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
472:           xp = (iw-0.5)*dx+(j-1)*dz*itb; uSW = u_back_arc*LidVelocity(xp,x_back_arc,back_arc);
473:         }
474:         dudxN = etaN * ( uNE - uNW )/dxp; dudxS = etaS * ( uSE - uSW )/dxp;
475:         dudzE = etaE * ( uNE - uSE )/dzp; dudzW = etaW * ( uNW - uSW )/dzp;
476:         /* ------ END ISOVISC BETA != 0 USE ONLY ------- */

478:         dudzN = etaN * ( x[j+1][i].u  - x[j][i].u   )/dzp;
479:         dudzS = etaS * ( x[j][i].u    - x[j-1][i].u )/dzp;
480:         dudxE = etaE * ( x[j][i+1].u  - x[j][i].u   )/dxp;
481:         dudxW = etaW * ( x[j][i].u    - x[j][i-1].u )/dxp;
482:         if (j==jc) {
483:           if (ibound==0) { /* apply isoviscous boundary condition */
484:             xp = iw*dx;
485:             dudzS = etaS * sb*(HorizVelocity(xp,eps,c,d)-HorizVelocity(xp,-eps,c,d))/eps/2.0;
486:           } else  /* force u=0 on the lid-wedge interface (off-grid point) */
487:             dudzS = etaS * ( 2.0*x[j][i].u - 2.0*x[j-1][i].u )/dzp;
488:         }

490:         f[j][i].u = -( pE - pW )/dxp                         /* X-MOMENTUM EQUATION*/
491:                     +( dudxE - dudxW )/dxp * (1.0+itb*itb)
492:                     +( dudzN - dudzS )/dzp * sbi*sbi
493:                     -( ( dudxN - dudxS )/dzp
494:                       +( dudzE - dudzW )/dxp ) * itb*sbi;

496:         if (ivisc>0) {
497:           f[j][i].u = f[j][i].u + ( dudxE - dudxW )/dxp
498:                                 + ( dwdxN - dwdxS )/dzp * sbi
499:                                 - ( dwdxE - dwdxW )/dxp * itb;
500:         }
501:       }

503:       /* Z-MOMENTUM/VELOCITY */
504:       if        (i<ic) {  /* within the slab */
505:         f[j][i].w = x[j][i].w - sb;
506:       } else if (j<jc) {  /* within the lid */
507:         f[j][i].w = x[j][i].w - 0.0;
508:       } else if (j==jlim) { /* On the OUTFLOW boundary */
509:         if (ibound==0) { /* isoviscous analytic soln */
510:           zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
511:           f[j][i].w = x[j][i].w - VertVelocity(xp,zp,c,d);
512:         } else { /* normal stress = 0 boundary condition */
513:           wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
514:           wW = WInterp(x,i-1,j-1,fr); if (i==ic) wW = sb;
515:           if (i==ilim) wE = ( x[j][i].w + x[j-1][i].w )/2.0;
516:           else wE = WInterp(x,i,j-1,fr);
517:           TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
518:           etaC = Viscosity(TC,ivisc);
519:           f[j][i].w = 2.0*etaC*( sbi*( wN - wS )/dzp
520:                                 -itb*( wE - wW )/dxp ) - pC;
521:         }
522:       } else if (i==ilim) { /* On the INFLOW boundary */
523:         if (ibound==0) { /* isoviscous analytic soln */
524:           zp = jw*dz; xp = (iw-0.5)*dx+zp*itb;
525:           f[j][i].w = x[j][i].w - VertVelocity(xp,zp,c,d);
526:         } else { /* shear stress = 0 boundary condition */
527:           uN = x[j+1][i-1].u; wE = x[j][i].w;
528:           uS = x[j][i-1].u;   wW = x[j][i-1].w;
529:           uW = UInterp(x,i-2,j,fr); uE = UInterp(x,i-1,j,fr);
530:           f[j][i].w =  sbi*( uN - uS )/dzp
531:                      - itb*( uE - uW )/dxp
532:                      +     ( wE - wW )/dxp;
533:           if (j==jlim-1) {
534:             f[j][i].w = x[j][i].w - x[j-1][i].w;
535:           }
536:         }
537:       } else {   /* Mantle wedge, vertical velocity */
538: 
539:         pE = PInterp(x,i,j,fr); pW = PInterp(x,i-1,j,fr); pS = x[j][i].p; pN = x[j+1][i].p;
540:               if ( (i==ic) && (ibound==0) ) { zp = jw*dz; xp = zp*itb; pW = Pressure(xp,zp,c,d); }
541: 
542:         TN = param->potentialT * x[j+1][i].T         * exp( (j+0.5)*dz*z_scale );
543:         TS = param->potentialT * x[j][i].T           * exp( (j-0.5)*dz*z_scale );
544:         TE = param->potentialT * TInterp(x,i,j,fr)   * exp(  j     *dz*z_scale );
545:         TW = param->potentialT * TInterp(x,i-1,j,fr) * exp(  j     *dz*z_scale );
546:         etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
547:         etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);

549:          /* ------ BGN VAR VISCOSITY USE ONLY ------- */
550:         dudzE = etaE * ( x[j+1][i].u   - x[j][i].u   )/dzp;
551:         dudzW = etaW * ( x[j+1][i-1].u - x[j][i-1].u )/dzp;
552:         uE = UInterp(x,i,j,fr);   uC = UInterp(x,i-1,j,fr);
553:         uW = UInterp(x,i-2,j,fr); if (i==ic) uW = 2.0*cb - uC;
554:         dudxE = etaE * ( uE - uC )/dxp;
555:         dudxW = etaW * ( uC - uW )/dxp;
556:         /* ------ END VAR VISCOSITY USE ONLY ------- */

558:         /* ------ BGN ISOVISC BETA != 0 USE ONLY ------- */
559:         wNE = WInterp(x,i,j,fr);   wSE = WInterp(x,i,j-1,fr);
560:         wNW = WInterp(x,i-1,j,fr); wSW = WInterp(x,i-1,j-1,fr);
561:         if (i==ic) { wNW = wSW = sb; }
562:         dwdzE = etaE * ( wNE - wSE )/dzp; dwdzW = etaW * ( wNW - wSW )/dzp;
563:         dwdxN = etaN * ( wNE - wNW )/dxp; dwdxS = etaS * ( wSE - wSW )/dxp;
564:         /* ------ END ISOVISC BETA != 0 USE ONLY ------- */

566:         dwdzN = etaN * ( x[j+1][i].w - x[j][i].w   )/dzp;
567:         dwdzS = etaS * ( x[j][i].w   - x[j-1][i].w )/dzp;
568:         dwdxE = etaE * ( x[j][i+1].w - x[j][i].w   )/dxp;
569:         dwdxW = etaW * ( x[j][i].w - x[j][i-1].w   )/dxp;
570:         if (i==ic) {
571:           if (ibound==0) { /* apply isoviscous boundary condition */
572:             zp = jw*dz; xp = itb*zp;
573:             dwdxW = etaW * ( VertVelocity(xp+eps,zp,c,d) - VertVelocity(xp,zp,c,d) )/eps;
574:           } else /*  force w=sin(beta) on the slab-wedge interface (off-grid point) */
575:             dwdxW = etaW * ( 2.0*x[j][i].w - 2.0*sb )/dxp;
576:         }

578:          /* Z-MOMENTUM */
579:         f[j][i].w =  ( pE - pW )/dxp * itb                 /* constant viscosity terms */
580:                     -( pN - pS )/dzp * sbi
581:                     +( dwdzN - dwdzS )/dzp * sbi*sbi
582:                     +( dwdxE - dwdxW )/dxp * (itb*itb+1.0)
583:                     -( ( dwdzE - dwdzW )/dxp
584:                       +( dwdxN - dwdxS )/dzp ) * itb*sbi;

586:         if (ivisc>0) {
587:           f[j][i].w = f[j][i].w + ( dwdxE - dwdxW )/dxp * itb*itb
588:                                 + ( dwdzN - dwdzS )/dzp * sbi*sbi
589:                                 -( ( dwdzE - dwdzW )/dxp
590:                                   +( dwdxN - dwdxS )/dzp ) * itb*sbi
591:                                 - ( dudxE - dudxW )/dxp * itb
592:                                 + ( dudzE - dudzW )/dxp * sbi;
593:         }
594:       }
595: 
596:       /* CONTINUITY/PRESSURE */
597:       if ( (j<jc) || (i<ic-1) ) { /* within slab or lid */
598:         f[j][i].p = x[j][i].p - 0.0;

600:       } else if ( (j==jlim)&&(i==ic-1) ) {
601:         f[j][i].p = x[j][i].p - 0.0;
602: 
603:       } else if ( (ibound==0) && ((i==ilim)||(j==jlim)) ) { /* isoviscous inflow/outflow BC */
604:         zp = (jw-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
605:         f[j][i].p = x[j][i].p - Pressure(xp,zp,c,d);

607:       } else if (i==ic-1) { /* just west of the slab-wedge interface constrain 
608:                                pressure using the x-momentum equation */

610:         pE = x[j][i+1].p; pW = x[j][i].p; /* pW IS THE UNKNOWN */
611: 
612:         TN = param->potentialT * TInterp(x,i,j,fr)   * exp(  j     *dz*z_scale );
613:         TS = param->potentialT * TInterp(x,i,j-1,fr) * exp( (j-1.0)*dz*z_scale );
614:         TE = param->potentialT * x[j][i+1].T         * exp( (j-0.5)*dz*z_scale );
615:         TW = param->potentialT * x[j][i].T           * exp( (j-0.5)*dz*z_scale );
616:         etaN = Viscosity(TN,ivisc); etaS = Viscosity(TS,ivisc);
617:         etaW = Viscosity(TW,ivisc); etaE = Viscosity(TE,ivisc);

619:         /* ------ BGN VAR VISCOSITY USE ONLY ------- */
620:         dwdxN = etaN * ( x[j][i+1].w   - (2*sb - x[j][i+1].w)   )/dxp;
621:         dwdxS = etaS * ( x[j-1][i+1].w - (2*sb - x[j-1][i+1].w) )/dxp;
622:         wE = WInterp(x,i+1,j-1,fr);   wC = sb; wW = 2*sb - wE;
623:         dwdxE = etaE * ( wE - wC )/dxp;
624:         dwdxW = etaW * ( wC - wW )/dxp;
625:         /* ------ END VAR VISCOSITY USE ONLY ------- */

627:         /* ------ BGN BETA != 0 USE ONLY ------- */
628:         uNE = UInterp(x,i,j,fr);   uNW = 2*cb - uNE;
629:         uSE = UInterp(x,i,j-1,fr); uSW = 2*cb - uSE;
630:         if (j==jc) uSE = 0.0;
631:         dudxN = etaN * ( uNE - uNW )/dxp; dudxS = etaS * ( uSE - uSW )/dxp;
632:         dudzE = etaE * ( uNE - uSE )/dzp; dudzW = 0.0;
633:         if (j==jc) dudxS = 0.0;
634:         /* ------ BGN BETA != 0 USE ONLY ------- */

636:         dudxE = etaE * ( x[j][i+1].u  - x[j][i].u   )/dxp;
637:         dudxW = 0.0;

639:         f[j][i].p = -( pE - pW )/dxp                                 /* X-MOMENTUM */
640:                     +( dudxE - dudxW )/dxp * (1.0+itb*itb)
641:                     -( ( dudxN - dudxS )/dzp
642:                       +( dudzE - dudzW )/dxp ) * itb*sbi;

644:         if (ivisc>0) {
645:           f[j][i].p = f[j][i].p + ( dudxE - dudxW )/dxp
646:                                 + ( dwdxN - dwdxS )/dzp * sbi
647:                                 - ( dwdxE - dwdxW )/dxp * itb;
648:         }
649:       } else { /* interior of the domain */
650:         uW = x[j][i-1].u; uE = x[j][i].u;
651:         wS = x[j-1][i].w; wN = x[j][i].w;
652:         wW = WInterp(x,i-1,j-1,fr); if (i==ic) wW = sb;
653:         if (i==ilim) wE = ( x[j][i].w + x[j-1][i].w )/2.0;
654:         else wE = WInterp(x,i,j-1,fr);

656:         f[j][i].p = ( uE - uW )/dxp   /* CONTINUITY ON PRESSURE POINTS */
657:                    -( wE - wW )/dxp * itb
658:                    +( wN - wS )/dzp * sbi;

660:         if ( (ivisc==10)&&(i<ilim)&&(j<jlim) ) {
661:           wW = x[j][i].w;   wN = WInterp(x,i,j,fr);
662:           wE = x[j][i+1].w; wS = WInterp(x,i,j-1,fr);
663:           uW = UInterp(x,i-1,j,fr);  uE = UInterp(x,i,j,fr);
664: 
665:           f[j][i].p = f[j][i].p +
666:                       ( uE - uW )/dxp -  /* CONTINUITY ON VACANT POINTS */
667:                       ( wE - wW )/dxp * itb +
668:                       ( wN - wS )/dzp * sbi;
669:         }
670: 
671:       }
672: 
673:       /* TEMPERATURE EQUATION */
674:       zp = (j-0.5)*dz; xp = (iw-0.5)*dx+zp*itb;
675:       if (i<=1) {    /* dirichlet on boundary along slab side */
676:         f[j][i].T = x[j][i].T - 1.0;
677:       } else if ( (j<=2) && (i<ic) ) {   /* slab inflow dirichlet */
678:         f[j][i].T = x[j][i].T - erf(-xp*sb*param->lid_depth/2.0/skt);
679:       } else if (j==0) {   /* force T=0 on surface */
680:         f[j][i].T = x[j][i].T + x[j+1][i].T;
681:       } else if (j>=jlim-1) { /* neumann on outflow boundary */
682:         if (x[j][i].w<0.0)
683:           f[j][i].T = x[j][i].T - 1.0;
684:         else
685:           f[j][i].T = x[j][i].T - x[j-1][i].T;
686:       } else if (i>=ilim-1) /* (xp>=20.0) */ {
687:         if (back_arc && (x[j][i].u>0))
688:           f[j][i].T = x[j][i].T - x[j][i-1].T;
689:         else
690:           f[j][i].T = x[j][i].T - PetscMin(zp,1.0);
691:       } else {                 /* interior of the domain */

693:         uW = x[j][i-1].u; uE = x[j][i].u;
694:         wS = x[j-1][i].w; wN = x[j][i].w;
695:         wE = WInterp(x,i,j-1,fr); wW = WInterp(x,i-1,j-1,fr);

697:         if (i==ic) { /* Just east of slab-wedge interface */
698:           if (j<jc) {
699:             wW = wE = wN = wS = 0.0;
700:             uW = uE = 0.0;
701:           } else {
702:             wW = sb;
703:           }
704:         }

706:         if (i==ic-1) wE = sb; /* Just west of slab-wedge interface */

708:         TNE =   TInterp(x,i,j,fr);              TNW =   TInterp(x,i-1,j,fr);
709:         TSE =   TInterp(x,i,j-1,fr);            TSW =   TInterp(x,i-1,j-1,fr);
710:         dTdxN = ( TNE - TNW )/dxp;
711:         dTdxS = ( TSE - TSW )/dxp;
712:         dTdzN = ( x[j+1][i].T - x[j][i].T   )/dzp;
713:         dTdzS = ( x[j][i].T   - x[j-1][i].T )/dzp;
714:         dTdxE = ( x[j][i+1].T - x[j][i].T   )/dxp;
715:         dTdxW = ( x[j][i].T   - x[j][i-1].T )/dxp;
716: 
717:         f[j][i].T = ( ( dTdzN - dTdzS )/dzp * sbi*sbi + /* diffusion term */
718:                       ( dTdxE - dTdxW )/dxp * (1.0+itb*itb) -
719:                       ( dTdxN - dTdxS )/dzp * 2.0*itb*sbi  )*dx*dz/pe;

721:         if ( (j<jc) && (i>=ic) ) { /* don't advect in lid */
722:             fE = fW = fN = fS = 0.0;

724:         } else if ( (ifromm==0) ||(i>=ilim-2)||(j>=jlim-2)||(i<=2) ) { /* finite volume advection */
725:           TN  = ( x[j][i].T + x[j+1][i].T )/2.0;  TS  = ( x[j][i].T + x[j-1][i].T )/2.0;
726:           TE  = ( x[j][i].T + x[j][i+1].T )/2.0;  TW  = ( x[j][i].T + x[j][i-1].T )/2.0;
727:           fN = wN*TN*dxp; fS = wS*TS*dxp;
728:           fE = ( uE*sb - wE*cb )*TE*dzp;
729:           fW = ( uW*sb - wW*cb )*TW*dzp;
730: 
731:         } else {         /* Fromm advection scheme */
732:           vN = wN; vS = wS; vE = uE*sb - wE*cb; vW = uW*sb - wW*cb;
733:           fE =     ( vE *(-x[j][i+2].T + 5.0*(x[j][i+1].T+x[j][i].T)-x[j][i-1].T)/8.0
734:                      - fabs(vE)*(-x[j][i+2].T + 3.0*(x[j][i+1].T-x[j][i].T)+x[j][i-1].T)/8.0 )*dzp;
735:           fW =     ( vW *(-x[j][i+1].T + 5.0*(x[j][i].T+x[j][i-1].T)-x[j][i-2].T)/8.0
736:                      - fabs(vW)*(-x[j][i+1].T + 3.0*(x[j][i].T-x[j][i-1].T)+x[j][i-2].T)/8.0 )*dzp;
737:           fN =     ( vN *(-x[j+2][i].T + 5.0*(x[j+1][i].T+x[j][i].T)-x[j-1][i].T)/8.0
738:                      - fabs(vN)*(-x[j+2][i].T + 3.0*(x[j+1][i].T-x[j][i].T)+x[j-1][i].T)/8.0 )*dxp;
739:           fS =     ( vS *(-x[j+1][i].T + 5.0*(x[j][i].T+x[j-1][i].T)-x[j-2][i].T)/8.0
740:                      - fabs(vS)*(-x[j+1][i].T + 3.0*(x[j][i].T-x[j-1][i].T)+x[j-2][i].T)/8.0 )*dxp;
741:         }
742: 
743:           f[j][i].T = f[j][i].T -
744:             ( fE - fW + fN - fS );
745:       }
746:     }
747:   }
748: 
749:   return(0);
750: }

752: /* ------------------------------------------------------------------- */
755: /* 
756:    FormInitialGuess - Forms initial approximation.

758:    Input Parameters:
759:    user - user-defined application context
760:    X - vector

762:    Output Parameter:
763:    X - vector
764:  */
765: /* ------------------------------------------------------------------- */
766: int FormInitialGuess(SNES snes,Vec X,void *ptr)
767: {
768:   DMMG      dmmg = (DMMG)ptr;
769:   AppCtx    *user = (AppCtx*)dmmg->user;
770:   int       ierr;

772:   VecCopy(user->Xold, X);

774:   /* calculate the residual on fine mesh, but only the first time this is called */
775:   if (user->fnorm_ini == 0.0) {
776:     SNESComputeFunction(snes,user->Xold,user->func);
777:     VecNorm(user->func,NORM_2,&user->fnorm_ini);
778:   }
779: 
780:   return 0;
781: }

783: /*--------------------------UTILITY FUNCTION BELOW THIS LINE-----------------------------*/

785: /*---------------------------------------------------------------------*/
788: int SetParams(Parameter *param, int mx, int mz)
789: /*---------------------------------------------------------------------*/
790: {
792:   PetscReal  SEC_PER_YR = 3600.00*24.00*356.2500;

794:   /* domain geometry */
795:   param->icorner       = (int)(mx/6+1);     /* gridpoints */
796:   param->jcorner       = (int)((mz-2)/12+1);/* gridpoints */
797:   param->width         = 1200.0;            /* km */
798:   param->depth         = 600.0;             /* km */
799:   param->slab_dip      = HALFPI;            /* 90 degrees */
800:   PetscOptionsGetInt(PETSC_NULL, "-icorner",&(param->icorner),PETSC_NULL);
801:   PetscOptionsGetInt(PETSC_NULL, "-jcorner",&(param->jcorner),PETSC_NULL);
802:   PetscOptionsGetReal(PETSC_NULL,"-width",&(param->width),PETSC_NULL);
803:   PetscOptionsGetReal(PETSC_NULL,"-depth",&(param->depth),PETSC_NULL);
804:   PetscOptionsGetReal(PETSC_NULL,"-slab_dip",&(param->slab_dip),PETSC_NULL);

806:   /* back-arc */
807:   param->back_arc      = PETSC_FALSE;       /* no back arc spreading */
808:   param->x_back_arc    = 0.0;               /* km */
809:   param->u_back_arc    = 1.0;               /* full spreading at velocity of slab */
810:   PetscOptionsHasName(PETSC_NULL,"-back_arc",&(param->back_arc));
811:   PetscOptionsGetReal(PETSC_NULL,"-back_arc",&(param->x_back_arc),PETSC_NULL);
812:   PetscOptionsGetReal(PETSC_NULL,"-back_arc_velocity",&(param->u_back_arc),PETSC_NULL);
813:   if (param->back_arc) {
814:     PetscPrintf(PETSC_COMM_WORLD,"Dist to back arc = %g km, ",param->x_back_arc);
815:     PetscPrintf(PETSC_COMM_WORLD,"Full spreading rate of back arc (scaled) = %g \n",param->u_back_arc);
816:   }

818:   /* physics parameters */
819:   param->slab_velocity = 5.0;               /* cm/yr */
820:   param->slab_age      = 50.0;              /* Ma */
821:   param->kappa         = 0.7272e-6;         /* m^2/sec */
822:   param->potentialT    = 1300.0;            /* degrees C */
823:   param->ivisc         = 0;                 /* 0=constant, 1=diffusion creep, 2=simple T dependent */
824:   PetscOptionsGetReal(PETSC_NULL,"-slab_velocity",&(param->slab_velocity),PETSC_NULL);
825:   PetscOptionsGetReal(PETSC_NULL,"-slab_age",&(param->slab_age),PETSC_NULL);
826:   PetscOptionsGetReal(PETSC_NULL,"-kappa",&(param->kappa),PETSC_NULL);
827:   PetscOptionsGetReal(PETSC_NULL,"-potentialT",&(param->potentialT),PETSC_NULL);
828:   PetscOptionsGetInt(PETSC_NULL, "-ivisc",&(param->ivisc),PETSC_NULL);

830:   /* boundaries */
831:   param->ibound = param->ivisc;       /* 0=isovisc analytic, 1,2,...= stress free */
832:   PetscOptionsGetInt(PETSC_NULL,"-ibound",&(param->ibound),PETSC_NULL);

834:   /* misc */
835:   param->ifromm = 1;                 /* advection scheme: 0=finite vol, 1=Fromm */
836:   PetscOptionsGetInt(PETSC_NULL,"-ifromm",&(param->ifromm),PETSC_NULL);

838:   /* unit conversions and derived parameters */
839:   param->lid_depth = (param->jcorner - 1) * param->depth/((double)(mz-2));
840:   PetscOptionsGetReal(PETSC_NULL,"-lid_depth",&(param->lid_depth),PETSC_NULL);
841:   param->slab_dip =     param->slab_dip*HALFPI/90.0;
842:   param->scaled_width = param->width/param->lid_depth;
843:   param->scaled_depth = param->depth/param->lid_depth;
844:   param->x_back_arc     = param->x_back_arc/param->lid_depth;
845:   param->peclet =  param->slab_velocity/100.0/SEC_PER_YR /* m/sec */
846:                  * param->lid_depth*1000.0               /* m */
847:                  / param->kappa;                         /* m^2/sec */
848:   PetscOptionsGetReal(PETSC_NULL,"-peclet",&(param->peclet),PETSC_NULL);
849:   PetscPrintf(PETSC_COMM_WORLD,"Lid depth = %g km, ",param->lid_depth);
850:   PetscPrintf(PETSC_COMM_WORLD,"Peclet number = %g\n",param->peclet);
851:   if ( (param->ibound==0) && (param->ivisc>0) )
852:     PetscPrintf(PETSC_COMM_WORLD,"Warning: isoviscous BC may be inconsistent w/ var viscosity!!\n");

854:   return 0;
855: }


858: /*---------------------------------------------------------------------*/
861: PetscScalar Viscosity(PetscScalar T, int iVisc)
862: /*---------------------------------------------------------------------*/
863: {
864:   PetscScalar  result;
865:   double       p1 = 1.32047792e-12, p2 = 335.0e3/8.314510;
866:   double       t1 = 1200.0;
867:   /*
868:     p1 = exp( -p2/(1473 K) ) so the cutoff for high viscosity 
869:     occurs at 1200 C.  Below this cutoff, all eta( T<1200 ) = 1.0;
870:     p2=Ea/R is from van Keken's subroutine
871:   */

873:   if (iVisc==0) {        /* constant viscosity */
874:     result = 1.0;
875:   } else if (iVisc==1) { /* diffusion creep rheology */
876:     result = p1*PetscExpScalar(p2/(T+273.0));
877:   } else if (iVisc==2) { /* ad hoc T-dependent rheology */
878:     result = t1/T;
879:   } else if (iVisc==3) {
880:     result = 1.0;
881:   }

883:   if (result<1.0)
884:     return result;
885:   else
886:     return 1.0;
887: }

889: /*---------------------------------------------------------------------*/
892: PassiveScalar LidVelocity(PassiveScalar x, PassiveScalar xBA, PetscTruth BA)
893: /*---------------------------------------------------------------------*/
894: {
895:   PassiveScalar localize = 10.0;

897:    if (BA)
898:      return ( 1.0 + tanh( (x - xBA)*localize ) )/2.0;
899:    else
900:      return 0.0;
901: }

903: /*---------------------------------------------------------------------*/
906: PetscScalar UInterp(Field **x, int i, int j, PetscScalar fr)
907: /*---------------------------------------------------------------------*/
908: {
909:   PetscScalar p,m;
910:   p = (1.0-fr)*x[j+1][i].u + fr*x[j+1][i+1].u;
911:   m = (1.0-fr)*x[j][i+1].u + fr*x[j][i].u;
912:   return (p + m)/2.0;
913: }

915: /*---------------------------------------------------------------------*/
918: PetscScalar WInterp(Field **x, int i, int j, PetscScalar fr)
919: /*---------------------------------------------------------------------*/
920: {
921:   PetscScalar p,m;
922:   p = (1.0-fr)*x[j+1][i].w + fr*x[j+1][i+1].w;
923:   m = (1.0-fr)*x[j][i+1].w + fr*x[j][i].w;
924:   return (p + m)/2.0;
925: }

927: /*---------------------------------------------------------------------*/
930: PetscScalar PInterp(Field **x, int i, int j, PetscScalar fr)
931: /*---------------------------------------------------------------------*/
932: {
933:   PetscScalar p,m;
934:   p = (1.0-fr)*x[j+1][i].p + fr*x[j+1][i+1].p;
935:   m = (1.0-fr)*x[j][i+1].p + fr*x[j][i].p;
936:   return (p + m)/2.0;
937: }

939: /*---------------------------------------------------------------------*/
942: PetscScalar TInterp(Field **x, int i, int j, PetscScalar fr)
943: /*---------------------------------------------------------------------*/
944: {
945:   PetscScalar p,m;
946:   p = (1.0-fr)*x[j+1][i].T + fr*x[j+1][i+1].T;
947:   m = (1.0-fr)*x[j][i+1].T + fr*x[j][i].T;
948:   return (p + m)/2.0;
949: }

951: /*---------------------------------------------------------------------*/
954: PetscScalar HorizVelocity(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
955: /*---------------------------------------------------------------------*/
956:  {
957:    PetscScalar r, st, ct, th;
958:    r = sqrt(x*x+z*z);
959:    st = z/r;  ct = x/r;  th = atan(z/x);
960:    return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
961: }

963: /*---------------------------------------------------------------------*/
966: PetscScalar VertVelocity(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
967: /*---------------------------------------------------------------------*/
968:  {
969:    PetscScalar r, st, ct, th;
970:    r = sqrt(x*x+z*z);
971:    st = z/r;  ct = x/r;  th = atan(z/x);
972:    return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
973: }

975: /*---------------------------------------------------------------------*/
978: PetscScalar Pressure(PetscScalar x, PetscScalar z, PetscScalar c, PetscScalar d)
979: /*---------------------------------------------------------------------*/
980: {
981:   PetscScalar r, st, ct;
982:   r = sqrt(x*x+z*z);
983:   st = z/r;  ct = x/r;
984:   return (-2.0*(c*ct-d*st)/r);
985: }

987: /*---------------------------------------------------------------------*/
990: void CalcJunk(double kappa,double slab_age,double beta,double *skt,double *cb,double *sb, double *itb,
991:               double *sbi)
992: /*---------------------------------------------------------------------*/
993: {
994:   PetscReal SEC_PER_YR = 3600.00*24.00*356.2500;

996:   *skt = sqrt(kappa*slab_age*SEC_PER_YR);
997:   *cb  = cos(beta); *sb = sin(beta);
998:   *itb = 1.0/tan(beta); *sbi = 1.0/(*sb);
999: }