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(¶m,&grid);
175: ReportParams(¶m,&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 = ¶m;
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 = ¶m;
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: DMMGSetFromOptions(dmmg);
227: DMMGSetInitialGuess(dmmg,FormInitialGuess);
228: SNESSetConvergenceTest(DMMGGetSNES(dmmg),SNESConverged_Interactive,(void*)user,PETSC_NULL);
229: PetscPushSignalHandler(InteractiveHandler,(void*)user);
230:
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Initialize and solve the nonlinear system
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: Initialize(dmmg);
235: UpdateSolution(dmmg,user,&nits);
236:
237: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: Output variables.
239: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240: DoOutput(dmmg,nits);
241:
242: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: Free work space.
244: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
245: VecDestroy(user->Xguess);
246: PetscFree(user);
247: DMMGDestroy(dmmg);
248:
249: PetscFinalize();
250: return 0;
251: }
253: /*=====================================================================
254: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
255: =====================================================================*/
257: /*---------------------------------------------------------------------*/
260: /* manages solve: adaptive continuation method */
261: PetscErrorCode UpdateSolution(DMMG *dmmg, AppCtx *user, PetscInt *nits)
262: {
263: SNES snes;
264: KSP ksp;
265: PC pc;
266: SNESConvergedReason reason;
267: Parameter *param = user->param;
268: PassiveScalar cont_incr=0.3;
269: PetscInt its;
270: PetscErrorCode ierr;
271: PetscTruth q = PETSC_FALSE;
274: snes = DMMGGetSNES(dmmg);
275: SNESGetKSP(snes,&ksp);
276: KSPGetPC(ksp, &pc);
277: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
279: *nits=0;
281: /* Isoviscous solve */
282: if (param->ivisc == VISC_CONST && !param->stop_solve) {
283: param->ivisc = VISC_CONST;
284: DMMGSolve(dmmg);
285: VecCopy(DMMGGetx(dmmg),user->Xguess);
286: SNESGetIterationNumber(snes, &its);
287: *nits +=its;
288: if (param->stop_solve) goto done;
289: }
291: /* Olivine diffusion creep */
292: if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
293: if (!q) PetscPrintf(PETSC_COMM_WORLD,"Computing Variable Viscosity Solution\n");
295: /* continuation method on viscosity cutoff */
296: for (param->continuation=0.0; ; param->continuation+=cont_incr) {
297: if (!q) PetscPrintf(PETSC_COMM_WORLD," Continuation parameter = %G\n", param->continuation);
299: /* solve the non-linear system */
300: DMMGSolve(dmmg);
301: SNESGetConvergedReason(snes,&reason);
302: SNESGetIterationNumber(snes,&its);
303: *nits += its;
304: if (!q) PetscPrintf(PETSC_COMM_WORLD," Newton iterations: %D, Cumulative: %D\n", its, *nits);
305: if (param->stop_solve) goto done;
307: if (reason<0) {
308: /* NOT converged */
309: cont_incr = -fabs(cont_incr)/2.0;
310: if (fabs(cont_incr)<0.01) goto done;
312: } else {
313: /* converged */
314: VecCopy(DMMGGetx(dmmg),user->Xguess);
315: if (param->continuation >= 1.0) goto done;
316: if (its<=3) {
317: cont_incr = 0.30001;
318: } else if (its<=8) {
319: cont_incr = 0.15001;
320: } else {
321: cont_incr = 0.10001;
322: }
323: if (param->continuation+cont_incr > 1.0) {
324: cont_incr = 1.0 - param->continuation;
325: }
326: } /* endif reason<0 */
327: }
328: }
329: done:
330: if (param->stop_solve && !q) PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");
331: if (reason<0 && !q) PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");
332: return(0);
333: }
335: /* ------------------------------------------------------------------- */
338: /* used by SNESSolve to get an initial guess for the solution X */
339: PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X)
340: /* ------------------------------------------------------------------- */
341: {
342: AppCtx *user = (AppCtx*)dmmg->user;
345: VecCopy(user->Xguess, X);
346: return 0;
347: }
349: /*=====================================================================
350: PHYSICS FUNCTIONS (compute the discrete residual)
351: =====================================================================*/
353: /*---------------------------------------------------------------------*/
356: /* main call-back function that computes the processor-local piece
357: of the residual */
358: PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
359: /*---------------------------------------------------------------------*/
360: {
361: AppCtx *user = (AppCtx*)ptr;
362: Parameter *param = user->param;
363: GridInfo *grid = user->grid;
364: PetscScalar mag_w, mag_u;
365: PetscInt i,j,mx,mz,ilim,jlim;
366: PetscInt is,ie,js,je,ivisc,ibound;
370: /* Define global and local grid parameters */
371: mx = info->mx; mz = info->my;
372: ilim = mx-1; jlim = mz-1;
373: is = info->xs; ie = info->xs+info->xm;
374: js = info->ys; je = info->ys+info->ym;
376: /* Define geometric and numeric parameters */
377: ivisc = param->ivisc; ibound = param->ibound;
379: for (j=js; j<je; j++) {
380: for (i=is; i<ie; i++) {
382: /************* X-MOMENTUM/VELOCITY *************/
383: if (i<j) {
384: f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
386: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
387: /* in the lithospheric lid */
388: f[j][i].u = x[j][i].u - 0.0;
390: } else if (i==ilim) {
391: /* on the right side boundary */
392: if (ibound==BC_ANALYTIC) {
393: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
394: } else {
395: f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
396: }
398: } else if (j==jlim) {
399: /* on the bottom boundary */
400: if (ibound==BC_ANALYTIC) {
401: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
402: } else if (ibound==BC_NOSTRESS) {
403: f[j][i].u = XMomentumResidual(x,i,j,user);
404: } else {
405: /* experimental boundary condition */
406: }
408: } else {
409: /* in the mantle wedge */
410: f[j][i].u = XMomentumResidual(x,i,j,user);
411: }
412:
413: /************* Z-MOMENTUM/VELOCITY *************/
414: if (i<=j) {
415: f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);
417: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
418: /* in the lithospheric lid */
419: f[j][i].w = x[j][i].w - 0.0;
421: } else if (j==jlim) {
422: /* on the bottom boundary */
423: if (ibound==BC_ANALYTIC) {
424: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
425: } else {
426: f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
427: }
429: } else if (i==ilim) {
430: /* on the right side boundary */
431: if (ibound==BC_ANALYTIC) {
432: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
433: } else if (ibound==BC_NOSTRESS) {
434: f[j][i].w = ZMomentumResidual(x,i,j,user);
435: } else {
436: /* experimental boundary condition */
437: }
439: } else {
440: /* in the mantle wedge */
441: f[j][i].w = ZMomentumResidual(x,i,j,user);
442: }
444: /************* CONTINUITY/PRESSURE *************/
445: if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
446: /* in the lid or slab */
447: f[j][i].p = x[j][i].p;
448:
449: } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
450: /* on an analytic boundary */
451: f[j][i].p = x[j][i].p - Pressure(i,j,user);
453: } else {
454: /* in the mantle wedge */
455: f[j][i].p = ContinuityResidual(x,i,j,user);
456: }
458: /************* TEMPERATURE *************/
459: if (j==0) {
460: /* on the surface */
461: f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(x[j][i].T,0.0);
463: } else if (i==0) {
464: /* slab inflow boundary */
465: f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user);
467: } else if (i==ilim) {
468: /* right side boundary */
469: mag_u = 1.0 - pow( (1.0-PetscMax(PetscMin(x[j][i-1].u/param->cb,1.0),0.0)), 5.0 );
470: 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);
472: } else if (j==jlim) {
473: /* bottom boundary */
474: mag_w = 1.0 - pow( (1.0-PetscMax(PetscMin(x[j-1][i].w/param->sb,1.0),0.0)), 5.0 );
475: f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);
477: } else {
478: /* in the mantle wedge */
479: f[j][i].T = EnergyResidual(x,i,j,user);
480: }
481: }
482: }
483: return(0);
484: }
486: /*---------------------------------------------------------------------*/
489: /* computes the residual of the x-component of eqn (1) above */
490: PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
491: /*---------------------------------------------------------------------*/
492: {
493: Parameter *param=user->param;
494: GridInfo *grid =user->grid;
495: PetscScalar dx = grid->dx, dz=grid->dz;
496: PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
497: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
498: PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
499: PetscInt jlim = grid->nj-1;
501: z_scale = param->z_scale;
503: if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
504: TS = param->potentialT * TInterp(x,i,j-1) * exp( (j-1.0)*dz*z_scale );
505: if (j==jlim) TN = TS;
506: else TN = param->potentialT * TInterp(x,i,j) * exp( j *dz*z_scale );
507: TW = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
508: TE = param->potentialT * x[j][i+1].T * exp( (j-0.5)*dz*z_scale );
509: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
510: epsN = CalcSecInv(x,i,j, CELL_CORNER,user);
511: epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
512: epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
513: epsW = CalcSecInv(x,i,j, CELL_CENTER,user);
514: }
515: }
516: etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
517: etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
518: etaW = Viscosity(TW,epsW,dz*j,param);
519: etaE = Viscosity(TE,epsE,dz*j,param);
521: dPdx = ( x[j][i+1].p - x[j][i].p )/dx;
522: if (j==jlim) dudzN = etaN * ( x[j][i].w - x[j][i+1].w )/dx;
523: else dudzN = etaN * ( x[j+1][i].u - x[j][i].u )/dz;
524: dudzS = etaS * ( x[j][i].u - x[j-1][i].u )/dz;
525: dudxE = etaE * ( x[j][i+1].u - x[j][i].u )/dx;
526: dudxW = etaW * ( x[j][i].u - x[j][i-1].u )/dx;
528: residual = -dPdx /* X-MOMENTUM EQUATION*/
529: +( dudxE - dudxW )/dx
530: +( dudzN - dudzS )/dz;
532: if ( param->ivisc!=VISC_CONST ) {
533: dwdxN = etaN * ( x[j][i+1].w - x[j][i].w )/dx;
534: dwdxS = etaS * ( x[j-1][i+1].w - x[j-1][i].w )/dx;
536: residual += ( dudxE - dudxW )/dx + ( dwdxN - dwdxS )/dz;
537: }
539: return residual;
540: }
542: /*---------------------------------------------------------------------*/
545: /* computes the residual of the z-component of eqn (1) above */
546: PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
547: /*---------------------------------------------------------------------*/
548: {
549: Parameter *param=user->param;
550: GridInfo *grid =user->grid;
551: PetscScalar dx = grid->dx, dz=grid->dz;
552: 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;
553: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
554: PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
555: PetscInt ilim = grid->ni-1;
557: /* geometric and other parameters */
558: z_scale = param->z_scale;
559:
560: /* viscosity */
561: if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
562: TN = param->potentialT * x[j+1][i].T * exp( (j+0.5)*dz*z_scale );
563: TS = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
564: TW = param->potentialT * TInterp(x,i-1,j) * exp( j *dz*z_scale );
565: if (i==ilim) TE = TW;
566: else TE = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
567: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
568: epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
569: epsS = CalcSecInv(x,i,j, CELL_CENTER,user);
570: epsE = CalcSecInv(x,i,j, CELL_CORNER,user);
571: epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
572: }
573: }
574: etaN = Viscosity(TN,epsN,dz*(j+1),param);
575: etaS = Viscosity(TS,epsS,dz*j,param);
576: etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
577: etaE = Viscosity(TE,epsE,dz*(j+0.5),param);
579: dPdz = ( x[j+1][i].p - x[j][i].p )/dz;
580: dwdzN = etaN * ( x[j+1][i].w - x[j][i].w )/dz;
581: dwdzS = etaS * ( x[j][i].w - x[j-1][i].w )/dz;
582: if (i==ilim) dwdxE = etaE * ( x[j][i].u - x[j+1][i].u )/dz;
583: else dwdxE = etaE * ( x[j][i+1].w - x[j][i].w )/dx;
584: dwdxW = 2.0*etaW * ( x[j][i].w - x[j][i-1].w )/dx;
585:
586: /* Z-MOMENTUM */
587: residual = -dPdz /* constant viscosity terms */
588: +( dwdzN - dwdzS )/dz
589: +( dwdxE - dwdxW )/dx;
591: if ( param->ivisc!=VISC_CONST ) {
592: dudzE = etaE * ( x[j+1][i].u - x[j][i].u )/dz;
593: dudzW = etaW * ( x[j+1][i-1].u - x[j][i-1].u )/dz;
595: residual += ( dwdzN - dwdzS )/dz + ( dudzE - dudzW )/dx;
596: }
598: return residual;
599: }
601: /*---------------------------------------------------------------------*/
604: /* computes the residual of eqn (2) above */
605: PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
606: /*---------------------------------------------------------------------*/
607: {
608: GridInfo *grid =user->grid;
609: PetscScalar uE,uW,wN,wS,dudx,dwdz;
611: uW = x[j][i-1].u; uE = x[j][i].u; dudx = ( uE - uW )/grid->dx;
612: wS = x[j-1][i].w; wN = x[j][i].w; dwdz = ( wN - wS )/grid->dz;
614: return dudx + dwdz;
615: }
617: /*---------------------------------------------------------------------*/
620: /* computes the residual of eqn (3) above */
621: PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
622: /*---------------------------------------------------------------------*/
623: {
624: Parameter *param=user->param;
625: GridInfo *grid =user->grid;
626: PetscScalar dx = grid->dx, dz=grid->dz;
627: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
628: PetscScalar TE, TN, TS, TW, residual;
629: PetscScalar uE,uW,wN,wS;
630: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;
632: dTdzN = ( x[j+1][i].T - x[j][i].T )/dz;
633: dTdzS = ( x[j][i].T - x[j-1][i].T )/dz;
634: dTdxE = ( x[j][i+1].T - x[j][i].T )/dx;
635: dTdxW = ( x[j][i].T - x[j][i-1].T )/dx;
636:
637: residual = ( ( dTdzN - dTdzS )/dz + /* diffusion term */
638: ( dTdxE - dTdxW )/dx )*dx*dz/param->peclet;
640: if (j<=jlid && i>=j) {
641: /* don't advect in the lid */
642: return residual;
644: } else if (i<j) {
645: /* beneath the slab sfc */
646: uW = uE = param->cb;
647: wS = wN = param->sb;
649: } else {
650: /* advect in the slab and wedge */
651: uW = x[j][i-1].u; uE = x[j][i].u;
652: wS = x[j-1][i].w; wN = x[j][i].w;
653: }
655: if ( param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1 ) {
656: /* finite volume advection */
657: TS = ( x[j][i].T + x[j-1][i].T )/2.0;
658: TN = ( x[j][i].T + x[j+1][i].T )/2.0;
659: TE = ( x[j][i].T + x[j][i+1].T )/2.0;
660: TW = ( x[j][i].T + x[j][i-1].T )/2.0;
661: fN = wN*TN*dx; fS = wS*TS*dx;
662: fE = uE*TE*dz; fW = uW*TW*dz;
663:
664: } else {
665: /* Fromm advection scheme */
666: 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
667: - 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;
668: 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
669: - 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;
670: 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
671: - 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;
672: 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
673: - 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;
674: }
675:
676: residual -= ( fE - fW + fN - fS );
678: return residual;
679: }
681: /*---------------------------------------------------------------------*/
684: /* computes the shear stress---used on the boundaries */
685: PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
686: /*---------------------------------------------------------------------*/
687: {
688: Parameter *param=user->param;
689: GridInfo *grid=user->grid;
690: PetscInt ilim=grid->ni-1, jlim=grid->nj-1;
691: PetscScalar uN, uS, wE, wW;
693: if ( j<=grid->jlid || i<j || i==ilim || j==jlim ) return EPS_ZERO;
695: if (ipos==CELL_CENTER) { /* on cell center */
697: wE = WInterp(x,i,j-1);
698: if (i==j) { wW = param->sb; uN = param->cb;}
699: else { wW = WInterp(x,i-1,j-1); uN = UInterp(x,i-1,j); }
700: if (j==grid->jlid+1) uS = 0.0;
701: else uS = UInterp(x,i-1,j-1);
703: } else { /* on cell corner */
705: uN = x[j+1][i].u; uS = x[j][i].u;
706: wW = x[j][i].w; wE = x[j][i+1].w;
708: }
710: return (uN-uS)/grid->dz + (wE-wW)/grid->dx;
711: }
713: /*---------------------------------------------------------------------*/
716: /* computes the normal stress---used on the boundaries */
717: PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
718: /*---------------------------------------------------------------------*/
719: {
720: Parameter *param=user->param;
721: GridInfo *grid =user->grid;
722: PetscScalar dx = grid->dx, dz=grid->dz;
723: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, ivisc;
724: PetscScalar epsC=0.0, etaC, TC, uE, uW, pC, z_scale;
725: if (i<j || j<=grid->jlid) return EPS_ZERO;
727: ivisc=param->ivisc; z_scale = param->z_scale;
729: if (ipos==CELL_CENTER) { /* on cell center */
731: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
732: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
733: etaC = Viscosity(TC,epsC,dz*j,param);
735: uW = x[j][i-1].u; uE = x[j][i].u;
736: pC = x[j][i].p;
738: } else { /* on cell corner */
739: if ( i==ilim || j==jlim ) return EPS_ZERO;
741: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
742: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
743: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
745: if (i==j) uW = param->sb;
746: else uW = UInterp(x,i-1,j);
747: uE = UInterp(x,i,j); pC = PInterp(x,i,j);
748: }
749:
750: return 2.0*etaC*(uE-uW)/dx - pC;
751: }
753: /*---------------------------------------------------------------------*/
756: /* computes the normal stress---used on the boundaries */
757: PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
758: /*---------------------------------------------------------------------*/
759: {
760: Parameter *param=user->param;
761: GridInfo *grid =user->grid;
762: PetscScalar dz=grid->dz;
763: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, ivisc;
764: PetscScalar epsC=0.0, etaC, TC;
765: PetscScalar pC, wN, wS, z_scale;
766: if (i<j || j<=grid->jlid) return EPS_ZERO;
768: ivisc=param->ivisc; z_scale = param->z_scale;
770: if (ipos==CELL_CENTER) { /* on cell center */
772: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
773: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
774: etaC = Viscosity(TC,epsC,dz*j,param);
775: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
777: } else { /* on cell corner */
778: if ( (i==ilim) || (j==jlim) ) return EPS_ZERO;
780: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
781: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
782: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
783: if (i==j) wN = param->sb;
784: else wN = WInterp(x,i,j);
785: wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
786: }
788: return 2.0*etaC*(wN-wS)/dz - pC;
789: }
791: /*---------------------------------------------------------------------*/
794: /* computes the second invariant of the strain rate tensor */
795: PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
796: /*---------------------------------------------------------------------*/
797: {
798: Parameter *param = user->param;
799: GridInfo *grid = user->grid;
800: PetscInt ilim=grid->ni-1, jlim=grid->nj-1;
801: PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
802: PetscScalar eps11, eps12, eps22;
804: if (i<j) return EPS_ZERO;
805: if (i==ilim) i--; if (j==jlim) j--;
807: if (ipos==CELL_CENTER) { /* on cell center */
808: if (j<=grid->jlid) return EPS_ZERO;
810: uE = x[j][i].u; uW = x[j][i-1].u;
811: wN = x[j][i].w; wS = x[j-1][i].w;
812: wE = WInterp(x,i,j-1);
813: if (i==j) { uN = param->cb; wW = param->sb; }
814: else { uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1); }
816: if (j==grid->jlid+1) uS = 0.0;
817: else uS = UInterp(x,i-1,j-1);
819: } else { /* on CELL_CORNER */
820: if (j<grid->jlid) return EPS_ZERO;
822: uN = x[j+1][i].u; uS = x[j][i].u;
823: wE = x[j][i+1].w; wW = x[j][i].w;
824: if (i==j) { wN = param->sb; uW = param->cb; }
825: else { wN = WInterp(x,i,j); uW = UInterp(x,i-1,j); }
827: if (j==grid->jlid) {
828: uE = 0.0; uW = 0.0;
829: uS = -uN;
830: wS = -wN;
831: } else {
832: uE = UInterp(x,i,j);
833: wS = WInterp(x,i,j-1);
834: }
835: }
837: eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz;
838: eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx);
840: return sqrt( 0.5*( eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22 ) );
841: }
843: /*---------------------------------------------------------------------*/
846: /* computes the shear viscosity */
847: PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PassiveScalar z,
848: Parameter *param)
849: /*---------------------------------------------------------------------*/
850: {
851: PetscScalar result=0.0;
852: ViscParam difn=param->diffusion, disl=param->dislocation;
853: PetscInt iVisc=param->ivisc;
854: double eps_scale=param->V/(param->L*1000.0);
855: double strain_power, v1, v2, P;
856: double rho_g = 32340.0, R=8.3144;
858: P = rho_g*(z*param->L*1000.0); /* Pa */
860: if (iVisc==VISC_CONST) {
861: /* constant viscosity */
862: return 1.0;
864: } else if (iVisc==VISC_DIFN) {
865: /* diffusion creep rheology */
866: result = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
868: } else if (iVisc==VISC_DISL) {
869: /* dislocation creep rheology */
870: strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
871: result = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
873: } else if (iVisc==VISC_FULL) {
874: /* dislocation/diffusion creep rheology */
875: strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
876: v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
877: v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
878: result = 1.0/(1.0/v1 + 1.0/v2);
879: }
881: /* max viscosity is param->eta0 */
882: result = PetscMin( result, 1.0 );
883: /* min viscosity is param->visc_cutoff */
884: result = PetscMax( result, param->visc_cutoff );
885: /* continuation method */
886: result = pow(result,param->continuation);
887: return result;
888: }
890: /*=====================================================================
891: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
892: =====================================================================*/
894: /*---------------------------------------------------------------------*/
897: /* initializes the problem parameters and checks for
898: command line changes */
899: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
900: /*---------------------------------------------------------------------*/
901: {
902: PetscErrorCode ierr, ierr_out=0;
903: PetscReal SEC_PER_YR = 3600.00*24.00*365.2500;
904: PetscReal PI = 3.14159265358979323846;
905: PetscReal alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
906:
907: /* domain geometry */
908: param->slab_dip = 45.0;
909: param->width = 320.0; /* km */
910: param->depth = 300.0; /* km */
911: param->lid_depth = 35.0; /* km */
912: param->fault_depth = 35.0; /* km */
913: PetscOptionsGetReal(PETSC_NULL,"-slab_dip",&(param->slab_dip),PETSC_NULL);
914: PetscOptionsGetReal(PETSC_NULL,"-width",&(param->width),PETSC_NULL);
915: PetscOptionsGetReal(PETSC_NULL,"-depth",&(param->depth),PETSC_NULL);
916: PetscOptionsGetReal(PETSC_NULL,"-lid_depth",&(param->lid_depth),PETSC_NULL);
917: PetscOptionsGetReal(PETSC_NULL,"-fault_depth",&(param->fault_depth),PETSC_NULL);
918: param->slab_dip = param->slab_dip*PI/180.0; /* radians */
920: /* grid information */
921: PetscOptionsGetInt(PETSC_NULL, "-jfault",&(grid->jfault),PETSC_NULL);
922: grid->ni = 82;
923: PetscOptionsGetInt(PETSC_NULL, "-ni",&(grid->ni),PETSC_NULL);
924: grid->dx = param->width/((double)(grid->ni-2)); /* km */
925: grid->dz = grid->dx*tan(param->slab_dip); /* km */
926: grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/
927: param->depth = grid->dz*(grid->nj-2); /* km */
928: grid->inose = 0; /* gridpoints*/
929: PetscOptionsGetInt(PETSC_NULL,"-inose",&(grid->inose),PETSC_NULL);
930: grid->periodic = DA_NONPERIODIC;
931: grid->stencil = DA_STENCIL_BOX;
932: grid->dof = 4;
933: grid->stencil_width = 2;
934: grid->mglevels = 1;
936: /* boundary conditions */
937: param->pv_analytic = PETSC_FALSE;
938: param->ibound = BC_NOSTRESS;
939: PetscOptionsGetInt(PETSC_NULL,"-ibound",&(param->ibound),PETSC_NULL);
941: /* physical constants */
942: param->slab_velocity = 5.0; /* cm/yr */
943: param->slab_age = 50.0; /* Ma */
944: param->lid_age = 50.0; /* Ma */
945: param->kappa = 0.7272e-6; /* m^2/sec */
946: param->potentialT = 1300.0; /* degrees C */
947: PetscOptionsGetReal(PETSC_NULL,"-slab_velocity",&(param->slab_velocity),PETSC_NULL);
948: PetscOptionsGetReal(PETSC_NULL,"-slab_age",&(param->slab_age),PETSC_NULL);
949: PetscOptionsGetReal(PETSC_NULL,"-lid_age",&(param->lid_age),PETSC_NULL);
950: PetscOptionsGetReal(PETSC_NULL,"-kappa",&(param->kappa),PETSC_NULL);
951: PetscOptionsGetReal(PETSC_NULL,"-potentialT",&(param->potentialT),PETSC_NULL);
953: /* viscosity */
954: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
955: param->eta0 = 1e24; /* Pa-s */
956: param->visc_cutoff = 0.0; /* factor of eta_0 */
957: param->continuation = 1.0;
958: /* constants for diffusion creep */
959: param->diffusion.A = 1.8e7; /* Pa-s */
960: param->diffusion.n = 1.0; /* dim'less */
961: param->diffusion.Estar = 375e3; /* J/mol */
962: param->diffusion.Vstar = 5e-6; /* m^3/mol */
963: /* constants for param->dislocationocation creep */
964: param->dislocation.A = 2.8969e4; /* Pa-s */
965: param->dislocation.n = 3.5; /* dim'less */
966: param->dislocation.Estar = 530e3; /* J/mol */
967: param->dislocation.Vstar = 14e-6; /* m^3/mol */
968: PetscOptionsGetInt(PETSC_NULL, "-ivisc",&(param->ivisc),PETSC_NULL);
969: PetscOptionsGetReal(PETSC_NULL,"-visc_cutoff",&(param->visc_cutoff),PETSC_NULL);
970: param->output_ivisc = param->ivisc;
971: PetscOptionsGetInt(PETSC_NULL,"-output_ivisc",&(param->output_ivisc),PETSC_NULL);
972: PetscOptionsGetReal(PETSC_NULL,"-vstar",&(param->dislocation.Vstar),PETSC_NULL);
974: /* output options */
975: param->quiet = PETSC_FALSE;
976: param->param_test = PETSC_FALSE;
977: PetscOptionsHasName(PETSC_NULL,"-quiet",&(param->quiet));
978: PetscOptionsHasName(PETSC_NULL,"-test",&(param->param_test));
979: PetscOptionsGetString(PETSC_NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));
981: /* advection */
982: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
983: PetscOptionsGetInt(PETSC_NULL,"-adv_scheme",&(param->adv_scheme),PETSC_NULL);
985: /* misc. flags */
986: param->stop_solve = PETSC_FALSE;
987: param->interrupted = PETSC_FALSE;
988: param->kspmon = PETSC_FALSE;
989: param->toggle_kspmon = PETSC_FALSE;
991: /* derived parameters for slab angle */
992: param->sb = sin(param->slab_dip);
993: param->cb = cos(param->slab_dip);
994: param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
995: param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);
997: /* length, velocity and time scale for non-dimensionalization */
998: param->L = PetscMin(param->width,param->depth); /* km */
999: param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */
1001: /* other unit conversions and derived parameters */
1002: param->scaled_width = param->width/param->L; /* dim'less */
1003: param->scaled_depth = param->depth/param->L; /* dim'less */
1004: param->lid_depth = param->lid_depth/param->L; /* dim'less */
1005: param->fault_depth = param->fault_depth/param->L; /* dim'less */
1006: grid->dx = grid->dx/param->L; /* dim'less */
1007: grid->dz = grid->dz/param->L; /* dim'less */
1008: grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */
1009: grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */
1010: param->lid_depth = grid->jlid*grid->dz; /* dim'less */
1011: param->fault_depth = grid->jfault*grid->dz; /* dim'less */
1012: grid->corner = grid->jlid+1; /* gridcells */
1013: param->peclet = param->V /* m/sec */
1014: * param->L*1000.0 /* m */
1015: / param->kappa; /* m^2/sec */
1016: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
1017: param->skt = sqrt(param->kappa*param->slab_age*SEC_PER_YR);
1018: PetscOptionsGetReal(PETSC_NULL,"-peclet",&(param->peclet),PETSC_NULL);
1019:
1020: return ierr_out;
1021: }
1023: /*---------------------------------------------------------------------*/
1026: /* prints a report of the problem parameters to stdout */
1027: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
1028: /*---------------------------------------------------------------------*/
1029: {
1030: PetscErrorCode ierr, ierr_out=0;
1031: char date[30];
1032: PetscReal PI = 3.14159265358979323846;
1034: PetscGetDate(date,30);
1036: if ( !(param->quiet) ) {
1037: PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
1038: PetscPrintf(PETSC_COMM_WORLD," %s\n",&(date[0]));
1040: PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
1041: PetscPrintf(PETSC_COMM_WORLD," Width = %G km, Depth = %G km\n",param->width,param->depth);
1042: PetscPrintf(PETSC_COMM_WORLD," Slab dip = %G degrees, Slab velocity = %G cm/yr\n",param->slab_dip*180.0/PI,param->slab_velocity);
1043: 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);
1045: PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
1046: 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);
1047: PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);
1048: PetscPrintf(PETSC_COMM_WORLD," Pe = %G\n",param->peclet);
1050: PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
1051: if (param->ivisc==VISC_CONST) {
1052: PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");
1053: if (param->pv_analytic)
1054: PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");
1055: } else if (param->ivisc==VISC_DIFN) {
1056: PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");
1057: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
1058: } else if (param->ivisc==VISC_DISL ) {
1059: PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");
1060: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
1061: } else if (param->ivisc==VISC_FULL ) {
1062: PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");
1063: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
1064: } else {
1065: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
1066: ierr_out=1;
1067: }
1069: PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
1070: if ( param->ibound==BC_ANALYTIC ) {
1071: PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");
1072: } else if ( param->ibound==BC_NOSTRESS ) {
1073: PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");
1074: } else if ( param->ibound==BC_EXPERMNT ) {
1075: PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");
1076: } else {
1077: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
1078: ierr_out=1;
1079: }
1081: if (param->output_to_file)
1082: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1083: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);
1084: #else
1085: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);
1086: #endif
1087: if ( param->output_ivisc != param->ivisc )
1088: PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);
1090: PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
1091: }
1092: if ( param->param_test ) PetscEnd();
1093: return ierr_out;
1094: }
1096: /* ------------------------------------------------------------------- */
1099: /* generates an inital guess using the analytic solution for isoviscous
1100: corner flow */
1101: PetscErrorCode Initialize(DMMG *dmmg)
1102: /* ------------------------------------------------------------------- */
1103: {
1104: AppCtx *user = (AppCtx*)dmmg[0]->user;
1105: Parameter *param = user->param;
1106: GridInfo *grid = user->grid;
1107: DA da;
1108: PetscInt i,j,is,js,im,jm;
1110: Field **x;
1112: /* Get the fine grid */
1113: da = (DA)(dmmg[0]->dm);
1114: DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1115: DAVecGetArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);
1117: /* Compute initial guess */
1118: for (j=js; j<js+jm; j++) {
1119: for (i=is; i<is+im; i++) {
1120: if (i<j) {
1121: x[j][i].u = param->cb;
1122: } else if (j<=grid->jlid) {
1123: x[j][i].u = 0.0;
1124: } else {
1125: x[j][i].u = HorizVelocity(i,j,user);
1126: }
1127: if (i<=j) {
1128: x[j][i].w = param->sb;
1129: } else if (j<=grid->jlid) {
1130: x[j][i].w = 0.0;
1131: } else {
1132: x[j][i].w = VertVelocity(i,j,user);
1133: }
1134: if (i<j || j<=grid->jlid) {
1135: x[j][i].p = 0.0;
1136: } else {
1137: x[j][i].p = Pressure(i,j,user);
1138: }
1139: x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1140: }
1141: }
1143: /* Restore x to Xguess */
1144: DAVecRestoreArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);
1146: return 0;
1147: }
1149: /*---------------------------------------------------------------------*/
1152: /* controls output to a file */
1153: PetscErrorCode DoOutput(DMMG *dmmg, PetscInt its)
1154: /*---------------------------------------------------------------------*/
1155: {
1156: AppCtx *user = (AppCtx*)dmmg[0]->user;
1157: Parameter *param = user->param;
1158: GridInfo *grid = user->grid;
1160: PetscMPIInt rank;
1161: PetscInt ivt=param->ivisc;
1162: PetscViewer viewer;
1163: Vec res, pars;
1164: MPI_Comm comm;
1166: param->ivisc = param->output_ivisc;
1168: /* compute final residual and final viscosity/strain rate fields */
1169: SNESGetFunction(DMMGGetSNES(dmmg), &res, PETSC_NULL, PETSC_NULL);
1170: ViscosityField(DMMGGetDMMG(dmmg), DMMGGetx(dmmg), ((AppCtx *)dmmg[0]->user)->Xguess);
1172: /* get the communicator and the rank of the processor */
1173: PetscObjectGetComm((PetscObject)DMMGGetSNES(dmmg), &comm);
1174: MPI_Comm_rank(comm, &rank);
1176: if (param->output_to_file) { /* send output to binary file */
1177: VecCreate(comm, &pars);
1178: if (rank == 0) { /* on processor 0 */
1179: VecSetSizes(pars, 20, PETSC_DETERMINE);
1180: VecSetFromOptions(pars);
1181: VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1182: VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1183: VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1184: VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1185: VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1186: VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1187: /* skipped 6 intentionally */
1188: VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1189: VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1190: VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1191: VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1192: VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1193: VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1194: VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1195: VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1196: VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1197: VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1198: } else { /* on some other processor */
1199: VecSetSizes(pars, 0, PETSC_DETERMINE);
1200: VecSetFromOptions(pars);
1201: }
1202: VecAssemblyBegin(pars); VecAssemblyEnd(pars);
1204: /* create viewer */
1205: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1206: PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1207: #else
1208: PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1209: #endif
1211: /* send vectors to viewer */
1212: PetscObjectSetName((PetscObject)res,"res");
1213: VecView(res,viewer);
1214: PetscObjectSetName((PetscObject)DMMGGetx(dmmg),"out");
1215: VecView(DMMGGetx(dmmg), viewer);
1216: PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1217: VecView(user->Xguess, viewer);
1218: StressField(dmmg); /* compute stress fields */
1219: PetscObjectSetName((PetscObject)(user->Xguess),"str");
1220: VecView(user->Xguess, viewer);
1221: PetscObjectSetName((PetscObject)pars,"par");
1222: VecView(pars, viewer);
1223:
1224: /* destroy viewer and vector */
1225: PetscViewerDestroy(viewer);
1226: VecDestroy(pars);
1227: }
1228:
1229: param->ivisc = ivt;
1230: return 0;
1231: }
1233: /* ------------------------------------------------------------------- */
1236: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1237: PetscErrorCode ViscosityField(DMMG dmmg, Vec X, Vec V)
1238: /* ------------------------------------------------------------------- */
1239: {
1240: DA da = (DA) dmmg->dm;
1241: AppCtx *user = (AppCtx *) dmmg->user;
1242: Parameter *param = user->param;
1243: GridInfo *grid = user->grid;
1244: Vec localX;
1245: Field **v, **x;
1246: PassiveReal eps, dx, dz, T, epsC, TC;
1247: PetscInt i,j,is,js,im,jm,ilim,jlim,ivt;
1251: ivt = param->ivisc;
1252: param->ivisc = param->output_ivisc;
1254: DACreateLocalVector(da, &localX);
1255: DAGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1256: DAGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1257: DAVecGetArray(da,localX,(void**)&x);
1258: DAVecGetArray(da,V,(void**)&v);
1260: /* Parameters */
1261: dx = grid->dx; dz = grid->dz;
1262: ilim = grid->ni-1; jlim = grid->nj-1;
1264: /* Compute real temperature, strain rate and viscosity */
1265: DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1266: for (j=js; j<js+jm; j++) {
1267: for (i=is; i<is+im; i++) {
1268: T = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*param->z_scale );
1269: if (i<ilim && j<jlim) {
1270: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*param->z_scale );
1271: } else {
1272: TC = T;
1273: }
1274: eps = CalcSecInv(x,i,j,CELL_CENTER,user);
1275: epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
1276: v[j][i].u = eps;
1277: v[j][i].w = epsC;
1278: v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1279: v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1280: }
1281: }
1282: DAVecRestoreArray(da,V,(void**)&v);
1283: DAVecRestoreArray(da,localX,(void**)&x);
1284: param->ivisc = ivt;
1285: return(0);
1286: }
1288: /* ------------------------------------------------------------------- */
1291: /* post-processing: compute stress everywhere */
1292: PetscErrorCode StressField(DMMG *dmmg)
1293: /* ------------------------------------------------------------------- */
1294: {
1295: AppCtx *user = (AppCtx*)dmmg[0]->user;
1296: PetscInt i,j,is,js,im,jm;
1298: DA da;
1299: Vec locVec;
1300: Field **x, **y;
1302: /* Get the fine grid of Xguess and X */
1303: da = (DA)(dmmg[0]->dm);
1304: DAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1305: DAVecGetArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);
1307: DACreateLocalVector(da, &locVec);
1308: DAGlobalToLocalBegin(da, DMMGGetx(dmmg), INSERT_VALUES, locVec);
1309: DAGlobalToLocalEnd(da, DMMGGetx(dmmg), INSERT_VALUES, locVec);
1310: DAVecGetArray(da,locVec,(void**)&y);
1312: /* Compute stress on the corner points */
1313: for (j=js; j<js+jm; j++) {
1314: for (i=is; i<is+im; i++) {
1315:
1316: x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1317: x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1318: x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1319: x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1320: }
1321: }
1323: /* Restore the fine grid of Xguess and X */
1324: DAVecRestoreArray(da,((AppCtx*)dmmg[0]->user)->Xguess,(void**)&x);
1325: DAVecRestoreArray(da,locVec,(void**)&y);
1327: return 0;
1328: }
1330: /*=====================================================================
1331: UTILITY FUNCTIONS
1332: =====================================================================*/
1334: /*---------------------------------------------------------------------*/
1337: /* returns the velocity of the subducting slab and handles fault nodes
1338: for BC */
1339: PassiveScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1340: /*---------------------------------------------------------------------*/
1341: {
1342: Parameter *param = user->param;
1343: GridInfo *grid = user->grid;
1345: if (c=='U' || c=='u') {
1346: if (i<j-1) {
1347: return param->cb;
1348: } else if (j<=grid->jfault) {
1349: return 0.0;
1350: } else
1351: return param->cb;
1353: } else {
1354: if (i<j) {
1355: return param->sb;
1356: } else if (j<=grid->jfault) {
1357: return 0.0;
1358: } else
1359: return param->sb;
1360: }
1361: }
1363: /*---------------------------------------------------------------------*/
1366: /* solution to diffusive half-space cooling model for BC */
1367: PassiveScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1368: /*---------------------------------------------------------------------*/
1369: {
1370: Parameter *param = user->param;
1371: PassiveScalar z;
1372: if (plate==PLATE_LID)
1373: z = (j-0.5)*user->grid->dz;
1374: else /* PLATE_SLAB */
1375: z = (j-0.5)*user->grid->dz*param->cb;
1376: #if defined (PETSC_HAVE_ERF)
1377: return erf(z*param->L/2.0/param->skt);
1378: #else
1379: SETERRQ(1,"erf() not available on this machine");
1380: #endif
1381: }
1383: /*---------------------------------------------------------------------*/
1386: PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
1387: /*---------------------------------------------------------------------*/
1388: {
1389: return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
1390: }
1392: /*---------------------------------------------------------------------*/
1395: PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
1396: /*---------------------------------------------------------------------*/
1397: {
1398: return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
1399: }
1401: /*---------------------------------------------------------------------*/
1404: PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
1405: /*---------------------------------------------------------------------*/
1406: {
1407: return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
1408: }
1410: /*---------------------------------------------------------------------*/
1413: PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
1414: /*---------------------------------------------------------------------*/
1415: {
1416: return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
1417: }
1419: /*---------------------------------------------------------------------*/
1422: /* isoviscous analytic solution for IC */
1423: PassiveScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
1424: /*---------------------------------------------------------------------*/
1425: {
1426: Parameter *param = user->param;
1427: GridInfo *grid = user->grid;
1428: PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
1429:
1430: x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
1431: r = sqrt(x*x+z*z); st = z/r; ct = x/r; th = atan(z/x);
1432: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
1433: }
1435: /*---------------------------------------------------------------------*/
1438: /* isoviscous analytic solution for IC */
1439: PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
1440: /*---------------------------------------------------------------------*/
1441: {
1442: Parameter *param = user->param;
1443: GridInfo *grid = user->grid;
1444: PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
1445:
1446: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz;
1447: r = sqrt(x*x+z*z); st = z/r; ct = x/r; th = atan(z/x);
1448: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
1449: }
1451: /*---------------------------------------------------------------------*/
1454: /* isoviscous analytic solution for IC */
1455: PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
1456: /*---------------------------------------------------------------------*/
1457: {
1458: Parameter *param = user->param;
1459: GridInfo *grid = user->grid;
1460: PetscScalar x, z, r, st, ct, c=param->c, d=param->d;
1462: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
1463: r = sqrt(x*x+z*z); st = z/r; ct = x/r;
1464: return (-2.0*(c*ct-d*st)/r);
1465: }
1467: /* ------------------------------------------------------------------- */
1470: /* utility function */
1471: PetscTruth OptionsHasName(const char pre[],const char name[])
1472: /* ------------------------------------------------------------------- */
1473: {
1474: PetscTruth retval;
1476: PetscOptionsHasName(pre,name,&retval);
1477: return retval;
1478: }
1480: /*=====================================================================
1481: INTERACTIVE SIGNAL HANDLING
1482: =====================================================================*/
1484: /* ------------------------------------------------------------------- */
1487: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal pnorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1488: /* ------------------------------------------------------------------- */
1489: {
1490: AppCtx *user = (AppCtx *) ctx;
1491: Parameter *param = user->param;
1492: KSP ksp;
1496: if (param->interrupted) {
1497: param->interrupted = PETSC_FALSE;
1498: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1499: *reason = SNES_CONVERGED_FNORM_ABS;
1500: return(0);
1501: } else if (param->toggle_kspmon) {
1502: param->toggle_kspmon = PETSC_FALSE;
1503: SNESGetKSP(snes, &ksp);
1504: if (param->kspmon) {
1505: KSPMonitorCancel(ksp);
1506: param->kspmon = PETSC_FALSE;
1507: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1508: } else {
1509: KSPMonitorSet(ksp,KSPMonitorSingularValue,PETSC_NULL,PETSC_NULL);
1510: param->kspmon = PETSC_TRUE;
1511: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1512: }
1513: }
1514: PetscFunctionReturn(SNESDefaultConverged(snes,it,xnorm,pnorm,fnorm,reason,ctx));
1515: }
1517: /* ------------------------------------------------------------------- */
1518: #include <signal.h>
1521: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1522: /* ------------------------------------------------------------------- */
1523: {
1524: AppCtx *user = (AppCtx *) ctx;
1525: Parameter *param = user->param;
1527: if (signum == SIGILL) {
1528: param->toggle_kspmon = PETSC_TRUE;
1529: } else if (signum == SIGCONT) {
1530: param->interrupted = PETSC_TRUE;
1531: } else if (signum == SIGURG) {
1532: param->stop_solve = PETSC_TRUE;
1533: }
1534: return 0;
1535: }