Actual source code: ex30.c

      1: static char help[] =
      2: "ex30: Steady-state 2D subduction flow, pressure and temperature solver.\n\
      3:        The flow is driven by the subducting slab.\n\
      4: ---------------------------------ex30 help---------------------------------\n\
      5:   -OPTION <DEFAULT> = (UNITS) DESCRIPTION.\n\n\
      6:   -width <320> = (km) width of domain.\n\
      7:   -depth <300> = (km) depth of domain.\n\
      8:   -slab_dip <45> = (degrees) dip angle of the slab (determines the grid aspect ratio).\n\
      9:   -lid_depth <35> = (km) depth of the static conductive lid.\n\
     10:   -fault_depth <35> = (km) depth of slab-wedge mechanical coupling\n\
     11:      ( fault dept >= lid depth ).\n\
     12: \n\
     13:   -ni <82> = grid cells in x-direction. (nj adjusts to accommodate\n\
     14:       the slab dip & depth). DO NOT USE -da_grid_x option!!!\n\
     15:   -ivisc <3> = rheology option.\n\
     16:       0 --- constant viscosity.\n\
     17:       1 --- olivine diffusion creep rheology (T&P-dependent, newtonian).\n\
     18:       2 --- olivine dislocation creep rheology (T&P-dependent, non-newtonian).\n\
     19:       3 --- Full mantle rheology, combination of 1 & 2.\n\
     20: \n\
     21:   -slab_velocity <5> = (cm/year) convergence rate of slab into subduction zone.\n\
     22:   -slab_age <50> = (million yrs) age of slab for thermal profile boundary condition.\n\
     23:   -lid_age <50> = (million yrs) age of lid for thermal profile boundary condition.\n\
     24: \n\
     25:   FOR OTHER PARAMETER OPTIONS AND THEIR DEFAULT VALUES, see SetParams() in ex30.c.\n\
     26: ---------------------------------ex30 help---------------------------------\n";

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

     34:     The problem is modeled by the partial differential equation system
     35:
     36:        (1)   -Grad(P) + Div[Eta (Grad(v) + Grad(v)^T)] = 0
     37:        (2)                                    Div(U,W) = 0
     38:        (3)             dT/dt + Div(vT) - 1/Pe Del^2(T) = 0
     39:        (4)   Eta(T,Eps_dot) = constant                             if ivisc==0
     40:                             = diffusion creep (T,P-dependent)      if ivisc==1
     41:                             = dislocation creep (T,P,v-dependent)  if ivisc==2
     42:                                   = mantle viscosity (difn & disl)       if ivisc==3

     44:     which is uniformly discretized on a staggered mesh:
     45:                       -------w_ij------
     46:                       |               |
     47:                   u_i-1j    P,T_ij   u_ij
     48:                       |               |
     49:                        ------w_ij-1-----

     51:   ------------------------------------------------------------------------- */

     53:  #include petscsnes.h
     54:  #include petscda.h
     55:  #include petscdmmg.h

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

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

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

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

     92: typedef struct { /* grid parameters */
     93:   DAPeriodicType periodic;
     94:   DAStencilType  stencil;
     95:   PetscInt       corner,ni,nj,jlid,jfault,inose;
     96:   PetscInt       dof,stencil_width,mglevels;
     97:   PassiveScalar  dx,dz;
     98: } GridInfo;

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

    106: /* Callback functions (static interface) */

    110: /* Main routines */

    117: /* Physics subroutines */

    128: /* Utilities for interpolation, ICs and BCs */

    139: /* Post-processing & misc */

    146: /*-----------------------------------------------------------------------*/
    149: int main(int argc,char **argv)
    150: /*-----------------------------------------------------------------------*/
    151: {
    152:   DMMG           *dmmg;               /* multilevel grid structure */
    153:   AppCtx         *user;               /* user-defined work context */
    154:   Parameter      param;
    155:   GridInfo       grid;
    156:   PetscInt       nits;
    158:   MPI_Comm       comm;
    159:   DA             da;

    161:   PetscInitialize(&argc,&argv,(char *)0,help);
    162:   PetscOptionsSetValue("-file","ex30_output");
    163:   PetscOptionsSetValue("-snes_monitor",PETSC_NULL);
    164:   PetscOptionsSetValue("-snes_max_it","20");
    165:   PetscOptionsSetValue("-ksp_max_it","1500");
    166:   PetscOptionsSetValue("-ksp_gmres_restart","300");
    167:   PetscOptionsInsert(&argc,&argv,PETSC_NULL);

    169:   comm = PETSC_COMM_WORLD;

    171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:      Set up the problem parameters.
    173:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    174:   SetParams(&param,&grid);
    175:   ReportParams(&param,&grid);

    177: #if 0
    178:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    179:      Create user context, set problem data, create vector data structures.
    180:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    181:   PetscMalloc(sizeof(AppCtx),&user);
    182:   user->param = &param;
    183:   user->grid  = &grid;

    185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    186:      Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
    187:      for principal unknowns (x) and governing residuals (f)
    188:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    189:   DMMGCreate(comm,grid.mglevels,user,&dmmg);
    190:   DACreate2d(comm,grid.periodic,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
    191:   DMMGSetDM(dmmg,(DM)da);
    192:   DADestroy(da);
    193:   DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
    194:   DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
    195:   DASetFieldName(DMMGGetDA(dmmg),2,"pressure");
    196:   DASetFieldName(DMMGGetDA(dmmg),3,"temperature");
    197:   VecDuplicate(dmmg[0]->x, &(user->Xguess));
    198: #else
    199:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    200:      Create distributed array multigrid object (DMMG) to manage parallel grid and vectors
    201:      for principal unknowns (x) and governing residuals (f)
    202:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    203:   DMMGCreate(comm,grid.mglevels,&user,&dmmg);
    204:   DACreate2d(comm,grid.periodic,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
    205:   DMMGSetDM(dmmg,(DM)da);
    206:   DADestroy(da);
    207:   DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");
    208:   DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");
    209:   DASetFieldName(DMMGGetDA(dmmg),2,"pressure");
    210:   DASetFieldName(DMMGGetDA(dmmg),3,"temperature");

    212:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    213:      Create user context, set problem data, create vector data structures.
    214:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    215:   PetscMalloc(sizeof(AppCtx),&user);
    216:   user->param   = &param;
    217:   user->grid    = &grid;
    218:   dmmg[0]->user = user;
    219:   VecDuplicate(dmmg[0]->x, &(user->Xguess));
    220: #endif

    222:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    223:      Set up the SNES solver with callback functions.
    224:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    225:   DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,0,0);
    226:   DMMGSetInitialGuess(dmmg,FormInitialGuess);
    227:   SNESSetConvergenceTest(DMMGGetSNES(dmmg),SNESConverged_Interactive,(void*)user);
    228:   PetscPushSignalHandler(InteractiveHandler,(void*)user);
    229:
    230:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    231:      Initialize and solve the nonlinear system
    232:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    233:   Initialize(dmmg);
    234:   UpdateSolution(dmmg,user,&nits);
    235:
    236:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    237:      Output variables.
    238:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    239:   DoOutput(dmmg,nits);
    240:
    241:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    242:      Free work space.
    243:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    244:   VecDestroy(user->Xguess);
    245:   PetscFree(user);
    246:   DMMGDestroy(dmmg);
    247:
    248:   PetscFinalize();
    249:   return 0;
    250: }

    252: /*=====================================================================
    253:   PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
    254:   =====================================================================*/

    256: /*---------------------------------------------------------------------*/
    259: /*  manages solve: adaptive continuation method  */
    260: PetscErrorCode UpdateSolution(DMMG *dmmg, AppCtx *user, PetscInt *nits)
    261: {
    262:   SNES                snes;
    263:   KSP                 ksp;
    264:   PC                  pc;
    265:   SNESConvergedReason reason;
    266:   Parameter           *param = user->param;
    267:   PassiveScalar       cont_incr=0.3;
    268:   PetscInt            its;
    269:   PetscErrorCode      ierr;
    270:   PetscTruth          q = PETSC_FALSE;

    273:   snes = DMMGGetSNES(dmmg);
    274:   SNESGetKSP(snes,&ksp);
    275:   KSPGetPC(ksp, &pc);
    276:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);

    278:   *nits=0;

    280:   /* Isoviscous solve */
    281:   if (param->ivisc == VISC_CONST && !param->stop_solve) {
    282:     param->ivisc = VISC_CONST;
    283:     DMMGSolve(dmmg);
    284:     VecCopy(DMMGGetx(dmmg),user->Xguess);
    285:     SNESGetIterationNumber(snes, &its);
    286:     *nits +=its;
    287:     if (param->stop_solve) goto done;
    288:   }

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

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

    298:       /* solve the non-linear system */
    299:       DMMGSolve(dmmg);
    300:       SNESGetConvergedReason(snes,&reason);
    301:       SNESGetIterationNumber(snes,&its);
    302:       *nits += its;
    303:       if (!q) PetscPrintf(PETSC_COMM_WORLD," Newton iterations: %D, Cumulative: %D\n", its, *nits);
    304:       if (param->stop_solve) goto done;

    306:       if (reason<0) {
    307:         /* NOT converged */
    308:         cont_incr = -fabs(cont_incr)/2.0;
    309:         if (fabs(cont_incr)<0.01) goto done;

    311:       } else {
    312:         /* converged */
    313:         VecCopy(DMMGGetx(dmmg),user->Xguess);
    314:         if (param->continuation >= 1.0) goto done;
    315:         if (its<=3) {
    316:           cont_incr = 0.30001;
    317:         } else if (its<=8) {
    318:           cont_incr = 0.15001;
    319:         } else {
    320:           cont_incr = 0.10001;
    321:         }
    322:         if (param->continuation+cont_incr > 1.0) {
    323:           cont_incr = 1.0 - param->continuation;
    324:         }
    325:       } /* endif reason<0 */
    326:     }
    327:   }
    328:   done:
    329:   if (param->stop_solve && !q) PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");
    330:   if (reason<0 && !q) PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");
    331:   return(0);
    332: }

    334: /* ------------------------------------------------------------------- */
    337: /*  used by SNESSolve to get an initial guess for the solution X */
    338: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
    339: /* ------------------------------------------------------------------- */
    340: {
    341:   AppCtx         *user = (AppCtx*)dmmg->user;

    344:   VecCopy(user->Xguess, X);
    345:   return 0;
    346: }

    348: /*=====================================================================
    349:   PHYSICS FUNCTIONS (compute the discrete residual)
    350:   =====================================================================*/

    352: /*---------------------------------------------------------------------*/
    355: /*  main call-back function that computes the processor-local piece
    356:     of the residual */
    357: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
    358: /*---------------------------------------------------------------------*/
    359: {
    360:   AppCtx         *user = (AppCtx*)ptr;
    361:   Parameter      *param = user->param;
    362:   GridInfo       *grid  = user->grid;
    363:   PetscScalar    mag_w, mag_u;
    364:   PetscInt       i,j,mx,mz,ilim,jlim;
    365:   PetscInt       is,ie,js,je,ivisc,ibound;

    369:   /* Define global and local grid parameters */
    370:   mx   = info->mx;     mz   = info->my;
    371:   ilim = mx-1;         jlim = mz-1;
    372:   is   = info->xs;     ie   = info->xs+info->xm;
    373:   js   = info->ys;     je   = info->ys+info->ym;

    375:   /* Define geometric and numeric parameters */
    376:   ivisc = param->ivisc;       ibound = param->ibound;

    378:   for (j=js; j<je; j++) {
    379:     for (i=is; i<ie; i++) {

    381:       /************* X-MOMENTUM/VELOCITY *************/
    382:       if (i<j) {
    383:           f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);

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

    389:       } else if (i==ilim) {
    390:         /* on the right side boundary */
    391:         if (ibound==BC_ANALYTIC) {
    392:           f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
    393:         } else {
    394:           f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
    395:         }

    397:       } else if (j==jlim) {
    398:         /* on the bottom boundary */
    399:         if (ibound==BC_ANALYTIC) {
    400:           f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
    401:         } else if (ibound==BC_NOSTRESS) {
    402:           f[j][i].u = XMomentumResidual(x,i,j,user);
    403:         } else {
    404:           /* experimental boundary condition */
    405:         }

    407:       } else {
    408:         /* in the mantle wedge */
    409:         f[j][i].u = XMomentumResidual(x,i,j,user);
    410:       }
    411:
    412:       /************* Z-MOMENTUM/VELOCITY *************/
    413:       if (i<=j) {
    414:         f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);

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

    420:       } else if (j==jlim) {
    421:         /* on the bottom boundary */
    422:         if (ibound==BC_ANALYTIC) {
    423:           f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
    424:         } else {
    425:           f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
    426:         }

    428:       } else if (i==ilim) {
    429:         /* on the right side boundary */
    430:         if (ibound==BC_ANALYTIC) {
    431:           f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
    432:         } else if (ibound==BC_NOSTRESS) {
    433:           f[j][i].w = ZMomentumResidual(x,i,j,user);
    434:         } else {
    435:           /* experimental boundary condition */
    436:         }

    438:       } else {
    439:         /* in the mantle wedge */
    440:         f[j][i].w =  ZMomentumResidual(x,i,j,user);
    441:       }

    443:       /************* CONTINUITY/PRESSURE *************/
    444:       if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
    445:         /* in the lid or slab */
    446:         f[j][i].p = x[j][i].p;
    447:
    448:       } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
    449:         /* on an analytic boundary */
    450:         f[j][i].p = x[j][i].p - Pressure(i,j,user);

    452:       } else {
    453:         /* in the mantle wedge */
    454:         f[j][i].p = ContinuityResidual(x,i,j,user);
    455:       }

    457:       /************* TEMPERATURE *************/
    458:       if (j==0) {
    459:         /* on the surface */
    460:         f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(x[j][i].T,0.0);

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

    466:       } else if (i==ilim) {
    467:         /* right side boundary */
    468:         mag_u = 1.0 - pow( (1.0-PetscMax(PetscMin(x[j][i-1].u/param->cb,1.0),0.0)), 5.0 );
    469:         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);

    471:       } else if (j==jlim) {
    472:         /* bottom boundary */
    473:         mag_w = 1.0 - pow( (1.0-PetscMax(PetscMin(x[j-1][i].w/param->sb,1.0),0.0)), 5.0 );
    474:         f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);

    476:       } else {
    477:         /* in the mantle wedge */
    478:         f[j][i].T = EnergyResidual(x,i,j,user);
    479:       }
    480:     }
    481:   }
    482:   return(0);
    483: }

    485: /*---------------------------------------------------------------------*/
    488: /*  computes the residual of the x-component of eqn (1) above */
    489: PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    490: /*---------------------------------------------------------------------*/
    491: {
    492:   Parameter      *param=user->param;
    493:   GridInfo       *grid =user->grid;
    494:   PetscScalar    dx = grid->dx, dz=grid->dz;
    495:   PetscScalar    etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
    496:   PetscScalar    TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
    497:   PetscScalar    dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
    498:   PetscInt            jlim = grid->nj-1;

    500:   z_scale = param->z_scale;

    502:   if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
    503:     TS = param->potentialT * TInterp(x,i,j-1) * exp( (j-1.0)*dz*z_scale );
    504:     if (j==jlim) TN = TS;
    505:     else         TN = param->potentialT * TInterp(x,i,j)   * exp(  j     *dz*z_scale );
    506:     TW = param->potentialT * x[j][i].T        * exp( (j-0.5)*dz*z_scale );
    507:     TE = param->potentialT * x[j][i+1].T      * exp( (j-0.5)*dz*z_scale );
    508:     if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
    509:       epsN = CalcSecInv(x,i,j,  CELL_CORNER,user);
    510:       epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
    511:       epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
    512:       epsW = CalcSecInv(x,i,j,  CELL_CENTER,user);
    513:     }
    514:   }
    515:   etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
    516:   etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
    517:   etaW = Viscosity(TW,epsW,dz*j,param);
    518:   etaE = Viscosity(TE,epsE,dz*j,param);

    520:   dPdx = ( x[j][i+1].p - x[j][i].p )/dx;
    521:   if (j==jlim) dudzN = etaN * ( x[j][i].w   - x[j][i+1].w )/dx;
    522:   else         dudzN = etaN * ( x[j+1][i].u - x[j][i].u   )/dz;
    523:   dudzS = etaS * ( x[j][i].u    - x[j-1][i].u )/dz;
    524:   dudxE = etaE * ( x[j][i+1].u  - x[j][i].u   )/dx;
    525:   dudxW = etaW * ( x[j][i].u    - x[j][i-1].u )/dx;

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

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

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

    538:   return residual;
    539: }

    541: /*---------------------------------------------------------------------*/
    544: /*  computes the residual of the z-component of eqn (1) above */
    545: PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    546: /*---------------------------------------------------------------------*/
    547: {
    548:   Parameter      *param=user->param;
    549:   GridInfo       *grid =user->grid;
    550:   PetscScalar    dx = grid->dx, dz=grid->dz;
    551:   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;
    552:   PetscScalar    TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
    553:   PetscScalar    dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
    554:   PetscInt            ilim = grid->ni-1;

    556:   /* geometric and other parameters */
    557:   z_scale = param->z_scale;
    558:
    559:   /* viscosity */
    560:   if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
    561:     TN = param->potentialT * x[j+1][i].T      * exp( (j+0.5)*dz*z_scale );
    562:     TS = param->potentialT * x[j][i].T        * exp( (j-0.5)*dz*z_scale );
    563:     TW = param->potentialT * TInterp(x,i-1,j) * exp(  j     *dz*z_scale );
    564:     if (i==ilim) TE = TW;
    565:     else         TE = param->potentialT * TInterp(x,i,j)   * exp(  j*dz*z_scale );
    566:     if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
    567:       epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
    568:       epsS = CalcSecInv(x,i,j,  CELL_CENTER,user);
    569:       epsE = CalcSecInv(x,i,j,  CELL_CORNER,user);
    570:       epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
    571:     }
    572:   }
    573:   etaN = Viscosity(TN,epsN,dz*(j+1),param);
    574:   etaS = Viscosity(TS,epsS,dz*j,param);
    575:   etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
    576:   etaE = Viscosity(TE,epsE,dz*(j+0.5),param);

    578:   dPdz = ( x[j+1][i].p - x[j][i].p )/dz;
    579:   dwdzN = etaN * ( x[j+1][i].w - x[j][i].w )/dz;
    580:   dwdzS = etaS * ( x[j][i].w - x[j-1][i].w )/dz;
    581:   if (i==ilim) dwdxE = etaE * ( x[j][i].u   - x[j+1][i].u )/dz;
    582:   else         dwdxE = etaE * ( x[j][i+1].w - x[j][i].w   )/dx;
    583:   dwdxW = 2.0*etaW * ( x[j][i].w - x[j][i-1].w )/dx;
    584:
    585:   /* Z-MOMENTUM */
    586:   residual  = -dPdz                /* constant viscosity terms */
    587:               +( dwdzN - dwdzS )/dz
    588:               +( dwdxE - dwdxW )/dx;

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

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

    597:   return residual;
    598: }

    600: /*---------------------------------------------------------------------*/
    603: /*  computes the residual of eqn (2) above */
    604: PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    605: /*---------------------------------------------------------------------*/
    606: {
    607:   GridInfo       *grid =user->grid;
    608:   PetscScalar    uE,uW,wN,wS,dudx,dwdz;

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

    613:   return dudx + dwdz;
    614: }

    616: /*---------------------------------------------------------------------*/
    619: /*  computes the residual of eqn (3) above */
    620: PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
    621: /*---------------------------------------------------------------------*/
    622: {
    623:   Parameter      *param=user->param;
    624:   GridInfo       *grid =user->grid;
    625:   PetscScalar    dx = grid->dx, dz=grid->dz;
    626:   PetscInt            ilim=grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
    627:   PetscScalar    TE, TN, TS, TW, residual;
    628:   PetscScalar    uE,uW,wN,wS;
    629:   PetscScalar    fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;

    631:   dTdzN = ( x[j+1][i].T - x[j][i].T   )/dz;
    632:   dTdzS = ( x[j][i].T   - x[j-1][i].T )/dz;
    633:   dTdxE = ( x[j][i+1].T - x[j][i].T   )/dx;
    634:   dTdxW = ( x[j][i].T   - x[j][i-1].T )/dx;
    635:
    636:   residual = ( ( dTdzN - dTdzS )/dz + /* diffusion term */
    637:                ( dTdxE - dTdxW )/dx  )*dx*dz/param->peclet;

    639:   if (j<=jlid && i>=j) {
    640:     /* don't advect in the lid */
    641:     return residual;

    643:   } else if (i<j) {
    644:     /* beneath the slab sfc */
    645:     uW = uE = param->cb;
    646:     wS = wN = param->sb;

    648:   } else {
    649:     /* advect in the slab and wedge */
    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:   }

    654:   if ( param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1 ) {
    655:     /* finite volume advection */
    656:     TS  = ( x[j][i].T + x[j-1][i].T )/2.0;
    657:     TN  = ( x[j][i].T + x[j+1][i].T )/2.0;
    658:     TE  = ( x[j][i].T + x[j][i+1].T )/2.0;
    659:     TW  = ( x[j][i].T + x[j][i-1].T )/2.0;
    660:     fN = wN*TN*dx; fS = wS*TS*dx;
    661:     fE = uE*TE*dz; fW = uW*TW*dz;
    662:
    663:   } else {
    664:     /* Fromm advection scheme */
    665:     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
    666:                - fabs(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;
    667:     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
    668:                - fabs(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;
    669:     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
    670:                - fabs(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;
    671:     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
    672:                - fabs(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;
    673:   }
    674:
    675:   residual -= ( fE - fW + fN - fS );

    677:   return residual;
    678: }

    680: /*---------------------------------------------------------------------*/
    683: /*  computes the shear stress---used on the boundaries */
    684: PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    685: /*---------------------------------------------------------------------*/
    686: {
    687:   Parameter    *param=user->param;
    688:   GridInfo     *grid=user->grid;
    689:   PetscInt          ilim=grid->ni-1, jlim=grid->nj-1;
    690:   PetscScalar  uN, uS, wE, wW;

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

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

    696:     wE = WInterp(x,i,j-1);
    697:     if (i==j) { wW = param->sb; uN = param->cb;}
    698:     else      { wW = WInterp(x,i-1,j-1); uN = UInterp(x,i-1,j); }
    699:     if (j==grid->jlid+1)  uS = 0.0;
    700:     else                  uS = UInterp(x,i-1,j-1);

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

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

    707:   }

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

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

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

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

    730:     TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
    731:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
    732:     etaC = Viscosity(TC,epsC,dz*j,param);

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

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

    740:     TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
    741:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
    742:     etaC = Viscosity(TC,epsC,dz*(j+0.5),param);

    744:     if (i==j) uW = param->sb;
    745:     else      uW = UInterp(x,i-1,j);
    746:     uE = UInterp(x,i,j); pC = PInterp(x,i,j);
    747:   }
    748:
    749:   return 2.0*etaC*(uE-uW)/dx - pC;
    750: }

    752: /*---------------------------------------------------------------------*/
    755: /*  computes the normal stress---used on the boundaries */
    756: PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    757: /*---------------------------------------------------------------------*/
    758: {
    759:   Parameter      *param=user->param;
    760:   GridInfo       *grid =user->grid;
    761:   PetscScalar    dz=grid->dz;
    762:   PetscInt            ilim=grid->ni-1, jlim=grid->nj-1, ivisc;
    763:   PetscScalar    epsC=0.0, etaC, TC;
    764:   PetscScalar    pC, wN, wS, z_scale;
    765:   if (i<j || j<=grid->jlid) return EPS_ZERO;

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

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

    771:     TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
    772:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
    773:     etaC = Viscosity(TC,epsC,dz*j,param);
    774:     wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;

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

    779:     TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
    780:     if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
    781:     etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
    782:     if (i==j) wN = param->sb;
    783:     else      wN = WInterp(x,i,j);
    784:     wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
    785:   }

    787:   return  2.0*etaC*(wN-wS)/dz - pC;
    788: }

    790: /*---------------------------------------------------------------------*/
    793: /*  computes the second invariant of the strain rate tensor */
    794: PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
    795: /*---------------------------------------------------------------------*/
    796: {
    797:   Parameter     *param = user->param;
    798:   GridInfo      *grid  = user->grid;
    799:   PetscInt           ilim=grid->ni-1, jlim=grid->nj-1;
    800:   PetscScalar   uN,uS,uE,uW,wN,wS,wE,wW;
    801:   PetscScalar   eps11, eps12, eps22;

    803:   if (i<j) return EPS_ZERO;
    804:   if (i==ilim) i--; if (j==jlim) j--;

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

    809:     uE = x[j][i].u; uW = x[j][i-1].u;
    810:     wN = x[j][i].w; wS = x[j-1][i].w;
    811:     wE = WInterp(x,i,j-1);
    812:     if (i==j) { uN = param->cb; wW = param->sb; }
    813:     else      { uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1); }

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

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

    821:     uN = x[j+1][i].u;  uS = x[j][i].u;
    822:     wE = x[j][i+1].w;  wW = x[j][i].w;
    823:     if (i==j) { wN = param->sb; uW = param->cb; }
    824:     else      { wN = WInterp(x,i,j); uW = UInterp(x,i-1,j); }

    826:     if (j==grid->jlid) {
    827:       uE = 0.0;  uW = 0.0;
    828:       uS = -uN;
    829:       wS = -wN;
    830:     } else {
    831:       uE = UInterp(x,i,j);
    832:       wS = WInterp(x,i,j-1);
    833:     }
    834:   }

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

    839:   return sqrt( 0.5*( eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22 ) );
    840: }

    842: /*---------------------------------------------------------------------*/
    845: /*  computes the shear viscosity */
    846: PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PassiveScalar z,
    847:                        Parameter *param)
    848: /*---------------------------------------------------------------------*/
    849: {
    850:   PetscScalar  result=0.0;
    851:   ViscParam    difn=param->diffusion, disl=param->dislocation;
    852:   PetscInt          iVisc=param->ivisc;
    853:   double       eps_scale=param->V/(param->L*1000.0);
    854:   double       strain_power, v1, v2, P;
    855:   double       rho_g = 32340.0, R=8.3144;

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

    859:   if        (iVisc==VISC_CONST) {
    860:     /* constant viscosity */
    861:     return 1.0;

    863:   } else if (iVisc==VISC_DIFN) {
    864:     /* diffusion creep rheology */
    865:     result = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;

    867:   } else if (iVisc==VISC_DISL) {
    868:     /* dislocation creep rheology */
    869:     strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
    870:     result = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;

    872:   } else if (iVisc==VISC_FULL) {
    873:     /* dislocation/diffusion creep rheology */
    874:     strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
    875:     v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
    876:     v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
    877:     result = 1.0/(1.0/v1 + 1.0/v2);
    878:   }

    880:   /* max viscosity is param->eta0 */
    881:   result = PetscMin( result, 1.0 );
    882:   /* min viscosity is param->visc_cutoff */
    883:   result = PetscMax( result, param->visc_cutoff );
    884:   /* continuation method */
    885:   result = pow(result,param->continuation);
    886:   return result;
    887: }

    889: /*=====================================================================
    890:   INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
    891:   =====================================================================*/

    893: /*---------------------------------------------------------------------*/
    896: /* initializes the problem parameters and checks for
    897:    command line changes */
    898: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
    899: /*---------------------------------------------------------------------*/
    900: {
    901:   PetscErrorCode ierr, ierr_out=0;
    902:   PetscReal      SEC_PER_YR = 3600.00*24.00*365.2500;
    903:   PetscReal      PI = 3.14159265358979323846;
    904:   PetscReal      alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
    905:
    906:   /* domain geometry */
    907:   param->slab_dip      = 45.0;
    908:   param->width         = 320.0;                                            /* km */
    909:   param->depth         = 300.0;                                            /* km */
    910:   param->lid_depth     = 35.0;                                             /* km */
    911:   param->fault_depth   = 35.0;                                             /* km */
    912:   PetscOptionsGetReal(PETSC_NULL,"-slab_dip",&(param->slab_dip),PETSC_NULL);
    913:   PetscOptionsGetReal(PETSC_NULL,"-width",&(param->width),PETSC_NULL);
    914:   PetscOptionsGetReal(PETSC_NULL,"-depth",&(param->depth),PETSC_NULL);
    915:   PetscOptionsGetReal(PETSC_NULL,"-lid_depth",&(param->lid_depth),PETSC_NULL);
    916:   PetscOptionsGetReal(PETSC_NULL,"-fault_depth",&(param->fault_depth),PETSC_NULL);
    917:   param->slab_dip      = param->slab_dip*PI/180.0;                    /* radians */

    919:   /* grid information */
    920:   PetscOptionsGetInt(PETSC_NULL, "-jfault",&(grid->jfault),PETSC_NULL);
    921:   grid->ni             = 82;
    922:   PetscOptionsGetInt(PETSC_NULL, "-ni",&(grid->ni),PETSC_NULL);
    923:   grid->dx             = param->width/((double)(grid->ni-2));              /* km */
    924:   grid->dz             = grid->dx*tan(param->slab_dip);                    /* km */
    925:   grid->nj             = (PetscInt)(param->depth/grid->dz + 3.0);        /* gridpoints*/
    926:   param->depth         = grid->dz*(grid->nj-2);                            /* km */
    927:   grid->inose          = 0;                                         /* gridpoints*/
    928:   PetscOptionsGetInt(PETSC_NULL,"-inose",&(grid->inose),PETSC_NULL);
    929:   grid->periodic       = DA_NONPERIODIC;
    930:   grid->stencil        = DA_STENCIL_BOX;
    931:   grid->dof            = 4;
    932:   grid->stencil_width  = 2;
    933:   grid->mglevels       = 1;

    935:   /* boundary conditions */
    936:   param->pv_analytic        = PETSC_FALSE;
    937:   param->ibound             = BC_NOSTRESS;
    938:   PetscOptionsGetInt(PETSC_NULL,"-ibound",&(param->ibound),PETSC_NULL);

    940:   /* physical constants */
    941:   param->slab_velocity = 5.0;               /* cm/yr */
    942:   param->slab_age      = 50.0;              /* Ma */
    943:   param->lid_age       = 50.0;              /* Ma */
    944:   param->kappa         = 0.7272e-6;         /* m^2/sec */
    945:   param->potentialT    = 1300.0;            /* degrees C */
    946:   PetscOptionsGetReal(PETSC_NULL,"-slab_velocity",&(param->slab_velocity),PETSC_NULL);
    947:   PetscOptionsGetReal(PETSC_NULL,"-slab_age",&(param->slab_age),PETSC_NULL);
    948:   PetscOptionsGetReal(PETSC_NULL,"-lid_age",&(param->lid_age),PETSC_NULL);
    949:   PetscOptionsGetReal(PETSC_NULL,"-kappa",&(param->kappa),PETSC_NULL);
    950:   PetscOptionsGetReal(PETSC_NULL,"-potentialT",&(param->potentialT),PETSC_NULL);

    952:   /* viscosity */
    953:   param->ivisc         = 3;                 /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
    954:   param->eta0          = 1e24;              /* Pa-s */
    955:   param->visc_cutoff   = 0.0;               /* factor of eta_0 */
    956:   param->continuation  = 1.0;
    957:   /* constants for diffusion creep */
    958:   param->diffusion.A       = 1.8e7;           /* Pa-s */
    959:   param->diffusion.n       = 1.0;             /* dim'less */
    960:   param->diffusion.Estar   = 375e3;           /* J/mol */
    961:   param->diffusion.Vstar   = 5e-6;            /* m^3/mol */
    962:   /* constants for param->dislocationocation creep */
    963:   param->dislocation.A     = 2.8969e4;        /* Pa-s */
    964:   param->dislocation.n     = 3.5;             /* dim'less */
    965:   param->dislocation.Estar = 530e3;           /* J/mol */
    966:   param->dislocation.Vstar = 14e-6;           /* m^3/mol */
    967:   PetscOptionsGetInt(PETSC_NULL, "-ivisc",&(param->ivisc),PETSC_NULL);
    968:   PetscOptionsGetReal(PETSC_NULL,"-visc_cutoff",&(param->visc_cutoff),PETSC_NULL);
    969:   param->output_ivisc  = param->ivisc;
    970:   PetscOptionsGetInt(PETSC_NULL,"-output_ivisc",&(param->output_ivisc),PETSC_NULL);
    971:   PetscOptionsGetReal(PETSC_NULL,"-vstar",&(param->dislocation.Vstar),PETSC_NULL);

    973:   /* output options */
    974:   param->quiet            = PETSC_FALSE;
    975:   param->param_test       = PETSC_FALSE;
    976:   PetscOptionsHasName(PETSC_NULL,"-quiet",&(param->quiet));
    977:   PetscOptionsHasName(PETSC_NULL,"-test",&(param->param_test));
    978:   PetscOptionsGetString(PETSC_NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));

    980:   /* advection */
    981:   param->adv_scheme       = ADVECT_FROMM;       /* advection scheme: 0=finite vol, 1=Fromm */
    982:   PetscOptionsGetInt(PETSC_NULL,"-adv_scheme",&(param->adv_scheme),PETSC_NULL);

    984:   /* misc. flags */
    985:   param->stop_solve          = PETSC_FALSE;
    986:   param->interrupted         = PETSC_FALSE;
    987:   param->kspmon              = PETSC_FALSE;
    988:   param->toggle_kspmon       = PETSC_FALSE;

    990:   /* derived parameters for slab angle */
    991:   param->sb  = sin(param->slab_dip);
    992:   param->cb  = cos(param->slab_dip);
    993:   param->c   =  param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
    994:   param->d   = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);

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

    1000:   /* other unit conversions and derived parameters */
    1001:   param->scaled_width  = param->width/param->L;                 /* dim'less */
    1002:   param->scaled_depth  = param->depth/param->L;                 /* dim'less */
    1003:   param->lid_depth     = param->lid_depth/param->L;             /* dim'less */
    1004:   param->fault_depth   = param->fault_depth/param->L;           /* dim'less */
    1005:   grid->dx             = grid->dx/param->L;                     /* dim'less */
    1006:   grid->dz             = grid->dz/param->L;                     /* dim'less */
    1007:   grid->jlid           = (PetscInt)(param->lid_depth/grid->dz);      /* gridcells */
    1008:   grid->jfault         = (PetscInt)(param->fault_depth/grid->dz);    /* gridcells */
    1009:   param->lid_depth     = grid->jlid*grid->dz;                   /* dim'less */
    1010:   param->fault_depth   = grid->jfault*grid->dz;                 /* dim'less */
    1011:   grid->corner         = grid->jlid+1;                          /* gridcells */
    1012:   param->peclet        = param->V                               /* m/sec */
    1013:                        * param->L*1000.0                        /* m */
    1014:                        / param->kappa;                          /* m^2/sec */
    1015:   param->z_scale       = param->L * alpha_g_on_cp_units_inverse_km;
    1016:   param->skt           = sqrt(param->kappa*param->slab_age*SEC_PER_YR);
    1017:   PetscOptionsGetReal(PETSC_NULL,"-peclet",&(param->peclet),PETSC_NULL);
    1018:
    1019:   return ierr_out;
    1020: }

    1022: /*---------------------------------------------------------------------*/
    1025: /*  prints a report of the problem parameters to stdout */
    1026: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
    1027: /*---------------------------------------------------------------------*/
    1028: {
    1029:   PetscErrorCode ierr, ierr_out=0;
    1030:   char           date[30];
    1031:   PetscReal      PI = 3.14159265358979323846;

    1033:   PetscGetDate(date,30);

    1035:   if ( !(param->quiet) ) {
    1036:     PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
    1037:     PetscPrintf(PETSC_COMM_WORLD,"                   %s\n",&(date[0]));

    1039:     PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
    1040:     PetscPrintf(PETSC_COMM_WORLD,"  Width = %g km,         Depth = %g km\n",param->width,param->depth);
    1041:     PetscPrintf(PETSC_COMM_WORLD,"  Slab dip = %g degrees,  Slab velocity = %g cm/yr\n",param->slab_dip*180.0/PI,param->slab_velocity);
    1042:     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);

    1044:     PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
    1045:     PetscPrintf(PETSC_COMM_WORLD,"  [ni,nj] = %D, %D       [dx,dz] = %g, %g km\n",grid->ni,grid->nj,grid->dx*param->L,grid->dz*param->L);
    1046:     PetscPrintf(PETSC_COMM_WORLD,"  jlid = %3D              jfault = %3D \n",grid->jlid,grid->jfault);
    1047:     PetscPrintf(PETSC_COMM_WORLD,"  Pe = %g\n",param->peclet);

    1049:     PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
    1050:     if (param->ivisc==VISC_CONST) {
    1051:       PetscPrintf(PETSC_COMM_WORLD,"                 Isoviscous \n");
    1052:       if (param->pv_analytic)
    1053:         PetscPrintf(PETSC_COMM_WORLD,"                          Pressure and Velocity prescribed! \n");
    1054:     } else if (param->ivisc==VISC_DIFN) {
    1055:       PetscPrintf(PETSC_COMM_WORLD,"                 Diffusion Creep (T-Dependent Newtonian) \n");
    1056:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
    1057:     } else if (param->ivisc==VISC_DISL ) {
    1058:       PetscPrintf(PETSC_COMM_WORLD,"                 Dislocation Creep (T-Dependent Non-Newtonian) \n");
    1059:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
    1060:     } else if (param->ivisc==VISC_FULL ) {
    1061:       PetscPrintf(PETSC_COMM_WORLD,"                 Full Rheology \n");
    1062:       PetscPrintf(PETSC_COMM_WORLD,"                          Viscosity range: %g--%g Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
    1063:     } else {
    1064:       PetscPrintf(PETSC_COMM_WORLD,"                 Invalid! \n");
    1065:       ierr_out=1;
    1066:     }

    1068:     PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
    1069:     if ( param->ibound==BC_ANALYTIC ) {
    1070:       PetscPrintf(PETSC_COMM_WORLD,"       Isoviscous Analytic Dirichlet \n");
    1071:     } else if ( param->ibound==BC_NOSTRESS ) {
    1072:       PetscPrintf(PETSC_COMM_WORLD,"       Stress-Free (normal & shear stress)\n");
    1073:     } else if ( param->ibound==BC_EXPERMNT ) {
    1074:       PetscPrintf(PETSC_COMM_WORLD,"       Experimental boundary condition \n");
    1075:     } else {
    1076:       PetscPrintf(PETSC_COMM_WORLD,"       Invalid! \n");
    1077:       ierr_out=1;
    1078:     }

    1080:     if (param->output_to_file)
    1081: #if defined(PETSC_HAVE_MATLAB)
    1082:       PetscPrintf(PETSC_COMM_WORLD,"Output Destination:       Mat file \"%s\"\n",param->filename);
    1083: #else
    1084:       PetscPrintf(PETSC_COMM_WORLD,"Output Destination:       PETSc binary file \"%s\"\n",param->filename);
    1085: #endif
    1086:     if ( param->output_ivisc != param->ivisc )
    1087:       PetscPrintf(PETSC_COMM_WORLD,"                          Output viscosity: -ivisc %D\n",param->output_ivisc);

    1089:     PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
    1090:   }
    1091:   if ( param->param_test ) PetscEnd();
    1092:   return ierr_out;
    1093: }

    1095: /* ------------------------------------------------------------------- */
    1098: /*  generates an inital guess using the analytic solution for isoviscous
    1099:     corner flow */
    1100: PetscErrorCode Initialize(DMMG *dmmg)
    1101: /* ------------------------------------------------------------------- */
    1102: {
    1103:   AppCtx         *user = (AppCtx*)dmmg[0]->user;
    1104:   Parameter      *param = user->param;
    1105:   GridInfo       *grid  = user->grid;
    1106:   DA             da;
    1107:   PetscInt       i,j,is,js,im,jm;
    1109:   Field          **x;

    1111:   /* Get the fine grid */
    1112:   da = (DA)(dmmg[0]->dm);
    1113:   DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
    1114:   DAVecGetArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);

    1116:   /* Compute initial guess */
    1117:   for (j=js; j<js+jm; j++) {
    1118:     for (i=is; i<is+im; i++) {
    1119:       if (i<j) {
    1120:         x[j][i].u  = param->cb;
    1121:       } else if (j<=grid->jlid) {
    1122:         x[j][i].u  = 0.0;
    1123:       } else {
    1124:         x[j][i].u  = HorizVelocity(i,j,user);
    1125:       }
    1126:       if (i<=j) {
    1127:         x[j][i].w = param->sb;
    1128:       } else if (j<=grid->jlid) {
    1129:         x[j][i].w = 0.0;
    1130:       } else {
    1131:         x[j][i].w = VertVelocity(i,j,user);
    1132:       }
    1133:       if (i<j || j<=grid->jlid) {
    1134:         x[j][i].p = 0.0;
    1135:       } else {
    1136:         x[j][i].p = Pressure(i,j,user);
    1137:       }
    1138:       x[j][i].T   = PetscMin(grid->dz*(j-0.5),1.0);
    1139:     }
    1140:   }

    1142:   /* Restore x to Xguess */
    1143:   DAVecRestoreArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);

    1145:   return 0;
    1146: }

    1148: /*---------------------------------------------------------------------*/
    1151: /*  controls output to a file */
    1152: PetscErrorCode DoOutput(DMMG *dmmg, PetscInt its)
    1153: /*---------------------------------------------------------------------*/
    1154: {
    1155:   AppCtx         *user = (AppCtx*)dmmg[0]->user;
    1156:   Parameter      *param = user->param;
    1157:   GridInfo       *grid  = user->grid;
    1159:   PetscMPIInt    rank;
    1160:   PetscInt       ivt=param->ivisc;
    1161:   PetscViewer    viewer;
    1162:   Vec            res, pars;
    1163:   MPI_Comm       comm;

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

    1167:   /* compute final residual and final viscosity/strain rate fields */
    1168:   SNESGetFunction(DMMGGetSNES(dmmg), &res, PETSC_NULL, PETSC_NULL);
    1169:   ViscosityField(DMMGGetDMMG(dmmg), DMMGGetx(dmmg), ((AppCtx *)dmmg[0]->user)->Xguess);

    1171:   /* get the communicator and the rank of the processor */
    1172:   PetscObjectGetComm((PetscObject)DMMGGetSNES(dmmg), &comm);
    1173:   MPI_Comm_rank(comm, &rank);

    1175:   if (param->output_to_file) { /* send output to binary file */
    1176:     VecCreate(comm, &pars);
    1177:     if (rank == 0) { /* on processor 0 */
    1178:       VecSetSizes(pars, 20, PETSC_DETERMINE);
    1179:       VecSetFromOptions(pars);
    1180:       VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
    1181:       VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
    1182:       VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
    1183:       VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
    1184:       VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
    1185:       VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
    1186:       /* skipped 6 intentionally */
    1187:       VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
    1188:       VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
    1189:       VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
    1190:       VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
    1191:       VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
    1192:       VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
    1193:       VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
    1194:       VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
    1195:       VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
    1196:       VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
    1197:     } else { /* on some other processor */
    1198:       VecSetSizes(pars, 0, PETSC_DETERMINE);
    1199:       VecSetFromOptions(pars);
    1200:     }
    1201:     VecAssemblyBegin(pars); VecAssemblyEnd(pars);

    1203:     /* create viewer */
    1204: #if defined(PETSC_HAVE_MATLAB)
    1205:     PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,PETSC_FILE_CREATE,&viewer);
    1206: #else
    1207:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,PETSC_FILE_CREATE,&viewer);
    1208: #endif

    1210:     /* send vectors to viewer */
    1211:     PetscObjectSetName((PetscObject)res,"res");
    1212:     VecView(res,viewer);
    1213:     PetscObjectSetName((PetscObject)DMMGGetx(dmmg),"out");
    1214:     VecView(DMMGGetx(dmmg), viewer);
    1215:     PetscObjectSetName((PetscObject)(user->Xguess),"aux");
    1216:     VecView(user->Xguess, viewer);
    1217:     StressField(dmmg); /* compute stress fields */
    1218:     PetscObjectSetName((PetscObject)(user->Xguess),"str");
    1219:     VecView(user->Xguess, viewer);
    1220:     PetscObjectSetName((PetscObject)pars,"par");
    1221:     VecView(pars, viewer);
    1222:
    1223:     /* destroy viewer and vector */
    1224:     PetscViewerDestroy(viewer);
    1225:     VecDestroy(pars);
    1226:   }
    1227:
    1228:   param->ivisc = ivt;
    1229:   return 0;
    1230: }

    1232: /* ------------------------------------------------------------------- */
    1235: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
    1236: PetscErrorCode ViscosityField(DMMG dmmg, Vec X, Vec V)
    1237: /* ------------------------------------------------------------------- */
    1238: {
    1239:   DA             da    = (DA) dmmg->dm;
    1240:   AppCtx         *user  = (AppCtx *) dmmg->user;
    1241:   Parameter      *param = user->param;
    1242:   GridInfo       *grid  = user->grid;
    1243:   Vec            localX;
    1244:   Field          **v, **x;
    1245:   PassiveReal    eps, dx, dz, T, epsC, TC;
    1246:   PetscInt       i,j,is,js,im,jm,ilim,jlim,ivt;

    1250:   ivt          = param->ivisc;
    1251:   param->ivisc = param->output_ivisc;

    1253:   DACreateLocalVector(da, &localX);
    1254:   DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
    1255:   DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
    1256:   DAVecGetArray(da,localX,(void**)&x);
    1257:   DAVecGetArray(da,V,(void**)&v);

    1259:   /* Parameters */
    1260:   dx   = grid->dx;   dz   = grid->dz;
    1261:   ilim = grid->ni-1; jlim = grid->nj-1;

    1263:   /* Compute real temperature, strain rate and viscosity */
    1264:   DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
    1265:   for (j=js; j<js+jm; j++) {
    1266:     for (i=is; i<is+im; i++) {
    1267:       T  = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*param->z_scale );
    1268:       if (i<ilim && j<jlim) {
    1269:         TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*param->z_scale );
    1270:       } else {
    1271:         TC = T;
    1272:       }
    1273:       eps  = CalcSecInv(x,i,j,CELL_CENTER,user);
    1274:       epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
    1275:       v[j][i].u = eps;
    1276:       v[j][i].w = epsC;
    1277:       v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
    1278:       v[j][i].T = Viscosity(TC,epsC,dz*j,param);
    1279:     }
    1280:   }
    1281:   DAVecRestoreArray(da,V,(void**)&v);
    1282:   DAVecRestoreArray(da,localX,(void**)&x);
    1283:   param->ivisc = ivt;
    1284:   return(0);
    1285: }

    1287: /* ------------------------------------------------------------------- */
    1290: /* post-processing: compute stress everywhere */
    1291: PetscErrorCode StressField(DMMG *dmmg)
    1292: /* ------------------------------------------------------------------- */
    1293: {
    1294:   AppCtx         *user = (AppCtx*)dmmg[0]->user;
    1295:   PetscInt       i,j,is,js,im,jm;
    1297:   DA             da;
    1298:   Vec            locVec;
    1299:   Field          **x, **y;

    1301:   /* Get the fine grid of Xguess and X */
    1302:   da = (DA)(dmmg[0]->dm);
    1303:   DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
    1304:   DAVecGetArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);

    1306:   DACreateLocalVector(da, &locVec);
    1307:   DAGlobalToLocalBegin(da, DMMGGetx(dmmg), INSERT_VALUES, locVec);
    1308:   DAGlobalToLocalEnd(da, DMMGGetx(dmmg), INSERT_VALUES, locVec);
    1309:   DAVecGetArray(da,locVec,(void**)&y);

    1311:   /* Compute stress on the corner points */
    1312:   for (j=js; j<js+jm; j++) {
    1313:      for (i=is; i<is+im; i++) {
    1314:
    1315:         x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
    1316:         x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
    1317:         x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
    1318:         x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
    1319:     }
    1320:   }

    1322:   /* Restore the fine grid of Xguess and X */
    1323:   DAVecRestoreArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);
    1324:   DAVecRestoreArray(da,locVec,(void**)&y);

    1326:   return 0;
    1327: }

    1329: /*=====================================================================
    1330:   UTILITY FUNCTIONS
    1331:   =====================================================================*/

    1333: /*---------------------------------------------------------------------*/
    1336: /* returns the velocity of the subducting slab and handles fault nodes
    1337:    for BC */
    1338: PassiveScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
    1339: /*---------------------------------------------------------------------*/
    1340: {
    1341:   Parameter     *param = user->param;
    1342:   GridInfo      *grid  = user->grid;

    1344:   if (c=='U' || c=='u') {
    1345:     if (i<j-1) {
    1346:       return param->cb;
    1347:     } else if (j<=grid->jfault) {
    1348:       return 0.0;
    1349:     } else
    1350:       return param->cb;

    1352:   } else {
    1353:     if (i<j) {
    1354:       return param->sb;
    1355:     } else if (j<=grid->jfault) {
    1356:       return 0.0;
    1357:     } else
    1358:       return param->sb;
    1359:   }
    1360: }

    1362: /*---------------------------------------------------------------------*/
    1365: /*  solution to diffusive half-space cooling model for BC */
    1366: PassiveScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
    1367: /*---------------------------------------------------------------------*/
    1368: {
    1369:   Parameter     *param = user->param;
    1370:   PassiveScalar z;
    1371:   if (plate==PLATE_LID)
    1372:     z = (j-0.5)*user->grid->dz;
    1373:   else /* PLATE_SLAB */
    1374:     z = (j-0.5)*user->grid->dz*param->cb;

    1376:   return erf(z*param->L/2.0/param->skt);
    1377: }

    1379: /*---------------------------------------------------------------------*/
    1382: PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
    1383: /*---------------------------------------------------------------------*/
    1384: {
    1385:   return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
    1386: }

    1388: /*---------------------------------------------------------------------*/
    1391: PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
    1392: /*---------------------------------------------------------------------*/
    1393: {
    1394:   return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
    1395: }

    1397: /*---------------------------------------------------------------------*/
    1400: PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
    1401: /*---------------------------------------------------------------------*/
    1402: {
    1403:   return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
    1404: }

    1406: /*---------------------------------------------------------------------*/
    1409: PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
    1410: /*---------------------------------------------------------------------*/
    1411: {
    1412:   return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
    1413: }

    1415: /*---------------------------------------------------------------------*/
    1418: /*  isoviscous analytic solution for IC */
    1419: PassiveScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
    1420: /*---------------------------------------------------------------------*/
    1421: {
    1422:   Parameter   *param = user->param;
    1423:   GridInfo    *grid  = user->grid;
    1424:   PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
    1425:
    1426:   x = (i - grid->jlid)*grid->dx;  z = (j - grid->jlid - 0.5)*grid->dz;
    1427:   r = sqrt(x*x+z*z); st = z/r;  ct = x/r;  th = atan(z/x);
    1428:   return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
    1429: }

    1431: /*---------------------------------------------------------------------*/
    1434: /*  isoviscous analytic solution for IC */
    1435: PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
    1436: /*---------------------------------------------------------------------*/
    1437: {
    1438:   Parameter   *param = user->param;
    1439:   GridInfo    *grid  = user->grid;
    1440:   PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
    1441:
    1442:   x = (i - grid->jlid - 0.5)*grid->dx;  z = (j - grid->jlid)*grid->dz;
    1443:   r = sqrt(x*x+z*z); st = z/r;  ct = x/r;  th = atan(z/x);
    1444:   return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
    1445: }

    1447: /*---------------------------------------------------------------------*/
    1450: /*  isoviscous analytic solution for IC */
    1451: PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
    1452: /*---------------------------------------------------------------------*/
    1453: {
    1454:   Parameter   *param = user->param;
    1455:   GridInfo    *grid  = user->grid;
    1456:   PetscScalar x, z, r, st, ct, c=param->c, d=param->d;

    1458:   x = (i - grid->jlid - 0.5)*grid->dx;  z = (j - grid->jlid - 0.5)*grid->dz;
    1459:   r = sqrt(x*x+z*z);  st = z/r;  ct = x/r;
    1460:   return (-2.0*(c*ct-d*st)/r);
    1461: }

    1463: /* ------------------------------------------------------------------- */
    1466: /*  utility function */
    1467: PetscTruth OptionsHasName(const char pre[],const char name[])
    1468: /* ------------------------------------------------------------------- */
    1469: {
    1470:   PetscTruth     retval;
    1472:   PetscOptionsHasName(pre,name,&retval);
    1473:   return retval;
    1474: }

    1476: /*=====================================================================
    1477:   INTERACTIVE SIGNAL HANDLING
    1478:   =====================================================================*/

    1480: /* ------------------------------------------------------------------- */
    1483: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
    1484: /* ------------------------------------------------------------------- */
    1485: {
    1486:   AppCtx        *user = (AppCtx *) ctx;
    1487:   Parameter     *param = user->param;
    1488:   KSP            ksp;

    1492:   if (param->interrupted) {
    1493:     param->interrupted = PETSC_FALSE;
    1494:     PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
    1495:     *reason = SNES_CONVERGED_FNORM_ABS;
    1496:     return(0);
    1497:   } else if (param->toggle_kspmon) {
    1498:     param->toggle_kspmon = PETSC_FALSE;
    1499:     SNESGetKSP(snes, &ksp);
    1500:     if (param->kspmon) {
    1501:       KSPClearMonitor(ksp);
    1502:       param->kspmon = PETSC_FALSE;
    1503:       PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
    1504:     } else {
    1505:       KSPSetMonitor(ksp,KSPSingularValueMonitor,PETSC_NULL,PETSC_NULL);
    1506:       param->kspmon = PETSC_TRUE;
    1507:       PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
    1508:     }
    1509:   }
    1510:   PetscFunctionReturn(SNESConverged_LS(snes, xnorm, pnorm, fnorm, reason, ctx));
    1511: }

    1513: /* ------------------------------------------------------------------- */
    1514: #include <signal.h>
    1517: PetscErrorCode InteractiveHandler(int signum, void *ctx)
    1518: /* ------------------------------------------------------------------- */
    1519: {
    1520:   AppCtx    *user = (AppCtx *) ctx;
    1521:   Parameter *param = user->param;

    1523:   if (signum == SIGILL) {
    1524:     param->toggle_kspmon = PETSC_TRUE;
    1525:   } else if (signum == SIGCONT) {
    1526:     param->interrupted = PETSC_TRUE;
    1527:   } else if (signum == SIGURG) {
    1528:     param->stop_solve = PETSC_TRUE;
    1529:   }
    1530:   return 0;
    1531: }