Actual source code: ex30.c
petsc-3.8.4 2018-03-24
1: static const char help[] = "Steady-state 2D subduction flow, pressure and temperature solver.\n\
2: The flow is driven by the subducting slab.\n\
3: ---------------------------------ex30 help---------------------------------\n\
4: -OPTION <DEFAULT> = (UNITS) DESCRIPTION.\n\n\
5: -width <320> = (km) width of domain.\n\
6: -depth <300> = (km) depth of domain.\n\
7: -slab_dip <45> = (degrees) dip angle of the slab (determines the grid aspect ratio).\n\
8: -lid_depth <35> = (km) depth of the static conductive lid.\n\
9: -fault_depth <35> = (km) depth of slab-wedge mechanical coupling\n\
10: (fault dept >= lid depth).\n\
11: \n\
12: -ni <82> = grid cells in x-direction. (nj adjusts to accommodate\n\
13: the slab dip & depth). DO NOT USE -da_grid_x option!!!\n\
14: -ivisc <3> = rheology option.\n\
15: 0 --- constant viscosity.\n\
16: 1 --- olivine diffusion creep rheology (T&P-dependent, newtonian).\n\
17: 2 --- olivine dislocation creep rheology (T&P-dependent, non-newtonian).\n\
18: 3 --- Full mantle rheology, combination of 1 & 2.\n\
19: \n\
20: -slab_velocity <5> = (cm/year) convergence rate of slab into subduction zone.\n\
21: -slab_age <50> = (million yrs) age of slab for thermal profile boundary condition.\n\
22: -lid_age <50> = (million yrs) age of lid for thermal profile boundary condition.\n\
23: \n\
24: FOR OTHER PARAMETER OPTIONS AND THEIR DEFAULT VALUES, see SetParams() in ex30.c.\n\
25: ---------------------------------ex30 help---------------------------------\n";
This PETSc 2.2.0 example by Richard F. Katz
http://www.ldeo.columbia.edu/~katz/
The problem is modeled by the partial differential equation system
\begin{eqnarray}
-\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\
\nabla \cdot v & = & 0 \\
dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\
\end{eqnarray}
\begin{eqnarray}
\eta(T,Eps\_dot) & = & \hbox{constant } \hbox{if ivisc} ==0 \\
& = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\
& = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\
& = & \hbox{mantle viscosity (difn and disl) } \hbox{if ivisc} ==3
\end{eqnarray}
which is uniformly discretized on a staggered mesh:
-------$w_{ij}$------
$u_{i-1j}$ $P,T_{ij}$ $u_{ij}$
------$w_{ij-1}$-----
55: #include <petscsnes.h>
56: #include <petscdm.h>
57: #include <petscdmda.h>
59: #define VISC_CONST 0
60: #define VISC_DIFN 1
61: #define VISC_DISL 2
62: #define VISC_FULL 3
63: #define CELL_CENTER 0
64: #define CELL_CORNER 1
65: #define BC_ANALYTIC 0
66: #define BC_NOSTRESS 1
67: #define BC_EXPERMNT 2
68: #define ADVECT_FV 0
69: #define ADVECT_FROMM 1
70: #define PLATE_SLAB 0
71: #define PLATE_LID 1
72: #define EPS_ZERO 0.00000001
74: typedef struct { /* holds the variables to be solved for */
75: PetscScalar u,w,p,T;
76: } Field;
78: typedef struct { /* parameters needed to compute viscosity */
79: PetscReal A,n,Estar,Vstar;
80: } ViscParam;
82: typedef struct { /* physical and miscelaneous parameters */
83: PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT;
84: PetscReal slab_dip, slab_age, slab_velocity, kappa, z_scale;
85: PetscReal c, d, sb, cb, skt, visc_cutoff, lid_age, eta0, continuation;
86: PetscReal L, V, lid_depth, fault_depth;
87: ViscParam diffusion, dislocation;
88: PetscInt ivisc, adv_scheme, ibound, output_ivisc;
89: PetscBool quiet, param_test, output_to_file, pv_analytic;
90: PetscBool interrupted, stop_solve, toggle_kspmon, kspmon;
91: char filename[PETSC_MAX_PATH_LEN];
92: } Parameter;
94: typedef struct { /* grid parameters */
95: DMBoundaryType bx,by;
96: DMDAStencilType stencil;
97: PetscInt corner,ni,nj,jlid,jfault,inose;
98: PetscInt dof,stencil_width,mglevels;
99: PetscReal dx,dz;
100: } GridInfo;
102: typedef struct { /* application context */
103: Vec x,Xguess;
104: Parameter *param;
105: GridInfo *grid;
106: } AppCtx;
108: /* Callback functions (static interface) */
109: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
111: /* Main routines */
112: extern PetscErrorCode SetParams(Parameter*, GridInfo*);
113: extern PetscErrorCode ReportParams(Parameter*, GridInfo*);
114: extern PetscErrorCode Initialize(DM);
115: extern PetscErrorCode UpdateSolution(SNES,AppCtx*, PetscInt*);
116: extern PetscErrorCode DoOutput(SNES,PetscInt);
118: /* Post-processing & misc */
119: extern PetscErrorCode ViscosityField(DM,Vec,Vec);
120: extern PetscErrorCode StressField(DM);
121: extern PetscErrorCode SNESConverged_Interactive(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason*, void*);
122: extern PetscErrorCode InteractiveHandler(int, void*);
123: extern PetscBool OptionsHasName(const char pre[],const char name[]);
125: /*-----------------------------------------------------------------------*/
126: int main(int argc,char **argv)
127: /*-----------------------------------------------------------------------*/
128: {
129: SNES snes;
130: AppCtx *user; /* user-defined work context */
131: Parameter param;
132: GridInfo grid;
133: PetscInt nits;
135: MPI_Comm comm;
136: DM da;
138: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
139: PetscOptionsSetValue(NULL,"-file","ex30_output");
140: PetscOptionsSetValue(NULL,"-snes_monitor_short",NULL);
141: PetscOptionsSetValue(NULL,"-snes_max_it","20");
142: PetscOptionsSetValue(NULL,"-ksp_max_it","1500");
143: PetscOptionsSetValue(NULL,"-ksp_gmres_restart","300");
144: PetscOptionsInsert(NULL,&argc,&argv,NULL);
146: comm = PETSC_COMM_WORLD;
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Set up the problem parameters.
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: SetParams(¶m,&grid);
152: ReportParams(¶m,&grid);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: SNESCreate(comm,&snes);
157: DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
158: DMSetFromOptions(da);
159: DMSetUp(da);
160: SNESSetDM(snes,da);
161: DMDASetFieldName(da,0,"x-velocity");
162: DMDASetFieldName(da,1,"y-velocity");
163: DMDASetFieldName(da,2,"pressure");
164: DMDASetFieldName(da,3,"temperature");
167: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: Create user context, set problem data, create vector data structures.
169: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscNew(&user);
171: user->param = ¶m;
172: user->grid = &grid;
173: DMSetApplicationContext(da,user);
174: DMCreateGlobalVector(da,&(user->Xguess));
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Set up the SNES solver with callback functions.
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,(void*)user);
181: SNESSetFromOptions(snes);
184: SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,NULL);
185: PetscPushSignalHandler(InteractiveHandler,(void*)user);
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Initialize and solve the nonlinear system
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: Initialize(da);
191: UpdateSolution(snes,user,&nits);
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Output variables.
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: DoOutput(snes,nits);
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Free work space.
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: VecDestroy(&user->Xguess);
202: VecDestroy(&user->x);
203: PetscFree(user);
204: SNESDestroy(&snes);
205: DMDestroy(&da);
206: PetscPopSignalHandler();
207: PetscFinalize();
208: return ierr;
209: }
211: /*=====================================================================
212: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
213: =====================================================================*/
215: /*---------------------------------------------------------------------*/
216: /* manages solve: adaptive continuation method */
217: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
218: {
219: KSP ksp;
220: PC pc;
221: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
222: Parameter *param = user->param;
223: PetscReal cont_incr=0.3;
224: PetscInt its;
225: PetscErrorCode ierr;
226: PetscBool q = PETSC_FALSE;
227: DM dm;
230: SNESGetDM(snes,&dm);
231: DMCreateGlobalVector(dm,&user->x);
232: SNESGetKSP(snes,&ksp);
233: KSPGetPC(ksp, &pc);
234: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
236: *nits=0;
238: /* Isoviscous solve */
239: if (param->ivisc == VISC_CONST && !param->stop_solve) {
240: 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", (double)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 = -PetscAbsReal(cont_incr)/2.0;
270: if (PetscAbsReal(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) cont_incr = 0.30001;
277: else if (its<=8) cont_incr = 0.15001;
278: else cont_incr = 0.10001;
280: if (param->continuation+cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
281: } /* endif reason<0 */
282: }
283: }
284: done:
285: if (param->stop_solve && !q) {PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");}
286: if (reason<0 && !q) {PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");}
287: return(0);
288: }
291: /*=====================================================================
292: PHYSICS FUNCTIONS (compute the discrete residual)
293: =====================================================================*/
296: /*---------------------------------------------------------------------*/
297: PETSC_STATIC_INLINE PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
298: /*---------------------------------------------------------------------*/
299: {
300: return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
301: }
303: /*---------------------------------------------------------------------*/
304: PETSC_STATIC_INLINE PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
305: /*---------------------------------------------------------------------*/
306: {
307: return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
308: }
310: /*---------------------------------------------------------------------*/
311: PETSC_STATIC_INLINE PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
312: /*---------------------------------------------------------------------*/
313: {
314: return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
315: }
317: /*---------------------------------------------------------------------*/
318: PETSC_STATIC_INLINE PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
319: /*---------------------------------------------------------------------*/
320: {
321: return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
322: }
324: /*---------------------------------------------------------------------*/
325: /* isoviscous analytic solution for IC */
326: PETSC_STATIC_INLINE PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
327: /*---------------------------------------------------------------------*/
328: {
329: Parameter *param = user->param;
330: GridInfo *grid = user->grid;
331: PetscScalar st, ct, th, c=param->c, d=param->d;
332: PetscReal x, z,r;
334: x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
335: r = PetscSqrtReal(x*x+z*z);
336: st = z/r;
337: ct = x/r;
338: th = PetscAtanReal(z/x);
339: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
340: }
342: /*---------------------------------------------------------------------*/
343: /* isoviscous analytic solution for IC */
344: PETSC_STATIC_INLINE PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
345: /*---------------------------------------------------------------------*/
346: {
347: Parameter *param = user->param;
348: GridInfo *grid = user->grid;
349: PetscScalar st, ct, th, c=param->c, d=param->d;
350: PetscReal x, z, r;
352: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz;
353: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r; th = PetscAtanReal(z/x);
354: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
355: }
357: /*---------------------------------------------------------------------*/
358: /* isoviscous analytic solution for IC */
359: PETSC_STATIC_INLINE PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
360: /*---------------------------------------------------------------------*/
361: {
362: Parameter *param = user->param;
363: GridInfo *grid = user->grid;
364: PetscScalar x, z, r, st, ct, c=param->c, d=param->d;
366: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
367: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r;
368: return (-2.0*(c*ct-d*st)/r);
369: }
371: /* computes the second invariant of the strain rate tensor */
372: PETSC_STATIC_INLINE PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
373: /*---------------------------------------------------------------------*/
374: {
375: Parameter *param = user->param;
376: GridInfo *grid = user->grid;
377: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
378: PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
379: PetscScalar eps11, eps12, eps22;
381: if (i<j) return EPS_ZERO;
382: if (i==ilim) i--;
383: if (j==jlim) j--;
385: if (ipos==CELL_CENTER) { /* on cell center */
386: if (j<=grid->jlid) return EPS_ZERO;
388: uE = x[j][i].u; uW = x[j][i-1].u;
389: wN = x[j][i].w; wS = x[j-1][i].w;
390: wE = WInterp(x,i,j-1);
391: if (i==j) {
392: uN = param->cb; wW = param->sb;
393: } else {
394: uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1);
395: }
397: if (j==grid->jlid+1) uS = 0.0;
398: else uS = UInterp(x,i-1,j-1);
400: } else { /* on CELL_CORNER */
401: if (j<grid->jlid) return EPS_ZERO;
403: uN = x[j+1][i].u; uS = x[j][i].u;
404: wE = x[j][i+1].w; wW = x[j][i].w;
405: if (i==j) {
406: wN = param->sb;
407: uW = param->cb;
408: } else {
409: wN = WInterp(x,i,j);
410: uW = UInterp(x,i-1,j);
411: }
413: if (j==grid->jlid) {
414: uE = 0.0; uW = 0.0;
415: uS = -uN;
416: wS = -wN;
417: } else {
418: uE = UInterp(x,i,j);
419: wS = WInterp(x,i,j-1);
420: }
421: }
423: eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz;
424: eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx);
426: return PetscSqrtReal(0.5*(eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22));
427: }
429: /*---------------------------------------------------------------------*/
430: /* computes the shear viscosity */
431: PETSC_STATIC_INLINE PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
432: /*---------------------------------------------------------------------*/
433: {
434: PetscReal result =0.0;
435: ViscParam difn =param->diffusion, disl=param->dislocation;
436: PetscInt iVisc =param->ivisc;
437: PetscScalar eps_scale=param->V/(param->L*1000.0);
438: PetscScalar strain_power, v1, v2, P;
439: PetscScalar rho_g = 32340.0, R=8.3144;
441: P = rho_g*(z*param->L*1000.0); /* Pa */
443: if (iVisc==VISC_CONST) {
444: /* constant viscosity */
445: return 1.0;
446: } else if (iVisc==VISC_DIFN) {
447: /* diffusion creep rheology */
448: result = PetscRealPart((difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0));
449: } else if (iVisc==VISC_DISL) {
450: /* dislocation creep rheology */
451: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
453: result = PetscRealPart(disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0);
454: } else if (iVisc==VISC_FULL) {
455: /* dislocation/diffusion creep rheology */
456: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
458: v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
459: v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
461: result = PetscRealPart(1.0/(1.0/v1 + 1.0/v2));
462: }
464: /* max viscosity is param->eta0 */
465: result = PetscMin(result, 1.0);
466: /* min viscosity is param->visc_cutoff */
467: result = PetscMax(result, param->visc_cutoff);
468: /* continuation method */
469: result = PetscPowReal(result,param->continuation);
470: return result;
471: }
473: /*---------------------------------------------------------------------*/
474: /* computes the residual of the x-component of eqn (1) above */
475: PETSC_STATIC_INLINE PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
476: /*---------------------------------------------------------------------*/
477: {
478: Parameter *param=user->param;
479: GridInfo *grid =user->grid;
480: PetscScalar dx = grid->dx, dz=grid->dz;
481: PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
482: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
483: PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
484: PetscInt jlim = grid->nj-1;
486: z_scale = param->z_scale;
488: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
489: TS = param->potentialT * TInterp(x,i,j-1) * PetscExpScalar((j-1.0)*dz*z_scale);
490: if (j==jlim) TN = TS;
491: else TN = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
492: TW = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
493: TE = param->potentialT * x[j][i+1].T * PetscExpScalar((j-0.5)*dz*z_scale);
494: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
495: epsN = CalcSecInv(x,i,j, CELL_CORNER,user);
496: epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
497: epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
498: epsW = CalcSecInv(x,i,j, CELL_CENTER,user);
499: }
500: }
501: etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
502: etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
503: etaW = Viscosity(TW,epsW,dz*j,param);
504: etaE = Viscosity(TE,epsE,dz*j,param);
506: dPdx = (x[j][i+1].p - x[j][i].p)/dx;
507: if (j==jlim) dudzN = etaN * (x[j][i].w - x[j][i+1].w)/dx;
508: else dudzN = etaN * (x[j+1][i].u - x[j][i].u) /dz;
509: dudzS = etaS * (x[j][i].u - x[j-1][i].u)/dz;
510: dudxE = etaE * (x[j][i+1].u - x[j][i].u) /dx;
511: dudxW = etaW * (x[j][i].u - x[j][i-1].u)/dx;
513: residual = -dPdx /* X-MOMENTUM EQUATION*/
514: +(dudxE - dudxW)/dx
515: +(dudzN - dudzS)/dz;
517: if (param->ivisc!=VISC_CONST) {
518: dwdxN = etaN * (x[j][i+1].w - x[j][i].w) /dx;
519: dwdxS = etaS * (x[j-1][i+1].w - x[j-1][i].w)/dx;
521: residual += (dudxE - dudxW)/dx + (dwdxN - dwdxS)/dz;
522: }
524: return residual;
525: }
527: /*---------------------------------------------------------------------*/
528: /* computes the residual of the z-component of eqn (1) above */
529: PETSC_STATIC_INLINE PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
530: /*---------------------------------------------------------------------*/
531: {
532: Parameter *param=user->param;
533: GridInfo *grid =user->grid;
534: PetscScalar dx = grid->dx, dz=grid->dz;
535: 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;
536: PetscScalar TE =0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
537: PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
538: PetscInt ilim = grid->ni-1;
540: /* geometric and other parameters */
541: z_scale = param->z_scale;
543: /* viscosity */
544: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
545: TN = param->potentialT * x[j+1][i].T * PetscExpScalar((j+0.5)*dz*z_scale);
546: TS = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
547: TW = param->potentialT * TInterp(x,i-1,j) * PetscExpScalar(j*dz*z_scale);
548: if (i==ilim) TE = TW;
549: else TE = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
550: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
551: epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
552: epsS = CalcSecInv(x,i,j, CELL_CENTER,user);
553: epsE = CalcSecInv(x,i,j, CELL_CORNER,user);
554: epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
555: }
556: }
557: etaN = Viscosity(TN,epsN,dz*(j+1.0),param);
558: etaS = Viscosity(TS,epsS,dz*(j+0.0),param);
559: etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
560: etaE = Viscosity(TE,epsE,dz*(j+0.5),param);
562: dPdz = (x[j+1][i].p - x[j][i].p)/dz;
563: dwdzN = etaN * (x[j+1][i].w - x[j][i].w)/dz;
564: dwdzS = etaS * (x[j][i].w - x[j-1][i].w)/dz;
565: if (i==ilim) dwdxE = etaE * (x[j][i].u - x[j+1][i].u)/dz;
566: else dwdxE = etaE * (x[j][i+1].w - x[j][i].w) /dx;
567: dwdxW = 2.0*etaW * (x[j][i].w - x[j][i-1].w)/dx;
569: /* Z-MOMENTUM */
570: residual = -dPdz /* constant viscosity terms */
571: +(dwdzN - dwdzS)/dz
572: +(dwdxE - dwdxW)/dx;
574: if (param->ivisc!=VISC_CONST) {
575: dudzE = etaE * (x[j+1][i].u - x[j][i].u)/dz;
576: dudzW = etaW * (x[j+1][i-1].u - x[j][i-1].u)/dz;
578: residual += (dwdzN - dwdzS)/dz + (dudzE - dudzW)/dx;
579: }
581: return residual;
582: }
584: /*---------------------------------------------------------------------*/
585: /* computes the residual of eqn (2) above */
586: PETSC_STATIC_INLINE PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
587: /*---------------------------------------------------------------------*/
588: {
589: GridInfo *grid =user->grid;
590: PetscScalar uE,uW,wN,wS,dudx,dwdz;
592: uW = x[j][i-1].u; uE = x[j][i].u; dudx = (uE - uW)/grid->dx;
593: wS = x[j-1][i].w; wN = x[j][i].w; dwdz = (wN - wS)/grid->dz;
595: return dudx + dwdz;
596: }
598: /*---------------------------------------------------------------------*/
599: /* computes the residual of eqn (3) above */
600: PETSC_STATIC_INLINE PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
601: /*---------------------------------------------------------------------*/
602: {
603: Parameter *param=user->param;
604: GridInfo *grid =user->grid;
605: PetscScalar dx = grid->dx, dz=grid->dz;
606: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
607: PetscScalar TE, TN, TS, TW, residual;
608: PetscScalar uE,uW,wN,wS;
609: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;
611: dTdzN = (x[j+1][i].T - x[j][i].T) /dz;
612: dTdzS = (x[j][i].T - x[j-1][i].T)/dz;
613: dTdxE = (x[j][i+1].T - x[j][i].T) /dx;
614: dTdxW = (x[j][i].T - x[j][i-1].T)/dx;
616: residual = ((dTdzN - dTdzS)/dz + /* diffusion term */
617: (dTdxE - dTdxW)/dx)*dx*dz/param->peclet;
619: if (j<=jlid && i>=j) {
620: /* don't advect in the lid */
621: return residual;
622: } else if (i<j) {
623: /* beneath the slab sfc */
624: uW = uE = param->cb;
625: wS = wN = param->sb;
626: } else {
627: /* advect in the slab and wedge */
628: uW = x[j][i-1].u; uE = x[j][i].u;
629: wS = x[j-1][i].w; wN = x[j][i].w;
630: }
632: if (param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1) {
633: /* finite volume advection */
634: TS = (x[j][i].T + x[j-1][i].T)/2.0;
635: TN = (x[j][i].T + x[j+1][i].T)/2.0;
636: TE = (x[j][i].T + x[j][i+1].T)/2.0;
637: TW = (x[j][i].T + x[j][i-1].T)/2.0;
638: fN = wN*TN*dx; fS = wS*TS*dx;
639: fE = uE*TE*dz; fW = uW*TW*dz;
641: } else {
642: /* Fromm advection scheme */
643: 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
644: - PetscAbsScalar(uE)*(-x[j][i+2].T + 3.0*(x[j][i+1].T-x[j][i].T)+x[j][i-1].T)/8.0)*dz;
645: 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
646: - PetscAbsScalar(uW)*(-x[j][i+1].T + 3.0*(x[j][i].T-x[j][i-1].T)+x[j][i-2].T)/8.0)*dz;
647: 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
648: - PetscAbsScalar(wN)*(-x[j+2][i].T + 3.0*(x[j+1][i].T-x[j][i].T)+x[j-1][i].T)/8.0)*dx;
649: 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
650: - PetscAbsScalar(wS)*(-x[j+1][i].T + 3.0*(x[j][i].T-x[j-1][i].T)+x[j-2][i].T)/8.0)*dx;
651: }
653: residual -= (fE - fW + fN - fS);
655: return residual;
656: }
658: /*---------------------------------------------------------------------*/
659: /* computes the shear stress---used on the boundaries */
660: PETSC_STATIC_INLINE PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
661: /*---------------------------------------------------------------------*/
662: {
663: Parameter *param=user->param;
664: GridInfo *grid =user->grid;
665: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
666: PetscScalar uN, uS, wE, wW;
668: if (j<=grid->jlid || i<j || i==ilim || j==jlim) return EPS_ZERO;
670: if (ipos==CELL_CENTER) { /* on cell center */
672: wE = WInterp(x,i,j-1);
673: if (i==j) {
674: wW = param->sb;
675: uN = param->cb;
676: } else {
677: wW = WInterp(x,i-1,j-1);
678: uN = UInterp(x,i-1,j);
679: }
680: if (j==grid->jlid+1) uS = 0.0;
681: else uS = UInterp(x,i-1,j-1);
683: } else { /* on cell corner */
685: uN = x[j+1][i].u; uS = x[j][i].u;
686: wW = x[j][i].w; wE = x[j][i+1].w;
688: }
690: return (uN-uS)/grid->dz + (wE-wW)/grid->dx;
691: }
693: /*---------------------------------------------------------------------*/
694: /* computes the normal stress---used on the boundaries */
695: PETSC_STATIC_INLINE PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
696: /*---------------------------------------------------------------------*/
697: {
698: Parameter *param=user->param;
699: GridInfo *grid =user->grid;
700: PetscScalar dx = grid->dx, dz=grid->dz;
701: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
702: PetscScalar epsC =0.0, etaC, TC, uE, uW, pC, z_scale;
703: if (i<j || j<=grid->jlid) return EPS_ZERO;
705: ivisc=param->ivisc; z_scale = param->z_scale;
707: if (ipos==CELL_CENTER) { /* on cell center */
709: TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
710: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
711: etaC = Viscosity(TC,epsC,dz*j,param);
713: uW = x[j][i-1].u; uE = x[j][i].u;
714: pC = x[j][i].p;
716: } else { /* on cell corner */
717: if (i==ilim || j==jlim) return EPS_ZERO;
719: TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
720: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
721: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
723: if (i==j) uW = param->sb;
724: else uW = UInterp(x,i-1,j);
725: uE = UInterp(x,i,j); pC = PInterp(x,i,j);
726: }
728: return 2.0*etaC*(uE-uW)/dx - pC;
729: }
731: /*---------------------------------------------------------------------*/
732: /* computes the normal stress---used on the boundaries */
733: PETSC_STATIC_INLINE PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
734: /*---------------------------------------------------------------------*/
735: {
736: Parameter *param=user->param;
737: GridInfo *grid =user->grid;
738: PetscScalar dz =grid->dz;
739: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
740: PetscScalar epsC =0.0, etaC, TC;
741: PetscScalar pC, wN, wS, z_scale;
742: if (i<j || j<=grid->jlid) return EPS_ZERO;
744: ivisc=param->ivisc; z_scale = param->z_scale;
746: if (ipos==CELL_CENTER) { /* on cell center */
748: TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
749: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
750: etaC = Viscosity(TC,epsC,dz*j,param);
751: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
753: } else { /* on cell corner */
754: if ((i==ilim) || (j==jlim)) return EPS_ZERO;
756: TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
757: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
758: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
759: if (i==j) wN = param->sb;
760: else wN = WInterp(x,i,j);
761: wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
762: }
764: return 2.0*etaC*(wN-wS)/dz - pC;
765: }
767: /*---------------------------------------------------------------------*/
769: /*=====================================================================
770: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
771: =====================================================================*/
773: /*---------------------------------------------------------------------*/
774: /* initializes the problem parameters and checks for
775: command line changes */
776: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
777: /*---------------------------------------------------------------------*/
778: {
779: PetscErrorCode ierr, ierr_out=0;
780: PetscReal SEC_PER_YR = 3600.00*24.00*365.2500;
781: PetscReal alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
783: /* domain geometry */
784: param->slab_dip = 45.0;
785: param->width = 320.0; /* km */
786: param->depth = 300.0; /* km */
787: param->lid_depth = 35.0; /* km */
788: param->fault_depth = 35.0; /* km */
790: PetscOptionsGetReal(NULL,NULL,"-slab_dip",&(param->slab_dip),NULL);
791: PetscOptionsGetReal(NULL,NULL,"-width",&(param->width),NULL);
792: PetscOptionsGetReal(NULL,NULL,"-depth",&(param->depth),NULL);
793: PetscOptionsGetReal(NULL,NULL,"-lid_depth",&(param->lid_depth),NULL);
794: PetscOptionsGetReal(NULL,NULL,"-fault_depth",&(param->fault_depth),NULL);
796: param->slab_dip = param->slab_dip*PETSC_PI/180.0; /* radians */
798: /* grid information */
799: PetscOptionsGetInt(NULL,NULL, "-jfault",&(grid->jfault),NULL);
800: grid->ni = 82;
801: PetscOptionsGetInt(NULL,NULL, "-ni",&(grid->ni),NULL);
803: grid->dx = param->width/((PetscReal)(grid->ni-2)); /* km */
804: grid->dz = grid->dx*PetscTanReal(param->slab_dip); /* km */
805: grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/
806: param->depth = grid->dz*(grid->nj-2); /* km */
807: grid->inose = 0; /* gridpoints*/
808: PetscOptionsGetInt(NULL,NULL,"-inose",&(grid->inose),NULL);
809: grid->bx = DM_BOUNDARY_NONE;
810: grid->by = DM_BOUNDARY_NONE;
811: grid->stencil = DMDA_STENCIL_BOX;
812: grid->dof = 4;
813: grid->stencil_width = 2;
814: grid->mglevels = 1;
816: /* boundary conditions */
817: param->pv_analytic = PETSC_FALSE;
818: param->ibound = BC_NOSTRESS;
819: PetscOptionsGetInt(NULL,NULL,"-ibound",&(param->ibound),NULL);
821: /* physical constants */
822: param->slab_velocity = 5.0; /* cm/yr */
823: param->slab_age = 50.0; /* Ma */
824: param->lid_age = 50.0; /* Ma */
825: param->kappa = 0.7272e-6; /* m^2/sec */
826: param->potentialT = 1300.0; /* degrees C */
828: PetscOptionsGetReal(NULL,NULL,"-slab_velocity",&(param->slab_velocity),NULL);
829: PetscOptionsGetReal(NULL,NULL,"-slab_age",&(param->slab_age),NULL);
830: PetscOptionsGetReal(NULL,NULL,"-lid_age",&(param->lid_age),NULL);
831: PetscOptionsGetReal(NULL,NULL,"-kappa",&(param->kappa),NULL);
832: PetscOptionsGetReal(NULL,NULL,"-potentialT",&(param->potentialT),NULL);
834: /* viscosity */
835: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
836: param->eta0 = 1e24; /* Pa-s */
837: param->visc_cutoff = 0.0; /* factor of eta_0 */
838: param->continuation = 1.0;
840: /* constants for diffusion creep */
841: param->diffusion.A = 1.8e7; /* Pa-s */
842: param->diffusion.n = 1.0; /* dim'less */
843: param->diffusion.Estar = 375e3; /* J/mol */
844: param->diffusion.Vstar = 5e-6; /* m^3/mol */
846: /* constants for param->dislocationocation creep */
847: param->dislocation.A = 2.8969e4; /* Pa-s */
848: param->dislocation.n = 3.5; /* dim'less */
849: param->dislocation.Estar = 530e3; /* J/mol */
850: param->dislocation.Vstar = 14e-6; /* m^3/mol */
852: PetscOptionsGetInt(NULL,NULL, "-ivisc",&(param->ivisc),NULL);
853: PetscOptionsGetReal(NULL,NULL,"-visc_cutoff",&(param->visc_cutoff),NULL);
855: param->output_ivisc = param->ivisc;
857: PetscOptionsGetInt(NULL,NULL,"-output_ivisc",&(param->output_ivisc),NULL);
858: PetscOptionsGetReal(NULL,NULL,"-vstar",&(param->dislocation.Vstar),NULL);
860: /* output options */
861: param->quiet = PETSC_FALSE;
862: param->param_test = PETSC_FALSE;
864: PetscOptionsHasName(NULL,NULL,"-quiet",&(param->quiet));
865: PetscOptionsHasName(NULL,NULL,"-test",&(param->param_test));
866: PetscOptionsGetString(NULL,NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));
868: /* advection */
869: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
871: PetscOptionsGetInt(NULL,NULL,"-adv_scheme",&(param->adv_scheme),NULL);
873: /* misc. flags */
874: param->stop_solve = PETSC_FALSE;
875: param->interrupted = PETSC_FALSE;
876: param->kspmon = PETSC_FALSE;
877: param->toggle_kspmon = PETSC_FALSE;
879: /* derived parameters for slab angle */
880: param->sb = PetscSinReal(param->slab_dip);
881: param->cb = PetscCosReal(param->slab_dip);
882: param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
883: param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);
885: /* length, velocity and time scale for non-dimensionalization */
886: param->L = PetscMin(param->width,param->depth); /* km */
887: param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */
889: /* other unit conversions and derived parameters */
890: param->scaled_width = param->width/param->L; /* dim'less */
891: param->scaled_depth = param->depth/param->L; /* dim'less */
892: param->lid_depth = param->lid_depth/param->L; /* dim'less */
893: param->fault_depth = param->fault_depth/param->L; /* dim'less */
894: grid->dx = grid->dx/param->L; /* dim'less */
895: grid->dz = grid->dz/param->L; /* dim'less */
896: grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */
897: grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */
898: param->lid_depth = grid->jlid*grid->dz; /* dim'less */
899: param->fault_depth = grid->jfault*grid->dz; /* dim'less */
900: grid->corner = grid->jlid+1; /* gridcells */
901: param->peclet = param->V /* m/sec */
902: * param->L*1000.0 /* m */
903: / param->kappa; /* m^2/sec */
904: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
905: param->skt = PetscSqrtReal(param->kappa*param->slab_age*SEC_PER_YR);
906: PetscOptionsGetReal(NULL,NULL,"-peclet",&(param->peclet),NULL);
908: return ierr_out;
909: }
911: /*---------------------------------------------------------------------*/
912: /* prints a report of the problem parameters to stdout */
913: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
914: /*---------------------------------------------------------------------*/
915: {
916: PetscErrorCode ierr, ierr_out=0;
917: char date[30];
919: PetscGetDate(date,30);
921: if (!(param->quiet)) {
922: PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
923: /* PetscPrintf(PETSC_COMM_WORLD," %s\n",&(date[0]));*/
925: PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
926: PetscPrintf(PETSC_COMM_WORLD," Width = %g km, Depth = %g km\n",(double)param->width,(double)param->depth);
927: PetscPrintf(PETSC_COMM_WORLD," Slab dip = %g degrees, Slab velocity = %g cm/yr\n",(double)(param->slab_dip*180.0/PETSC_PI),(double)param->slab_velocity);
928: 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);
930: PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
931: PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %D, %D [dx,dz] = %g, %g km\n",grid->ni,grid->nj,(double)grid->dx*param->L,(double)(grid->dz*param->L));
932: PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);
933: PetscPrintf(PETSC_COMM_WORLD," Pe = %g\n",(double)param->peclet);
935: PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
936: if (param->ivisc==VISC_CONST) {
937: PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");
938: if (param->pv_analytic) {
939: PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");
940: }
941: } else if (param->ivisc==VISC_DIFN) {
942: PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");
943: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
944: } else if (param->ivisc==VISC_DISL) {
945: PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");
946: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
947: } else if (param->ivisc==VISC_FULL) {
948: PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");
949: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
950: } else {
951: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
952: ierr_out = 1;
953: }
955: PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
956: if (param->ibound==BC_ANALYTIC) {
957: PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");
958: } else if (param->ibound==BC_NOSTRESS) {
959: PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");
960: } else if (param->ibound==BC_EXPERMNT) {
961: PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");
962: } else {
963: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
964: ierr_out = 1;
965: }
967: if (param->output_to_file)
968: #if defined(PETSC_HAVE_MATLAB_ENGINE)
969: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);
970: #else
971: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);
972: #endif
973: if (param->output_ivisc != param->ivisc) {
974: PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);
975: }
977: PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
978: }
979: if (param->param_test) PetscEnd();
980: return ierr_out;
981: }
983: /* ------------------------------------------------------------------- */
984: /* generates an inital guess using the analytic solution for isoviscous
985: corner flow */
986: PetscErrorCode Initialize(DM da)
987: /* ------------------------------------------------------------------- */
988: {
989: AppCtx *user;
990: Parameter *param;
991: GridInfo *grid;
992: PetscInt i,j,is,js,im,jm;
994: Field **x;
995: Vec Xguess;
997: /* Get the fine grid */
998: DMGetApplicationContext(da,&user);
999: Xguess = user->Xguess;
1000: param = user->param;
1001: grid = user->grid;
1002: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1003: DMDAVecGetArray(da,Xguess,(void**)&x);
1005: /* Compute initial guess */
1006: for (j=js; j<js+jm; j++) {
1007: for (i=is; i<is+im; i++) {
1008: if (i<j) x[j][i].u = param->cb;
1009: else if (j<=grid->jlid) x[j][i].u = 0.0;
1010: else x[j][i].u = HorizVelocity(i,j,user);
1012: if (i<=j) x[j][i].w = param->sb;
1013: else if (j<=grid->jlid) x[j][i].w = 0.0;
1014: else x[j][i].w = VertVelocity(i,j,user);
1016: if (i<j || j<=grid->jlid) x[j][i].p = 0.0;
1017: else x[j][i].p = Pressure(i,j,user);
1019: x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1020: }
1021: }
1023: /* Restore x to Xguess */
1024: DMDAVecRestoreArray(da,Xguess,(void**)&x);
1026: return 0;
1027: }
1029: /*---------------------------------------------------------------------*/
1030: /* controls output to a file */
1031: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1032: /*---------------------------------------------------------------------*/
1033: {
1034: AppCtx *user;
1035: Parameter *param;
1036: GridInfo *grid;
1037: PetscInt ivt;
1039: PetscMPIInt rank;
1040: PetscViewer viewer;
1041: Vec res, pars;
1042: MPI_Comm comm;
1043: DM da;
1045: SNESGetDM(snes,&da);
1046: DMGetApplicationContext(da,&user);
1047: param = user->param;
1048: grid = user->grid;
1049: ivt = param->ivisc;
1051: param->ivisc = param->output_ivisc;
1053: /* compute final residual and final viscosity/strain rate fields */
1054: SNESGetFunction(snes, &res, NULL, NULL);
1055: ViscosityField(da, user->x, user->Xguess);
1057: /* get the communicator and the rank of the processor */
1058: PetscObjectGetComm((PetscObject)snes, &comm);
1059: MPI_Comm_rank(comm, &rank);
1061: if (param->output_to_file) { /* send output to binary file */
1062: VecCreate(comm, &pars);
1063: if (!rank) { /* on processor 0 */
1064: VecSetSizes(pars, 20, PETSC_DETERMINE);
1065: VecSetFromOptions(pars);
1066: VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1067: VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1068: VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1069: VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1070: VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1071: VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1072: /* skipped 6 intentionally */
1073: VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1074: VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1075: VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1076: VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1077: VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1078: VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1079: VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1080: VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1081: VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1082: VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1083: } else { /* on some other processor */
1084: VecSetSizes(pars, 0, PETSC_DETERMINE);
1085: VecSetFromOptions(pars);
1086: }
1087: VecAssemblyBegin(pars); VecAssemblyEnd(pars);
1089: /* create viewer */
1090: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1091: PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1092: #else
1093: PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1094: #endif
1096: /* send vectors to viewer */
1097: PetscObjectSetName((PetscObject)res,"res");
1098: VecView(res,viewer);
1099: PetscObjectSetName((PetscObject)user->x,"out");
1100: VecView(user->x, viewer);
1101: PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1102: VecView(user->Xguess, viewer);
1103: StressField(da); /* compute stress fields */
1104: PetscObjectSetName((PetscObject)(user->Xguess),"str");
1105: VecView(user->Xguess, viewer);
1106: PetscObjectSetName((PetscObject)pars,"par");
1107: VecView(pars, viewer);
1109: /* destroy viewer and vector */
1110: PetscViewerDestroy(&viewer);
1111: VecDestroy(&pars);
1112: }
1114: param->ivisc = ivt;
1115: return 0;
1116: }
1118: /* ------------------------------------------------------------------- */
1119: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1120: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1121: /* ------------------------------------------------------------------- */
1122: {
1123: AppCtx *user;
1124: Parameter *param;
1125: GridInfo *grid;
1126: Vec localX;
1127: Field **v, **x;
1128: PetscReal eps, /* dx,*/ dz, T, epsC, TC;
1129: PetscInt i,j,is,js,im,jm,ilim,jlim,ivt;
1133: DMGetApplicationContext(da,&user);
1134: param = user->param;
1135: grid = user->grid;
1136: ivt = param->ivisc;
1137: param->ivisc = param->output_ivisc;
1139: DMGetLocalVector(da, &localX);
1140: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1141: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1142: DMDAVecGetArray(da,localX,(void**)&x);
1143: DMDAVecGetArray(da,V,(void**)&v);
1145: /* Parameters */
1146: /* dx = grid->dx; */ dz = grid->dz;
1148: ilim = grid->ni-1; jlim = grid->nj-1;
1150: /* Compute real temperature, strain rate and viscosity */
1151: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1152: for (j=js; j<js+jm; j++) {
1153: for (i=is; i<is+im; i++) {
1154: T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*param->z_scale));
1155: if (i<ilim && j<jlim) {
1156: TC = PetscRealPart(param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*param->z_scale));
1157: } else {
1158: TC = T;
1159: }
1160: eps = PetscRealPart((CalcSecInv(x,i,j,CELL_CENTER,user)));
1161: epsC = PetscRealPart(CalcSecInv(x,i,j,CELL_CORNER,user));
1163: v[j][i].u = eps;
1164: v[j][i].w = epsC;
1165: v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1166: v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1167: }
1168: }
1169: DMDAVecRestoreArray(da,V,(void**)&v);
1170: DMDAVecRestoreArray(da,localX,(void**)&x);
1171: DMRestoreLocalVector(da, &localX);
1173: param->ivisc = ivt;
1174: return(0);
1175: }
1177: /* ------------------------------------------------------------------- */
1178: /* post-processing: compute stress everywhere */
1179: PetscErrorCode StressField(DM da)
1180: /* ------------------------------------------------------------------- */
1181: {
1182: AppCtx *user;
1183: PetscInt i,j,is,js,im,jm;
1185: Vec locVec;
1186: Field **x, **y;
1188: DMGetApplicationContext(da,&user);
1190: /* Get the fine grid of Xguess and X */
1191: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1192: DMDAVecGetArray(da,user->Xguess,(void**)&x);
1194: DMGetLocalVector(da, &locVec);
1195: DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);
1196: DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);
1197: DMDAVecGetArray(da,locVec,(void**)&y);
1199: /* Compute stress on the corner points */
1200: for (j=js; j<js+jm; j++) {
1201: for (i=is; i<is+im; i++) {
1202: x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1203: x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1204: x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1205: x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1206: }
1207: }
1209: /* Restore the fine grid of Xguess and X */
1210: DMDAVecRestoreArray(da,user->Xguess,(void**)&x);
1211: DMDAVecRestoreArray(da,locVec,(void**)&y);
1212: DMRestoreLocalVector(da, &locVec);
1213: return 0;
1214: }
1216: /*=====================================================================
1217: UTILITY FUNCTIONS
1218: =====================================================================*/
1220: /*---------------------------------------------------------------------*/
1221: /* returns the velocity of the subducting slab and handles fault nodes
1222: for BC */
1223: PETSC_STATIC_INLINE PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1224: /*---------------------------------------------------------------------*/
1225: {
1226: Parameter *param = user->param;
1227: GridInfo *grid = user->grid;
1229: if (c=='U' || c=='u') {
1230: if (i<j-1) return param->cb;
1231: else if (j<=grid->jfault) return 0.0;
1232: else return param->cb;
1234: } else {
1235: if (i<j) return param->sb;
1236: else if (j<=grid->jfault) return 0.0;
1237: else return param->sb;
1238: }
1239: }
1241: /*---------------------------------------------------------------------*/
1242: /* solution to diffusive half-space cooling model for BC */
1243: PETSC_STATIC_INLINE PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1244: /*---------------------------------------------------------------------*/
1245: {
1246: Parameter *param = user->param;
1247: PetscScalar z;
1248: if (plate==PLATE_LID) z = (j-0.5)*user->grid->dz;
1249: else z = (j-0.5)*user->grid->dz*param->cb; /* PLATE_SLAB */
1250: #if defined(PETSC_HAVE_ERF)
1251: return (PetscReal)(erf((double)PetscRealPart(z*param->L/2.0/param->skt)));
1252: #else
1253: SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
1254: #endif
1255: }
1258: /* ------------------------------------------------------------------- */
1259: /* utility function */
1260: PetscBool OptionsHasName(const char pre[],const char name[])
1261: /* ------------------------------------------------------------------- */
1262: {
1263: PetscBool retval;
1265: PetscOptionsHasName(NULL,pre,name,&retval);CHKERRABORT(PETSC_COMM_WORLD,ierr);
1266: return retval;
1267: }
1269: /*=====================================================================
1270: INTERACTIVE SIGNAL HANDLING
1271: =====================================================================*/
1273: /* ------------------------------------------------------------------- */
1274: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1275: /* ------------------------------------------------------------------- */
1276: {
1277: AppCtx *user = (AppCtx*) ctx;
1278: Parameter *param = user->param;
1279: KSP ksp;
1283: if (param->interrupted) {
1284: param->interrupted = PETSC_FALSE;
1285: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1286: *reason = SNES_CONVERGED_FNORM_ABS;
1287: return(0);
1288: } else if (param->toggle_kspmon) {
1289: param->toggle_kspmon = PETSC_FALSE;
1291: SNESGetKSP(snes, &ksp);
1293: if (param->kspmon) {
1294: KSPMonitorCancel(ksp);
1296: param->kspmon = PETSC_FALSE;
1297: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1298: } else {
1299: PetscViewerAndFormat *vf;
1300: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
1301: KSPMonitorSet(ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSingularValue,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
1303: param->kspmon = PETSC_TRUE;
1304: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1305: }
1306: }
1307: PetscFunctionReturn(SNESConvergedDefault(snes,it,xnorm,snorm,fnorm,reason,ctx));
1308: }
1310: /* ------------------------------------------------------------------- */
1311: #include <signal.h>
1312: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1313: /* ------------------------------------------------------------------- */
1314: {
1315: AppCtx *user = (AppCtx*) ctx;
1316: Parameter *param = user->param;
1318: if (signum == SIGILL) {
1319: param->toggle_kspmon = PETSC_TRUE;
1320: #if !defined(PETSC_MISSING_SIGCONT)
1321: } else if (signum == SIGCONT) {
1322: param->interrupted = PETSC_TRUE;
1323: #endif
1324: #if !defined(PETSC_MISSING_SIGURG)
1325: } else if (signum == SIGURG) {
1326: param->stop_solve = PETSC_TRUE;
1327: #endif
1328: }
1329: return 0;
1330: }
1332: /*---------------------------------------------------------------------*/
1333: /* main call-back function that computes the processor-local piece
1334: of the residual */
1335: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
1336: /*---------------------------------------------------------------------*/
1337: {
1338: AppCtx *user = (AppCtx*)ptr;
1339: Parameter *param = user->param;
1340: GridInfo *grid = user->grid;
1341: PetscScalar mag_w, mag_u;
1342: PetscInt i,j,mx,mz,ilim,jlim;
1343: PetscInt is,ie,js,je,ibound; /* ,ivisc */
1346: /* Define global and local grid parameters */
1347: mx = info->mx; mz = info->my;
1348: ilim = mx-1; jlim = mz-1;
1349: is = info->xs; ie = info->xs+info->xm;
1350: js = info->ys; je = info->ys+info->ym;
1352: /* Define geometric and numeric parameters */
1353: /* ivisc = param->ivisc; */ ibound = param->ibound;
1355: for (j=js; j<je; j++) {
1356: for (i=is; i<ie; i++) {
1358: /************* X-MOMENTUM/VELOCITY *************/
1359: if (i<j) f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
1360: else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1361: /* in the lithospheric lid */
1362: f[j][i].u = x[j][i].u - 0.0;
1363: } else if (i==ilim) {
1364: /* on the right side boundary */
1365: if (ibound==BC_ANALYTIC) {
1366: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1367: } else {
1368: f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1369: }
1371: } else if (j==jlim) {
1372: /* on the bottom boundary */
1373: if (ibound==BC_ANALYTIC) {
1374: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1375: } else if (ibound==BC_NOSTRESS) {
1376: f[j][i].u = XMomentumResidual(x,i,j,user);
1377: } else {
1378: /* experimental boundary condition */
1379: }
1381: } else {
1382: /* in the mantle wedge */
1383: f[j][i].u = XMomentumResidual(x,i,j,user);
1384: }
1386: /************* Z-MOMENTUM/VELOCITY *************/
1387: if (i<=j) {
1388: f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);
1390: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1391: /* in the lithospheric lid */
1392: f[j][i].w = x[j][i].w - 0.0;
1394: } else if (j==jlim) {
1395: /* on the bottom boundary */
1396: if (ibound==BC_ANALYTIC) {
1397: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1398: } else {
1399: f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1400: }
1402: } else if (i==ilim) {
1403: /* on the right side boundary */
1404: if (ibound==BC_ANALYTIC) {
1405: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1406: } else if (ibound==BC_NOSTRESS) {
1407: f[j][i].w = ZMomentumResidual(x,i,j,user);
1408: } else {
1409: /* experimental boundary condition */
1410: }
1412: } else {
1413: /* in the mantle wedge */
1414: f[j][i].w = ZMomentumResidual(x,i,j,user);
1415: }
1417: /************* CONTINUITY/PRESSURE *************/
1418: if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1419: /* in the lid or slab */
1420: f[j][i].p = x[j][i].p;
1422: } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
1423: /* on an analytic boundary */
1424: f[j][i].p = x[j][i].p - Pressure(i,j,user);
1426: } else {
1427: /* in the mantle wedge */
1428: f[j][i].p = ContinuityResidual(x,i,j,user);
1429: }
1431: /************* TEMPERATURE *************/
1432: if (j==0) {
1433: /* on the surface */
1434: f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0);
1436: } else if (i==0) {
1437: /* slab inflow boundary */
1438: f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user);
1440: } else if (i==ilim) {
1441: /* right side boundary */
1442: mag_u = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j][i-1].u)/param->cb,1.0),0.0)), 5);
1443: 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);
1445: } else if (j==jlim) {
1446: /* bottom boundary */
1447: mag_w = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j-1][i].w)/param->sb,1.0),0.0)), 5);
1448: f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);
1450: } else {
1451: /* in the mantle wedge */
1452: f[j][i].T = EnergyResidual(x,i,j,user);
1453: }
1454: }
1455: }
1456: return(0);
1457: }