Actual source code: ex30.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\
2: The flow is driven by the subducting slab.\n\
3: ---------------------------------ex30 help---------------------------------\n\
4: -OPTION <DEFAULT> = (UNITS) DESCRIPTION.\n\n\
5: -width <320> = (km) width of domain.\n\
6: -depth <300> = (km) depth of domain.\n\
7: -slab_dip <45> = (degrees) dip angle of the slab (determines the grid aspect ratio).\n\
8: -lid_depth <35> = (km) depth of the static conductive lid.\n\
9: -fault_depth <35> = (km) depth of slab-wedge mechanical coupling\n\
10: ( fault dept >= lid depth ).\n\
11: \n\
12: -ni <82> = grid cells in x-direction. (nj adjusts to accommodate\n\
13: the slab dip & depth). DO NOT USE -da_grid_x option!!!\n\
14: -ivisc <3> = rheology option.\n\
15: 0 --- constant viscosity.\n\
16: 1 --- olivine diffusion creep rheology (T&P-dependent, newtonian).\n\
17: 2 --- olivine dislocation creep rheology (T&P-dependent, non-newtonian).\n\
18: 3 --- Full mantle rheology, combination of 1 & 2.\n\
19: \n\
20: -slab_velocity <5> = (cm/year) convergence rate of slab into subduction zone.\n\
21: -slab_age <50> = (million yrs) age of slab for thermal profile boundary condition.\n\
22: -lid_age <50> = (million yrs) age of lid for thermal profile boundary condition.\n\
23: \n\
24: FOR OTHER PARAMETER OPTIONS AND THEIR DEFAULT VALUES, see SetParams() in ex30.c.\n\
25: ---------------------------------ex30 help---------------------------------\n";
28: /*F-----------------------------------------------------------------------
This PETSc 2.2.0 example by Richard F. Katz
http://www.ldeo.columbia.edu/~katz/
The problem is modeled by the partial differential equation system
\begin{eqnarray}
-\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\
\nabla \cdot v & = & 0 \\
dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\
\end{eqnarray}
\begin{eqnarray}
\eta(T,Eps\_dot) & = & \hbox{constant } \hbox{if ivisc} ==0 \\
& = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\
& = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\
& = & \hbox{mantle viscosity (difn and disl) } \hbox{if ivisc} ==3
\end{eqnarray}
which is uniformly discretized on a staggered mesh:
-------$w_{ij}$------
$u_{i-1j}$ $P,T_{ij}$ $u_{ij}$
------$w_{ij-1}$-----
53: ------------------------------------------------------------------------F*/
55: #include <petscsnes.h>
56: #include <petscdmda.h>
58: #define VISC_CONST 0
59: #define VISC_DIFN 1
60: #define VISC_DISL 2
61: #define VISC_FULL 3
62: #define CELL_CENTER 0
63: #define CELL_CORNER 1
64: #define BC_ANALYTIC 0
65: #define BC_NOSTRESS 1
66: #define BC_EXPERMNT 2
67: #define ADVECT_FV 0
68: #define ADVECT_FROMM 1
69: #define PLATE_SLAB 0
70: #define PLATE_LID 1
71: #define EPS_ZERO 0.00000001
73: typedef struct { /* holds the variables to be solved for */
74: PetscScalar u,w,p,T;
75: } Field;
77: typedef struct { /* parameters needed to compute viscosity */
78: PetscReal A,n,Estar,Vstar;
79: } ViscParam;
81: typedef struct { /* physical and miscelaneous parameters */
82: PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT;
83: PetscReal slab_dip, slab_age, slab_velocity, kappa, z_scale;
84: PetscReal c, d, sb, cb, skt, visc_cutoff, lid_age, eta0, continuation;
85: PetscReal L, V, lid_depth, fault_depth;
86: ViscParam diffusion, dislocation;
87: PetscInt ivisc, adv_scheme, ibound, output_ivisc;
88: PetscBool quiet, param_test, output_to_file, pv_analytic;
89: PetscBool interrupted, stop_solve, toggle_kspmon, kspmon;
90: char filename[PETSC_MAX_PATH_LEN];
91: } Parameter;
93: typedef struct { /* grid parameters */
94: DMDABoundaryType bx,by;
95: DMDAStencilType stencil;
96: PetscInt corner,ni,nj,jlid,jfault,inose;
97: PetscInt dof,stencil_width,mglevels;
98: PetscReal dx,dz;
99: } GridInfo;
101: typedef struct { /* application context */
102: Vec x,Xguess;
103: Parameter *param;
104: GridInfo *grid;
105: } AppCtx;
107: /* Callback functions (static interface) */
108: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
110: /* Main routines */
111: extern PetscErrorCode SetParams(Parameter*, GridInfo *);
112: extern PetscErrorCode ReportParams(Parameter *, GridInfo *);
113: extern PetscErrorCode Initialize(DM);
114: extern PetscErrorCode UpdateSolution(SNES,AppCtx*, PetscInt*);
115: extern PetscErrorCode DoOutput(SNES,PetscInt);
117: /* Post-processing & misc */
118: extern PetscErrorCode ViscosityField(DM,Vec,Vec);
119: extern PetscErrorCode StressField(DM);
120: extern PetscErrorCode SNESConverged_Interactive(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason *, void *);
121: extern PetscErrorCode InteractiveHandler(int, void *);
122: extern PetscBool OptionsHasName(const char pre[],const char name[]);
124: /*-----------------------------------------------------------------------*/
127: int main(int argc,char **argv)
128: /*-----------------------------------------------------------------------*/
129: {
130: SNES snes;
131: AppCtx *user; /* user-defined work context */
132: Parameter param;
133: GridInfo grid;
134: PetscInt nits;
136: MPI_Comm comm;
137: DM da;
139: PetscInitialize(&argc,&argv,(char *)0,help);
140: PetscOptionsSetValue("-file","ex30_output");
141: PetscOptionsSetValue("-snes_monitor_short",PETSC_NULL);
142: PetscOptionsSetValue("-snes_max_it","20");
143: PetscOptionsSetValue("-ksp_max_it","1500");
144: PetscOptionsSetValue("-ksp_gmres_restart","300");
145: PetscOptionsInsert(&argc,&argv,PETSC_NULL);
147: comm = PETSC_COMM_WORLD;
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Set up the problem parameters.
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: SetParams(¶m,&grid);
153: ReportParams(¶m,&grid);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: SNESCreate(comm,&snes);
158: DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
159: SNESSetDM(snes,da);
160: DMDASetFieldName(da,0,"x-velocity");
161: DMDASetFieldName(da,1,"y-velocity");
162: DMDASetFieldName(da,2,"pressure");
163: DMDASetFieldName(da,3,"temperature");
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Create user context, set problem data, create vector data structures.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: PetscMalloc(sizeof(AppCtx),&user);
170: user->param = ¶m;
171: user->grid = &grid;
172: DMSetApplicationContext(da,user);
173: DMCreateGlobalVector(da,&(user->Xguess));
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Set up the SNES solver with callback functions.
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: DMDASetLocalFunction(da,(DMDALocalFunction1)FormFunctionLocal);
180: SNESSetFromOptions(snes);
183: SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,PETSC_NULL);
184: PetscPushSignalHandler(InteractiveHandler,(void*)user);
185:
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Initialize and solve the nonlinear system
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: Initialize(da);
190: UpdateSolution(snes,user,&nits);
191:
192: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: Output variables.
194: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195: DoOutput(snes,nits);
196:
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Free work space.
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200: VecDestroy(&user->Xguess);
201: VecDestroy(&user->x);
202: PetscFree(user);
203: SNESDestroy(&snes);
204: DMDestroy(&da);
205: PetscPopSignalHandler();
206: PetscFinalize();
207: return 0;
208: }
210: /*=====================================================================
211: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
212: =====================================================================*/
214: /*---------------------------------------------------------------------*/
217: /* manages solve: adaptive continuation method */
218: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
219: {
220: KSP ksp;
221: PC pc;
222: SNESConvergedReason reason;
223: Parameter *param = user->param;
224: PetscReal cont_incr=0.3;
225: PetscInt its;
226: PetscErrorCode ierr;
227: PetscBool q = PETSC_FALSE;
228: DM dm;
231: SNESGetDM(snes,&dm);
232: DMCreateGlobalVector(dm,&user->x);
233: SNESGetKSP(snes,&ksp);
234: KSPGetPC(ksp, &pc);
235: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
237: *nits=0;
239: /* Isoviscous solve */
240: if (param->ivisc == VISC_CONST && !param->stop_solve) {
241: param->ivisc = VISC_CONST;
242: SNESSolve(snes,0,user->x);
243: SNESGetConvergedReason(snes,&reason);
244: SNESGetIterationNumber(snes,&its);
245: *nits +=its;
246: VecCopy(user->x,user->Xguess);
247: if (param->stop_solve) goto done;
248: }
250: /* Olivine diffusion creep */
251: if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
252: if (!q) PetscPrintf(PETSC_COMM_WORLD,"Computing Variable Viscosity Solution\n");
254: /* continuation method on viscosity cutoff */
255: for (param->continuation=0.0; ; param->continuation+=cont_incr) {
256: if (!q) PetscPrintf(PETSC_COMM_WORLD," Continuation parameter = %G\n", param->continuation);
258: /* solve the non-linear system */
259: VecCopy(user->Xguess,user->x);
260: SNESSolve(snes,0,user->x);
261: SNESGetConvergedReason(snes,&reason);
262: SNESGetIterationNumber(snes,&its);
263: *nits += its;
264: if (!q) PetscPrintf(PETSC_COMM_WORLD," SNES iterations: %D, Cumulative: %D\n", its, *nits);
265: if (param->stop_solve) goto done;
267: if (reason<0) {
268: /* NOT converged */
269: cont_incr = -fabs(cont_incr)/2.0;
270: if (fabs(cont_incr)<0.01) goto done;
272: } else {
273: /* converged */
274: VecCopy(user->x,user->Xguess);
275: if (param->continuation >= 1.0) goto done;
276: if (its<=3) {
277: cont_incr = 0.30001;
278: } else if (its<=8) {
279: cont_incr = 0.15001;
280: } else {
281: cont_incr = 0.10001;
282: }
283: if (param->continuation+cont_incr > 1.0) {
284: cont_incr = 1.0 - param->continuation;
285: }
286: } /* endif reason<0 */
287: }
288: }
289: done:
290: if (param->stop_solve && !q) PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");
291: if (reason<0 && !q) PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");
292: return(0);
293: }
296: /*=====================================================================
297: PHYSICS FUNCTIONS (compute the discrete residual)
298: =====================================================================*/
301: /*---------------------------------------------------------------------*/
304: PETSC_STATIC_INLINE PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
305: /*---------------------------------------------------------------------*/
306: {
307: return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
308: }
310: /*---------------------------------------------------------------------*/
313: PETSC_STATIC_INLINE PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
314: /*---------------------------------------------------------------------*/
315: {
316: return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
317: }
319: /*---------------------------------------------------------------------*/
322: PETSC_STATIC_INLINE PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
323: /*---------------------------------------------------------------------*/
324: {
325: return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
326: }
328: /*---------------------------------------------------------------------*/
331: PETSC_STATIC_INLINE PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
332: /*---------------------------------------------------------------------*/
333: {
334: return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
335: }
337: /*---------------------------------------------------------------------*/
340: /* isoviscous analytic solution for IC */
341: PETSC_STATIC_INLINE PassiveScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
342: /*---------------------------------------------------------------------*/
343: {
344: Parameter *param = user->param;
345: GridInfo *grid = user->grid;
346: PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
347:
348: x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
349: r = sqrt(x*x+z*z); st = z/r; ct = x/r; th = atan(z/x);
350: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
351: }
353: /*---------------------------------------------------------------------*/
356: /* isoviscous analytic solution for IC */
357: PETSC_STATIC_INLINE PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
358: /*---------------------------------------------------------------------*/
359: {
360: Parameter *param = user->param;
361: GridInfo *grid = user->grid;
362: PetscScalar x, z, r, st, ct, th, c=param->c, d=param->d;
363:
364: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz;
365: r = sqrt(x*x+z*z); st = z/r; ct = x/r; th = atan(z/x);
366: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
367: }
369: /*---------------------------------------------------------------------*/
372: /* isoviscous analytic solution for IC */
373: PETSC_STATIC_INLINE PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
374: /*---------------------------------------------------------------------*/
375: {
376: Parameter *param = user->param;
377: GridInfo *grid = user->grid;
378: PetscScalar x, z, r, st, ct, c=param->c, d=param->d;
380: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
381: r = sqrt(x*x+z*z); st = z/r; ct = x/r;
382: return (-2.0*(c*ct-d*st)/r);
383: }
387: /* computes the second invariant of the strain rate tensor */
388: PETSC_STATIC_INLINE PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
389: /*---------------------------------------------------------------------*/
390: {
391: Parameter *param = user->param;
392: GridInfo *grid = user->grid;
393: PetscInt ilim=grid->ni-1, jlim=grid->nj-1;
394: PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
395: PetscScalar eps11, eps12, eps22;
397: if (i<j) return EPS_ZERO;
398: if (i==ilim) i--; if (j==jlim) j--;
400: if (ipos==CELL_CENTER) { /* on cell center */
401: if (j<=grid->jlid) return EPS_ZERO;
403: uE = x[j][i].u; uW = x[j][i-1].u;
404: wN = x[j][i].w; wS = x[j-1][i].w;
405: wE = WInterp(x,i,j-1);
406: if (i==j) { uN = param->cb; wW = param->sb; }
407: else { uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1); }
409: if (j==grid->jlid+1) uS = 0.0;
410: else uS = UInterp(x,i-1,j-1);
412: } else { /* on CELL_CORNER */
413: if (j<grid->jlid) return EPS_ZERO;
415: uN = x[j+1][i].u; uS = x[j][i].u;
416: wE = x[j][i+1].w; wW = x[j][i].w;
417: if (i==j) { wN = param->sb; uW = param->cb; }
418: else { wN = WInterp(x,i,j); uW = UInterp(x,i-1,j); }
420: if (j==grid->jlid) {
421: uE = 0.0; uW = 0.0;
422: uS = -uN;
423: wS = -wN;
424: } else {
425: uE = UInterp(x,i,j);
426: wS = WInterp(x,i,j-1);
427: }
428: }
430: eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz;
431: eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx);
433: return sqrt( 0.5*( eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22 ) );
434: }
436: /*---------------------------------------------------------------------*/
439: /* computes the shear viscosity */
440: PETSC_STATIC_INLINE PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PassiveScalar z, Parameter *param)
441: /*---------------------------------------------------------------------*/
442: {
443: PetscReal result=0.0;
444: ViscParam difn=param->diffusion, disl=param->dislocation;
445: PetscInt iVisc=param->ivisc;
446: double eps_scale=param->V/(param->L*1000.0);
447: double strain_power, v1, v2, P;
448: double rho_g = 32340.0, R=8.3144;
450: P = rho_g*(z*param->L*1000.0); /* Pa */
452: if (iVisc==VISC_CONST) {
453: /* constant viscosity */
454: return 1.0;
456: } else if (iVisc==VISC_DIFN) {
457: /* diffusion creep rheology */
458: result = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
460: } else if (iVisc==VISC_DISL) {
461: /* dislocation creep rheology */
462: strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
463: result = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
465: } else if (iVisc==VISC_FULL) {
466: /* dislocation/diffusion creep rheology */
467: strain_power = pow( eps*eps_scale, (1.0-disl.n)/disl.n );
468: v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
469: v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
470: result = 1.0/(1.0/v1 + 1.0/v2);
471: }
473: /* max viscosity is param->eta0 */
474: result = PetscMin( result, 1.0 );
475: /* min viscosity is param->visc_cutoff */
476: result = PetscMax( result, param->visc_cutoff );
477: /* continuation method */
478: result = pow(result,param->continuation);
479: return result;
480: }
482: /*---------------------------------------------------------------------*/
485: /* computes the residual of the x-component of eqn (1) above */
486: PETSC_STATIC_INLINE PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
487: /*---------------------------------------------------------------------*/
488: {
489: Parameter *param=user->param;
490: GridInfo *grid =user->grid;
491: PetscScalar dx = grid->dx, dz=grid->dz;
492: PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
493: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
494: PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
495: PetscInt jlim = grid->nj-1;
497: z_scale = param->z_scale;
499: if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
500: TS = param->potentialT * TInterp(x,i,j-1) * exp( (j-1.0)*dz*z_scale );
501: if (j==jlim) TN = TS;
502: else TN = param->potentialT * TInterp(x,i,j) * exp( j *dz*z_scale );
503: TW = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
504: TE = param->potentialT * x[j][i+1].T * exp( (j-0.5)*dz*z_scale );
505: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
506: epsN = CalcSecInv(x,i,j, CELL_CORNER,user);
507: epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
508: epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
509: epsW = CalcSecInv(x,i,j, CELL_CENTER,user);
510: }
511: }
512: etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
513: etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
514: etaW = Viscosity(TW,epsW,dz*j,param);
515: etaE = Viscosity(TE,epsE,dz*j,param);
517: dPdx = ( x[j][i+1].p - x[j][i].p )/dx;
518: if (j==jlim) dudzN = etaN * ( x[j][i].w - x[j][i+1].w )/dx;
519: else dudzN = etaN * ( x[j+1][i].u - x[j][i].u )/dz;
520: dudzS = etaS * ( x[j][i].u - x[j-1][i].u )/dz;
521: dudxE = etaE * ( x[j][i+1].u - x[j][i].u )/dx;
522: dudxW = etaW * ( x[j][i].u - x[j][i-1].u )/dx;
524: residual = -dPdx /* X-MOMENTUM EQUATION*/
525: +( dudxE - dudxW )/dx
526: +( dudzN - dudzS )/dz;
528: if ( param->ivisc!=VISC_CONST ) {
529: dwdxN = etaN * ( x[j][i+1].w - x[j][i].w )/dx;
530: dwdxS = etaS * ( x[j-1][i+1].w - x[j-1][i].w )/dx;
532: residual += ( dudxE - dudxW )/dx + ( dwdxN - dwdxS )/dz;
533: }
535: return residual;
536: }
538: /*---------------------------------------------------------------------*/
541: /* computes the residual of the z-component of eqn (1) above */
542: PETSC_STATIC_INLINE PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
543: /*---------------------------------------------------------------------*/
544: {
545: Parameter *param=user->param;
546: GridInfo *grid =user->grid;
547: PetscScalar dx = grid->dx, dz=grid->dz;
548: 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;
549: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
550: PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
551: PetscInt ilim = grid->ni-1;
553: /* geometric and other parameters */
554: z_scale = param->z_scale;
555:
556: /* viscosity */
557: if ( param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL ) { /* viscosity is T-dependent */
558: TN = param->potentialT * x[j+1][i].T * exp( (j+0.5)*dz*z_scale );
559: TS = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
560: TW = param->potentialT * TInterp(x,i-1,j) * exp( j *dz*z_scale );
561: if (i==ilim) TE = TW;
562: else TE = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
563: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
564: epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
565: epsS = CalcSecInv(x,i,j, CELL_CENTER,user);
566: epsE = CalcSecInv(x,i,j, CELL_CORNER,user);
567: epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
568: }
569: }
570: etaN = Viscosity(TN,epsN,dz*(j+1),param);
571: etaS = Viscosity(TS,epsS,dz*j,param);
572: etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
573: etaE = Viscosity(TE,epsE,dz*(j+0.5),param);
575: dPdz = ( x[j+1][i].p - x[j][i].p )/dz;
576: dwdzN = etaN * ( x[j+1][i].w - x[j][i].w )/dz;
577: dwdzS = etaS * ( x[j][i].w - x[j-1][i].w )/dz;
578: if (i==ilim) dwdxE = etaE * ( x[j][i].u - x[j+1][i].u )/dz;
579: else dwdxE = etaE * ( x[j][i+1].w - x[j][i].w )/dx;
580: dwdxW = 2.0*etaW * ( x[j][i].w - x[j][i-1].w )/dx;
581:
582: /* Z-MOMENTUM */
583: residual = -dPdz /* constant viscosity terms */
584: +( dwdzN - dwdzS )/dz
585: +( dwdxE - dwdxW )/dx;
587: if ( param->ivisc!=VISC_CONST ) {
588: dudzE = etaE * ( x[j+1][i].u - x[j][i].u )/dz;
589: dudzW = etaW * ( x[j+1][i-1].u - x[j][i-1].u )/dz;
591: residual += ( dwdzN - dwdzS )/dz + ( dudzE - dudzW )/dx;
592: }
594: return residual;
595: }
597: /*---------------------------------------------------------------------*/
600: /* computes the residual of eqn (2) above */
601: PETSC_STATIC_INLINE PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
602: /*---------------------------------------------------------------------*/
603: {
604: GridInfo *grid =user->grid;
605: PetscScalar uE,uW,wN,wS,dudx,dwdz;
607: uW = x[j][i-1].u; uE = x[j][i].u; dudx = ( uE - uW )/grid->dx;
608: wS = x[j-1][i].w; wN = x[j][i].w; dwdz = ( wN - wS )/grid->dz;
610: return dudx + dwdz;
611: }
613: /*---------------------------------------------------------------------*/
616: /* computes the residual of eqn (3) above */
617: PETSC_STATIC_INLINE PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
618: /*---------------------------------------------------------------------*/
619: {
620: Parameter *param=user->param;
621: GridInfo *grid =user->grid;
622: PetscScalar dx = grid->dx, dz=grid->dz;
623: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
624: PetscScalar TE, TN, TS, TW, residual;
625: PetscScalar uE,uW,wN,wS;
626: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;
628: dTdzN = ( x[j+1][i].T - x[j][i].T )/dz;
629: dTdzS = ( x[j][i].T - x[j-1][i].T )/dz;
630: dTdxE = ( x[j][i+1].T - x[j][i].T )/dx;
631: dTdxW = ( x[j][i].T - x[j][i-1].T )/dx;
632:
633: residual = ( ( dTdzN - dTdzS )/dz + /* diffusion term */
634: ( dTdxE - dTdxW )/dx )*dx*dz/param->peclet;
636: if (j<=jlid && i>=j) {
637: /* don't advect in the lid */
638: return residual;
640: } else if (i<j) {
641: /* beneath the slab sfc */
642: uW = uE = param->cb;
643: wS = wN = param->sb;
645: } else {
646: /* advect in the slab and wedge */
647: uW = x[j][i-1].u; uE = x[j][i].u;
648: wS = x[j-1][i].w; wN = x[j][i].w;
649: }
651: if ( param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1 ) {
652: /* finite volume advection */
653: TS = ( x[j][i].T + x[j-1][i].T )/2.0;
654: TN = ( x[j][i].T + x[j+1][i].T )/2.0;
655: TE = ( x[j][i].T + x[j][i+1].T )/2.0;
656: TW = ( x[j][i].T + x[j][i-1].T )/2.0;
657: fN = wN*TN*dx; fS = wS*TS*dx;
658: fE = uE*TE*dz; fW = uW*TW*dz;
659:
660: } else {
661: /* Fromm advection scheme */
662: 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
663: - 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;
664: 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
665: - 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;
666: 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
667: - 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;
668: 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
669: - 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;
670: }
671:
672: residual -= ( fE - fW + fN - fS );
674: return residual;
675: }
677: /*---------------------------------------------------------------------*/
680: /* computes the shear stress---used on the boundaries */
681: PETSC_STATIC_INLINE PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
682: /*---------------------------------------------------------------------*/
683: {
684: Parameter *param=user->param;
685: GridInfo *grid=user->grid;
686: PetscInt ilim=grid->ni-1, jlim=grid->nj-1;
687: PetscScalar uN, uS, wE, wW;
689: if ( j<=grid->jlid || i<j || i==ilim || j==jlim ) return EPS_ZERO;
691: if (ipos==CELL_CENTER) { /* on cell center */
693: wE = WInterp(x,i,j-1);
694: if (i==j) { wW = param->sb; uN = param->cb;}
695: else { wW = WInterp(x,i-1,j-1); uN = UInterp(x,i-1,j); }
696: if (j==grid->jlid+1) uS = 0.0;
697: else uS = UInterp(x,i-1,j-1);
699: } else { /* on cell corner */
701: uN = x[j+1][i].u; uS = x[j][i].u;
702: wW = x[j][i].w; wE = x[j][i+1].w;
704: }
706: return (uN-uS)/grid->dz + (wE-wW)/grid->dx;
707: }
709: /*---------------------------------------------------------------------*/
712: /* computes the normal stress---used on the boundaries */
713: PETSC_STATIC_INLINE PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
714: /*---------------------------------------------------------------------*/
715: {
716: Parameter *param=user->param;
717: GridInfo *grid =user->grid;
718: PetscScalar dx = grid->dx, dz=grid->dz;
719: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, ivisc;
720: PetscScalar epsC=0.0, etaC, TC, uE, uW, pC, z_scale;
721: if (i<j || j<=grid->jlid) return EPS_ZERO;
723: ivisc=param->ivisc; z_scale = param->z_scale;
725: if (ipos==CELL_CENTER) { /* on cell center */
727: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
728: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
729: etaC = Viscosity(TC,epsC,dz*j,param);
731: uW = x[j][i-1].u; uE = x[j][i].u;
732: pC = x[j][i].p;
734: } else { /* on cell corner */
735: if ( i==ilim || j==jlim ) return EPS_ZERO;
737: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
738: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
739: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
741: if (i==j) uW = param->sb;
742: else uW = UInterp(x,i-1,j);
743: uE = UInterp(x,i,j); pC = PInterp(x,i,j);
744: }
745:
746: return 2.0*etaC*(uE-uW)/dx - pC;
747: }
749: /*---------------------------------------------------------------------*/
752: /* computes the normal stress---used on the boundaries */
753: PETSC_STATIC_INLINE PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
754: /*---------------------------------------------------------------------*/
755: {
756: Parameter *param=user->param;
757: GridInfo *grid =user->grid;
758: PetscScalar dz=grid->dz;
759: PetscInt ilim=grid->ni-1, jlim=grid->nj-1, ivisc;
760: PetscScalar epsC=0.0, etaC, TC;
761: PetscScalar pC, wN, wS, z_scale;
762: if (i<j || j<=grid->jlid) return EPS_ZERO;
764: ivisc=param->ivisc; z_scale = param->z_scale;
766: if (ipos==CELL_CENTER) { /* on cell center */
768: TC = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*z_scale );
769: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
770: etaC = Viscosity(TC,epsC,dz*j,param);
771: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
773: } else { /* on cell corner */
774: if ( (i==ilim) || (j==jlim) ) return EPS_ZERO;
776: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*z_scale );
777: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
778: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
779: if (i==j) wN = param->sb;
780: else wN = WInterp(x,i,j);
781: wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
782: }
784: return 2.0*etaC*(wN-wS)/dz - pC;
785: }
787: /*---------------------------------------------------------------------*/
789: /*=====================================================================
790: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
791: =====================================================================*/
793: /*---------------------------------------------------------------------*/
796: /* initializes the problem parameters and checks for
797: command line changes */
798: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
799: /*---------------------------------------------------------------------*/
800: {
801: PetscErrorCode ierr, ierr_out=0;
802: PetscReal SEC_PER_YR = 3600.00*24.00*365.2500;
803: PetscReal PI = 3.14159265358979323846;
804: PetscReal alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
805:
806: /* domain geometry */
807: param->slab_dip = 45.0;
808: param->width = 320.0; /* km */
809: param->depth = 300.0; /* km */
810: param->lid_depth = 35.0; /* km */
811: param->fault_depth = 35.0; /* km */
812: PetscOptionsGetReal(PETSC_NULL,"-slab_dip",&(param->slab_dip),PETSC_NULL);
813: PetscOptionsGetReal(PETSC_NULL,"-width",&(param->width),PETSC_NULL);
814: PetscOptionsGetReal(PETSC_NULL,"-depth",&(param->depth),PETSC_NULL);
815: PetscOptionsGetReal(PETSC_NULL,"-lid_depth",&(param->lid_depth),PETSC_NULL);
816: PetscOptionsGetReal(PETSC_NULL,"-fault_depth",&(param->fault_depth),PETSC_NULL);
817: param->slab_dip = param->slab_dip*PI/180.0; /* radians */
819: /* grid information */
820: PetscOptionsGetInt(PETSC_NULL, "-jfault",&(grid->jfault),PETSC_NULL);
821: grid->ni = 82;
822: PetscOptionsGetInt(PETSC_NULL, "-ni",&(grid->ni),PETSC_NULL);
823: grid->dx = param->width/((double)(grid->ni-2)); /* km */
824: grid->dz = grid->dx*tan(param->slab_dip); /* km */
825: grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/
826: param->depth = grid->dz*(grid->nj-2); /* km */
827: grid->inose = 0; /* gridpoints*/
828: PetscOptionsGetInt(PETSC_NULL,"-inose",&(grid->inose),PETSC_NULL);
829: grid->bx = DMDA_BOUNDARY_NONE;
830: grid->by = DMDA_BOUNDARY_NONE;
831: grid->stencil = DMDA_STENCIL_BOX;
832: grid->dof = 4;
833: grid->stencil_width = 2;
834: grid->mglevels = 1;
836: /* boundary conditions */
837: param->pv_analytic = PETSC_FALSE;
838: param->ibound = BC_NOSTRESS;
839: PetscOptionsGetInt(PETSC_NULL,"-ibound",&(param->ibound),PETSC_NULL);
841: /* physical constants */
842: param->slab_velocity = 5.0; /* cm/yr */
843: param->slab_age = 50.0; /* Ma */
844: param->lid_age = 50.0; /* Ma */
845: param->kappa = 0.7272e-6; /* m^2/sec */
846: param->potentialT = 1300.0; /* degrees C */
847: PetscOptionsGetReal(PETSC_NULL,"-slab_velocity",&(param->slab_velocity),PETSC_NULL);
848: PetscOptionsGetReal(PETSC_NULL,"-slab_age",&(param->slab_age),PETSC_NULL);
849: PetscOptionsGetReal(PETSC_NULL,"-lid_age",&(param->lid_age),PETSC_NULL);
850: PetscOptionsGetReal(PETSC_NULL,"-kappa",&(param->kappa),PETSC_NULL);
851: PetscOptionsGetReal(PETSC_NULL,"-potentialT",&(param->potentialT),PETSC_NULL);
853: /* viscosity */
854: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
855: param->eta0 = 1e24; /* Pa-s */
856: param->visc_cutoff = 0.0; /* factor of eta_0 */
857: param->continuation = 1.0;
858: /* constants for diffusion creep */
859: param->diffusion.A = 1.8e7; /* Pa-s */
860: param->diffusion.n = 1.0; /* dim'less */
861: param->diffusion.Estar = 375e3; /* J/mol */
862: param->diffusion.Vstar = 5e-6; /* m^3/mol */
863: /* constants for param->dislocationocation creep */
864: param->dislocation.A = 2.8969e4; /* Pa-s */
865: param->dislocation.n = 3.5; /* dim'less */
866: param->dislocation.Estar = 530e3; /* J/mol */
867: param->dislocation.Vstar = 14e-6; /* m^3/mol */
868: PetscOptionsGetInt(PETSC_NULL, "-ivisc",&(param->ivisc),PETSC_NULL);
869: PetscOptionsGetReal(PETSC_NULL,"-visc_cutoff",&(param->visc_cutoff),PETSC_NULL);
870: param->output_ivisc = param->ivisc;
871: PetscOptionsGetInt(PETSC_NULL,"-output_ivisc",&(param->output_ivisc),PETSC_NULL);
872: PetscOptionsGetReal(PETSC_NULL,"-vstar",&(param->dislocation.Vstar),PETSC_NULL);
874: /* output options */
875: param->quiet = PETSC_FALSE;
876: param->param_test = PETSC_FALSE;
877: PetscOptionsHasName(PETSC_NULL,"-quiet",&(param->quiet));
878: PetscOptionsHasName(PETSC_NULL,"-test",&(param->param_test));
879: PetscOptionsGetString(PETSC_NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));
881: /* advection */
882: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
883: PetscOptionsGetInt(PETSC_NULL,"-adv_scheme",&(param->adv_scheme),PETSC_NULL);
885: /* misc. flags */
886: param->stop_solve = PETSC_FALSE;
887: param->interrupted = PETSC_FALSE;
888: param->kspmon = PETSC_FALSE;
889: param->toggle_kspmon = PETSC_FALSE;
891: /* derived parameters for slab angle */
892: param->sb = sin(param->slab_dip);
893: param->cb = cos(param->slab_dip);
894: param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
895: param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);
897: /* length, velocity and time scale for non-dimensionalization */
898: param->L = PetscMin(param->width,param->depth); /* km */
899: param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */
901: /* other unit conversions and derived parameters */
902: param->scaled_width = param->width/param->L; /* dim'less */
903: param->scaled_depth = param->depth/param->L; /* dim'less */
904: param->lid_depth = param->lid_depth/param->L; /* dim'less */
905: param->fault_depth = param->fault_depth/param->L; /* dim'less */
906: grid->dx = grid->dx/param->L; /* dim'less */
907: grid->dz = grid->dz/param->L; /* dim'less */
908: grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */
909: grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */
910: param->lid_depth = grid->jlid*grid->dz; /* dim'less */
911: param->fault_depth = grid->jfault*grid->dz; /* dim'less */
912: grid->corner = grid->jlid+1; /* gridcells */
913: param->peclet = param->V /* m/sec */
914: * param->L*1000.0 /* m */
915: / param->kappa; /* m^2/sec */
916: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
917: param->skt = sqrt(param->kappa*param->slab_age*SEC_PER_YR);
918: PetscOptionsGetReal(PETSC_NULL,"-peclet",&(param->peclet),PETSC_NULL);
919:
920: return ierr_out;
921: }
923: /*---------------------------------------------------------------------*/
926: /* prints a report of the problem parameters to stdout */
927: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
928: /*---------------------------------------------------------------------*/
929: {
930: PetscErrorCode ierr, ierr_out=0;
931: char date[30];
932: PetscReal PI = 3.14159265358979323846;
934: PetscGetDate(date,30);
936: if ( !(param->quiet) ) {
937: PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
938: /* PetscPrintf(PETSC_COMM_WORLD," %s\n",&(date[0]));*/
940: PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
941: PetscPrintf(PETSC_COMM_WORLD," Width = %G km, Depth = %G km\n",param->width,param->depth);
942: PetscPrintf(PETSC_COMM_WORLD," Slab dip = %G degrees, Slab velocity = %G cm/yr\n",param->slab_dip*180.0/PI,param->slab_velocity);
943: 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);
945: PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
946: 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);
947: PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);
948: PetscPrintf(PETSC_COMM_WORLD," Pe = %G\n",param->peclet);
950: PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
951: if (param->ivisc==VISC_CONST) {
952: PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");
953: if (param->pv_analytic)
954: PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");
955: } else if (param->ivisc==VISC_DIFN) {
956: PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");
957: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
958: } else if (param->ivisc==VISC_DISL ) {
959: PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");
960: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
961: } else if (param->ivisc==VISC_FULL ) {
962: PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");
963: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
964: } else {
965: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
966: ierr_out=1;
967: }
969: PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
970: if ( param->ibound==BC_ANALYTIC ) {
971: PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");
972: } else if ( param->ibound==BC_NOSTRESS ) {
973: PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");
974: } else if ( param->ibound==BC_EXPERMNT ) {
975: PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");
976: } else {
977: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
978: ierr_out=1;
979: }
981: if (param->output_to_file)
982: #if defined(PETSC_HAVE_MATLAB_ENGINE)
983: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);
984: #else
985: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);
986: #endif
987: if ( param->output_ivisc != param->ivisc )
988: PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);
990: PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
991: }
992: if ( param->param_test ) PetscEnd();
993: return ierr_out;
994: }
996: /* ------------------------------------------------------------------- */
999: /* generates an inital guess using the analytic solution for isoviscous
1000: corner flow */
1001: PetscErrorCode Initialize(DM da)
1002: /* ------------------------------------------------------------------- */
1003: {
1004: AppCtx *user;
1005: Parameter *param;
1006: GridInfo *grid;
1007: PetscInt i,j,is,js,im,jm;
1009: Field **x;
1010: Vec Xguess;
1012: /* Get the fine grid */
1013: DMGetApplicationContext(da,&user);
1014: Xguess = user->Xguess;
1015: param = user->param;
1016: grid = user->grid;
1017: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1018: DMDAVecGetArray(da,Xguess,(void**)&x);
1020: /* Compute initial guess */
1021: for (j=js; j<js+jm; j++) {
1022: for (i=is; i<is+im; i++) {
1023: if (i<j) {
1024: x[j][i].u = param->cb;
1025: } else if (j<=grid->jlid) {
1026: x[j][i].u = 0.0;
1027: } else {
1028: x[j][i].u = HorizVelocity(i,j,user);
1029: }
1030: if (i<=j) {
1031: x[j][i].w = param->sb;
1032: } else if (j<=grid->jlid) {
1033: x[j][i].w = 0.0;
1034: } else {
1035: x[j][i].w = VertVelocity(i,j,user);
1036: }
1037: if (i<j || j<=grid->jlid) {
1038: x[j][i].p = 0.0;
1039: } else {
1040: x[j][i].p = Pressure(i,j,user);
1041: }
1042: x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1043: }
1044: }
1046: /* Restore x to Xguess */
1047: DMDAVecRestoreArray(da,Xguess,(void**)&x);
1049: return 0;
1050: }
1052: /*---------------------------------------------------------------------*/
1055: /* controls output to a file */
1056: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1057: /*---------------------------------------------------------------------*/
1058: {
1059: AppCtx *user;
1060: Parameter *param;
1061: GridInfo *grid;
1062: PetscInt ivt;
1064: PetscMPIInt rank;
1065: PetscViewer viewer;
1066: Vec res, pars;
1067: MPI_Comm comm;
1068: DM da;
1070: SNESGetDM(snes,&da);
1071: DMGetApplicationContext(da,&user);
1072: param = user->param;
1073: grid = user->grid;
1074: ivt = param->ivisc;
1076: param->ivisc = param->output_ivisc;
1078: /* compute final residual and final viscosity/strain rate fields */
1079: SNESGetFunction(snes, &res, PETSC_NULL, PETSC_NULL);
1080: ViscosityField(da, user->x, user->Xguess);
1082: /* get the communicator and the rank of the processor */
1083: PetscObjectGetComm((PetscObject)snes, &comm);
1084: MPI_Comm_rank(comm, &rank);
1086: if (param->output_to_file) { /* send output to binary file */
1087: VecCreate(comm, &pars);
1088: if (rank == 0) { /* on processor 0 */
1089: VecSetSizes(pars, 20, PETSC_DETERMINE);
1090: VecSetFromOptions(pars);
1091: VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1092: VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1093: VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1094: VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1095: VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1096: VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1097: /* skipped 6 intentionally */
1098: VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1099: VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1100: VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1101: VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1102: VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1103: VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1104: VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1105: VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1106: VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1107: VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1108: } else { /* on some other processor */
1109: VecSetSizes(pars, 0, PETSC_DETERMINE);
1110: VecSetFromOptions(pars);
1111: }
1112: VecAssemblyBegin(pars); VecAssemblyEnd(pars);
1114: /* create viewer */
1115: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1116: PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1117: #else
1118: PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1119: #endif
1121: /* send vectors to viewer */
1122: PetscObjectSetName((PetscObject)res,"res");
1123: VecView(res,viewer);
1124: PetscObjectSetName((PetscObject)user->x,"out");
1125: VecView(user->x, viewer);
1126: PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1127: VecView(user->Xguess, viewer);
1128: StressField(da); /* compute stress fields */
1129: PetscObjectSetName((PetscObject)(user->Xguess),"str");
1130: VecView(user->Xguess, viewer);
1131: PetscObjectSetName((PetscObject)pars,"par");
1132: VecView(pars, viewer);
1133:
1134: /* destroy viewer and vector */
1135: PetscViewerDestroy(&viewer);
1136: VecDestroy(&pars);
1137: }
1138:
1139: param->ivisc = ivt;
1140: return 0;
1141: }
1143: /* ------------------------------------------------------------------- */
1146: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1147: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1148: /* ------------------------------------------------------------------- */
1149: {
1150: AppCtx *user;
1151: Parameter *param;
1152: GridInfo *grid;
1153: Vec localX;
1154: Field **v, **x;
1155: PassiveReal eps, /* dx,*/ dz, T, epsC, TC;
1156: PetscInt i,j,is,js,im,jm,ilim,jlim,ivt;
1160: DMGetApplicationContext(da,&user);
1161: param = user->param;
1162: grid = user->grid;
1163: ivt = param->ivisc;
1164: param->ivisc = param->output_ivisc;
1166: DMGetLocalVector(da, &localX);
1167: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1168: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1169: DMDAVecGetArray(da,localX,(void**)&x);
1170: DMDAVecGetArray(da,V,(void**)&v);
1172: /* Parameters */
1173: /* dx = grid->dx; */ dz = grid->dz;
1174: ilim = grid->ni-1; jlim = grid->nj-1;
1176: /* Compute real temperature, strain rate and viscosity */
1177: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1178: for (j=js; j<js+jm; j++) {
1179: for (i=is; i<is+im; i++) {
1180: T = param->potentialT * x[j][i].T * exp( (j-0.5)*dz*param->z_scale );
1181: if (i<ilim && j<jlim) {
1182: TC = param->potentialT * TInterp(x,i,j) * exp( j*dz*param->z_scale );
1183: } else {
1184: TC = T;
1185: }
1186: eps = CalcSecInv(x,i,j,CELL_CENTER,user);
1187: epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
1188: v[j][i].u = eps;
1189: v[j][i].w = epsC;
1190: v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1191: v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1192: }
1193: }
1194: DMDAVecRestoreArray(da,V,(void**)&v);
1195: DMDAVecRestoreArray(da,localX,(void**)&x);
1196: DMRestoreLocalVector(da, &localX);
1197: param->ivisc = ivt;
1198: return(0);
1199: }
1201: /* ------------------------------------------------------------------- */
1204: /* post-processing: compute stress everywhere */
1205: PetscErrorCode StressField(DM da)
1206: /* ------------------------------------------------------------------- */
1207: {
1208: AppCtx *user;
1209: PetscInt i,j,is,js,im,jm;
1211: Vec locVec;
1212: Field **x, **y;
1214: DMGetApplicationContext(da,&user);
1216: /* Get the fine grid of Xguess and X */
1217: DMDAGetCorners(da,&is,&js,PETSC_NULL,&im,&jm,PETSC_NULL);
1218: DMDAVecGetArray(da,user->Xguess,(void**)&x);
1220: DMGetLocalVector(da, &locVec);
1221: DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);
1222: DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);
1223: DMDAVecGetArray(da,locVec,(void**)&y);
1225: /* Compute stress on the corner points */
1226: for (j=js; j<js+jm; j++) {
1227: for (i=is; i<is+im; i++) {
1228:
1229: x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1230: x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1231: x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1232: x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1233: }
1234: }
1236: /* Restore the fine grid of Xguess and X */
1237: DMDAVecRestoreArray(da,user->Xguess,(void**)&x);
1238: DMDAVecRestoreArray(da,locVec,(void**)&y);
1239: DMRestoreLocalVector(da, &locVec);
1240: return 0;
1241: }
1243: /*=====================================================================
1244: UTILITY FUNCTIONS
1245: =====================================================================*/
1247: /*---------------------------------------------------------------------*/
1250: /* returns the velocity of the subducting slab and handles fault nodes
1251: for BC */
1252: PETSC_STATIC_INLINE PassiveScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1253: /*---------------------------------------------------------------------*/
1254: {
1255: Parameter *param = user->param;
1256: GridInfo *grid = user->grid;
1258: if (c=='U' || c=='u') {
1259: if (i<j-1) {
1260: return param->cb;
1261: } else if (j<=grid->jfault) {
1262: return 0.0;
1263: } else
1264: return param->cb;
1266: } else {
1267: if (i<j) {
1268: return param->sb;
1269: } else if (j<=grid->jfault) {
1270: return 0.0;
1271: } else
1272: return param->sb;
1273: }
1274: }
1276: /*---------------------------------------------------------------------*/
1279: /* solution to diffusive half-space cooling model for BC */
1280: PETSC_STATIC_INLINE PassiveScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1281: /*---------------------------------------------------------------------*/
1282: {
1283: Parameter *param = user->param;
1284: PassiveScalar z;
1285: if (plate==PLATE_LID)
1286: z = (j-0.5)*user->grid->dz;
1287: else /* PLATE_SLAB */
1288: z = (j-0.5)*user->grid->dz*param->cb;
1289: #if defined (PETSC_HAVE_ERF)
1290: return erf(z*param->L/2.0/param->skt);
1291: #else
1292: SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
1293: #endif
1294: }
1297: /* ------------------------------------------------------------------- */
1300: /* utility function */
1301: PetscBool OptionsHasName(const char pre[],const char name[])
1302: /* ------------------------------------------------------------------- */
1303: {
1304: PetscBool retval;
1306: PetscOptionsHasName(pre,name,&retval);CHKERRABORT(PETSC_COMM_WORLD,ierr);
1307: return retval;
1308: }
1310: /*=====================================================================
1311: INTERACTIVE SIGNAL HANDLING
1312: =====================================================================*/
1314: /* ------------------------------------------------------------------- */
1317: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1318: /* ------------------------------------------------------------------- */
1319: {
1320: AppCtx *user = (AppCtx *) ctx;
1321: Parameter *param = user->param;
1322: KSP ksp;
1326: if (param->interrupted) {
1327: param->interrupted = PETSC_FALSE;
1328: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1329: *reason = SNES_CONVERGED_FNORM_ABS;
1330: return(0);
1331: } else if (param->toggle_kspmon) {
1332: param->toggle_kspmon = PETSC_FALSE;
1333: SNESGetKSP(snes, &ksp);
1334: if (param->kspmon) {
1335: KSPMonitorCancel(ksp);
1336: param->kspmon = PETSC_FALSE;
1337: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1338: } else {
1339: KSPMonitorSet(ksp,KSPMonitorSingularValue,PETSC_NULL,PETSC_NULL);
1340: param->kspmon = PETSC_TRUE;
1341: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1342: }
1343: }
1344: PetscFunctionReturn(SNESDefaultConverged(snes,it,xnorm,snorm,fnorm,reason,ctx));
1345: }
1347: /* ------------------------------------------------------------------- */
1348: #include <signal.h>
1351: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1352: /* ------------------------------------------------------------------- */
1353: {
1354: AppCtx *user = (AppCtx *) ctx;
1355: Parameter *param = user->param;
1357: if (signum == SIGILL) {
1358: param->toggle_kspmon = PETSC_TRUE;
1359: #if !defined(PETSC_MISSING_SIGCONT)
1360: } else if (signum == SIGCONT) {
1361: param->interrupted = PETSC_TRUE;
1362: #endif
1363: #if !defined(PETSC_MISSING_SIGURG)
1364: } else if (signum == SIGURG) {
1365: param->stop_solve = PETSC_TRUE;
1366: #endif
1367: }
1368: return 0;
1369: }
1371: /*---------------------------------------------------------------------*/
1374: /* main call-back function that computes the processor-local piece
1375: of the residual */
1376: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
1377: /*---------------------------------------------------------------------*/
1378: {
1379: AppCtx *user = (AppCtx*)ptr;
1380: Parameter *param = user->param;
1381: GridInfo *grid = user->grid;
1382: PetscScalar mag_w, mag_u;
1383: PetscInt i,j,mx,mz,ilim,jlim;
1384: PetscInt is,ie,js,je,ibound; /* ,ivisc */
1388: /* Define global and local grid parameters */
1389: mx = info->mx; mz = info->my;
1390: ilim = mx-1; jlim = mz-1;
1391: is = info->xs; ie = info->xs+info->xm;
1392: js = info->ys; je = info->ys+info->ym;
1394: /* Define geometric and numeric parameters */
1395: /* ivisc = param->ivisc; */ ibound = param->ibound;
1397: for (j=js; j<je; j++) {
1398: for (i=is; i<ie; i++) {
1400: /************* X-MOMENTUM/VELOCITY *************/
1401: if (i<j) {
1402: f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
1404: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1405: /* in the lithospheric lid */
1406: f[j][i].u = x[j][i].u - 0.0;
1408: } else if (i==ilim) {
1409: /* on the right side boundary */
1410: if (ibound==BC_ANALYTIC) {
1411: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1412: } else {
1413: f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1414: }
1416: } else if (j==jlim) {
1417: /* on the bottom boundary */
1418: if (ibound==BC_ANALYTIC) {
1419: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1420: } else if (ibound==BC_NOSTRESS) {
1421: f[j][i].u = XMomentumResidual(x,i,j,user);
1422: } else {
1423: /* experimental boundary condition */
1424: }
1426: } else {
1427: /* in the mantle wedge */
1428: f[j][i].u = XMomentumResidual(x,i,j,user);
1429: }
1430:
1431: /************* Z-MOMENTUM/VELOCITY *************/
1432: if (i<=j) {
1433: f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);
1435: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1436: /* in the lithospheric lid */
1437: f[j][i].w = x[j][i].w - 0.0;
1439: } else if (j==jlim) {
1440: /* on the bottom boundary */
1441: if (ibound==BC_ANALYTIC) {
1442: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1443: } else {
1444: f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1445: }
1447: } else if (i==ilim) {
1448: /* on the right side boundary */
1449: if (ibound==BC_ANALYTIC) {
1450: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1451: } else if (ibound==BC_NOSTRESS) {
1452: f[j][i].w = ZMomentumResidual(x,i,j,user);
1453: } else {
1454: /* experimental boundary condition */
1455: }
1457: } else {
1458: /* in the mantle wedge */
1459: f[j][i].w = ZMomentumResidual(x,i,j,user);
1460: }
1462: /************* CONTINUITY/PRESSURE *************/
1463: if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1464: /* in the lid or slab */
1465: f[j][i].p = x[j][i].p;
1466:
1467: } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
1468: /* on an analytic boundary */
1469: f[j][i].p = x[j][i].p - Pressure(i,j,user);
1471: } else {
1472: /* in the mantle wedge */
1473: f[j][i].p = ContinuityResidual(x,i,j,user);
1474: }
1476: /************* TEMPERATURE *************/
1477: if (j==0) {
1478: /* on the surface */
1479: f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0);
1481: } else if (i==0) {
1482: /* slab inflow boundary */
1483: f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user);
1485: } else if (i==ilim) {
1486: /* right side boundary */
1487: mag_u = 1.0 - pow( (1.0-PetscMax(PetscMin(PetscRealPart(x[j][i-1].u)/param->cb,1.0),0.0)), 5.0 );
1488: 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);
1490: } else if (j==jlim) {
1491: /* bottom boundary */
1492: mag_w = 1.0 - pow( (1.0-PetscMax(PetscMin(PetscRealPart(x[j-1][i].w)/param->sb,1.0),0.0)), 5.0 );
1493: f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);
1495: } else {
1496: /* in the mantle wedge */
1497: f[j][i].T = EnergyResidual(x,i,j,user);
1498: }
1499: }
1500: }
1501: return(0);
1502: }