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";
27: /*F-----------------------------------------------------------------------
29: This PETSc 2.2.0 example by Richard F. Katz
30: http://www.ldeo.columbia.edu/~katz/
32: The problem is modeled by the partial differential equation system
34: \begin{eqnarray}
35: -\nabla P + \nabla \cdot [\eta (\nabla v + \nabla v^T)] & = & 0 \\
36: \nabla \cdot v & = & 0 \\
37: dT/dt + \nabla \cdot (vT) - 1/Pe \triangle^2(T) & = & 0 \\
38: \end{eqnarray}
40: \begin{eqnarray}
41: \eta(T,Eps\_dot) & = & \hbox{constant } \hbox{if ivisc} ==0 \\
42: & = & \hbox{diffusion creep (T,P-dependent) } \hbox{if ivisc} ==1 \\
43: & = & \hbox{dislocation creep (T,P,v-dependent)} \hbox{if ivisc} ==2 \\
44: & = & \hbox{mantle viscosity (difn and disl) } \hbox{if ivisc} ==3
45: \end{eqnarray}
47: which is uniformly discretized on a staggered mesh:
48: -------$w_{ij}$------
49: $u_{i-1j}$ $P,T_{ij}$ $u_{ij}$
50: ------$w_{ij-1}$-----
52: ------------------------------------------------------------------------F*/
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 miscellaneous 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: int main(int argc, char **argv)
124: {
125: SNES snes;
126: AppCtx *user; /* user-defined work context */
127: Parameter param;
128: GridInfo grid;
129: PetscInt nits;
130: MPI_Comm comm;
131: DM da;
133: PetscFunctionBeginUser;
134: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
135: PetscCall(PetscOptionsSetValue(NULL, "-file", "ex30_output"));
136: PetscCall(PetscOptionsSetValue(NULL, "-snes_monitor_short", NULL));
137: PetscCall(PetscOptionsSetValue(NULL, "-snes_max_it", "20"));
138: PetscCall(PetscOptionsSetValue(NULL, "-ksp_max_it", "1500"));
139: PetscCall(PetscOptionsSetValue(NULL, "-ksp_gmres_restart", "300"));
140: PetscCall(PetscOptionsInsert(NULL, &argc, &argv, NULL));
142: comm = PETSC_COMM_WORLD;
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Set up the problem parameters.
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: PetscCall(SetParams(¶m, &grid));
148: PetscCall(ReportParams(¶m, &grid));
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: PetscCall(SNESCreate(comm, &snes));
153: PetscCall(DMDACreate2d(comm, grid.bx, grid.by, grid.stencil, grid.ni, grid.nj, PETSC_DECIDE, PETSC_DECIDE, grid.dof, grid.stencil_width, 0, 0, &da));
154: PetscCall(DMSetFromOptions(da));
155: PetscCall(DMSetUp(da));
156: PetscCall(SNESSetDM(snes, da));
157: PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
158: PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
159: PetscCall(DMDASetFieldName(da, 2, "pressure"));
160: PetscCall(DMDASetFieldName(da, 3, "temperature"));
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Create user context, set problem data, create vector data structures.
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: PetscCall(PetscNew(&user));
166: user->param = ¶m;
167: user->grid = &grid;
168: PetscCall(DMSetApplicationContext(da, user));
169: PetscCall(DMCreateGlobalVector(da, &user->Xguess));
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Set up the SNES solver with callback functions.
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: PetscCall(DMDASNESSetFunctionLocal(da, INSERT_VALUES, (PetscErrorCode(*)(DMDALocalInfo *, void *, void *, void *))FormFunctionLocal, (void *)user));
175: PetscCall(SNESSetFromOptions(snes));
177: PetscCall(SNESSetConvergenceTest(snes, SNESConverged_Interactive, (void *)user, NULL));
178: PetscCall(PetscPushSignalHandler(InteractiveHandler, (void *)user));
180: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: Initialize and solve the nonlinear system
182: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: PetscCall(Initialize(da));
184: PetscCall(UpdateSolution(snes, user, &nits));
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Output variables.
188: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: PetscCall(DoOutput(snes, nits));
191: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: Free work space.
193: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194: PetscCall(VecDestroy(&user->Xguess));
195: PetscCall(VecDestroy(&user->x));
196: PetscCall(PetscFree(user));
197: PetscCall(SNESDestroy(&snes));
198: PetscCall(DMDestroy(&da));
199: PetscCall(PetscPopSignalHandler());
200: PetscCall(PetscFinalize());
201: return 0;
202: }
204: /*=====================================================================
205: PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
206: =====================================================================*/
208: /* manages solve: adaptive continuation method */
209: PetscErrorCode UpdateSolution(SNES snes, AppCtx *user, PetscInt *nits)
210: {
211: KSP ksp;
212: PC pc;
213: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
214: Parameter *param = user->param;
215: PetscReal cont_incr = 0.3;
216: PetscInt its;
217: PetscBool q = PETSC_FALSE;
218: DM dm;
220: PetscFunctionBeginUser;
221: PetscCall(SNESGetDM(snes, &dm));
222: PetscCall(DMCreateGlobalVector(dm, &user->x));
223: PetscCall(SNESGetKSP(snes, &ksp));
224: PetscCall(KSPGetPC(ksp, &pc));
225: PetscCall(KSPSetComputeSingularValues(ksp, PETSC_TRUE));
227: *nits = 0;
229: /* Isoviscous solve */
230: if (param->ivisc == VISC_CONST && !param->stop_solve) {
231: param->ivisc = VISC_CONST;
233: PetscCall(SNESSolve(snes, 0, user->x));
234: PetscCall(SNESGetConvergedReason(snes, &reason));
235: PetscCall(SNESGetIterationNumber(snes, &its));
236: *nits += its;
237: PetscCall(VecCopy(user->x, user->Xguess));
238: if (param->stop_solve) goto done;
239: }
241: /* Olivine diffusion creep */
242: if (param->ivisc >= VISC_DIFN && !param->stop_solve) {
243: if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Computing Variable Viscosity Solution\n"));
245: /* continuation method on viscosity cutoff */
246: for (param->continuation = 0.0;; param->continuation += cont_incr) {
247: if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Continuation parameter = %g\n", (double)param->continuation));
249: /* solve the non-linear system */
250: PetscCall(VecCopy(user->Xguess, user->x));
251: PetscCall(SNESSolve(snes, 0, user->x));
252: PetscCall(SNESGetConvergedReason(snes, &reason));
253: PetscCall(SNESGetIterationNumber(snes, &its));
254: *nits += its;
255: if (!q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " SNES iterations: %" PetscInt_FMT ", Cumulative: %" PetscInt_FMT "\n", its, *nits));
256: if (param->stop_solve) goto done;
258: if (reason < 0) {
259: /* NOT converged */
260: cont_incr = -PetscAbsReal(cont_incr) / 2.0;
261: if (PetscAbsReal(cont_incr) < 0.01) goto done;
263: } else {
264: /* converged */
265: PetscCall(VecCopy(user->x, user->Xguess));
266: if (param->continuation >= 1.0) goto done;
267: if (its <= 3) cont_incr = 0.30001;
268: else if (its <= 8) cont_incr = 0.15001;
269: else cont_incr = 0.10001;
271: if (param->continuation + cont_incr > 1.0) cont_incr = 1.0 - param->continuation;
272: } /* endif reason<0 */
273: }
274: }
275: done:
276: if (param->stop_solve && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: stopping solve.\n"));
277: if (reason < 0 && !q) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "FAILED TO CONVERGE: stopping solve.\n"));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*=====================================================================
282: PHYSICS FUNCTIONS (compute the discrete residual)
283: =====================================================================*/
285: static inline PetscScalar UInterp(Field **x, PetscInt i, PetscInt j)
286: {
287: return 0.25 * (x[j][i].u + x[j + 1][i].u + x[j][i + 1].u + x[j + 1][i + 1].u);
288: }
290: static inline PetscScalar WInterp(Field **x, PetscInt i, PetscInt j)
291: {
292: return 0.25 * (x[j][i].w + x[j + 1][i].w + x[j][i + 1].w + x[j + 1][i + 1].w);
293: }
295: static inline PetscScalar PInterp(Field **x, PetscInt i, PetscInt j)
296: {
297: return 0.25 * (x[j][i].p + x[j + 1][i].p + x[j][i + 1].p + x[j + 1][i + 1].p);
298: }
300: static inline PetscScalar TInterp(Field **x, PetscInt i, PetscInt j)
301: {
302: return 0.25 * (x[j][i].T + x[j + 1][i].T + x[j][i + 1].T + x[j + 1][i + 1].T);
303: }
305: /* isoviscous analytic solution for IC */
306: static inline PetscScalar HorizVelocity(PetscInt i, PetscInt j, AppCtx *user)
307: {
308: Parameter *param = user->param;
309: GridInfo *grid = user->grid;
310: PetscScalar st, ct, th, c = param->c, d = param->d;
311: PetscReal x, z, r;
313: x = (i - grid->jlid) * grid->dx;
314: z = (j - grid->jlid - 0.5) * grid->dz;
315: r = PetscSqrtReal(x * x + z * z);
316: st = z / r;
317: ct = x / r;
318: th = PetscAtanReal(z / x);
319: return ct * (c * th * st + d * (st + th * ct)) + st * (c * (st - th * ct) + d * th * st);
320: }
322: /* isoviscous analytic solution for IC */
323: static inline PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
325: {
326: Parameter *param = user->param;
327: GridInfo *grid = user->grid;
328: PetscScalar st, ct, th, c = param->c, d = param->d;
329: PetscReal x, z, r;
331: x = (i - grid->jlid - 0.5) * grid->dx;
332: z = (j - grid->jlid) * grid->dz;
333: r = PetscSqrtReal(x * x + z * z);
334: st = z / r;
335: ct = x / r;
336: th = PetscAtanReal(z / x);
337: return st * (c * th * st + d * (st + th * ct)) - ct * (c * (st - th * ct) + d * th * st);
338: }
340: /* isoviscous analytic solution for IC */
341: static inline PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
342: {
343: Parameter *param = user->param;
344: GridInfo *grid = user->grid;
345: PetscScalar x, z, r, st, ct, c = param->c, d = param->d;
347: x = (i - grid->jlid - 0.5) * grid->dx;
348: z = (j - grid->jlid - 0.5) * grid->dz;
349: r = PetscSqrtReal(x * x + z * z);
350: st = z / r;
351: ct = x / r;
352: return -2.0 * (c * ct - d * st) / r;
353: }
355: /* computes the second invariant of the strain rate tensor */
356: static inline PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
357: {
358: Parameter *param = user->param;
359: GridInfo *grid = user->grid;
360: PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1;
361: PetscScalar uN, uS, uE, uW, wN, wS, wE, wW;
362: PetscScalar eps11, eps12, eps22;
364: if (i < j) return EPS_ZERO;
365: if (i == ilim) i--;
366: if (j == jlim) j--;
368: if (ipos == CELL_CENTER) { /* on cell center */
369: if (j <= grid->jlid) return EPS_ZERO;
371: uE = x[j][i].u;
372: uW = x[j][i - 1].u;
373: wN = x[j][i].w;
374: wS = x[j - 1][i].w;
375: wE = WInterp(x, i, j - 1);
376: if (i == j) {
377: uN = param->cb;
378: wW = param->sb;
379: } else {
380: uN = UInterp(x, i - 1, j);
381: wW = WInterp(x, i - 1, j - 1);
382: }
384: if (j == grid->jlid + 1) uS = 0.0;
385: else uS = UInterp(x, i - 1, j - 1);
387: } else { /* on CELL_CORNER */
388: if (j < grid->jlid) return EPS_ZERO;
390: uN = x[j + 1][i].u;
391: uS = x[j][i].u;
392: wE = x[j][i + 1].w;
393: wW = x[j][i].w;
394: if (i == j) {
395: wN = param->sb;
396: uW = param->cb;
397: } else {
398: wN = WInterp(x, i, j);
399: uW = UInterp(x, i - 1, j);
400: }
402: if (j == grid->jlid) {
403: uE = 0.0;
404: uW = 0.0;
405: uS = -uN;
406: wS = -wN;
407: } else {
408: uE = UInterp(x, i, j);
409: wS = WInterp(x, i, j - 1);
410: }
411: }
413: eps11 = (uE - uW) / grid->dx;
414: 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: /* computes the shear viscosity */
421: static inline PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
422: {
423: PetscReal result = 0.0;
424: ViscParam difn = param->diffusion, disl = param->dislocation;
425: PetscInt iVisc = param->ivisc;
426: PetscScalar eps_scale = param->V / (param->L * 1000.0);
427: PetscScalar strain_power, v1, v2, P;
428: PetscScalar rho_g = 32340.0, R = 8.3144;
430: P = rho_g * (z * param->L * 1000.0); /* Pa */
432: if (iVisc == VISC_CONST) {
433: /* constant viscosity */
434: return 1.0;
435: } else if (iVisc == VISC_DIFN) {
436: /* diffusion creep rheology */
437: result = PetscRealPart(difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0);
438: } else if (iVisc == VISC_DISL) {
439: /* dislocation creep rheology */
440: strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);
442: result = PetscRealPart(disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0);
443: } else if (iVisc == VISC_FULL) {
444: /* dislocation/diffusion creep rheology */
445: strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);
447: v1 = difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0;
448: v2 = disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0;
450: result = PetscRealPart(1.0 / (1.0 / v1 + 1.0 / v2));
451: }
453: /* max viscosity is param->eta0 */
454: result = PetscMin(result, 1.0);
455: /* min viscosity is param->visc_cutoff */
456: result = PetscMax(result, param->visc_cutoff);
457: /* continuation method */
458: result = PetscPowReal(result, param->continuation);
459: return result;
460: }
462: /* computes the residual of the x-component of eqn (1) above */
463: static inline PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
464: {
465: Parameter *param = user->param;
466: GridInfo *grid = user->grid;
467: PetscScalar dx = grid->dx, dz = grid->dz;
468: PetscScalar etaN, etaS, etaE, etaW, epsN = 0.0, epsS = 0.0, epsE = 0.0, epsW = 0.0;
469: PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdx, residual, z_scale;
470: PetscScalar dudxW, dudxE, dudzN, dudzS, dwdxN, dwdxS;
471: PetscInt jlim = grid->nj - 1;
473: z_scale = param->z_scale;
475: if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
476: TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale);
477: if (j == jlim) TN = TS;
478: else TN = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
479: TW = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
480: TE = param->potentialT * x[j][i + 1].T * PetscExpScalar((j - 0.5) * dz * z_scale);
481: if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
482: epsN = CalcSecInv(x, i, j, CELL_CORNER, user);
483: epsS = CalcSecInv(x, i, j - 1, CELL_CORNER, user);
484: epsE = CalcSecInv(x, i + 1, j, CELL_CENTER, user);
485: epsW = CalcSecInv(x, i, j, CELL_CENTER, user);
486: }
487: }
488: etaN = Viscosity(TN, epsN, dz * (j + 0.5), param);
489: etaS = Viscosity(TS, epsS, dz * (j - 0.5), param);
490: etaW = Viscosity(TW, epsW, dz * j, param);
491: etaE = Viscosity(TE, epsE, dz * j, param);
493: dPdx = (x[j][i + 1].p - x[j][i].p) / dx;
494: if (j == jlim) dudzN = etaN * (x[j][i].w - x[j][i + 1].w) / dx;
495: else dudzN = etaN * (x[j + 1][i].u - x[j][i].u) / dz;
496: dudzS = etaS * (x[j][i].u - x[j - 1][i].u) / dz;
497: dudxE = etaE * (x[j][i + 1].u - x[j][i].u) / dx;
498: dudxW = etaW * (x[j][i].u - x[j][i - 1].u) / dx;
500: residual = -dPdx /* X-MOMENTUM EQUATION*/
501: + (dudxE - dudxW) / dx + (dudzN - dudzS) / dz;
503: if (param->ivisc != VISC_CONST) {
504: dwdxN = etaN * (x[j][i + 1].w - x[j][i].w) / dx;
505: dwdxS = etaS * (x[j - 1][i + 1].w - x[j - 1][i].w) / dx;
507: residual += (dudxE - dudxW) / dx + (dwdxN - dwdxS) / dz;
508: }
510: return residual;
511: }
513: /* computes the residual of the z-component of eqn (1) above */
514: static inline PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
516: {
517: Parameter *param = user->param;
518: GridInfo *grid = user->grid;
519: PetscScalar dx = grid->dx, dz = grid->dz;
520: 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;
521: PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdz, residual, z_scale;
522: PetscScalar dudzE, dudzW, dwdxW, dwdxE, dwdzN, dwdzS;
523: PetscInt ilim = grid->ni - 1;
525: /* geometric and other parameters */
526: z_scale = param->z_scale;
528: /* viscosity */
529: if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
530: TN = param->potentialT * x[j + 1][i].T * PetscExpScalar((j + 0.5) * dz * z_scale);
531: TS = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
532: TW = param->potentialT * TInterp(x, i - 1, j) * PetscExpScalar(j * dz * z_scale);
533: if (i == ilim) TE = TW;
534: else TE = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
535: if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
536: epsN = CalcSecInv(x, i, j + 1, CELL_CENTER, user);
537: epsS = CalcSecInv(x, i, j, CELL_CENTER, user);
538: epsE = CalcSecInv(x, i, j, CELL_CORNER, user);
539: epsW = CalcSecInv(x, i - 1, j, CELL_CORNER, user);
540: }
541: }
542: etaN = Viscosity(TN, epsN, dz * (j + 1.0), param);
543: etaS = Viscosity(TS, epsS, dz * (j + 0.0), param);
544: etaW = Viscosity(TW, epsW, dz * (j + 0.5), param);
545: etaE = Viscosity(TE, epsE, dz * (j + 0.5), param);
547: dPdz = (x[j + 1][i].p - x[j][i].p) / dz;
548: dwdzN = etaN * (x[j + 1][i].w - x[j][i].w) / dz;
549: dwdzS = etaS * (x[j][i].w - x[j - 1][i].w) / dz;
550: if (i == ilim) dwdxE = etaE * (x[j][i].u - x[j + 1][i].u) / dz;
551: else dwdxE = etaE * (x[j][i + 1].w - x[j][i].w) / dx;
552: dwdxW = 2.0 * etaW * (x[j][i].w - x[j][i - 1].w) / dx;
554: /* Z-MOMENTUM */
555: residual = -dPdz /* constant viscosity terms */
556: + (dwdzN - dwdzS) / dz + (dwdxE - dwdxW) / dx;
558: if (param->ivisc != VISC_CONST) {
559: dudzE = etaE * (x[j + 1][i].u - x[j][i].u) / dz;
560: dudzW = etaW * (x[j + 1][i - 1].u - x[j][i - 1].u) / dz;
562: residual += (dwdzN - dwdzS) / dz + (dudzE - dudzW) / dx;
563: }
565: return residual;
566: }
568: /* computes the residual of eqn (2) above */
569: static inline PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
570: {
571: GridInfo *grid = user->grid;
572: PetscScalar uE, uW, wN, wS, dudx, dwdz;
574: uW = x[j][i - 1].u;
575: uE = x[j][i].u;
576: dudx = (uE - uW) / grid->dx;
577: wS = x[j - 1][i].w;
578: wN = x[j][i].w;
579: dwdz = (wN - wS) / grid->dz;
581: return dudx + dwdz;
582: }
584: /* computes the residual of eqn (3) above */
585: static inline PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
586: {
587: Parameter *param = user->param;
588: GridInfo *grid = user->grid;
589: PetscScalar dx = grid->dx, dz = grid->dz;
590: PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, jlid = grid->jlid;
591: PetscScalar TE, TN, TS, TW, residual;
592: PetscScalar uE, uW, wN, wS;
593: PetscScalar fN, fS, fE, fW, dTdxW, dTdxE, dTdzN, dTdzS;
595: dTdzN = (x[j + 1][i].T - x[j][i].T) / dz;
596: dTdzS = (x[j][i].T - x[j - 1][i].T) / dz;
597: dTdxE = (x[j][i + 1].T - x[j][i].T) / dx;
598: dTdxW = (x[j][i].T - x[j][i - 1].T) / dx;
600: residual = ((dTdzN - dTdzS) / dz + /* diffusion term */
601: (dTdxE - dTdxW) / dx) *
602: dx * dz / param->peclet;
604: if (j <= jlid && i >= j) {
605: /* don't advect in the lid */
606: return residual;
607: } else if (i < j) {
608: /* beneath the slab sfc */
609: uW = uE = param->cb;
610: wS = wN = param->sb;
611: } else {
612: /* advect in the slab and wedge */
613: uW = x[j][i - 1].u;
614: uE = x[j][i].u;
615: wS = x[j - 1][i].w;
616: wN = x[j][i].w;
617: }
619: if (param->adv_scheme == ADVECT_FV || i == ilim - 1 || j == jlim - 1 || i == 1 || j == 1) {
620: /* finite volume advection */
621: TS = (x[j][i].T + x[j - 1][i].T) / 2.0;
622: TN = (x[j][i].T + x[j + 1][i].T) / 2.0;
623: TE = (x[j][i].T + x[j][i + 1].T) / 2.0;
624: TW = (x[j][i].T + x[j][i - 1].T) / 2.0;
625: fN = wN * TN * dx;
626: fS = wS * TS * dx;
627: fE = uE * TE * dz;
628: fW = uW * TW * dz;
630: } else {
631: /* Fromm advection scheme */
632: 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 - 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;
633: 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 - 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;
634: 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 - 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;
635: 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 - 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;
636: }
638: residual -= (fE - fW + fN - fS);
640: return residual;
641: }
643: /* computes the shear stress---used on the boundaries */
644: static inline PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
645: {
646: Parameter *param = user->param;
647: GridInfo *grid = user->grid;
648: PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1;
649: PetscScalar uN, uS, wE, wW;
651: if (j <= grid->jlid || i < j || i == ilim || j == jlim) return EPS_ZERO;
653: if (ipos == CELL_CENTER) { /* on cell center */
655: wE = WInterp(x, i, j - 1);
656: if (i == j) {
657: wW = param->sb;
658: uN = param->cb;
659: } else {
660: wW = WInterp(x, i - 1, j - 1);
661: uN = UInterp(x, i - 1, j);
662: }
663: if (j == grid->jlid + 1) uS = 0.0;
664: else uS = UInterp(x, i - 1, j - 1);
666: } else { /* on cell corner */
668: uN = x[j + 1][i].u;
669: uS = x[j][i].u;
670: wW = x[j][i].w;
671: wE = x[j][i + 1].w;
672: }
674: return (uN - uS) / grid->dz + (wE - wW) / grid->dx;
675: }
677: /* computes the normal stress---used on the boundaries */
678: static inline PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
679: {
680: Parameter *param = user->param;
681: GridInfo *grid = user->grid;
682: PetscScalar dx = grid->dx, dz = grid->dz;
683: PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
684: PetscScalar epsC = 0.0, etaC, TC, uE, uW, pC, z_scale;
685: if (i < j || j <= grid->jlid) return EPS_ZERO;
687: ivisc = param->ivisc;
688: z_scale = param->z_scale;
690: if (ipos == CELL_CENTER) { /* on cell center */
692: TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
693: if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
694: etaC = Viscosity(TC, epsC, dz * j, param);
696: uW = x[j][i - 1].u;
697: uE = x[j][i].u;
698: pC = x[j][i].p;
700: } else { /* on cell corner */
701: if (i == ilim || j == jlim) return EPS_ZERO;
703: TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
704: if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
705: etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);
707: if (i == j) uW = param->sb;
708: else uW = UInterp(x, i - 1, j);
709: uE = UInterp(x, i, j);
710: pC = PInterp(x, i, j);
711: }
713: return 2.0 * etaC * (uE - uW) / dx - pC;
714: }
716: /* computes the normal stress---used on the boundaries */
717: static inline PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
718: {
719: Parameter *param = user->param;
720: GridInfo *grid = user->grid;
721: PetscScalar dz = grid->dz;
722: PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
723: PetscScalar epsC = 0.0, etaC, TC;
724: PetscScalar pC, wN, wS, z_scale;
725: if (i < j || j <= grid->jlid) return EPS_ZERO;
727: ivisc = param->ivisc;
728: z_scale = param->z_scale;
730: if (ipos == CELL_CENTER) { /* on cell center */
732: TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
733: if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
734: etaC = Viscosity(TC, epsC, dz * j, param);
735: wN = x[j][i].w;
736: wS = x[j - 1][i].w;
737: pC = x[j][i].p;
739: } else { /* on cell corner */
740: if ((i == ilim) || (j == jlim)) return EPS_ZERO;
742: TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
743: if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
744: etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);
745: if (i == j) wN = param->sb;
746: else wN = WInterp(x, i, j);
747: wS = WInterp(x, i, j - 1);
748: pC = PInterp(x, i, j);
749: }
751: return 2.0 * etaC * (wN - wS) / dz - pC;
752: }
754: /*=====================================================================
755: INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
756: =====================================================================*/
758: /* initializes the problem parameters and checks for
759: command line changes */
760: PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
761: {
762: PetscReal SEC_PER_YR = 3600.00 * 24.00 * 365.2500;
763: PetscReal alpha_g_on_cp_units_inverse_km = 4.0e-5 * 9.8;
765: PetscFunctionBeginUser;
766: /* domain geometry */
767: param->slab_dip = 45.0;
768: param->width = 320.0; /* km */
769: param->depth = 300.0; /* km */
770: param->lid_depth = 35.0; /* km */
771: param->fault_depth = 35.0; /* km */
773: PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_dip", ¶m->slab_dip, NULL));
774: PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", ¶m->width, NULL));
775: PetscCall(PetscOptionsGetReal(NULL, NULL, "-depth", ¶m->depth, NULL));
776: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_depth", ¶m->lid_depth, NULL));
777: PetscCall(PetscOptionsGetReal(NULL, NULL, "-fault_depth", ¶m->fault_depth, NULL));
779: param->slab_dip = param->slab_dip * PETSC_PI / 180.0; /* radians */
781: /* grid information */
782: PetscCall(PetscOptionsGetInt(NULL, NULL, "-jfault", &grid->jfault, NULL));
783: grid->ni = 82;
784: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ni", &grid->ni, NULL));
786: grid->dx = param->width / ((PetscReal)(grid->ni - 2)); /* km */
787: grid->dz = grid->dx * PetscTanReal(param->slab_dip); /* km */
788: grid->nj = (PetscInt)(param->depth / grid->dz + 3.0); /* gridpoints*/
789: param->depth = grid->dz * (grid->nj - 2); /* km */
790: grid->inose = 0; /* gridpoints*/
791: PetscCall(PetscOptionsGetInt(NULL, NULL, "-inose", &grid->inose, NULL));
792: grid->bx = DM_BOUNDARY_NONE;
793: grid->by = DM_BOUNDARY_NONE;
794: grid->stencil = DMDA_STENCIL_BOX;
795: grid->dof = 4;
796: grid->stencil_width = 2;
797: grid->mglevels = 1;
799: /* boundary conditions */
800: param->pv_analytic = PETSC_FALSE;
801: param->ibound = BC_NOSTRESS;
802: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ibound", ¶m->ibound, NULL));
804: /* physical constants */
805: param->slab_velocity = 5.0; /* cm/yr */
806: param->slab_age = 50.0; /* Ma */
807: param->lid_age = 50.0; /* Ma */
808: param->kappa = 0.7272e-6; /* m^2/sec */
809: param->potentialT = 1300.0; /* degrees C */
811: PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_velocity", ¶m->slab_velocity, NULL));
812: PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_age", ¶m->slab_age, NULL));
813: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_age", ¶m->lid_age, NULL));
814: PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", ¶m->kappa, NULL));
815: PetscCall(PetscOptionsGetReal(NULL, NULL, "-potentialT", ¶m->potentialT, NULL));
817: /* viscosity */
818: param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
819: param->eta0 = 1e24; /* Pa-s */
820: param->visc_cutoff = 0.0; /* factor of eta_0 */
821: param->continuation = 1.0;
823: /* constants for diffusion creep */
824: param->diffusion.A = 1.8e7; /* Pa-s */
825: param->diffusion.n = 1.0; /* dim'less */
826: param->diffusion.Estar = 375e3; /* J/mol */
827: param->diffusion.Vstar = 5e-6; /* m^3/mol */
829: /* constants for param->dislocationocation creep */
830: param->dislocation.A = 2.8969e4; /* Pa-s */
831: param->dislocation.n = 3.5; /* dim'less */
832: param->dislocation.Estar = 530e3; /* J/mol */
833: param->dislocation.Vstar = 14e-6; /* m^3/mol */
835: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ivisc", ¶m->ivisc, NULL));
836: PetscCall(PetscOptionsGetReal(NULL, NULL, "-visc_cutoff", ¶m->visc_cutoff, NULL));
838: param->output_ivisc = param->ivisc;
840: PetscCall(PetscOptionsGetInt(NULL, NULL, "-output_ivisc", ¶m->output_ivisc, NULL));
841: PetscCall(PetscOptionsGetReal(NULL, NULL, "-vstar", ¶m->dislocation.Vstar, NULL));
843: /* output options */
844: param->quiet = PETSC_FALSE;
845: param->param_test = PETSC_FALSE;
847: PetscCall(PetscOptionsHasName(NULL, NULL, "-quiet", ¶m->quiet));
848: PetscCall(PetscOptionsHasName(NULL, NULL, "-test", ¶m->param_test));
849: PetscCall(PetscOptionsGetString(NULL, NULL, "-file", param->filename, sizeof(param->filename), ¶m->output_to_file));
851: /* advection */
852: param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
854: PetscCall(PetscOptionsGetInt(NULL, NULL, "-adv_scheme", ¶m->adv_scheme, NULL));
856: /* misc. flags */
857: param->stop_solve = PETSC_FALSE;
858: param->interrupted = PETSC_FALSE;
859: param->kspmon = PETSC_FALSE;
860: param->toggle_kspmon = PETSC_FALSE;
862: /* derived parameters for slab angle */
863: param->sb = PetscSinReal(param->slab_dip);
864: param->cb = PetscCosReal(param->slab_dip);
865: param->c = param->slab_dip * param->sb / (param->slab_dip * param->slab_dip - param->sb * param->sb);
866: param->d = (param->slab_dip * param->cb - param->sb) / (param->slab_dip * param->slab_dip - param->sb * param->sb);
868: /* length, velocity and time scale for non-dimensionalization */
869: param->L = PetscMin(param->width, param->depth); /* km */
870: param->V = param->slab_velocity / 100.0 / SEC_PER_YR; /* m/sec */
872: /* other unit conversions and derived parameters */
873: param->scaled_width = param->width / param->L; /* dim'less */
874: param->scaled_depth = param->depth / param->L; /* dim'less */
875: param->lid_depth = param->lid_depth / param->L; /* dim'less */
876: param->fault_depth = param->fault_depth / param->L; /* dim'less */
877: grid->dx = grid->dx / param->L; /* dim'less */
878: grid->dz = grid->dz / param->L; /* dim'less */
879: grid->jlid = (PetscInt)(param->lid_depth / grid->dz); /* gridcells */
880: grid->jfault = (PetscInt)(param->fault_depth / grid->dz); /* gridcells */
881: param->lid_depth = grid->jlid * grid->dz; /* dim'less */
882: param->fault_depth = grid->jfault * grid->dz; /* dim'less */
883: grid->corner = grid->jlid + 1; /* gridcells */
884: param->peclet = param->V /* m/sec */
885: * param->L * 1000.0 /* m */
886: / param->kappa; /* m^2/sec */
887: param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
888: param->skt = PetscSqrtReal(param->kappa * param->slab_age * SEC_PER_YR);
889: PetscCall(PetscOptionsGetReal(NULL, NULL, "-peclet", ¶m->peclet, NULL));
890: PetscFunctionReturn(PETSC_SUCCESS);
891: }
893: /* prints a report of the problem parameters to stdout */
894: PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
895: {
896: char date[30];
898: PetscFunctionBeginUser;
899: PetscCall(PetscGetDate(date, 30));
901: if (!param->quiet) {
902: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------BEGIN ex30 PARAM REPORT-------------------\n"));
903: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Domain: \n"));
904: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Width = %g km, Depth = %g km\n", (double)param->width, (double)param->depth));
905: PetscCall(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));
906: PetscCall(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)));
908: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nGrid: \n"));
909: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " [ni,nj] = %" PetscInt_FMT ", %" PetscInt_FMT " [dx,dz] = %g, %g km\n", grid->ni, grid->nj, (double)(grid->dx * param->L), (double)(grid->dz * param->L)));
910: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " jlid = %3" PetscInt_FMT " jfault = %3" PetscInt_FMT " \n", grid->jlid, grid->jfault));
911: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pe = %g\n", (double)param->peclet));
913: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nRheology:"));
914: if (param->ivisc == VISC_CONST) {
915: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Isoviscous \n"));
916: if (param->pv_analytic) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pressure and Velocity prescribed! \n"));
917: } else if (param->ivisc == VISC_DIFN) {
918: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Diffusion Creep (T-Dependent Newtonian) \n"));
919: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
920: } else if (param->ivisc == VISC_DISL) {
921: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Dislocation Creep (T-Dependent Non-Newtonian) \n"));
922: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
923: } else if (param->ivisc == VISC_FULL) {
924: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Full Rheology \n"));
925: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
926: } else {
927: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Invalid! \n"));
928: PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
929: }
931: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Boundary condition:"));
932: if (param->ibound == BC_ANALYTIC) {
933: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Isoviscous Analytic Dirichlet \n"));
934: } else if (param->ibound == BC_NOSTRESS) {
935: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Stress-Free (normal & shear stress)\n"));
936: } else if (param->ibound == BC_EXPERMNT) {
937: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Experimental boundary condition \n"));
938: } else {
939: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Invalid! \n"));
940: PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
941: }
943: if (param->output_to_file) {
944: #if defined(PETSC_HAVE_MATLAB)
945: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: Mat file \"%s\"\n", param->filename));
946: #else
947: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: PETSc binary file \"%s\"\n", param->filename));
948: #endif
949: }
950: if (param->output_ivisc != param->ivisc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Output viscosity: -ivisc %" PetscInt_FMT "\n", param->output_ivisc));
952: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------END ex30 PARAM REPORT---------------------\n"));
953: }
954: if (param->param_test) PetscCall(PetscEnd());
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: /* ------------------------------------------------------------------- */
959: /* generates an initial guess using the analytic solution for isoviscous
960: corner flow */
961: PetscErrorCode Initialize(DM da)
962: /* ------------------------------------------------------------------- */
963: {
964: AppCtx *user;
965: Parameter *param;
966: GridInfo *grid;
967: PetscInt i, j, is, js, im, jm;
968: Field **x;
969: Vec Xguess;
971: PetscFunctionBeginUser;
972: /* Get the fine grid */
973: PetscCall(DMGetApplicationContext(da, &user));
974: Xguess = user->Xguess;
975: param = user->param;
976: grid = user->grid;
977: PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
978: PetscCall(DMDAVecGetArray(da, Xguess, (void **)&x));
980: /* Compute initial guess */
981: for (j = js; j < js + jm; j++) {
982: for (i = is; i < is + im; i++) {
983: if (i < j) x[j][i].u = param->cb;
984: else if (j <= grid->jlid) x[j][i].u = 0.0;
985: else x[j][i].u = HorizVelocity(i, j, user);
987: if (i <= j) x[j][i].w = param->sb;
988: else if (j <= grid->jlid) x[j][i].w = 0.0;
989: else x[j][i].w = VertVelocity(i, j, user);
991: if (i < j || j <= grid->jlid) x[j][i].p = 0.0;
992: else x[j][i].p = Pressure(i, j, user);
994: x[j][i].T = PetscMin(grid->dz * (j - 0.5), 1.0);
995: }
996: }
998: /* Restore x to Xguess */
999: PetscCall(DMDAVecRestoreArray(da, Xguess, (void **)&x));
1000: PetscFunctionReturn(PETSC_SUCCESS);
1001: }
1003: /* controls output to a file */
1004: PetscErrorCode DoOutput(SNES snes, PetscInt its)
1005: {
1006: AppCtx *user;
1007: Parameter *param;
1008: GridInfo *grid;
1009: PetscInt ivt;
1010: PetscMPIInt rank;
1011: PetscViewer viewer;
1012: Vec res, pars;
1013: MPI_Comm comm;
1014: DM da;
1016: PetscFunctionBeginUser;
1017: PetscCall(SNESGetDM(snes, &da));
1018: PetscCall(DMGetApplicationContext(da, &user));
1019: param = user->param;
1020: grid = user->grid;
1021: ivt = param->ivisc;
1023: param->ivisc = param->output_ivisc;
1025: /* compute final residual and final viscosity/strain rate fields */
1026: PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
1027: PetscCall(ViscosityField(da, user->x, user->Xguess));
1029: /* get the communicator and the rank of the processor */
1030: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1031: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1033: if (param->output_to_file) { /* send output to binary file */
1034: PetscCall(VecCreate(comm, &pars));
1035: if (rank == 0) { /* on processor 0 */
1036: PetscCall(VecSetSizes(pars, 20, PETSC_DETERMINE));
1037: PetscCall(VecSetFromOptions(pars));
1038: PetscCall(VecSetValue(pars, 0, (PetscScalar)grid->ni, INSERT_VALUES));
1039: PetscCall(VecSetValue(pars, 1, (PetscScalar)grid->nj, INSERT_VALUES));
1040: PetscCall(VecSetValue(pars, 2, (PetscScalar)grid->dx, INSERT_VALUES));
1041: PetscCall(VecSetValue(pars, 3, (PetscScalar)grid->dz, INSERT_VALUES));
1042: PetscCall(VecSetValue(pars, 4, (PetscScalar)param->L, INSERT_VALUES));
1043: PetscCall(VecSetValue(pars, 5, (PetscScalar)param->V, INSERT_VALUES));
1044: /* skipped 6 intentionally */
1045: PetscCall(VecSetValue(pars, 7, (PetscScalar)param->slab_dip, INSERT_VALUES));
1046: PetscCall(VecSetValue(pars, 8, (PetscScalar)grid->jlid, INSERT_VALUES));
1047: PetscCall(VecSetValue(pars, 9, (PetscScalar)param->lid_depth, INSERT_VALUES));
1048: PetscCall(VecSetValue(pars, 10, (PetscScalar)grid->jfault, INSERT_VALUES));
1049: PetscCall(VecSetValue(pars, 11, (PetscScalar)param->fault_depth, INSERT_VALUES));
1050: PetscCall(VecSetValue(pars, 12, (PetscScalar)param->potentialT, INSERT_VALUES));
1051: PetscCall(VecSetValue(pars, 13, (PetscScalar)param->ivisc, INSERT_VALUES));
1052: PetscCall(VecSetValue(pars, 14, (PetscScalar)param->visc_cutoff, INSERT_VALUES));
1053: PetscCall(VecSetValue(pars, 15, (PetscScalar)param->ibound, INSERT_VALUES));
1054: PetscCall(VecSetValue(pars, 16, (PetscScalar)(its), INSERT_VALUES));
1055: } else { /* on some other processor */
1056: PetscCall(VecSetSizes(pars, 0, PETSC_DETERMINE));
1057: PetscCall(VecSetFromOptions(pars));
1058: }
1059: PetscCall(VecAssemblyBegin(pars));
1060: PetscCall(VecAssemblyEnd(pars));
1062: /* create viewer */
1063: #if defined(PETSC_HAVE_MATLAB)
1064: PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
1065: #else
1066: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
1067: #endif
1069: /* send vectors to viewer */
1070: PetscCall(PetscObjectSetName((PetscObject)res, "res"));
1071: PetscCall(VecView(res, viewer));
1072: PetscCall(PetscObjectSetName((PetscObject)user->x, "out"));
1073: PetscCall(VecView(user->x, viewer));
1074: PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "aux"));
1075: PetscCall(VecView(user->Xguess, viewer));
1076: PetscCall(StressField(da)); /* compute stress fields */
1077: PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "str"));
1078: PetscCall(VecView(user->Xguess, viewer));
1079: PetscCall(PetscObjectSetName((PetscObject)pars, "par"));
1080: PetscCall(VecView(pars, viewer));
1082: /* destroy viewer and vector */
1083: PetscCall(PetscViewerDestroy(&viewer));
1084: PetscCall(VecDestroy(&pars));
1085: }
1087: param->ivisc = ivt;
1088: PetscFunctionReturn(PETSC_SUCCESS);
1089: }
1091: /* ------------------------------------------------------------------- */
1092: /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
1093: PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1094: /* ------------------------------------------------------------------- */
1095: {
1096: AppCtx *user;
1097: Parameter *param;
1098: GridInfo *grid;
1099: Vec localX;
1100: Field **v, **x;
1101: PetscReal eps, /* dx,*/ dz, T, epsC, TC;
1102: PetscInt i, j, is, js, im, jm, ilim, jlim, ivt;
1104: PetscFunctionBeginUser;
1105: PetscCall(DMGetApplicationContext(da, &user));
1106: param = user->param;
1107: grid = user->grid;
1108: ivt = param->ivisc;
1109: param->ivisc = param->output_ivisc;
1111: PetscCall(DMGetLocalVector(da, &localX));
1112: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
1113: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
1114: PetscCall(DMDAVecGetArray(da, localX, (void **)&x));
1115: PetscCall(DMDAVecGetArray(da, V, (void **)&v));
1117: /* Parameters */
1118: /* dx = grid->dx; */ dz = grid->dz;
1120: ilim = grid->ni - 1;
1121: jlim = grid->nj - 1;
1123: /* Compute real temperature, strain rate and viscosity */
1124: PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
1125: for (j = js; j < js + jm; j++) {
1126: for (i = is; i < is + im; i++) {
1127: T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale));
1128: if (i < ilim && j < jlim) {
1129: TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale));
1130: } else {
1131: TC = T;
1132: }
1133: eps = PetscRealPart(CalcSecInv(x, i, j, CELL_CENTER, user));
1134: epsC = PetscRealPart(CalcSecInv(x, i, j, CELL_CORNER, user));
1136: v[j][i].u = eps;
1137: v[j][i].w = epsC;
1138: v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param);
1139: v[j][i].T = Viscosity(TC, epsC, dz * j, param);
1140: }
1141: }
1142: PetscCall(DMDAVecRestoreArray(da, V, (void **)&v));
1143: PetscCall(DMDAVecRestoreArray(da, localX, (void **)&x));
1144: PetscCall(DMRestoreLocalVector(da, &localX));
1146: param->ivisc = ivt;
1147: PetscFunctionReturn(PETSC_SUCCESS);
1148: }
1150: /* ------------------------------------------------------------------- */
1151: /* post-processing: compute stress everywhere */
1152: PetscErrorCode StressField(DM da)
1153: /* ------------------------------------------------------------------- */
1154: {
1155: AppCtx *user;
1156: PetscInt i, j, is, js, im, jm;
1157: Vec locVec;
1158: Field **x, **y;
1160: PetscFunctionBeginUser;
1161: PetscCall(DMGetApplicationContext(da, &user));
1163: /* Get the fine grid of Xguess and X */
1164: PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
1165: PetscCall(DMDAVecGetArray(da, user->Xguess, (void **)&x));
1167: PetscCall(DMGetLocalVector(da, &locVec));
1168: PetscCall(DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec));
1169: PetscCall(DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec));
1170: PetscCall(DMDAVecGetArray(da, locVec, (void **)&y));
1172: /* Compute stress on the corner points */
1173: for (j = js; j < js + jm; j++) {
1174: for (i = is; i < is + im; i++) {
1175: x[j][i].u = ShearStress(y, i, j, CELL_CENTER, user);
1176: x[j][i].w = ShearStress(y, i, j, CELL_CORNER, user);
1177: x[j][i].p = XNormalStress(y, i, j, CELL_CENTER, user);
1178: x[j][i].T = ZNormalStress(y, i, j, CELL_CENTER, user);
1179: }
1180: }
1182: /* Restore the fine grid of Xguess and X */
1183: PetscCall(DMDAVecRestoreArray(da, user->Xguess, (void **)&x));
1184: PetscCall(DMDAVecRestoreArray(da, locVec, (void **)&y));
1185: PetscCall(DMRestoreLocalVector(da, &locVec));
1186: PetscFunctionReturn(PETSC_SUCCESS);
1187: }
1189: /*=====================================================================
1190: UTILITY FUNCTIONS
1191: =====================================================================*/
1193: /* returns the velocity of the subducting slab and handles fault nodes for BC */
1194: static inline PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1195: {
1196: Parameter *param = user->param;
1197: GridInfo *grid = user->grid;
1199: if (c == 'U' || c == 'u') {
1200: if (i < j - 1) return param->cb;
1201: else if (j <= grid->jfault) return 0.0;
1202: else return param->cb;
1204: } else {
1205: if (i < j) return param->sb;
1206: else if (j <= grid->jfault) return 0.0;
1207: else return param->sb;
1208: }
1209: }
1211: /* solution to diffusive half-space cooling model for BC */
1212: static inline PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1213: {
1214: Parameter *param = user->param;
1215: PetscScalar z;
1216: if (plate == PLATE_LID) z = (j - 0.5) * user->grid->dz;
1217: else z = (j - 0.5) * user->grid->dz * param->cb; /* PLATE_SLAB */
1218: #if defined(PETSC_HAVE_ERF)
1219: return (PetscReal)(erf((double)PetscRealPart(z * param->L / 2.0 / param->skt)));
1220: #else
1221: (*PetscErrorPrintf)("erf() not available on this machine\n");
1222: MPI_Abort(PETSC_COMM_SELF, 1);
1223: #endif
1224: }
1226: /*=====================================================================
1227: INTERACTIVE SIGNAL HANDLING
1228: =====================================================================*/
1230: /* ------------------------------------------------------------------- */
1231: PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it, PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, void *ctx)
1232: /* ------------------------------------------------------------------- */
1233: {
1234: AppCtx *user = (AppCtx *)ctx;
1235: Parameter *param = user->param;
1236: KSP ksp;
1238: PetscFunctionBeginUser;
1239: if (param->interrupted) {
1240: param->interrupted = PETSC_FALSE;
1241: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: exiting SNES solve. \n"));
1242: *reason = SNES_CONVERGED_FNORM_ABS;
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: } else if (param->toggle_kspmon) {
1245: param->toggle_kspmon = PETSC_FALSE;
1247: PetscCall(SNESGetKSP(snes, &ksp));
1249: if (param->kspmon) {
1250: PetscCall(KSPMonitorCancel(ksp));
1252: param->kspmon = PETSC_FALSE;
1253: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: deactivating ksp singular value monitor. \n"));
1254: } else {
1255: PetscViewerAndFormat *vf;
1256: PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
1257: PetscCall(KSPMonitorSet(ksp, (PetscErrorCode(*)(KSP, PetscInt, PetscReal, void *))KSPMonitorSingularValue, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy));
1259: param->kspmon = PETSC_TRUE;
1260: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: activating ksp singular value monitor. \n"));
1261: }
1262: }
1263: PetscCall(SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason, ctx));
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: /* ------------------------------------------------------------------- */
1268: #include <signal.h>
1269: PetscErrorCode InteractiveHandler(int signum, void *ctx)
1270: /* ------------------------------------------------------------------- */
1271: {
1272: AppCtx *user = (AppCtx *)ctx;
1273: Parameter *param = user->param;
1275: if (signum == SIGILL) {
1276: param->toggle_kspmon = PETSC_TRUE;
1277: #if !defined(PETSC_MISSING_SIGCONT)
1278: } else if (signum == SIGCONT) {
1279: param->interrupted = PETSC_TRUE;
1280: #endif
1281: #if !defined(PETSC_MISSING_SIGURG)
1282: } else if (signum == SIGURG) {
1283: param->stop_solve = PETSC_TRUE;
1284: #endif
1285: }
1286: return PETSC_SUCCESS;
1287: }
1289: /* main call-back function that computes the processor-local piece of the residual */
1290: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr)
1291: {
1292: AppCtx *user = (AppCtx *)ptr;
1293: Parameter *param = user->param;
1294: GridInfo *grid = user->grid;
1295: PetscScalar mag_w, mag_u;
1296: PetscInt i, j, mx, mz, ilim, jlim;
1297: PetscInt is, ie, js, je, ibound; /* ,ivisc */
1299: PetscFunctionBeginUser;
1300: /* Define global and local grid parameters */
1301: mx = info->mx;
1302: mz = info->my;
1303: ilim = mx - 1;
1304: jlim = mz - 1;
1305: is = info->xs;
1306: ie = info->xs + info->xm;
1307: js = info->ys;
1308: je = info->ys + info->ym;
1310: /* Define geometric and numeric parameters */
1311: /* ivisc = param->ivisc; */ ibound = param->ibound;
1313: for (j = js; j < je; j++) {
1314: for (i = is; i < ie; i++) {
1315: /************* X-MOMENTUM/VELOCITY *************/
1316: if (i < j) f[j][i].u = x[j][i].u - SlabVel('U', i, j, user);
1317: else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1318: /* in the lithospheric lid */
1319: f[j][i].u = x[j][i].u - 0.0;
1320: } else if (i == ilim) {
1321: /* on the right side boundary */
1322: if (ibound == BC_ANALYTIC) {
1323: f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
1324: } else {
1325: f[j][i].u = XNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
1326: }
1328: } else if (j == jlim) {
1329: /* on the bottom boundary */
1330: if (ibound == BC_ANALYTIC) {
1331: f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
1332: } else if (ibound == BC_NOSTRESS) {
1333: f[j][i].u = XMomentumResidual(x, i, j, user);
1334: } else {
1335: /* experimental boundary condition */
1336: }
1338: } else {
1339: /* in the mantle wedge */
1340: f[j][i].u = XMomentumResidual(x, i, j, user);
1341: }
1343: /************* Z-MOMENTUM/VELOCITY *************/
1344: if (i <= j) {
1345: f[j][i].w = x[j][i].w - SlabVel('W', i, j, user);
1347: } else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1348: /* in the lithospheric lid */
1349: f[j][i].w = x[j][i].w - 0.0;
1351: } else if (j == jlim) {
1352: /* on the bottom boundary */
1353: if (ibound == BC_ANALYTIC) {
1354: f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
1355: } else {
1356: f[j][i].w = ZNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
1357: }
1359: } else if (i == ilim) {
1360: /* on the right side boundary */
1361: if (ibound == BC_ANALYTIC) {
1362: f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
1363: } else if (ibound == BC_NOSTRESS) {
1364: f[j][i].w = ZMomentumResidual(x, i, j, user);
1365: } else {
1366: /* experimental boundary condition */
1367: }
1369: } else {
1370: /* in the mantle wedge */
1371: f[j][i].w = ZMomentumResidual(x, i, j, user);
1372: }
1374: /************* CONTINUITY/PRESSURE *************/
1375: if (i < j || j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1376: /* in the lid or slab */
1377: f[j][i].p = x[j][i].p;
1379: } else if ((i == ilim || j == jlim) && ibound == BC_ANALYTIC) {
1380: /* on an analytic boundary */
1381: f[j][i].p = x[j][i].p - Pressure(i, j, user);
1383: } else {
1384: /* in the mantle wedge */
1385: f[j][i].p = ContinuityResidual(x, i, j, user);
1386: }
1388: /************* TEMPERATURE *************/
1389: if (j == 0) {
1390: /* on the surface */
1391: f[j][i].T = x[j][i].T + x[j + 1][i].T + PetscMax(PetscRealPart(x[j][i].T), 0.0);
1393: } else if (i == 0) {
1394: /* slab inflow boundary */
1395: f[j][i].T = x[j][i].T - PlateModel(j, PLATE_SLAB, user);
1397: } else if (i == ilim) {
1398: /* right side boundary */
1399: mag_u = 1.0 - PetscPowRealInt((1.0 - PetscMax(PetscMin(PetscRealPart(x[j][i - 1].u) / param->cb, 1.0), 0.0)), 5);
1400: 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);
1402: } else if (j == jlim) {
1403: /* bottom boundary */
1404: mag_w = 1.0 - PetscPowRealInt((1.0 - PetscMax(PetscMin(PetscRealPart(x[j - 1][i].w) / param->sb, 1.0), 0.0)), 5);
1405: f[j][i].T = x[j][i].T - mag_w * x[j - 1][i - 1].T - (1.0 - mag_w);
1407: } else {
1408: /* in the mantle wedge */
1409: f[j][i].T = EnergyResidual(x, i, j, user);
1410: }
1411: }
1412: }
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: /*TEST
1418: build:
1419: requires: !complex erf
1421: test:
1422: args: -ni 18 -fp_trap 0
1423: filter: grep -v Destination
1424: requires: !single
1426: TEST*/