Actual source code: ex30.c
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}$-----
54: #include <petscsnes.h>
55: #include <petscdm.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: DMBoundaryType 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*);
123: /*-----------------------------------------------------------------------*/
124: int main(int argc,char **argv)
125: /*-----------------------------------------------------------------------*/
126: {
127: SNES snes;
128: AppCtx *user; /* user-defined work context */
129: Parameter param;
130: GridInfo grid;
131: PetscInt nits;
132: MPI_Comm comm;
133: DM da;
135: PetscInitialize(&argc,&argv,(char*)0,help);
136: PetscOptionsSetValue(NULL,"-file","ex30_output");
137: PetscOptionsSetValue(NULL,"-snes_monitor_short",NULL);
138: PetscOptionsSetValue(NULL,"-snes_max_it","20");
139: PetscOptionsSetValue(NULL,"-ksp_max_it","1500");
140: PetscOptionsSetValue(NULL,"-ksp_gmres_restart","300");
141: PetscOptionsInsert(NULL,&argc,&argv,NULL);
143: comm = PETSC_COMM_WORLD;
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Set up the problem parameters.
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: SetParams(¶m,&grid);
149: ReportParams(¶m,&grid);
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: SNESCreate(comm,&snes);
154: DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
155: DMSetFromOptions(da);
156: DMSetUp(da);
157: SNESSetDM(snes,da);
158: DMDASetFieldName(da,0,"x-velocity");
159: DMDASetFieldName(da,1,"y-velocity");
160: DMDASetFieldName(da,2,"pressure");
161: DMDASetFieldName(da,3,"temperature");
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Create user context, set problem data, create vector data structures.
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: PetscNew(&user);
167: user->param = ¶m;
168: user->grid = &grid;
169: DMSetApplicationContext(da,user);
170: DMCreateGlobalVector(da,&(user->Xguess));
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Set up the SNES solver with callback functions.
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,(void*)user);
176: SNESSetFromOptions(snes);
178: SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,NULL);
179: PetscPushSignalHandler(InteractiveHandler,(void*)user);
181: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
182: Initialize and solve the nonlinear system
183: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184: Initialize(da);
185: UpdateSolution(snes,user,&nits);
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Output variables.
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: DoOutput(snes,nits);
192: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: Free work space.
194: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195: VecDestroy(&user->Xguess);
196: VecDestroy(&user->x);
197: PetscFree(user);
198: SNESDestroy(&snes);
199: DMDestroy(&da);
200: PetscPopSignalHandler();
201: PetscFinalize();
202: return 0;
203: }
205: /*=====================================================================
206: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
207: =====================================================================*/
209: /*---------------------------------------------------------------------*/
210: /* manages solve: adaptive continuation method */
211: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
212: {
213: KSP ksp;
214: PC pc;
215: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
216: Parameter *param = user->param;
217: PetscReal cont_incr=0.3;
218: PetscInt its;
219: PetscBool q = PETSC_FALSE;
220: DM dm;
223: SNESGetDM(snes,&dm);
224: DMCreateGlobalVector(dm,&user->x);
225: SNESGetKSP(snes,&ksp);
226: KSPGetPC(ksp, &pc);
227: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
229: *nits=0;
231: /* Isoviscous solve */
232: if (param->ivisc == VISC_CONST && !param->stop_solve) {
233: param->ivisc = VISC_CONST;
235: SNESSolve(snes,0,user->x);
236: SNESGetConvergedReason(snes,&reason);
237: SNESGetIterationNumber(snes,&its);
238: *nits += its;
239: VecCopy(user->x,user->Xguess);
240: if (param->stop_solve) goto done;
241: }
243: /* Olivine diffusion creep */
244: if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
245: if (!q) PetscPrintf(PETSC_COMM_WORLD,"Computing Variable Viscosity Solution\n");
247: /* continuation method on viscosity cutoff */
248: for (param->continuation=0.0;; param->continuation+=cont_incr) {
249: if (!q) PetscPrintf(PETSC_COMM_WORLD," Continuation parameter = %g\n", (double)param->continuation);
251: /* solve the non-linear system */
252: VecCopy(user->Xguess,user->x);
253: SNESSolve(snes,0,user->x);
254: SNESGetConvergedReason(snes,&reason);
255: SNESGetIterationNumber(snes,&its);
256: *nits += its;
257: if (!q) PetscPrintf(PETSC_COMM_WORLD," SNES iterations: %D, Cumulative: %D\n", its, *nits);
258: if (param->stop_solve) goto done;
260: if (reason<0) {
261: /* NOT converged */
262: cont_incr = -PetscAbsReal(cont_incr)/2.0;
263: if (PetscAbsReal(cont_incr)<0.01) goto done;
265: } else {
266: /* converged */
267: VecCopy(user->x,user->Xguess);
268: if (param->continuation >= 1.0) goto done;
269: if (its<=3) cont_incr = 0.30001;
270: else if (its<=8) cont_incr = 0.15001;
271: else cont_incr = 0.10001;
273: if (param->continuation+cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
274: } /* endif reason<0 */
275: }
276: }
277: done:
278: if (param->stop_solve && !q) PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");
279: if (reason<0 && !q) PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");
280: return 0;
281: }
283: /*=====================================================================
284: PHYSICS FUNCTIONS (compute the discrete residual)
285: =====================================================================*/
287: /*---------------------------------------------------------------------*/
288: static inline PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
289: /*---------------------------------------------------------------------*/
290: {
291: return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
292: }
294: /*---------------------------------------------------------------------*/
295: static inline PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
296: /*---------------------------------------------------------------------*/
297: {
298: return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
299: }
301: /*---------------------------------------------------------------------*/
302: static inline PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
303: /*---------------------------------------------------------------------*/
304: {
305: return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
306: }
308: /*---------------------------------------------------------------------*/
309: static inline PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
310: /*---------------------------------------------------------------------*/
311: {
312: return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
313: }
315: /*---------------------------------------------------------------------*/
316: /* isoviscous analytic solution for IC */
317: static inline PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
318: /*---------------------------------------------------------------------*/
319: {
320: Parameter *param = user->param;
321: GridInfo *grid = user->grid;
322: PetscScalar st, ct, th, c=param->c, d=param->d;
323: PetscReal x, z,r;
325: x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
326: r = PetscSqrtReal(x*x+z*z);
327: st = z/r;
328: ct = x/r;
329: th = PetscAtanReal(z/x);
330: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
331: }
333: /*---------------------------------------------------------------------*/
334: /* isoviscous analytic solution for IC */
335: static inline PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
336: /*---------------------------------------------------------------------*/
337: {
338: Parameter *param = user->param;
339: GridInfo *grid = user->grid;
340: PetscScalar st, ct, th, c=param->c, d=param->d;
341: PetscReal x, z, r;
343: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz;
344: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r; th = PetscAtanReal(z/x);
345: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
346: }
348: /*---------------------------------------------------------------------*/
349: /* isoviscous analytic solution for IC */
350: static inline PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
351: /*---------------------------------------------------------------------*/
352: {
353: Parameter *param = user->param;
354: GridInfo *grid = user->grid;
355: PetscScalar x, z, r, st, ct, c=param->c, d=param->d;
357: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
358: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r;
359: return (-2.0*(c*ct-d*st)/r);
360: }
362: /* computes the second invariant of the strain rate tensor */
363: static inline PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
364: /*---------------------------------------------------------------------*/
365: {
366: Parameter *param = user->param;
367: GridInfo *grid = user->grid;
368: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
369: PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
370: PetscScalar eps11, eps12, eps22;
372: if (i<j) return EPS_ZERO;
373: if (i==ilim) i--;
374: if (j==jlim) j--;
376: if (ipos==CELL_CENTER) { /* on cell center */
377: if (j<=grid->jlid) return EPS_ZERO;
379: uE = x[j][i].u; uW = x[j][i-1].u;
380: wN = x[j][i].w; wS = x[j-1][i].w;
381: wE = WInterp(x,i,j-1);
382: if (i==j) {
383: uN = param->cb; wW = param->sb;
384: } else {
385: uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1);
386: }
388: if (j==grid->jlid+1) uS = 0.0;
389: else uS = UInterp(x,i-1,j-1);
391: } else { /* on CELL_CORNER */
392: if (j<grid->jlid) return EPS_ZERO;
394: uN = x[j+1][i].u; uS = x[j][i].u;
395: wE = x[j][i+1].w; wW = x[j][i].w;
396: if (i==j) {
397: wN = param->sb;
398: uW = param->cb;
399: } else {
400: wN = WInterp(x,i,j);
401: uW = UInterp(x,i-1,j);
402: }
404: if (j==grid->jlid) {
405: uE = 0.0; uW = 0.0;
406: uS = -uN;
407: wS = -wN;
408: } else {
409: uE = UInterp(x,i,j);
410: wS = WInterp(x,i,j-1);
411: }
412: }
414: eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz;
415: eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx);
417: return PetscSqrtReal(0.5*(eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22));
418: }
420: /*---------------------------------------------------------------------*/
421: /* computes the shear viscosity */
422: static inline PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
423: /*---------------------------------------------------------------------*/
424: {
425: PetscReal result =0.0;
426: ViscParam difn =param->diffusion, disl=param->dislocation;
427: PetscInt iVisc =param->ivisc;
428: PetscScalar eps_scale=param->V/(param->L*1000.0);
429: PetscScalar strain_power, v1, v2, P;
430: PetscScalar rho_g = 32340.0, R=8.3144;
432: P = rho_g*(z*param->L*1000.0); /* Pa */
434: if (iVisc==VISC_CONST) {
435: /* constant viscosity */
436: return 1.0;
437: } else if (iVisc==VISC_DIFN) {
438: /* diffusion creep rheology */
439: result = PetscRealPart((difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0));
440: } else if (iVisc==VISC_DISL) {
441: /* dislocation creep rheology */
442: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
444: result = PetscRealPart(disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0);
445: } else if (iVisc==VISC_FULL) {
446: /* dislocation/diffusion creep rheology */
447: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
449: v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
450: v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
452: result = PetscRealPart(1.0/(1.0/v1 + 1.0/v2));
453: }
455: /* max viscosity is param->eta0 */
456: result = PetscMin(result, 1.0);
457: /* min viscosity is param->visc_cutoff */
458: result = PetscMax(result, param->visc_cutoff);
459: /* continuation method */
460: result = PetscPowReal(result,param->continuation);
461: return result;
462: }
464: /*---------------------------------------------------------------------*/
465: /* computes the residual of the x-component of eqn (1) above */
466: static inline PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
467: /*---------------------------------------------------------------------*/
468: {
469: Parameter *param=user->param;
470: GridInfo *grid =user->grid;
471: PetscScalar dx = grid->dx, dz=grid->dz;
472: PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
473: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
474: PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
475: PetscInt jlim = grid->nj-1;
477: z_scale = param->z_scale;
479: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
480: TS = param->potentialT * TInterp(x,i,j-1) * PetscExpScalar((j-1.0)*dz*z_scale);
481: if (j==jlim) TN = TS;
482: else TN = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
483: TW = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
484: TE = param->potentialT * x[j][i+1].T * PetscExpScalar((j-0.5)*dz*z_scale);
485: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
486: epsN = CalcSecInv(x,i,j, CELL_CORNER,user);
487: epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
488: epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
489: epsW = CalcSecInv(x,i,j, CELL_CENTER,user);
490: }
491: }
492: etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
493: etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
494: etaW = Viscosity(TW,epsW,dz*j,param);
495: etaE = Viscosity(TE,epsE,dz*j,param);
497: dPdx = (x[j][i+1].p - x[j][i].p)/dx;
498: if (j==jlim) dudzN = etaN * (x[j][i].w - x[j][i+1].w)/dx;
499: else dudzN = etaN * (x[j+1][i].u - x[j][i].u) /dz;
500: dudzS = etaS * (x[j][i].u - x[j-1][i].u)/dz;
501: dudxE = etaE * (x[j][i+1].u - x[j][i].u) /dx;
502: dudxW = etaW * (x[j][i].u - x[j][i-1].u)/dx;
504: residual = -dPdx /* X-MOMENTUM EQUATION*/
505: +(dudxE - dudxW)/dx
506: +(dudzN - dudzS)/dz;
508: if (param->ivisc!=VISC_CONST) {
509: dwdxN = etaN * (x[j][i+1].w - x[j][i].w) /dx;
510: dwdxS = etaS * (x[j-1][i+1].w - x[j-1][i].w)/dx;
512: residual += (dudxE - dudxW)/dx + (dwdxN - dwdxS)/dz;
513: }
515: return residual;
516: }
518: /*---------------------------------------------------------------------*/
519: /* computes the residual of the z-component of eqn (1) above */
520: static inline PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
521: /*---------------------------------------------------------------------*/
522: {
523: Parameter *param=user->param;
524: GridInfo *grid =user->grid;
525: PetscScalar dx = grid->dx, dz=grid->dz;
526: 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;
527: PetscScalar TE =0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
528: PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
529: PetscInt ilim = grid->ni-1;
531: /* geometric and other parameters */
532: z_scale = param->z_scale;
534: /* viscosity */
535: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
536: TN = param->potentialT * x[j+1][i].T * PetscExpScalar((j+0.5)*dz*z_scale);
537: TS = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
538: TW = param->potentialT * TInterp(x,i-1,j) * PetscExpScalar(j*dz*z_scale);
539: if (i==ilim) TE = TW;
540: else TE = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
541: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
542: epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
543: epsS = CalcSecInv(x,i,j, CELL_CENTER,user);
544: epsE = CalcSecInv(x,i,j, CELL_CORNER,user);
545: epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
546: }
547: }
548: etaN = Viscosity(TN,epsN,dz*(j+1.0),param);
549: etaS = Viscosity(TS,epsS,dz*(j+0.0),param);
550: etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
551: etaE = Viscosity(TE,epsE,dz*(j+0.5),param);
553: dPdz = (x[j+1][i].p - x[j][i].p)/dz;
554: dwdzN = etaN * (x[j+1][i].w - x[j][i].w)/dz;
555: dwdzS = etaS * (x[j][i].w - x[j-1][i].w)/dz;
556: if (i==ilim) dwdxE = etaE * (x[j][i].u - x[j+1][i].u)/dz;
557: else dwdxE = etaE * (x[j][i+1].w - x[j][i].w) /dx;
558: dwdxW = 2.0*etaW * (x[j][i].w - x[j][i-1].w)/dx;
560: /* Z-MOMENTUM */
561: residual = -dPdz /* constant viscosity terms */
562: +(dwdzN - dwdzS)/dz
563: +(dwdxE - dwdxW)/dx;
565: if (param->ivisc!=VISC_CONST) {
566: dudzE = etaE * (x[j+1][i].u - x[j][i].u)/dz;
567: dudzW = etaW * (x[j+1][i-1].u - x[j][i-1].u)/dz;
569: residual += (dwdzN - dwdzS)/dz + (dudzE - dudzW)/dx;
570: }
572: return residual;
573: }
575: /*---------------------------------------------------------------------*/
576: /* computes the residual of eqn (2) above */
577: static inline PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
578: /*---------------------------------------------------------------------*/
579: {
580: GridInfo *grid =user->grid;
581: PetscScalar uE,uW,wN,wS,dudx,dwdz;
583: uW = x[j][i-1].u; uE = x[j][i].u; dudx = (uE - uW)/grid->dx;
584: wS = x[j-1][i].w; wN = x[j][i].w; dwdz = (wN - wS)/grid->dz;
586: return dudx + dwdz;
587: }
589: /*---------------------------------------------------------------------*/
590: /* computes the residual of eqn (3) above */
591: static inline PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
592: /*---------------------------------------------------------------------*/
593: {
594: Parameter *param=user->param;
595: GridInfo *grid =user->grid;
596: PetscScalar dx = grid->dx, dz=grid->dz;
597: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
598: PetscScalar TE, TN, TS, TW, residual;
599: PetscScalar uE,uW,wN,wS;
600: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;
602: dTdzN = (x[j+1][i].T - x[j][i].T) /dz;
603: dTdzS = (x[j][i].T - x[j-1][i].T)/dz;
604: dTdxE = (x[j][i+1].T - x[j][i].T) /dx;
605: dTdxW = (x[j][i].T - x[j][i-1].T)/dx;
607: residual = ((dTdzN - dTdzS)/dz + /* diffusion term */
608: (dTdxE - dTdxW)/dx)*dx*dz/param->peclet;
610: if (j<=jlid && i>=j) {
611: /* don't advect in the lid */
612: return residual;
613: } else if (i<j) {
614: /* beneath the slab sfc */
615: uW = uE = param->cb;
616: wS = wN = param->sb;
617: } else {
618: /* advect in the slab and wedge */
619: uW = x[j][i-1].u; uE = x[j][i].u;
620: wS = x[j-1][i].w; wN = x[j][i].w;
621: }
623: if (param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1) {
624: /* finite volume advection */
625: TS = (x[j][i].T + x[j-1][i].T)/2.0;
626: TN = (x[j][i].T + x[j+1][i].T)/2.0;
627: TE = (x[j][i].T + x[j][i+1].T)/2.0;
628: TW = (x[j][i].T + x[j][i-1].T)/2.0;
629: fN = wN*TN*dx; fS = wS*TS*dx;
630: fE = uE*TE*dz; fW = uW*TW*dz;
632: } else {
633: /* Fromm advection scheme */
634: 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
635: - 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;
636: 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
637: - 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;
638: 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
639: - 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;
640: 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
641: - 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;
642: }
644: residual -= (fE - fW + fN - fS);
646: return residual;
647: }
649: /*---------------------------------------------------------------------*/
650: /* computes the shear stress---used on the boundaries */
651: static inline PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
652: /*---------------------------------------------------------------------*/
653: {
654: Parameter *param=user->param;
655: GridInfo *grid =user->grid;
656: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
657: PetscScalar uN, uS, wE, wW;
659: if (j<=grid->jlid || i<j || i==ilim || j==jlim) return EPS_ZERO;
661: if (ipos==CELL_CENTER) { /* on cell center */
663: wE = WInterp(x,i,j-1);
664: if (i==j) {
665: wW = param->sb;
666: uN = param->cb;
667: } else {
668: wW = WInterp(x,i-1,j-1);
669: uN = UInterp(x,i-1,j);
670: }
671: if (j==grid->jlid+1) uS = 0.0;
672: else uS = UInterp(x,i-1,j-1);
674: } else { /* on cell corner */
676: uN = x[j+1][i].u; uS = x[j][i].u;
677: wW = x[j][i].w; wE = x[j][i+1].w;
679: }
681: return (uN-uS)/grid->dz + (wE-wW)/grid->dx;
682: }
684: /*---------------------------------------------------------------------*/
685: /* computes the normal stress---used on the boundaries */
686: static inline PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
687: /*---------------------------------------------------------------------*/
688: {
689: Parameter *param=user->param;
690: GridInfo *grid =user->grid;
691: PetscScalar dx = grid->dx, dz=grid->dz;
692: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
693: PetscScalar epsC =0.0, etaC, TC, uE, uW, pC, z_scale;
694: if (i<j || j<=grid->jlid) return EPS_ZERO;
696: ivisc=param->ivisc; z_scale = param->z_scale;
698: if (ipos==CELL_CENTER) { /* on cell center */
700: TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
701: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
702: etaC = Viscosity(TC,epsC,dz*j,param);
704: uW = x[j][i-1].u; uE = x[j][i].u;
705: pC = x[j][i].p;
707: } else { /* on cell corner */
708: if (i==ilim || j==jlim) return EPS_ZERO;
710: TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
711: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
712: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
714: if (i==j) uW = param->sb;
715: else uW = UInterp(x,i-1,j);
716: uE = UInterp(x,i,j); pC = PInterp(x,i,j);
717: }
719: return 2.0*etaC*(uE-uW)/dx - pC;
720: }
722: /*---------------------------------------------------------------------*/
723: /* computes the normal stress---used on the boundaries */
724: static inline PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
725: /*---------------------------------------------------------------------*/
726: {
727: Parameter *param=user->param;
728: GridInfo *grid =user->grid;
729: PetscScalar dz =grid->dz;
730: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
731: PetscScalar epsC =0.0, etaC, TC;
732: PetscScalar pC, wN, wS, z_scale;
733: if (i<j || j<=grid->jlid) return EPS_ZERO;
735: ivisc=param->ivisc; z_scale = param->z_scale;
737: if (ipos==CELL_CENTER) { /* on cell center */
739: TC = param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*z_scale);
740: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
741: etaC = Viscosity(TC,epsC,dz*j,param);
742: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
744: } else { /* on cell corner */
745: if ((i==ilim) || (j==jlim)) return EPS_ZERO;
747: TC = param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*z_scale);
748: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
749: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
750: if (i==j) wN = param->sb;
751: else wN = WInterp(x,i,j);
752: wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
753: }
755: return 2.0*etaC*(wN-wS)/dz - pC;
756: }
758: /*---------------------------------------------------------------------*/
760: /*=====================================================================
761: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
762: =====================================================================*/
764: /*---------------------------------------------------------------------*/
765: /* initializes the problem parameters and checks for
766: command line changes */
767: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
768: /*---------------------------------------------------------------------*/
769: {
770: PetscReal SEC_PER_YR = 3600.00*24.00*365.2500;
771: PetscReal alpha_g_on_cp_units_inverse_km = 4.0e-5*9.8;
773: /* domain geometry */
774: param->slab_dip = 45.0;
775: param->width = 320.0; /* km */
776: param->depth = 300.0; /* km */
777: param->lid_depth = 35.0; /* km */
778: param->fault_depth = 35.0; /* km */
780: PetscOptionsGetReal(NULL,NULL,"-slab_dip",&(param->slab_dip),NULL);
781: PetscOptionsGetReal(NULL,NULL,"-width",&(param->width),NULL);
782: PetscOptionsGetReal(NULL,NULL,"-depth",&(param->depth),NULL);
783: PetscOptionsGetReal(NULL,NULL,"-lid_depth",&(param->lid_depth),NULL);
784: PetscOptionsGetReal(NULL,NULL,"-fault_depth",&(param->fault_depth),NULL);
786: param->slab_dip = param->slab_dip*PETSC_PI/180.0; /* radians */
788: /* grid information */
789: PetscOptionsGetInt(NULL,NULL, "-jfault",&(grid->jfault),NULL);
790: grid->ni = 82;
791: PetscOptionsGetInt(NULL,NULL, "-ni",&(grid->ni),NULL);
793: grid->dx = param->width/((PetscReal)(grid->ni-2)); /* km */
794: grid->dz = grid->dx*PetscTanReal(param->slab_dip); /* km */
795: grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/
796: param->depth = grid->dz*(grid->nj-2); /* km */
797: grid->inose = 0; /* gridpoints*/
798: PetscOptionsGetInt(NULL,NULL,"-inose",&(grid->inose),NULL);
799: grid->bx = DM_BOUNDARY_NONE;
800: grid->by = DM_BOUNDARY_NONE;
801: grid->stencil = DMDA_STENCIL_BOX;
802: grid->dof = 4;
803: grid->stencil_width = 2;
804: grid->mglevels = 1;
806: /* boundary conditions */
807: param->pv_analytic = PETSC_FALSE;
808: param->ibound = BC_NOSTRESS;
809: PetscOptionsGetInt(NULL,NULL,"-ibound",&(param->ibound),NULL);
811: /* physical constants */
812: param->slab_velocity = 5.0; /* cm/yr */
813: param->slab_age = 50.0; /* Ma */
814: param->lid_age = 50.0; /* Ma */
815: param->kappa = 0.7272e-6; /* m^2/sec */
816: param->potentialT = 1300.0; /* degrees C */
818: PetscOptionsGetReal(NULL,NULL,"-slab_velocity",&(param->slab_velocity),NULL);
819: PetscOptionsGetReal(NULL,NULL,"-slab_age",&(param->slab_age),NULL);
820: PetscOptionsGetReal(NULL,NULL,"-lid_age",&(param->lid_age),NULL);
821: PetscOptionsGetReal(NULL,NULL,"-kappa",&(param->kappa),NULL);
822: PetscOptionsGetReal(NULL,NULL,"-potentialT",&(param->potentialT),NULL);
824: /* viscosity */
825: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
826: param->eta0 = 1e24; /* Pa-s */
827: param->visc_cutoff = 0.0; /* factor of eta_0 */
828: param->continuation = 1.0;
830: /* constants for diffusion creep */
831: param->diffusion.A = 1.8e7; /* Pa-s */
832: param->diffusion.n = 1.0; /* dim'less */
833: param->diffusion.Estar = 375e3; /* J/mol */
834: param->diffusion.Vstar = 5e-6; /* m^3/mol */
836: /* constants for param->dislocationocation creep */
837: param->dislocation.A = 2.8969e4; /* Pa-s */
838: param->dislocation.n = 3.5; /* dim'less */
839: param->dislocation.Estar = 530e3; /* J/mol */
840: param->dislocation.Vstar = 14e-6; /* m^3/mol */
842: PetscOptionsGetInt(NULL,NULL, "-ivisc",&(param->ivisc),NULL);
843: PetscOptionsGetReal(NULL,NULL,"-visc_cutoff",&(param->visc_cutoff),NULL);
845: param->output_ivisc = param->ivisc;
847: PetscOptionsGetInt(NULL,NULL,"-output_ivisc",&(param->output_ivisc),NULL);
848: PetscOptionsGetReal(NULL,NULL,"-vstar",&(param->dislocation.Vstar),NULL);
850: /* output options */
851: param->quiet = PETSC_FALSE;
852: param->param_test = PETSC_FALSE;
854: PetscOptionsHasName(NULL,NULL,"-quiet",&(param->quiet));
855: PetscOptionsHasName(NULL,NULL,"-test",&(param->param_test));
856: PetscOptionsGetString(NULL,NULL,"-file",param->filename,sizeof(param->filename),&(param->output_to_file));
858: /* advection */
859: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
861: PetscOptionsGetInt(NULL,NULL,"-adv_scheme",&(param->adv_scheme),NULL);
863: /* misc. flags */
864: param->stop_solve = PETSC_FALSE;
865: param->interrupted = PETSC_FALSE;
866: param->kspmon = PETSC_FALSE;
867: param->toggle_kspmon = PETSC_FALSE;
869: /* derived parameters for slab angle */
870: param->sb = PetscSinReal(param->slab_dip);
871: param->cb = PetscCosReal(param->slab_dip);
872: param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
873: param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);
875: /* length, velocity and time scale for non-dimensionalization */
876: param->L = PetscMin(param->width,param->depth); /* km */
877: param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */
879: /* other unit conversions and derived parameters */
880: param->scaled_width = param->width/param->L; /* dim'less */
881: param->scaled_depth = param->depth/param->L; /* dim'less */
882: param->lid_depth = param->lid_depth/param->L; /* dim'less */
883: param->fault_depth = param->fault_depth/param->L; /* dim'less */
884: grid->dx = grid->dx/param->L; /* dim'less */
885: grid->dz = grid->dz/param->L; /* dim'less */
886: grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */
887: grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */
888: param->lid_depth = grid->jlid*grid->dz; /* dim'less */
889: param->fault_depth = grid->jfault*grid->dz; /* dim'less */
890: grid->corner = grid->jlid+1; /* gridcells */
891: param->peclet = param->V /* m/sec */
892: * param->L*1000.0 /* m */
893: / param->kappa; /* m^2/sec */
894: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
895: param->skt = PetscSqrtReal(param->kappa*param->slab_age*SEC_PER_YR);
896: PetscOptionsGetReal(NULL,NULL,"-peclet",&(param->peclet),NULL);
898: return 0;
899: }
901: /*---------------------------------------------------------------------*/
902: /* prints a report of the problem parameters to stdout */
903: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
904: /*---------------------------------------------------------------------*/
905: {
906: char date[30];
908: PetscGetDate(date,30);
910: if (!(param->quiet)) {
911: PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
912: PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
913: PetscPrintf(PETSC_COMM_WORLD," Width = %g km, Depth = %g km\n",(double)param->width,(double)param->depth);
914: 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);
915: PetscPrintf(PETSC_COMM_WORLD," Lid depth = %5.2f km, Fault depth = %5.2f km\n",(double)(param->lid_depth*param->L),(double)(param->fault_depth*param->L));
917: PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
918: 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));
919: PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);
920: PetscPrintf(PETSC_COMM_WORLD," Pe = %g\n",(double)param->peclet);
922: PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
923: if (param->ivisc==VISC_CONST) {
924: PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");
925: if (param->pv_analytic) {
926: PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");
927: }
928: } else if (param->ivisc==VISC_DIFN) {
929: PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");
930: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
931: } else if (param->ivisc==VISC_DISL) {
932: PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");
933: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
934: } else if (param->ivisc==VISC_FULL) {
935: PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");
936: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %g--%g Pa-sec \n",(double)param->eta0,(double)(param->visc_cutoff*param->eta0));
937: } else {
938: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
939: return 1;
940: }
942: PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
943: if (param->ibound==BC_ANALYTIC) {
944: PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");
945: } else if (param->ibound==BC_NOSTRESS) {
946: PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");
947: } else if (param->ibound==BC_EXPERMNT) {
948: PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");
949: } else {
950: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
951: return 1;
952: }
954: if (param->output_to_file) {
955: #if defined(PETSC_HAVE_MATLAB_ENGINE)
956: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);
957: #else
958: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);
959: #endif
960: }
961: if (param->output_ivisc != param->ivisc) {
962: PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);
963: }
965: PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
966: }
967: if (param->param_test) PetscEnd();
968: return 0;
969: }
971: /* ------------------------------------------------------------------- */
972: /* generates an initial guess using the analytic solution for isoviscous
973: corner flow */
974: PetscErrorCode Initialize(DM da)
975: /* ------------------------------------------------------------------- */
976: {
977: AppCtx *user;
978: Parameter *param;
979: GridInfo *grid;
980: PetscInt i,j,is,js,im,jm;
981: Field **x;
982: Vec Xguess;
984: /* Get the fine grid */
985: DMGetApplicationContext(da,&user);
986: Xguess = user->Xguess;
987: param = user->param;
988: grid = user->grid;
989: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
990: DMDAVecGetArray(da,Xguess,(void**)&x);
992: /* Compute initial guess */
993: for (j=js; j<js+jm; j++) {
994: for (i=is; i<is+im; i++) {
995: if (i<j) x[j][i].u = param->cb;
996: else if (j<=grid->jlid) x[j][i].u = 0.0;
997: else x[j][i].u = HorizVelocity(i,j,user);
999: if (i<=j) x[j][i].w = param->sb;
1000: else if (j<=grid->jlid) x[j][i].w = 0.0;
1001: else x[j][i].w = VertVelocity(i,j,user);
1003: if (i<j || j<=grid->jlid) x[j][i].p = 0.0;
1004: else x[j][i].p = Pressure(i,j,user);
1006: x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1007: }
1008: }
1010: /* Restore x to Xguess */
1011: DMDAVecRestoreArray(da,Xguess,(void**)&x);
1013: return 0;
1014: }
1016: /*---------------------------------------------------------------------*/
1017: /* controls output to a file */
1018: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1019: /*---------------------------------------------------------------------*/
1020: {
1021: AppCtx *user;
1022: Parameter *param;
1023: GridInfo *grid;
1024: PetscInt ivt;
1025: PetscMPIInt rank;
1026: PetscViewer viewer;
1027: Vec res, pars;
1028: MPI_Comm comm;
1029: DM da;
1031: SNESGetDM(snes,&da);
1032: DMGetApplicationContext(da,&user);
1033: param = user->param;
1034: grid = user->grid;
1035: ivt = param->ivisc;
1037: param->ivisc = param->output_ivisc;
1039: /* compute final residual and final viscosity/strain rate fields */
1040: SNESGetFunction(snes, &res, NULL, NULL);
1041: ViscosityField(da, user->x, user->Xguess);
1043: /* get the communicator and the rank of the processor */
1044: PetscObjectGetComm((PetscObject)snes, &comm);
1045: MPI_Comm_rank(comm, &rank);
1047: if (param->output_to_file) { /* send output to binary file */
1048: VecCreate(comm, &pars);
1049: if (rank == 0) { /* on processor 0 */
1050: VecSetSizes(pars, 20, PETSC_DETERMINE);
1051: VecSetFromOptions(pars);
1052: VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1053: VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1054: VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1055: VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1056: VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1057: VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1058: /* skipped 6 intentionally */
1059: VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1060: VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1061: VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1062: VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1063: VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1064: VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1065: VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1066: VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1067: VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1068: VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1069: } else { /* on some other processor */
1070: VecSetSizes(pars, 0, PETSC_DETERMINE);
1071: VecSetFromOptions(pars);
1072: }
1073: VecAssemblyBegin(pars)); PetscCall(VecAssemblyEnd(pars);
1075: /* create viewer */
1076: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1077: PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1078: #else
1079: PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1080: #endif
1082: /* send vectors to viewer */
1083: PetscObjectSetName((PetscObject)res,"res");
1084: VecView(res,viewer);
1085: PetscObjectSetName((PetscObject)user->x,"out");
1086: VecView(user->x, viewer);
1087: PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1088: VecView(user->Xguess, viewer);
1089: StressField(da); /* compute stress fields */
1090: PetscObjectSetName((PetscObject)(user->Xguess),"str");
1091: VecView(user->Xguess, viewer);
1092: PetscObjectSetName((PetscObject)pars,"par");
1093: VecView(pars, viewer);
1095: /* destroy viewer and vector */
1096: PetscViewerDestroy(&viewer);
1097: VecDestroy(&pars);
1098: }
1100: param->ivisc = ivt;
1101: return 0;
1102: }
1104: /* ------------------------------------------------------------------- */
1105: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1106: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1107: /* ------------------------------------------------------------------- */
1108: {
1109: AppCtx *user;
1110: Parameter *param;
1111: GridInfo *grid;
1112: Vec localX;
1113: Field **v, **x;
1114: PetscReal eps, /* dx,*/ dz, T, epsC, TC;
1115: PetscInt i,j,is,js,im,jm,ilim,jlim,ivt;
1118: DMGetApplicationContext(da,&user);
1119: param = user->param;
1120: grid = user->grid;
1121: ivt = param->ivisc;
1122: param->ivisc = param->output_ivisc;
1124: DMGetLocalVector(da, &localX);
1125: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1126: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1127: DMDAVecGetArray(da,localX,(void**)&x);
1128: DMDAVecGetArray(da,V,(void**)&v);
1130: /* Parameters */
1131: /* dx = grid->dx; */ dz = grid->dz;
1133: ilim = grid->ni-1; jlim = grid->nj-1;
1135: /* Compute real temperature, strain rate and viscosity */
1136: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1137: for (j=js; j<js+jm; j++) {
1138: for (i=is; i<is+im; i++) {
1139: T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j-0.5)*dz*param->z_scale));
1140: if (i<ilim && j<jlim) {
1141: TC = PetscRealPart(param->potentialT * TInterp(x,i,j) * PetscExpScalar(j*dz*param->z_scale));
1142: } else {
1143: TC = T;
1144: }
1145: eps = PetscRealPart((CalcSecInv(x,i,j,CELL_CENTER,user)));
1146: epsC = PetscRealPart(CalcSecInv(x,i,j,CELL_CORNER,user));
1148: v[j][i].u = eps;
1149: v[j][i].w = epsC;
1150: v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1151: v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1152: }
1153: }
1154: DMDAVecRestoreArray(da,V,(void**)&v);
1155: DMDAVecRestoreArray(da,localX,(void**)&x);
1156: DMRestoreLocalVector(da, &localX);
1158: param->ivisc = ivt;
1159: return 0;
1160: }
1162: /* ------------------------------------------------------------------- */
1163: /* post-processing: compute stress everywhere */
1164: PetscErrorCode StressField(DM da)
1165: /* ------------------------------------------------------------------- */
1166: {
1167: AppCtx *user;
1168: PetscInt i,j,is,js,im,jm;
1169: Vec locVec;
1170: Field **x, **y;
1172: DMGetApplicationContext(da,&user);
1174: /* Get the fine grid of Xguess and X */
1175: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1176: DMDAVecGetArray(da,user->Xguess,(void**)&x);
1178: DMGetLocalVector(da, &locVec);
1179: DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);
1180: DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);
1181: DMDAVecGetArray(da,locVec,(void**)&y);
1183: /* Compute stress on the corner points */
1184: for (j=js; j<js+jm; j++) {
1185: for (i=is; i<is+im; i++) {
1186: x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1187: x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1188: x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1189: x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1190: }
1191: }
1193: /* Restore the fine grid of Xguess and X */
1194: DMDAVecRestoreArray(da,user->Xguess,(void**)&x);
1195: DMDAVecRestoreArray(da,locVec,(void**)&y);
1196: DMRestoreLocalVector(da, &locVec);
1197: return 0;
1198: }
1200: /*=====================================================================
1201: UTILITY FUNCTIONS
1202: =====================================================================*/
1204: /*---------------------------------------------------------------------*/
1205: /* returns the velocity of the subducting slab and handles fault nodes
1206: for BC */
1207: static inline PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1208: /*---------------------------------------------------------------------*/
1209: {
1210: Parameter *param = user->param;
1211: GridInfo *grid = user->grid;
1213: if (c=='U' || c=='u') {
1214: if (i<j-1) return param->cb;
1215: else if (j<=grid->jfault) return 0.0;
1216: else return param->cb;
1218: } else {
1219: if (i<j) return param->sb;
1220: else if (j<=grid->jfault) return 0.0;
1221: else return param->sb;
1222: }
1223: }
1225: /*---------------------------------------------------------------------*/
1226: /* solution to diffusive half-space cooling model for BC */
1227: static inline PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1228: /*---------------------------------------------------------------------*/
1229: {
1230: Parameter *param = user->param;
1231: PetscScalar z;
1232: if (plate==PLATE_LID) z = (j-0.5)*user->grid->dz;
1233: else z = (j-0.5)*user->grid->dz*param->cb; /* PLATE_SLAB */
1234: #if defined(PETSC_HAVE_ERF)
1235: return (PetscReal)(erf((double)PetscRealPart(z*param->L/2.0/param->skt)));
1236: #else
1237: (*PetscErrorPrintf)("erf() not available on this machine\n");
1238: MPI_Abort(PETSC_COMM_SELF,1);
1239: #endif
1240: }
1242: /*=====================================================================
1243: INTERACTIVE SIGNAL HANDLING
1244: =====================================================================*/
1246: /* ------------------------------------------------------------------- */
1247: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1248: /* ------------------------------------------------------------------- */
1249: {
1250: AppCtx *user = (AppCtx*) ctx;
1251: Parameter *param = user->param;
1252: KSP ksp;
1255: if (param->interrupted) {
1256: param->interrupted = PETSC_FALSE;
1257: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1258: *reason = SNES_CONVERGED_FNORM_ABS;
1259: return 0;
1260: } else if (param->toggle_kspmon) {
1261: param->toggle_kspmon = PETSC_FALSE;
1263: SNESGetKSP(snes, &ksp);
1265: if (param->kspmon) {
1266: KSPMonitorCancel(ksp);
1268: param->kspmon = PETSC_FALSE;
1269: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1270: } else {
1271: PetscViewerAndFormat *vf;
1272: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);
1273: KSPMonitorSet(ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))KSPMonitorSingularValue,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);
1275: param->kspmon = PETSC_TRUE;
1276: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1277: }
1278: }
1279: PetscFunctionReturn(SNESConvergedDefault(snes,it,xnorm,snorm,fnorm,reason,ctx));
1280: }
1282: /* ------------------------------------------------------------------- */
1283: #include <signal.h>
1284: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1285: /* ------------------------------------------------------------------- */
1286: {
1287: AppCtx *user = (AppCtx*) ctx;
1288: Parameter *param = user->param;
1290: if (signum == SIGILL) {
1291: param->toggle_kspmon = PETSC_TRUE;
1292: #if !defined(PETSC_MISSING_SIGCONT)
1293: } else if (signum == SIGCONT) {
1294: param->interrupted = PETSC_TRUE;
1295: #endif
1296: #if !defined(PETSC_MISSING_SIGURG)
1297: } else if (signum == SIGURG) {
1298: param->stop_solve = PETSC_TRUE;
1299: #endif
1300: }
1301: return 0;
1302: }
1304: /*---------------------------------------------------------------------*/
1305: /* main call-back function that computes the processor-local piece
1306: of the residual */
1307: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
1308: /*---------------------------------------------------------------------*/
1309: {
1310: AppCtx *user = (AppCtx*)ptr;
1311: Parameter *param = user->param;
1312: GridInfo *grid = user->grid;
1313: PetscScalar mag_w, mag_u;
1314: PetscInt i,j,mx,mz,ilim,jlim;
1315: PetscInt is,ie,js,je,ibound; /* ,ivisc */
1318: /* Define global and local grid parameters */
1319: mx = info->mx; mz = info->my;
1320: ilim = mx-1; jlim = mz-1;
1321: is = info->xs; ie = info->xs+info->xm;
1322: js = info->ys; je = info->ys+info->ym;
1324: /* Define geometric and numeric parameters */
1325: /* ivisc = param->ivisc; */ ibound = param->ibound;
1327: for (j=js; j<je; j++) {
1328: for (i=is; i<ie; i++) {
1330: /************* X-MOMENTUM/VELOCITY *************/
1331: if (i<j) f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
1332: else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1333: /* in the lithospheric lid */
1334: f[j][i].u = x[j][i].u - 0.0;
1335: } else if (i==ilim) {
1336: /* on the right side boundary */
1337: if (ibound==BC_ANALYTIC) {
1338: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1339: } else {
1340: f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1341: }
1343: } else if (j==jlim) {
1344: /* on the bottom boundary */
1345: if (ibound==BC_ANALYTIC) {
1346: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1347: } else if (ibound==BC_NOSTRESS) {
1348: f[j][i].u = XMomentumResidual(x,i,j,user);
1349: } else {
1350: /* experimental boundary condition */
1351: }
1353: } else {
1354: /* in the mantle wedge */
1355: f[j][i].u = XMomentumResidual(x,i,j,user);
1356: }
1358: /************* Z-MOMENTUM/VELOCITY *************/
1359: if (i<=j) {
1360: f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);
1362: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1363: /* in the lithospheric lid */
1364: f[j][i].w = x[j][i].w - 0.0;
1366: } else if (j==jlim) {
1367: /* on the bottom boundary */
1368: if (ibound==BC_ANALYTIC) {
1369: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1370: } else {
1371: f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1372: }
1374: } else if (i==ilim) {
1375: /* on the right side boundary */
1376: if (ibound==BC_ANALYTIC) {
1377: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1378: } else if (ibound==BC_NOSTRESS) {
1379: f[j][i].w = ZMomentumResidual(x,i,j,user);
1380: } else {
1381: /* experimental boundary condition */
1382: }
1384: } else {
1385: /* in the mantle wedge */
1386: f[j][i].w = ZMomentumResidual(x,i,j,user);
1387: }
1389: /************* CONTINUITY/PRESSURE *************/
1390: if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1391: /* in the lid or slab */
1392: f[j][i].p = x[j][i].p;
1394: } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
1395: /* on an analytic boundary */
1396: f[j][i].p = x[j][i].p - Pressure(i,j,user);
1398: } else {
1399: /* in the mantle wedge */
1400: f[j][i].p = ContinuityResidual(x,i,j,user);
1401: }
1403: /************* TEMPERATURE *************/
1404: if (j==0) {
1405: /* on the surface */
1406: f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0);
1408: } else if (i==0) {
1409: /* slab inflow boundary */
1410: f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user);
1412: } else if (i==ilim) {
1413: /* right side boundary */
1414: mag_u = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j][i-1].u)/param->cb,1.0),0.0)), 5);
1415: 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);
1417: } else if (j==jlim) {
1418: /* bottom boundary */
1419: mag_w = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j-1][i].w)/param->sb,1.0),0.0)), 5);
1420: f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);
1422: } else {
1423: /* in the mantle wedge */
1424: f[j][i].T = EnergyResidual(x,i,j,user);
1425: }
1426: }
1427: }
1428: return 0;
1429: }
1431: /*TEST
1433: build:
1434: requires: !complex erf
1436: test:
1437: args: -ni 18
1438: filter: grep -v Destination
1439: requires: !single
1441: TEST*/