Actual source code: ex30.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
58: #define VISC_CONST 0
59: #define VISC_DIFN 1
60: #define VISC_DISL 2
61: #define VISC_FULL 3
62: #define CELL_CENTER 0
63: #define CELL_CORNER 1
64: #define BC_ANALYTIC 0
65: #define BC_NOSTRESS 1
66: #define BC_EXPERMNT 2
67: #define ADVECT_FV 0
68: #define ADVECT_FROMM 1
69: #define PLATE_SLAB 0
70: #define PLATE_LID 1
71: #define EPS_ZERO 0.00000001
73: typedef struct { /* holds the variables to be solved for */
74: PetscScalar u,w,p,T;
75: } Field;
77: typedef struct { /* parameters needed to compute viscosity */
78: PetscReal A,n,Estar,Vstar;
79: } ViscParam;
81: typedef struct { /* physical and miscelaneous parameters */
82: PetscReal width, depth, scaled_width, scaled_depth, peclet, potentialT;
83: PetscReal slab_dip, slab_age, slab_velocity, kappa, z_scale;
84: PetscReal c, d, sb, cb, skt, visc_cutoff, lid_age, eta0, continuation;
85: PetscReal L, V, lid_depth, fault_depth;
86: ViscParam diffusion, dislocation;
87: PetscInt ivisc, adv_scheme, ibound, output_ivisc;
88: PetscBool quiet, param_test, output_to_file, pv_analytic;
89: PetscBool interrupted, stop_solve, toggle_kspmon, kspmon;
90: char filename[PETSC_MAX_PATH_LEN];
91: } Parameter;
93: typedef struct { /* grid parameters */
94: DMDABoundaryType bx,by;
95: DMDAStencilType stencil;
96: PetscInt corner,ni,nj,jlid,jfault,inose;
97: PetscInt dof,stencil_width,mglevels;
98: PetscReal dx,dz;
99: } GridInfo;
101: typedef struct { /* application context */
102: Vec x,Xguess;
103: Parameter *param;
104: GridInfo *grid;
105: } AppCtx;
107: /* Callback functions (static interface) */
108: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
110: /* Main routines */
111: extern PetscErrorCode SetParams(Parameter*, GridInfo*);
112: extern PetscErrorCode ReportParams(Parameter*, GridInfo*);
113: extern PetscErrorCode Initialize(DM);
114: extern PetscErrorCode UpdateSolution(SNES,AppCtx*, PetscInt*);
115: extern PetscErrorCode DoOutput(SNES,PetscInt);
117: /* Post-processing & misc */
118: extern PetscErrorCode ViscosityField(DM,Vec,Vec);
119: extern PetscErrorCode StressField(DM);
120: extern PetscErrorCode SNESConverged_Interactive(SNES, PetscInt, PetscReal, PetscReal, PetscReal, SNESConvergedReason*, void*);
121: extern PetscErrorCode InteractiveHandler(int, void*);
122: extern PetscBool OptionsHasName(const char pre[],const char name[]);
124: /*-----------------------------------------------------------------------*/
127: int main(int argc,char **argv)
128: /*-----------------------------------------------------------------------*/
129: {
130: SNES snes;
131: AppCtx *user; /* user-defined work context */
132: Parameter param;
133: GridInfo grid;
134: PetscInt nits;
136: MPI_Comm comm;
137: DM da;
139: PetscInitialize(&argc,&argv,(char*)0,help);
140: PetscOptionsSetValue("-file","ex30_output");
141: PetscOptionsSetValue("-snes_monitor_short",NULL);
142: PetscOptionsSetValue("-snes_max_it","20");
143: PetscOptionsSetValue("-ksp_max_it","1500");
144: PetscOptionsSetValue("-ksp_gmres_restart","300");
145: PetscOptionsInsert(&argc,&argv,NULL);
147: comm = PETSC_COMM_WORLD;
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Set up the problem parameters.
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: SetParams(¶m,&grid);
153: ReportParams(¶m,&grid);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: SNESCreate(comm,&snes);
158: DMDACreate2d(comm,grid.bx,grid.by,grid.stencil,grid.ni,grid.nj,PETSC_DECIDE,PETSC_DECIDE,grid.dof,grid.stencil_width,0,0,&da);
159: SNESSetDM(snes,da);
160: DMDASetFieldName(da,0,"x-velocity");
161: DMDASetFieldName(da,1,"y-velocity");
162: DMDASetFieldName(da,2,"pressure");
163: DMDASetFieldName(da,3,"temperature");
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Create user context, set problem data, create vector data structures.
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: PetscMalloc(sizeof(AppCtx),&user);
170: user->param = ¶m;
171: user->grid = &grid;
172: DMSetApplicationContext(da,user);
173: DMCreateGlobalVector(da,&(user->Xguess));
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Set up the SNES solver with callback functions.
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,(void*)user);
180: SNESSetFromOptions(snes);
183: SNESSetConvergenceTest(snes,SNESConverged_Interactive,(void*)user,NULL);
184: PetscPushSignalHandler(InteractiveHandler,(void*)user);
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Initialize and solve the nonlinear system
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: Initialize(da);
190: UpdateSolution(snes,user,&nits);
192: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: Output variables.
194: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195: DoOutput(snes,nits);
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Free work space.
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200: VecDestroy(&user->Xguess);
201: VecDestroy(&user->x);
202: PetscFree(user);
203: SNESDestroy(&snes);
204: DMDestroy(&da);
205: PetscPopSignalHandler();
206: PetscFinalize();
207: return 0;
208: }
210: /*=====================================================================
211: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
212: =====================================================================*/
214: /*---------------------------------------------------------------------*/
217: /* manages solve: adaptive continuation method */
218: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
219: {
220: KSP ksp;
221: PC pc;
222: SNESConvergedReason reason;
223: Parameter *param = user->param;
224: PetscReal cont_incr=0.3;
225: PetscInt its;
226: PetscErrorCode ierr;
227: PetscBool q = PETSC_FALSE;
228: DM dm;
231: SNESGetDM(snes,&dm);
232: DMCreateGlobalVector(dm,&user->x);
233: SNESGetKSP(snes,&ksp);
234: KSPGetPC(ksp, &pc);
235: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
237: *nits=0;
239: /* Isoviscous solve */
240: if (param->ivisc == VISC_CONST && !param->stop_solve) {
241: param->ivisc = VISC_CONST;
243: SNESSolve(snes,0,user->x);
244: SNESGetConvergedReason(snes,&reason);
245: SNESGetIterationNumber(snes,&its);
246: *nits += its;
247: VecCopy(user->x,user->Xguess);
248: if (param->stop_solve) goto done;
249: }
251: /* Olivine diffusion creep */
252: if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
253: if (!q) PetscPrintf(PETSC_COMM_WORLD,"Computing Variable Viscosity Solution\n");
255: /* continuation method on viscosity cutoff */
256: for (param->continuation=0.0;; param->continuation+=cont_incr) {
257: if (!q) PetscPrintf(PETSC_COMM_WORLD," Continuation parameter = %G\n", param->continuation);
259: /* solve the non-linear system */
260: VecCopy(user->Xguess,user->x);
261: SNESSolve(snes,0,user->x);
262: SNESGetConvergedReason(snes,&reason);
263: SNESGetIterationNumber(snes,&its);
264: *nits += its;
265: if (!q) PetscPrintf(PETSC_COMM_WORLD," SNES iterations: %D, Cumulative: %D\n", its, *nits);
266: if (param->stop_solve) goto done;
268: if (reason<0) {
269: /* NOT converged */
270: cont_incr = -fabs(cont_incr)/2.0;
271: if (fabs(cont_incr)<0.01) goto done;
273: } else {
274: /* converged */
275: VecCopy(user->x,user->Xguess);
276: if (param->continuation >= 1.0) goto done;
277: if (its<=3) cont_incr = 0.30001;
278: else if (its<=8) cont_incr = 0.15001;
279: else cont_incr = 0.10001;
281: if (param->continuation+cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
282: } /* endif reason<0 */
283: }
284: }
285: done:
286: if (param->stop_solve && !q) PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: stopping solve.\n");
287: if (reason<0 && !q) PetscPrintf(PETSC_COMM_WORLD,"FAILED TO CONVERGE: stopping solve.\n");
288: return(0);
289: }
292: /*=====================================================================
293: PHYSICS FUNCTIONS (compute the discrete residual)
294: =====================================================================*/
297: /*---------------------------------------------------------------------*/
300: PETSC_STATIC_INLINE PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
301: /*---------------------------------------------------------------------*/
302: {
303: return 0.25*(x[j][i].u+x[j+1][i].u+x[j][i+1].u+x[j+1][i+1].u);
304: }
306: /*---------------------------------------------------------------------*/
309: PETSC_STATIC_INLINE PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
310: /*---------------------------------------------------------------------*/
311: {
312: return 0.25*(x[j][i].w+x[j+1][i].w+x[j][i+1].w+x[j+1][i+1].w);
313: }
315: /*---------------------------------------------------------------------*/
318: PETSC_STATIC_INLINE PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
319: /*---------------------------------------------------------------------*/
320: {
321: return 0.25*(x[j][i].p+x[j+1][i].p+x[j][i+1].p+x[j+1][i+1].p);
322: }
324: /*---------------------------------------------------------------------*/
327: PETSC_STATIC_INLINE PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
328: /*---------------------------------------------------------------------*/
329: {
330: return 0.25*(x[j][i].T+x[j+1][i].T+x[j][i+1].T+x[j+1][i+1].T);
331: }
333: /*---------------------------------------------------------------------*/
336: /* isoviscous analytic solution for IC */
337: PETSC_STATIC_INLINE PassiveScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
338: /*---------------------------------------------------------------------*/
339: {
340: Parameter *param = user->param;
341: GridInfo *grid = user->grid;
342: PetscScalar st, ct, th, c=param->c, d=param->d;
343: PetscReal x, z,r;
345: x = (i - grid->jlid)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
346: r = PetscSqrtReal(x*x+z*z);
347: st = z/r;
348: ct = x/r;
349: th = atan(z/x);
350: return ct*(c*th*st+d*(st+th*ct)) + st*(c*(st-th*ct)+d*th*st);
351: }
353: /*---------------------------------------------------------------------*/
356: /* isoviscous analytic solution for IC */
357: PETSC_STATIC_INLINE PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
358: /*---------------------------------------------------------------------*/
359: {
360: Parameter *param = user->param;
361: GridInfo *grid = user->grid;
362: PetscScalar st, ct, th, c=param->c, d=param->d;
363: PetscReal x, z, r;
365: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid)*grid->dz;
366: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r; th = atan(z/x);
367: return st*(c*th*st+d*(st+th*ct)) - ct*(c*(st-th*ct)+d*th*st);
368: }
370: /*---------------------------------------------------------------------*/
373: /* isoviscous analytic solution for IC */
374: PETSC_STATIC_INLINE PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
375: /*---------------------------------------------------------------------*/
376: {
377: Parameter *param = user->param;
378: GridInfo *grid = user->grid;
379: PetscScalar x, z, r, st, ct, c=param->c, d=param->d;
381: x = (i - grid->jlid - 0.5)*grid->dx; z = (j - grid->jlid - 0.5)*grid->dz;
382: r = PetscSqrtReal(x*x+z*z); st = z/r; ct = x/r;
383: return (-2.0*(c*ct-d*st)/r);
384: }
388: /* computes the second invariant of the strain rate tensor */
389: PETSC_STATIC_INLINE PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
390: /*---------------------------------------------------------------------*/
391: {
392: Parameter *param = user->param;
393: GridInfo *grid = user->grid;
394: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
395: PetscScalar uN,uS,uE,uW,wN,wS,wE,wW;
396: PetscScalar eps11, eps12, eps22;
398: if (i<j) return EPS_ZERO;
399: if (i==ilim) i--;
400: if (j==jlim) j--;
402: if (ipos==CELL_CENTER) { /* on cell center */
403: if (j<=grid->jlid) return EPS_ZERO;
405: uE = x[j][i].u; uW = x[j][i-1].u;
406: wN = x[j][i].w; wS = x[j-1][i].w;
407: wE = WInterp(x,i,j-1);
408: if (i==j) {
409: uN = param->cb; wW = param->sb;
410: } else {
411: uN = UInterp(x,i-1,j); wW = WInterp(x,i-1,j-1);
412: }
414: if (j==grid->jlid+1) uS = 0.0;
415: else uS = UInterp(x,i-1,j-1);
417: } else { /* on CELL_CORNER */
418: if (j<grid->jlid) return EPS_ZERO;
420: uN = x[j+1][i].u; uS = x[j][i].u;
421: wE = x[j][i+1].w; wW = x[j][i].w;
422: if (i==j) {
423: wN = param->sb;
424: uW = param->cb;
425: } else {
426: wN = WInterp(x,i,j);
427: uW = UInterp(x,i-1,j);
428: }
430: if (j==grid->jlid) {
431: uE = 0.0; uW = 0.0;
432: uS = -uN;
433: wS = -wN;
434: } else {
435: uE = UInterp(x,i,j);
436: wS = WInterp(x,i,j-1);
437: }
438: }
440: eps11 = (uE-uW)/grid->dx; eps22 = (wN-wS)/grid->dz;
441: eps12 = 0.5*((uN-uS)/grid->dz + (wE-wW)/grid->dx);
443: return PetscSqrtReal(0.5*(eps11*eps11 + 2.0*eps12*eps12 + eps22*eps22));
444: }
446: /*---------------------------------------------------------------------*/
449: /* computes the shear viscosity */
450: PETSC_STATIC_INLINE PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PassiveScalar z, Parameter *param)
451: /*---------------------------------------------------------------------*/
452: {
453: PetscReal result =0.0;
454: ViscParam difn =param->diffusion, disl=param->dislocation;
455: PetscInt iVisc =param->ivisc;
456: PetscScalar eps_scale=param->V/(param->L*1000.0);
457: PetscScalar strain_power, v1, v2, P;
458: PetscScalar rho_g = 32340.0, R=8.3144;
460: P = rho_g*(z*param->L*1000.0); /* Pa */
462: if (iVisc==VISC_CONST) {
463: /* constant viscosity */
464: return 1.0;
465: } else if (iVisc==VISC_DIFN) {
466: /* diffusion creep rheology */
467: result = PetscRealPart((difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0));
468: } else if (iVisc==VISC_DISL) {
469: /* dislocation creep rheology */
470: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
472: result = PetscRealPart(disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0);
473: } else if (iVisc==VISC_FULL) {
474: /* dislocation/diffusion creep rheology */
475: strain_power = PetscPowScalar(eps*eps_scale, (1.0-disl.n)/disl.n);
477: v1 = difn.A*PetscExpScalar((difn.Estar + P*difn.Vstar)/R/(T+273.0))/param->eta0;
478: v2 = disl.A*PetscExpScalar((disl.Estar + P*disl.Vstar)/disl.n/R/(T+273.0))*strain_power/param->eta0;
480: result = PetscRealPart(1.0/(1.0/v1 + 1.0/v2));
481: }
483: /* max viscosity is param->eta0 */
484: result = PetscMin(result, 1.0);
485: /* min viscosity is param->visc_cutoff */
486: result = PetscMax(result, param->visc_cutoff);
487: /* continuation method */
488: result = PetscPowReal(result,param->continuation);
489: return result;
490: }
492: /*---------------------------------------------------------------------*/
495: /* computes the residual of the x-component of eqn (1) above */
496: PETSC_STATIC_INLINE PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
497: /*---------------------------------------------------------------------*/
498: {
499: Parameter *param=user->param;
500: GridInfo *grid =user->grid;
501: PetscScalar dx = grid->dx, dz=grid->dz;
502: PetscScalar etaN,etaS,etaE,etaW,epsN=0.0,epsS=0.0,epsE=0.0,epsW=0.0;
503: PetscScalar TE=0.0,TN=0.0,TS=0.0,TW=0.0, dPdx, residual, z_scale;
504: PetscScalar dudxW,dudxE,dudzN,dudzS,dwdxN,dwdxS;
505: PetscInt jlim = grid->nj-1;
507: z_scale = param->z_scale;
509: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
510: TS = param->potentialT * TInterp(x,i,j-1) * exp((j-1.0)*dz*z_scale);
511: if (j==jlim) TN = TS;
512: else TN = param->potentialT * TInterp(x,i,j) * exp(j*dz*z_scale);
513: TW = param->potentialT * x[j][i].T * exp((j-0.5)*dz*z_scale);
514: TE = param->potentialT * x[j][i+1].T * exp((j-0.5)*dz*z_scale);
515: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
516: epsN = CalcSecInv(x,i,j, CELL_CORNER,user);
517: epsS = CalcSecInv(x,i,j-1,CELL_CORNER,user);
518: epsE = CalcSecInv(x,i+1,j,CELL_CENTER,user);
519: epsW = CalcSecInv(x,i,j, CELL_CENTER,user);
520: }
521: }
522: etaN = Viscosity(TN,epsN,dz*(j+0.5),param);
523: etaS = Viscosity(TS,epsS,dz*(j-0.5),param);
524: etaW = Viscosity(TW,epsW,dz*j,param);
525: etaE = Viscosity(TE,epsE,dz*j,param);
527: dPdx = (x[j][i+1].p - x[j][i].p)/dx;
528: if (j==jlim) dudzN = etaN * (x[j][i].w - x[j][i+1].w)/dx;
529: else dudzN = etaN * (x[j+1][i].u - x[j][i].u) /dz;
530: dudzS = etaS * (x[j][i].u - x[j-1][i].u)/dz;
531: dudxE = etaE * (x[j][i+1].u - x[j][i].u) /dx;
532: dudxW = etaW * (x[j][i].u - x[j][i-1].u)/dx;
534: residual = -dPdx /* X-MOMENTUM EQUATION*/
535: +(dudxE - dudxW)/dx
536: +(dudzN - dudzS)/dz;
538: if (param->ivisc!=VISC_CONST) {
539: dwdxN = etaN * (x[j][i+1].w - x[j][i].w) /dx;
540: dwdxS = etaS * (x[j-1][i+1].w - x[j-1][i].w)/dx;
542: residual += (dudxE - dudxW)/dx + (dwdxN - dwdxS)/dz;
543: }
545: return residual;
546: }
548: /*---------------------------------------------------------------------*/
551: /* computes the residual of the z-component of eqn (1) above */
552: PETSC_STATIC_INLINE PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
553: /*---------------------------------------------------------------------*/
554: {
555: Parameter *param=user->param;
556: GridInfo *grid =user->grid;
557: PetscScalar dx = grid->dx, dz=grid->dz;
558: 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;
559: PetscScalar TE =0.0,TN=0.0,TS=0.0,TW=0.0, dPdz, residual,z_scale;
560: PetscScalar dudzE,dudzW,dwdxW,dwdxE,dwdzN,dwdzS;
561: PetscInt ilim = grid->ni-1;
563: /* geometric and other parameters */
564: z_scale = param->z_scale;
566: /* viscosity */
567: if (param->ivisc==VISC_DIFN || param->ivisc>=VISC_DISL) { /* viscosity is T-dependent */
568: TN = param->potentialT * x[j+1][i].T * exp((j+0.5)*dz*z_scale);
569: TS = param->potentialT * x[j][i].T * exp((j-0.5)*dz*z_scale);
570: TW = param->potentialT * TInterp(x,i-1,j) * exp(j*dz*z_scale);
571: if (i==ilim) TE = TW;
572: else TE = param->potentialT * TInterp(x,i,j) * exp(j*dz*z_scale);
573: if (param->ivisc>=VISC_DISL) { /* olivine dislocation creep */
574: epsN = CalcSecInv(x,i,j+1,CELL_CENTER,user);
575: epsS = CalcSecInv(x,i,j, CELL_CENTER,user);
576: epsE = CalcSecInv(x,i,j, CELL_CORNER,user);
577: epsW = CalcSecInv(x,i-1,j,CELL_CORNER,user);
578: }
579: }
580: etaN = Viscosity(TN,epsN,dz*(j+1.0),param);
581: etaS = Viscosity(TS,epsS,dz*(j+0.0),param);
582: etaW = Viscosity(TW,epsW,dz*(j+0.5),param);
583: etaE = Viscosity(TE,epsE,dz*(j+0.5),param);
585: dPdz = (x[j+1][i].p - x[j][i].p)/dz;
586: dwdzN = etaN * (x[j+1][i].w - x[j][i].w)/dz;
587: dwdzS = etaS * (x[j][i].w - x[j-1][i].w)/dz;
588: if (i==ilim) dwdxE = etaE * (x[j][i].u - x[j+1][i].u)/dz;
589: else dwdxE = etaE * (x[j][i+1].w - x[j][i].w) /dx;
590: dwdxW = 2.0*etaW * (x[j][i].w - x[j][i-1].w)/dx;
592: /* Z-MOMENTUM */
593: residual = -dPdz /* constant viscosity terms */
594: +(dwdzN - dwdzS)/dz
595: +(dwdxE - dwdxW)/dx;
597: if (param->ivisc!=VISC_CONST) {
598: dudzE = etaE * (x[j+1][i].u - x[j][i].u)/dz;
599: dudzW = etaW * (x[j+1][i-1].u - x[j][i-1].u)/dz;
601: residual += (dwdzN - dwdzS)/dz + (dudzE - dudzW)/dx;
602: }
604: return residual;
605: }
607: /*---------------------------------------------------------------------*/
610: /* computes the residual of eqn (2) above */
611: PETSC_STATIC_INLINE PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
612: /*---------------------------------------------------------------------*/
613: {
614: GridInfo *grid =user->grid;
615: PetscScalar uE,uW,wN,wS,dudx,dwdz;
617: uW = x[j][i-1].u; uE = x[j][i].u; dudx = (uE - uW)/grid->dx;
618: wS = x[j-1][i].w; wN = x[j][i].w; dwdz = (wN - wS)/grid->dz;
620: return dudx + dwdz;
621: }
623: /*---------------------------------------------------------------------*/
626: /* computes the residual of eqn (3) above */
627: PETSC_STATIC_INLINE PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
628: /*---------------------------------------------------------------------*/
629: {
630: Parameter *param=user->param;
631: GridInfo *grid =user->grid;
632: PetscScalar dx = grid->dx, dz=grid->dz;
633: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, jlid=grid->jlid;
634: PetscScalar TE, TN, TS, TW, residual;
635: PetscScalar uE,uW,wN,wS;
636: PetscScalar fN,fS,fE,fW,dTdxW,dTdxE,dTdzN,dTdzS;
638: dTdzN = (x[j+1][i].T - x[j][i].T) /dz;
639: dTdzS = (x[j][i].T - x[j-1][i].T)/dz;
640: dTdxE = (x[j][i+1].T - x[j][i].T) /dx;
641: dTdxW = (x[j][i].T - x[j][i-1].T)/dx;
643: residual = ((dTdzN - dTdzS)/dz + /* diffusion term */
644: (dTdxE - dTdxW)/dx)*dx*dz/param->peclet;
646: if (j<=jlid && i>=j) {
647: /* don't advect in the lid */
648: return residual;
649: } else if (i<j) {
650: /* beneath the slab sfc */
651: uW = uE = param->cb;
652: wS = wN = param->sb;
653: } else {
654: /* advect in the slab and wedge */
655: uW = x[j][i-1].u; uE = x[j][i].u;
656: wS = x[j-1][i].w; wN = x[j][i].w;
657: }
659: if (param->adv_scheme==ADVECT_FV || i==ilim-1 || j==jlim-1 || i==1 || j==1) {
660: /* finite volume advection */
661: TS = (x[j][i].T + x[j-1][i].T)/2.0;
662: TN = (x[j][i].T + x[j+1][i].T)/2.0;
663: TE = (x[j][i].T + x[j][i+1].T)/2.0;
664: TW = (x[j][i].T + x[j][i-1].T)/2.0;
665: fN = wN*TN*dx; fS = wS*TS*dx;
666: fE = uE*TE*dz; fW = uW*TW*dz;
668: } else {
669: /* Fromm advection scheme */
670: 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
671: - 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;
672: 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
673: - 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;
674: 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
675: - 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;
676: 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
677: - 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;
678: }
680: residual -= (fE - fW + fN - fS);
682: return residual;
683: }
685: /*---------------------------------------------------------------------*/
688: /* computes the shear stress---used on the boundaries */
689: PETSC_STATIC_INLINE PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
690: /*---------------------------------------------------------------------*/
691: {
692: Parameter *param=user->param;
693: GridInfo *grid =user->grid;
694: PetscInt ilim =grid->ni-1, jlim=grid->nj-1;
695: PetscScalar uN, uS, wE, wW;
697: if (j<=grid->jlid || i<j || i==ilim || j==jlim) return EPS_ZERO;
699: if (ipos==CELL_CENTER) { /* on cell center */
701: wE = WInterp(x,i,j-1);
702: if (i==j) {
703: wW = param->sb;
704: uN = param->cb;
705: } else {
706: wW = WInterp(x,i-1,j-1);
707: uN = UInterp(x,i-1,j);
708: }
709: if (j==grid->jlid+1) uS = 0.0;
710: else uS = UInterp(x,i-1,j-1);
712: } else { /* on cell corner */
714: uN = x[j+1][i].u; uS = x[j][i].u;
715: wW = x[j][i].w; wE = x[j][i+1].w;
717: }
719: return (uN-uS)/grid->dz + (wE-wW)/grid->dx;
720: }
722: /*---------------------------------------------------------------------*/
725: /* computes the normal stress---used on the boundaries */
726: PETSC_STATIC_INLINE PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
727: /*---------------------------------------------------------------------*/
728: {
729: Parameter *param=user->param;
730: GridInfo *grid =user->grid;
731: PetscScalar dx = grid->dx, dz=grid->dz;
732: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
733: PetscScalar epsC =0.0, etaC, TC, uE, uW, pC, z_scale;
734: if (i<j || j<=grid->jlid) return EPS_ZERO;
736: ivisc=param->ivisc; z_scale = param->z_scale;
738: if (ipos==CELL_CENTER) { /* on cell center */
740: TC = param->potentialT * x[j][i].T * exp((j-0.5)*dz*z_scale);
741: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
742: etaC = Viscosity(TC,epsC,dz*j,param);
744: uW = x[j][i-1].u; uE = x[j][i].u;
745: pC = x[j][i].p;
747: } else { /* on cell corner */
748: if (i==ilim || j==jlim) return EPS_ZERO;
750: TC = param->potentialT * TInterp(x,i,j) * exp(j*dz*z_scale);
751: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
752: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
754: if (i==j) uW = param->sb;
755: else uW = UInterp(x,i-1,j);
756: uE = UInterp(x,i,j); pC = PInterp(x,i,j);
757: }
759: return 2.0*etaC*(uE-uW)/dx - pC;
760: }
762: /*---------------------------------------------------------------------*/
765: /* computes the normal stress---used on the boundaries */
766: PETSC_STATIC_INLINE PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
767: /*---------------------------------------------------------------------*/
768: {
769: Parameter *param=user->param;
770: GridInfo *grid =user->grid;
771: PetscScalar dz =grid->dz;
772: PetscInt ilim =grid->ni-1, jlim=grid->nj-1, ivisc;
773: PetscScalar epsC =0.0, etaC, TC;
774: PetscScalar pC, wN, wS, z_scale;
775: if (i<j || j<=grid->jlid) return EPS_ZERO;
777: ivisc=param->ivisc; z_scale = param->z_scale;
779: if (ipos==CELL_CENTER) { /* on cell center */
781: TC = param->potentialT * x[j][i].T * exp((j-0.5)*dz*z_scale);
782: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CENTER,user);
783: etaC = Viscosity(TC,epsC,dz*j,param);
784: wN = x[j][i].w; wS = x[j-1][i].w; pC = x[j][i].p;
786: } else { /* on cell corner */
787: if ((i==ilim) || (j==jlim)) return EPS_ZERO;
789: TC = param->potentialT * TInterp(x,i,j) * exp(j*dz*z_scale);
790: if (ivisc>=VISC_DISL) epsC = CalcSecInv(x,i,j,CELL_CORNER,user);
791: etaC = Viscosity(TC,epsC,dz*(j+0.5),param);
792: if (i==j) wN = param->sb;
793: else wN = WInterp(x,i,j);
794: wS = WInterp(x,i,j-1); pC = PInterp(x,i,j);
795: }
797: return 2.0*etaC*(wN-wS)/dz - pC;
798: }
800: /*---------------------------------------------------------------------*/
802: /*=====================================================================
803: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
804: =====================================================================*/
806: /*---------------------------------------------------------------------*/
809: /* initializes the problem parameters and checks for
810: command line changes */
811: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
812: /*---------------------------------------------------------------------*/
813: {
814: PetscErrorCode ierr, ierr_out=0;
815: PetscReal SEC_PER_YR = 3600.00*24.00*365.2500;
816: PetscReal PI = 3.14159265358979323846;
817: PetscReal alpha_g_on_cp_units_inverse_km=4.0e-5*9.8;
819: /* domain geometry */
820: param->slab_dip = 45.0;
821: param->width = 320.0; /* km */
822: param->depth = 300.0; /* km */
823: param->lid_depth = 35.0; /* km */
824: param->fault_depth = 35.0; /* km */
826: PetscOptionsGetReal(NULL,"-slab_dip",&(param->slab_dip),NULL);
827: PetscOptionsGetReal(NULL,"-width",&(param->width),NULL);
828: PetscOptionsGetReal(NULL,"-depth",&(param->depth),NULL);
829: PetscOptionsGetReal(NULL,"-lid_depth",&(param->lid_depth),NULL);
830: PetscOptionsGetReal(NULL,"-fault_depth",&(param->fault_depth),NULL);
832: param->slab_dip = param->slab_dip*PI/180.0; /* radians */
834: /* grid information */
835: PetscOptionsGetInt(NULL, "-jfault",&(grid->jfault),NULL);
836: grid->ni = 82;
837: PetscOptionsGetInt(NULL, "-ni",&(grid->ni),NULL);
839: grid->dx = param->width/((double)(grid->ni-2)); /* km */
840: grid->dz = grid->dx*tan(param->slab_dip); /* km */
841: grid->nj = (PetscInt)(param->depth/grid->dz + 3.0); /* gridpoints*/
842: param->depth = grid->dz*(grid->nj-2); /* km */
843: grid->inose = 0; /* gridpoints*/
844: PetscOptionsGetInt(NULL,"-inose",&(grid->inose),NULL);
845: grid->bx = DMDA_BOUNDARY_NONE;
846: grid->by = DMDA_BOUNDARY_NONE;
847: grid->stencil = DMDA_STENCIL_BOX;
848: grid->dof = 4;
849: grid->stencil_width = 2;
850: grid->mglevels = 1;
852: /* boundary conditions */
853: param->pv_analytic = PETSC_FALSE;
854: param->ibound = BC_NOSTRESS;
855: PetscOptionsGetInt(NULL,"-ibound",&(param->ibound),NULL);
857: /* physical constants */
858: param->slab_velocity = 5.0; /* cm/yr */
859: param->slab_age = 50.0; /* Ma */
860: param->lid_age = 50.0; /* Ma */
861: param->kappa = 0.7272e-6; /* m^2/sec */
862: param->potentialT = 1300.0; /* degrees C */
864: PetscOptionsGetReal(NULL,"-slab_velocity",&(param->slab_velocity),NULL);
865: PetscOptionsGetReal(NULL,"-slab_age",&(param->slab_age),NULL);
866: PetscOptionsGetReal(NULL,"-lid_age",&(param->lid_age),NULL);
867: PetscOptionsGetReal(NULL,"-kappa",&(param->kappa),NULL);
868: PetscOptionsGetReal(NULL,"-potentialT",&(param->potentialT),NULL);
870: /* viscosity */
871: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
872: param->eta0 = 1e24; /* Pa-s */
873: param->visc_cutoff = 0.0; /* factor of eta_0 */
874: param->continuation = 1.0;
876: /* constants for diffusion creep */
877: param->diffusion.A = 1.8e7; /* Pa-s */
878: param->diffusion.n = 1.0; /* dim'less */
879: param->diffusion.Estar = 375e3; /* J/mol */
880: param->diffusion.Vstar = 5e-6; /* m^3/mol */
882: /* constants for param->dislocationocation creep */
883: param->dislocation.A = 2.8969e4; /* Pa-s */
884: param->dislocation.n = 3.5; /* dim'less */
885: param->dislocation.Estar = 530e3; /* J/mol */
886: param->dislocation.Vstar = 14e-6; /* m^3/mol */
888: PetscOptionsGetInt(NULL, "-ivisc",&(param->ivisc),NULL);
889: PetscOptionsGetReal(NULL,"-visc_cutoff",&(param->visc_cutoff),NULL);
891: param->output_ivisc = param->ivisc;
893: PetscOptionsGetInt(NULL,"-output_ivisc",&(param->output_ivisc),NULL);
894: PetscOptionsGetReal(NULL,"-vstar",&(param->dislocation.Vstar),NULL);
896: /* output options */
897: param->quiet = PETSC_FALSE;
898: param->param_test = PETSC_FALSE;
900: PetscOptionsHasName(NULL,"-quiet",&(param->quiet));
901: PetscOptionsHasName(NULL,"-test",&(param->param_test));
902: PetscOptionsGetString(NULL,"-file",param->filename,PETSC_MAX_PATH_LEN,&(param->output_to_file));
904: /* advection */
905: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
907: PetscOptionsGetInt(NULL,"-adv_scheme",&(param->adv_scheme),NULL);
909: /* misc. flags */
910: param->stop_solve = PETSC_FALSE;
911: param->interrupted = PETSC_FALSE;
912: param->kspmon = PETSC_FALSE;
913: param->toggle_kspmon = PETSC_FALSE;
915: /* derived parameters for slab angle */
916: param->sb = sin(param->slab_dip);
917: param->cb = cos(param->slab_dip);
918: param->c = param->slab_dip*param->sb/(param->slab_dip*param->slab_dip-param->sb*param->sb);
919: param->d = (param->slab_dip*param->cb-param->sb)/(param->slab_dip*param->slab_dip-param->sb*param->sb);
921: /* length, velocity and time scale for non-dimensionalization */
922: param->L = PetscMin(param->width,param->depth); /* km */
923: param->V = param->slab_velocity/100.0/SEC_PER_YR; /* m/sec */
925: /* other unit conversions and derived parameters */
926: param->scaled_width = param->width/param->L; /* dim'less */
927: param->scaled_depth = param->depth/param->L; /* dim'less */
928: param->lid_depth = param->lid_depth/param->L; /* dim'less */
929: param->fault_depth = param->fault_depth/param->L; /* dim'less */
930: grid->dx = grid->dx/param->L; /* dim'less */
931: grid->dz = grid->dz/param->L; /* dim'less */
932: grid->jlid = (PetscInt)(param->lid_depth/grid->dz); /* gridcells */
933: grid->jfault = (PetscInt)(param->fault_depth/grid->dz); /* gridcells */
934: param->lid_depth = grid->jlid*grid->dz; /* dim'less */
935: param->fault_depth = grid->jfault*grid->dz; /* dim'less */
936: grid->corner = grid->jlid+1; /* gridcells */
937: param->peclet = param->V /* m/sec */
938: * param->L*1000.0 /* m */
939: / param->kappa; /* m^2/sec */
940: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
941: param->skt = PetscSqrtReal(param->kappa*param->slab_age*SEC_PER_YR);
942: PetscOptionsGetReal(NULL,"-peclet",&(param->peclet),NULL);
944: return ierr_out;
945: }
947: /*---------------------------------------------------------------------*/
950: /* prints a report of the problem parameters to stdout */
951: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
952: /*---------------------------------------------------------------------*/
953: {
954: PetscErrorCode ierr, ierr_out=0;
955: char date[30];
956: PetscReal PI = 3.14159265358979323846;
958: PetscGetDate(date,30);
960: if (!(param->quiet)) {
961: PetscPrintf(PETSC_COMM_WORLD,"---------------------BEGIN ex30 PARAM REPORT-------------------\n");
962: /* PetscPrintf(PETSC_COMM_WORLD," %s\n",&(date[0]));*/
964: PetscPrintf(PETSC_COMM_WORLD,"Domain: \n");
965: PetscPrintf(PETSC_COMM_WORLD," Width = %G km, Depth = %G km\n",param->width,param->depth);
966: PetscPrintf(PETSC_COMM_WORLD," Slab dip = %G degrees, Slab velocity = %G cm/yr\n",param->slab_dip*180.0/PI,param->slab_velocity);
967: 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);
969: PetscPrintf(PETSC_COMM_WORLD,"\nGrid: \n");
970: PetscPrintf(PETSC_COMM_WORLD," [ni,nj] = %D, %D [dx,dz] = %G, %G km\n",grid->ni,grid->nj,grid->dx*param->L,grid->dz*param->L);
971: PetscPrintf(PETSC_COMM_WORLD," jlid = %3D jfault = %3D \n",grid->jlid,grid->jfault);
972: PetscPrintf(PETSC_COMM_WORLD," Pe = %G\n",param->peclet);
974: PetscPrintf(PETSC_COMM_WORLD,"\nRheology:");
975: if (param->ivisc==VISC_CONST) {
976: PetscPrintf(PETSC_COMM_WORLD," Isoviscous \n");
977: if (param->pv_analytic) {
978: PetscPrintf(PETSC_COMM_WORLD," Pressure and Velocity prescribed! \n");
979: }
980: } else if (param->ivisc==VISC_DIFN) {
981: PetscPrintf(PETSC_COMM_WORLD," Diffusion Creep (T-Dependent Newtonian) \n");
982: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
983: } else if (param->ivisc==VISC_DISL) {
984: PetscPrintf(PETSC_COMM_WORLD," Dislocation Creep (T-Dependent Non-Newtonian) \n");
985: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
986: } else if (param->ivisc==VISC_FULL) {
987: PetscPrintf(PETSC_COMM_WORLD," Full Rheology \n");
988: PetscPrintf(PETSC_COMM_WORLD," Viscosity range: %G--%G Pa-sec \n",param->eta0,param->visc_cutoff*param->eta0);
989: } else {
990: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
991: ierr_out = 1;
992: }
994: PetscPrintf(PETSC_COMM_WORLD,"Boundary condition:");
995: if (param->ibound==BC_ANALYTIC) {
996: PetscPrintf(PETSC_COMM_WORLD," Isoviscous Analytic Dirichlet \n");
997: } else if (param->ibound==BC_NOSTRESS) {
998: PetscPrintf(PETSC_COMM_WORLD," Stress-Free (normal & shear stress)\n");
999: } else if (param->ibound==BC_EXPERMNT) {
1000: PetscPrintf(PETSC_COMM_WORLD," Experimental boundary condition \n");
1001: } else {
1002: PetscPrintf(PETSC_COMM_WORLD," Invalid! \n");
1003: ierr_out = 1;
1004: }
1006: if (param->output_to_file)
1007: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1008: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: Mat file \"%s\"\n",param->filename);
1009: #else
1010: PetscPrintf(PETSC_COMM_WORLD,"Output Destination: PETSc binary file \"%s\"\n",param->filename);
1011: #endif
1012: if (param->output_ivisc != param->ivisc) {
1013: PetscPrintf(PETSC_COMM_WORLD," Output viscosity: -ivisc %D\n",param->output_ivisc);
1014: }
1016: PetscPrintf(PETSC_COMM_WORLD,"---------------------END ex30 PARAM REPORT---------------------\n");
1017: }
1018: if (param->param_test) PetscEnd();
1019: return ierr_out;
1020: }
1022: /* ------------------------------------------------------------------- */
1025: /* generates an inital guess using the analytic solution for isoviscous
1026: corner flow */
1027: PetscErrorCode Initialize(DM da)
1028: /* ------------------------------------------------------------------- */
1029: {
1030: AppCtx *user;
1031: Parameter *param;
1032: GridInfo *grid;
1033: PetscInt i,j,is,js,im,jm;
1035: Field **x;
1036: Vec Xguess;
1038: /* Get the fine grid */
1039: DMGetApplicationContext(da,&user);
1040: Xguess = user->Xguess;
1041: param = user->param;
1042: grid = user->grid;
1043: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1044: DMDAVecGetArray(da,Xguess,(void**)&x);
1046: /* Compute initial guess */
1047: for (j=js; j<js+jm; j++) {
1048: for (i=is; i<is+im; i++) {
1049: if (i<j) x[j][i].u = param->cb;
1050: else if (j<=grid->jlid) x[j][i].u = 0.0;
1051: else x[j][i].u = HorizVelocity(i,j,user);
1053: if (i<=j) x[j][i].w = param->sb;
1054: else if (j<=grid->jlid) x[j][i].w = 0.0;
1055: else x[j][i].w = VertVelocity(i,j,user);
1057: if (i<j || j<=grid->jlid) x[j][i].p = 0.0;
1058: else x[j][i].p = Pressure(i,j,user);
1060: x[j][i].T = PetscMin(grid->dz*(j-0.5),1.0);
1061: }
1062: }
1064: /* Restore x to Xguess */
1065: DMDAVecRestoreArray(da,Xguess,(void**)&x);
1067: return 0;
1068: }
1070: /*---------------------------------------------------------------------*/
1073: /* controls output to a file */
1074: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1075: /*---------------------------------------------------------------------*/
1076: {
1077: AppCtx *user;
1078: Parameter *param;
1079: GridInfo *grid;
1080: PetscInt ivt;
1082: PetscMPIInt rank;
1083: PetscViewer viewer;
1084: Vec res, pars;
1085: MPI_Comm comm;
1086: DM da;
1088: SNESGetDM(snes,&da);
1089: DMGetApplicationContext(da,&user);
1090: param = user->param;
1091: grid = user->grid;
1092: ivt = param->ivisc;
1094: param->ivisc = param->output_ivisc;
1096: /* compute final residual and final viscosity/strain rate fields */
1097: SNESGetFunction(snes, &res, NULL, NULL);
1098: ViscosityField(da, user->x, user->Xguess);
1100: /* get the communicator and the rank of the processor */
1101: PetscObjectGetComm((PetscObject)snes, &comm);
1102: MPI_Comm_rank(comm, &rank);
1104: if (param->output_to_file) { /* send output to binary file */
1105: VecCreate(comm, &pars);
1106: if (!rank) { /* on processor 0 */
1107: VecSetSizes(pars, 20, PETSC_DETERMINE);
1108: VecSetFromOptions(pars);
1109: VecSetValue(pars,0, (PetscScalar)(grid->ni),INSERT_VALUES);
1110: VecSetValue(pars,1, (PetscScalar)(grid->nj),INSERT_VALUES);
1111: VecSetValue(pars,2, (PetscScalar)(grid->dx),INSERT_VALUES);
1112: VecSetValue(pars,3, (PetscScalar)(grid->dz),INSERT_VALUES);
1113: VecSetValue(pars,4, (PetscScalar)(param->L),INSERT_VALUES);
1114: VecSetValue(pars,5, (PetscScalar)(param->V),INSERT_VALUES);
1115: /* skipped 6 intentionally */
1116: VecSetValue(pars,7, (PetscScalar)(param->slab_dip),INSERT_VALUES);
1117: VecSetValue(pars,8, (PetscScalar)(grid->jlid),INSERT_VALUES);
1118: VecSetValue(pars,9, (PetscScalar)(param->lid_depth),INSERT_VALUES);
1119: VecSetValue(pars,10,(PetscScalar)(grid->jfault),INSERT_VALUES);
1120: VecSetValue(pars,11,(PetscScalar)(param->fault_depth),INSERT_VALUES);
1121: VecSetValue(pars,12,(PetscScalar)(param->potentialT),INSERT_VALUES);
1122: VecSetValue(pars,13,(PetscScalar)(param->ivisc),INSERT_VALUES);
1123: VecSetValue(pars,14,(PetscScalar)(param->visc_cutoff),INSERT_VALUES);
1124: VecSetValue(pars,15,(PetscScalar)(param->ibound),INSERT_VALUES);
1125: VecSetValue(pars,16,(PetscScalar)(its),INSERT_VALUES);
1126: } else { /* on some other processor */
1127: VecSetSizes(pars, 0, PETSC_DETERMINE);
1128: VecSetFromOptions(pars);
1129: }
1130: VecAssemblyBegin(pars); VecAssemblyEnd(pars);
1132: /* create viewer */
1133: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1134: PetscViewerMatlabOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1135: #else
1136: PetscViewerBinaryOpen(PETSC_COMM_WORLD,param->filename,FILE_MODE_WRITE,&viewer);
1137: #endif
1139: /* send vectors to viewer */
1140: PetscObjectSetName((PetscObject)res,"res");
1141: VecView(res,viewer);
1142: PetscObjectSetName((PetscObject)user->x,"out");
1143: VecView(user->x, viewer);
1144: PetscObjectSetName((PetscObject)(user->Xguess),"aux");
1145: VecView(user->Xguess, viewer);
1146: StressField(da); /* compute stress fields */
1147: PetscObjectSetName((PetscObject)(user->Xguess),"str");
1148: VecView(user->Xguess, viewer);
1149: PetscObjectSetName((PetscObject)pars,"par");
1150: VecView(pars, viewer);
1152: /* destroy viewer and vector */
1153: PetscViewerDestroy(&viewer);
1154: VecDestroy(&pars);
1155: }
1157: param->ivisc = ivt;
1158: return 0;
1159: }
1161: /* ------------------------------------------------------------------- */
1164: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1165: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1166: /* ------------------------------------------------------------------- */
1167: {
1168: AppCtx *user;
1169: Parameter *param;
1170: GridInfo *grid;
1171: Vec localX;
1172: Field **v, **x;
1173: PassiveReal eps, /* dx,*/ dz, T, epsC, TC;
1174: PetscInt i,j,is,js,im,jm,ilim,jlim,ivt;
1178: DMGetApplicationContext(da,&user);
1179: param = user->param;
1180: grid = user->grid;
1181: ivt = param->ivisc;
1182: param->ivisc = param->output_ivisc;
1184: DMGetLocalVector(da, &localX);
1185: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
1186: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
1187: DMDAVecGetArray(da,localX,(void**)&x);
1188: DMDAVecGetArray(da,V,(void**)&v);
1190: /* Parameters */
1191: /* dx = grid->dx; */ dz = grid->dz;
1193: ilim = grid->ni-1; jlim = grid->nj-1;
1195: /* Compute real temperature, strain rate and viscosity */
1196: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1197: for (j=js; j<js+jm; j++) {
1198: for (i=is; i<is+im; i++) {
1199: T = PetscRealPart(param->potentialT * x[j][i].T * exp((j-0.5)*dz*param->z_scale));
1200: if (i<ilim && j<jlim) {
1201: TC = PetscRealPart(param->potentialT * TInterp(x,i,j) * exp(j*dz*param->z_scale));
1202: } else {
1203: TC = T;
1204: }
1205: eps = PetscRealPart((CalcSecInv(x,i,j,CELL_CENTER,user)));
1206: epsC = PetscRealPart(CalcSecInv(x,i,j,CELL_CORNER,user));
1208: v[j][i].u = eps;
1209: v[j][i].w = epsC;
1210: v[j][i].p = Viscosity(T,eps,dz*(j-0.5),param);
1211: v[j][i].T = Viscosity(TC,epsC,dz*j,param);
1212: }
1213: }
1214: DMDAVecRestoreArray(da,V,(void**)&v);
1215: DMDAVecRestoreArray(da,localX,(void**)&x);
1216: DMRestoreLocalVector(da, &localX);
1218: param->ivisc = ivt;
1219: return(0);
1220: }
1222: /* ------------------------------------------------------------------- */
1225: /* post-processing: compute stress everywhere */
1226: PetscErrorCode StressField(DM da)
1227: /* ------------------------------------------------------------------- */
1228: {
1229: AppCtx *user;
1230: PetscInt i,j,is,js,im,jm;
1232: Vec locVec;
1233: Field **x, **y;
1235: DMGetApplicationContext(da,&user);
1237: /* Get the fine grid of Xguess and X */
1238: DMDAGetCorners(da,&is,&js,NULL,&im,&jm,NULL);
1239: DMDAVecGetArray(da,user->Xguess,(void**)&x);
1241: DMGetLocalVector(da, &locVec);
1242: DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec);
1243: DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec);
1244: DMDAVecGetArray(da,locVec,(void**)&y);
1246: /* Compute stress on the corner points */
1247: for (j=js; j<js+jm; j++) {
1248: for (i=is; i<is+im; i++) {
1249: x[j][i].u = ShearStress(y,i,j,CELL_CENTER,user);
1250: x[j][i].w = ShearStress(y,i,j,CELL_CORNER,user);
1251: x[j][i].p = XNormalStress(y,i,j,CELL_CENTER,user);
1252: x[j][i].T = ZNormalStress(y,i,j,CELL_CENTER,user);
1253: }
1254: }
1256: /* Restore the fine grid of Xguess and X */
1257: DMDAVecRestoreArray(da,user->Xguess,(void**)&x);
1258: DMDAVecRestoreArray(da,locVec,(void**)&y);
1259: DMRestoreLocalVector(da, &locVec);
1260: return 0;
1261: }
1263: /*=====================================================================
1264: UTILITY FUNCTIONS
1265: =====================================================================*/
1267: /*---------------------------------------------------------------------*/
1270: /* returns the velocity of the subducting slab and handles fault nodes
1271: for BC */
1272: PETSC_STATIC_INLINE PassiveScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1273: /*---------------------------------------------------------------------*/
1274: {
1275: Parameter *param = user->param;
1276: GridInfo *grid = user->grid;
1278: if (c=='U' || c=='u') {
1279: if (i<j-1) return param->cb;
1280: else if (j<=grid->jfault) return 0.0;
1281: else return param->cb;
1283: } else {
1284: if (i<j) return param->sb;
1285: else if (j<=grid->jfault) return 0.0;
1286: else return param->sb;
1287: }
1288: }
1290: /*---------------------------------------------------------------------*/
1293: /* solution to diffusive half-space cooling model for BC */
1294: PETSC_STATIC_INLINE PassiveScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1295: /*---------------------------------------------------------------------*/
1296: {
1297: Parameter *param = user->param;
1298: PassiveScalar z;
1299: if (plate==PLATE_LID) z = (j-0.5)*user->grid->dz;
1300: else z = (j-0.5)*user->grid->dz*param->cb; /* PLATE_SLAB */
1301: #if defined(PETSC_HAVE_ERF)
1302: return (erf(PetscRealPart(z*param->L/2.0/param->skt)));
1303: #else
1304: SETERRQ(PETSC_COMM_SELF,1,"erf() not available on this machine");
1305: #endif
1306: }
1309: /* ------------------------------------------------------------------- */
1312: /* utility function */
1313: PetscBool OptionsHasName(const char pre[],const char name[])
1314: /* ------------------------------------------------------------------- */
1315: {
1316: PetscBool retval;
1318: PetscOptionsHasName(pre,name,&retval);CHKERRABORT(PETSC_COMM_WORLD,ierr);
1319: return retval;
1320: }
1322: /*=====================================================================
1323: INTERACTIVE SIGNAL HANDLING
1324: =====================================================================*/
1326: /* ------------------------------------------------------------------- */
1329: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it,PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1330: /* ------------------------------------------------------------------- */
1331: {
1332: AppCtx *user = (AppCtx*) ctx;
1333: Parameter *param = user->param;
1334: KSP ksp;
1338: if (param->interrupted) {
1339: param->interrupted = PETSC_FALSE;
1340: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: exiting SNES solve. \n");
1341: *reason = SNES_CONVERGED_FNORM_ABS;
1342: return(0);
1343: } else if (param->toggle_kspmon) {
1344: param->toggle_kspmon = PETSC_FALSE;
1346: SNESGetKSP(snes, &ksp);
1348: if (param->kspmon) {
1349: KSPMonitorCancel(ksp);
1351: param->kspmon = PETSC_FALSE;
1352: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: deactivating ksp singular value monitor. \n");
1353: } else {
1354: KSPMonitorSet(ksp,KSPMonitorSingularValue,NULL,NULL);
1356: param->kspmon = PETSC_TRUE;
1357: PetscPrintf(PETSC_COMM_WORLD,"USER SIGNAL: activating ksp singular value monitor. \n");
1358: }
1359: }
1360: PetscFunctionReturn(SNESConvergedDefault(snes,it,xnorm,snorm,fnorm,reason,ctx));
1361: }
1363: /* ------------------------------------------------------------------- */
1364: #include <signal.h>
1367: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1368: /* ------------------------------------------------------------------- */
1369: {
1370: AppCtx *user = (AppCtx*) ctx;
1371: Parameter *param = user->param;
1373: if (signum == SIGILL) {
1374: param->toggle_kspmon = PETSC_TRUE;
1375: #if !defined(PETSC_MISSING_SIGCONT)
1376: } else if (signum == SIGCONT) {
1377: param->interrupted = PETSC_TRUE;
1378: #endif
1379: #if !defined(PETSC_MISSING_SIGURG)
1380: } else if (signum == SIGURG) {
1381: param->stop_solve = PETSC_TRUE;
1382: #endif
1383: }
1384: return 0;
1385: }
1387: /*---------------------------------------------------------------------*/
1390: /* main call-back function that computes the processor-local piece
1391: of the residual */
1392: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,Field **x,Field **f,void *ptr)
1393: /*---------------------------------------------------------------------*/
1394: {
1395: AppCtx *user = (AppCtx*)ptr;
1396: Parameter *param = user->param;
1397: GridInfo *grid = user->grid;
1398: PetscScalar mag_w, mag_u;
1399: PetscInt i,j,mx,mz,ilim,jlim;
1400: PetscInt is,ie,js,je,ibound; /* ,ivisc */
1403: /* Define global and local grid parameters */
1404: mx = info->mx; mz = info->my;
1405: ilim = mx-1; jlim = mz-1;
1406: is = info->xs; ie = info->xs+info->xm;
1407: js = info->ys; je = info->ys+info->ym;
1409: /* Define geometric and numeric parameters */
1410: /* ivisc = param->ivisc; */ ibound = param->ibound;
1412: for (j=js; j<je; j++) {
1413: for (i=is; i<ie; i++) {
1415: /************* X-MOMENTUM/VELOCITY *************/
1416: if (i<j) f[j][i].u = x[j][i].u - SlabVel('U',i,j,user);
1417: else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1418: /* in the lithospheric lid */
1419: f[j][i].u = x[j][i].u - 0.0;
1420: } else if (i==ilim) {
1421: /* on the right side boundary */
1422: if (ibound==BC_ANALYTIC) {
1423: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1424: } else {
1425: f[j][i].u = XNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1426: }
1428: } else if (j==jlim) {
1429: /* on the bottom boundary */
1430: if (ibound==BC_ANALYTIC) {
1431: f[j][i].u = x[j][i].u - HorizVelocity(i,j,user);
1432: } else if (ibound==BC_NOSTRESS) {
1433: f[j][i].u = XMomentumResidual(x,i,j,user);
1434: } else {
1435: /* experimental boundary condition */
1436: }
1438: } else {
1439: /* in the mantle wedge */
1440: f[j][i].u = XMomentumResidual(x,i,j,user);
1441: }
1443: /************* Z-MOMENTUM/VELOCITY *************/
1444: if (i<=j) {
1445: f[j][i].w = x[j][i].w - SlabVel('W',i,j,user);
1447: } else if (j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1448: /* in the lithospheric lid */
1449: f[j][i].w = x[j][i].w - 0.0;
1451: } else if (j==jlim) {
1452: /* on the bottom boundary */
1453: if (ibound==BC_ANALYTIC) {
1454: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1455: } else {
1456: f[j][i].w = ZNormalStress(x,i,j,CELL_CENTER,user) - EPS_ZERO;
1457: }
1459: } else if (i==ilim) {
1460: /* on the right side boundary */
1461: if (ibound==BC_ANALYTIC) {
1462: f[j][i].w = x[j][i].w - VertVelocity(i,j,user);
1463: } else if (ibound==BC_NOSTRESS) {
1464: f[j][i].w = ZMomentumResidual(x,i,j,user);
1465: } else {
1466: /* experimental boundary condition */
1467: }
1469: } else {
1470: /* in the mantle wedge */
1471: f[j][i].w = ZMomentumResidual(x,i,j,user);
1472: }
1474: /************* CONTINUITY/PRESSURE *************/
1475: if (i<j || j<=grid->jlid || (j<grid->corner+grid->inose && i<grid->corner+grid->inose)) {
1476: /* in the lid or slab */
1477: f[j][i].p = x[j][i].p;
1479: } else if ((i==ilim || j==jlim) && ibound==BC_ANALYTIC) {
1480: /* on an analytic boundary */
1481: f[j][i].p = x[j][i].p - Pressure(i,j,user);
1483: } else {
1484: /* in the mantle wedge */
1485: f[j][i].p = ContinuityResidual(x,i,j,user);
1486: }
1488: /************* TEMPERATURE *************/
1489: if (j==0) {
1490: /* on the surface */
1491: f[j][i].T = x[j][i].T + x[j+1][i].T + PetscMax(PetscRealPart(x[j][i].T),0.0);
1493: } else if (i==0) {
1494: /* slab inflow boundary */
1495: f[j][i].T = x[j][i].T - PlateModel(j,PLATE_SLAB,user);
1497: } else if (i==ilim) {
1498: /* right side boundary */
1499: mag_u = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j][i-1].u)/param->cb,1.0),0.0)), 5);
1500: 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);
1502: } else if (j==jlim) {
1503: /* bottom boundary */
1504: mag_w = 1.0 - PetscPowRealInt((1.0-PetscMax(PetscMin(PetscRealPart(x[j-1][i].w)/param->sb,1.0),0.0)), 5);
1505: f[j][i].T = x[j][i].T - mag_w*x[j-1][i-1].T - (1.0-mag_w);
1507: } else {
1508: /* in the mantle wedge */
1509: f[j][i].T = EnergyResidual(x,i,j,user);
1510: }
1511: }
1512: }
1513: return(0);
1514: }