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";
26
27 /*F-----------------------------------------------------------------------
28
29 This PETSc 2.2.0 example by Richard F. Katz
30 http://www.ldeo.columbia.edu/~katz/
31
32 The problem is modeled by the partial differential equation system
33
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}
39
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}
46
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}$-----
51
52 ------------------------------------------------------------------------F*/
53
54 #include <petscsnes.h>
55 #include <petscdm.h>
56 #include <petscdmda.h>
57
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
72
73 typedef struct { /* holds the variables to be solved for */
74 PetscScalar u, w, p, T;
75 } Field;
76
77 typedef struct { /* parameters needed to compute viscosity */
78 PetscReal A, n, Estar, Vstar;
79 } ViscParam;
80
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;
92
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;
100
101 typedef struct { /* application context */
102 Vec x, Xguess;
103 Parameter *param;
104 GridInfo *grid;
105 } AppCtx;
106
107 /* Callback functions (static interface) */
108 extern PetscErrorCode FormFunctionLocal(DMDALocalInfo *, Field **, Field **, void *);
109
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);
116
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
main(int argc,char ** argv)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;
132
133 PetscFunctionBeginUser;
134 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
141
142 comm = PETSC_COMM_WORLD;
143
144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145 Set up the problem parameters.
146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147 PetscCall(SetParams(¶m, &grid));
148 PetscCall(ReportParams(¶m, &grid));
149
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"));
161
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));
170
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));
176
177 PetscCall(SNESSetConvergenceTest(snes, SNESConverged_Interactive, (void *)user, NULL));
178 PetscCall(PetscPushSignalHandler(InteractiveHandler, (void *)user));
179
180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181 Initialize and solve the nonlinear system
182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183 PetscCall(Initialize(da));
184 PetscCall(UpdateSolution(snes, user, &nits));
185
186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187 Output variables.
188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189 PetscCall(DoOutput(snes, nits));
190
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 }
203
204 /*=====================================================================
205 PETSc INTERACTION FUNCTIONS (initialize & call SNESSolve)
206 =====================================================================*/
207
208 /* manages solve: adaptive continuation method */
UpdateSolution(SNES snes,AppCtx * user,PetscInt * nits)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;
219
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));
226
227 *nits = 0;
228
229 /* Isoviscous solve */
230 if (param->ivisc == VISC_CONST && !param->stop_solve) {
231 param->ivisc = VISC_CONST;
232
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 }
240
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"));
244
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));
248
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;
257
258 if (reason < 0) {
259 /* NOT converged */
260 cont_incr = -PetscAbsReal(cont_incr) / 2.0;
261 if (PetscAbsReal(cont_incr) < 0.01) goto done;
262
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;
270
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 }
280
281 /*=====================================================================
282 PHYSICS FUNCTIONS (compute the discrete residual)
283 =====================================================================*/
284
UInterp(Field ** x,PetscInt i,PetscInt j)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 }
289
WInterp(Field ** x,PetscInt i,PetscInt j)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 }
294
PInterp(Field ** x,PetscInt i,PetscInt j)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 }
299
TInterp(Field ** x,PetscInt i,PetscInt j)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 }
304
305 /* isoviscous analytic solution for IC */
HorizVelocity(PetscInt i,PetscInt j,AppCtx * user)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;
312
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 }
321
322 /* isoviscous analytic solution for IC */
VertVelocity(PetscInt i,PetscInt j,AppCtx * user)323 static inline PetscScalar VertVelocity(PetscInt i, PetscInt j, AppCtx *user)
324 {
325 Parameter *param = user->param;
326 GridInfo *grid = user->grid;
327 PetscScalar st, ct, th, c = param->c, d = param->d;
328 PetscReal x, z, r;
329
330 x = (i - grid->jlid - 0.5) * grid->dx;
331 z = (j - grid->jlid) * grid->dz;
332 r = PetscSqrtReal(x * x + z * z);
333 st = z / r;
334 ct = x / r;
335 th = PetscAtanReal(z / x);
336 return st * (c * th * st + d * (st + th * ct)) - ct * (c * (st - th * ct) + d * th * st);
337 }
338
339 /* isoviscous analytic solution for IC */
Pressure(PetscInt i,PetscInt j,AppCtx * user)340 static inline PetscScalar Pressure(PetscInt i, PetscInt j, AppCtx *user)
341 {
342 Parameter *param = user->param;
343 GridInfo *grid = user->grid;
344 PetscScalar x, z, r, st, ct, c = param->c, d = param->d;
345
346 x = (i - grid->jlid - 0.5) * grid->dx;
347 z = (j - grid->jlid - 0.5) * grid->dz;
348 r = PetscSqrtReal(x * x + z * z);
349 st = z / r;
350 ct = x / r;
351 return -2.0 * (c * ct - d * st) / r;
352 }
353
354 /* computes the second invariant of the strain rate tensor */
CalcSecInv(Field ** x,PetscInt i,PetscInt j,PetscInt ipos,AppCtx * user)355 static inline PetscScalar CalcSecInv(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
356 {
357 Parameter *param = user->param;
358 GridInfo *grid = user->grid;
359 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1;
360 PetscScalar uN, uS, uE, uW, wN, wS, wE, wW;
361 PetscScalar eps11, eps12, eps22;
362
363 if (i < j) return EPS_ZERO;
364 if (i == ilim) i--;
365 if (j == jlim) j--;
366
367 if (ipos == CELL_CENTER) { /* on cell center */
368 if (j <= grid->jlid) return EPS_ZERO;
369
370 uE = x[j][i].u;
371 uW = x[j][i - 1].u;
372 wN = x[j][i].w;
373 wS = x[j - 1][i].w;
374 wE = WInterp(x, i, j - 1);
375 if (i == j) {
376 uN = param->cb;
377 wW = param->sb;
378 } else {
379 uN = UInterp(x, i - 1, j);
380 wW = WInterp(x, i - 1, j - 1);
381 }
382
383 if (j == grid->jlid + 1) uS = 0.0;
384 else uS = UInterp(x, i - 1, j - 1);
385
386 } else { /* on CELL_CORNER */
387 if (j < grid->jlid) return EPS_ZERO;
388
389 uN = x[j + 1][i].u;
390 uS = x[j][i].u;
391 wE = x[j][i + 1].w;
392 wW = x[j][i].w;
393 if (i == j) {
394 wN = param->sb;
395 uW = param->cb;
396 } else {
397 wN = WInterp(x, i, j);
398 uW = UInterp(x, i - 1, j);
399 }
400
401 if (j == grid->jlid) {
402 uE = 0.0;
403 uW = 0.0;
404 uS = -uN;
405 wS = -wN;
406 } else {
407 uE = UInterp(x, i, j);
408 wS = WInterp(x, i, j - 1);
409 }
410 }
411
412 eps11 = (uE - uW) / grid->dx;
413 eps22 = (wN - wS) / grid->dz;
414 eps12 = 0.5 * ((uN - uS) / grid->dz + (wE - wW) / grid->dx);
415
416 return PetscSqrtReal(0.5 * (eps11 * eps11 + 2.0 * eps12 * eps12 + eps22 * eps22));
417 }
418
419 /* computes the shear viscosity */
Viscosity(PetscScalar T,PetscScalar eps,PetscScalar z,Parameter * param)420 static inline PetscScalar Viscosity(PetscScalar T, PetscScalar eps, PetscScalar z, Parameter *param)
421 {
422 PetscReal result = 0.0;
423 ViscParam difn = param->diffusion, disl = param->dislocation;
424 PetscInt iVisc = param->ivisc;
425 PetscScalar eps_scale = param->V / (param->L * 1000.0);
426 PetscScalar strain_power, v1, v2, P;
427 PetscScalar rho_g = 32340.0, R = 8.3144;
428
429 P = rho_g * (z * param->L * 1000.0); /* Pa */
430
431 if (iVisc == VISC_CONST) {
432 /* constant viscosity */
433 return 1.0;
434 } else if (iVisc == VISC_DIFN) {
435 /* diffusion creep rheology */
436 result = PetscRealPart(difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0);
437 } else if (iVisc == VISC_DISL) {
438 /* dislocation creep rheology */
439 strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);
440
441 result = PetscRealPart(disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0);
442 } else if (iVisc == VISC_FULL) {
443 /* dislocation/diffusion creep rheology */
444 strain_power = PetscPowScalar(eps * eps_scale, (1.0 - disl.n) / disl.n);
445
446 v1 = difn.A * PetscExpScalar((difn.Estar + P * difn.Vstar) / R / (T + 273.0)) / param->eta0;
447 v2 = disl.A * PetscExpScalar((disl.Estar + P * disl.Vstar) / disl.n / R / (T + 273.0)) * strain_power / param->eta0;
448
449 result = PetscRealPart(1.0 / (1.0 / v1 + 1.0 / v2));
450 }
451
452 /* max viscosity is param->eta0 */
453 result = PetscMin(result, 1.0);
454 /* min viscosity is param->visc_cutoff */
455 result = PetscMax(result, param->visc_cutoff);
456 /* continuation method */
457 result = PetscPowReal(result, param->continuation);
458 return result;
459 }
460
461 /* computes the residual of the x-component of eqn (1) above */
XMomentumResidual(Field ** x,PetscInt i,PetscInt j,AppCtx * user)462 static inline PetscScalar XMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
463 {
464 Parameter *param = user->param;
465 GridInfo *grid = user->grid;
466 PetscScalar dx = grid->dx, dz = grid->dz;
467 PetscScalar etaN, etaS, etaE, etaW, epsN = 0.0, epsS = 0.0, epsE = 0.0, epsW = 0.0;
468 PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdx, residual, z_scale;
469 PetscScalar dudxW, dudxE, dudzN, dudzS, dwdxN, dwdxS;
470 PetscInt jlim = grid->nj - 1;
471
472 z_scale = param->z_scale;
473
474 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
475 TS = param->potentialT * TInterp(x, i, j - 1) * PetscExpScalar((j - 1.0) * dz * z_scale);
476 if (j == jlim) TN = TS;
477 else TN = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
478 TW = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
479 TE = param->potentialT * x[j][i + 1].T * PetscExpScalar((j - 0.5) * dz * z_scale);
480 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
481 epsN = CalcSecInv(x, i, j, CELL_CORNER, user);
482 epsS = CalcSecInv(x, i, j - 1, CELL_CORNER, user);
483 epsE = CalcSecInv(x, i + 1, j, CELL_CENTER, user);
484 epsW = CalcSecInv(x, i, j, CELL_CENTER, user);
485 }
486 }
487 etaN = Viscosity(TN, epsN, dz * (j + 0.5), param);
488 etaS = Viscosity(TS, epsS, dz * (j - 0.5), param);
489 etaW = Viscosity(TW, epsW, dz * j, param);
490 etaE = Viscosity(TE, epsE, dz * j, param);
491
492 dPdx = (x[j][i + 1].p - x[j][i].p) / dx;
493 if (j == jlim) dudzN = etaN * (x[j][i].w - x[j][i + 1].w) / dx;
494 else dudzN = etaN * (x[j + 1][i].u - x[j][i].u) / dz;
495 dudzS = etaS * (x[j][i].u - x[j - 1][i].u) / dz;
496 dudxE = etaE * (x[j][i + 1].u - x[j][i].u) / dx;
497 dudxW = etaW * (x[j][i].u - x[j][i - 1].u) / dx;
498
499 residual = -dPdx /* X-MOMENTUM EQUATION*/
500 + (dudxE - dudxW) / dx + (dudzN - dudzS) / dz;
501
502 if (param->ivisc != VISC_CONST) {
503 dwdxN = etaN * (x[j][i + 1].w - x[j][i].w) / dx;
504 dwdxS = etaS * (x[j - 1][i + 1].w - x[j - 1][i].w) / dx;
505
506 residual += (dudxE - dudxW) / dx + (dwdxN - dwdxS) / dz;
507 }
508
509 return residual;
510 }
511
512 /* computes the residual of the z-component of eqn (1) above */
ZMomentumResidual(Field ** x,PetscInt i,PetscInt j,AppCtx * user)513 static inline PetscScalar ZMomentumResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
514 {
515 Parameter *param = user->param;
516 GridInfo *grid = user->grid;
517 PetscScalar dx = grid->dx, dz = grid->dz;
518 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;
519 PetscScalar TE = 0.0, TN = 0.0, TS = 0.0, TW = 0.0, dPdz, residual, z_scale;
520 PetscScalar dudzE, dudzW, dwdxW, dwdxE, dwdzN, dwdzS;
521 PetscInt ilim = grid->ni - 1;
522
523 /* geometric and other parameters */
524 z_scale = param->z_scale;
525
526 /* viscosity */
527 if (param->ivisc == VISC_DIFN || param->ivisc >= VISC_DISL) { /* viscosity is T-dependent */
528 TN = param->potentialT * x[j + 1][i].T * PetscExpScalar((j + 0.5) * dz * z_scale);
529 TS = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
530 TW = param->potentialT * TInterp(x, i - 1, j) * PetscExpScalar(j * dz * z_scale);
531 if (i == ilim) TE = TW;
532 else TE = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
533 if (param->ivisc >= VISC_DISL) { /* olivine dislocation creep */
534 epsN = CalcSecInv(x, i, j + 1, CELL_CENTER, user);
535 epsS = CalcSecInv(x, i, j, CELL_CENTER, user);
536 epsE = CalcSecInv(x, i, j, CELL_CORNER, user);
537 epsW = CalcSecInv(x, i - 1, j, CELL_CORNER, user);
538 }
539 }
540 etaN = Viscosity(TN, epsN, dz * (j + 1.0), param);
541 etaS = Viscosity(TS, epsS, dz * (j + 0.0), param);
542 etaW = Viscosity(TW, epsW, dz * (j + 0.5), param);
543 etaE = Viscosity(TE, epsE, dz * (j + 0.5), param);
544
545 dPdz = (x[j + 1][i].p - x[j][i].p) / dz;
546 dwdzN = etaN * (x[j + 1][i].w - x[j][i].w) / dz;
547 dwdzS = etaS * (x[j][i].w - x[j - 1][i].w) / dz;
548 if (i == ilim) dwdxE = etaE * (x[j][i].u - x[j + 1][i].u) / dz;
549 else dwdxE = etaE * (x[j][i + 1].w - x[j][i].w) / dx;
550 dwdxW = 2.0 * etaW * (x[j][i].w - x[j][i - 1].w) / dx;
551
552 /* Z-MOMENTUM */
553 residual = -dPdz /* constant viscosity terms */
554 + (dwdzN - dwdzS) / dz + (dwdxE - dwdxW) / dx;
555
556 if (param->ivisc != VISC_CONST) {
557 dudzE = etaE * (x[j + 1][i].u - x[j][i].u) / dz;
558 dudzW = etaW * (x[j + 1][i - 1].u - x[j][i - 1].u) / dz;
559
560 residual += (dwdzN - dwdzS) / dz + (dudzE - dudzW) / dx;
561 }
562
563 return residual;
564 }
565
566 /* computes the residual of eqn (2) above */
ContinuityResidual(Field ** x,PetscInt i,PetscInt j,AppCtx * user)567 static inline PetscScalar ContinuityResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
568 {
569 GridInfo *grid = user->grid;
570 PetscScalar uE, uW, wN, wS, dudx, dwdz;
571
572 uW = x[j][i - 1].u;
573 uE = x[j][i].u;
574 dudx = (uE - uW) / grid->dx;
575 wS = x[j - 1][i].w;
576 wN = x[j][i].w;
577 dwdz = (wN - wS) / grid->dz;
578
579 return dudx + dwdz;
580 }
581
582 /* computes the residual of eqn (3) above */
EnergyResidual(Field ** x,PetscInt i,PetscInt j,AppCtx * user)583 static inline PetscScalar EnergyResidual(Field **x, PetscInt i, PetscInt j, AppCtx *user)
584 {
585 Parameter *param = user->param;
586 GridInfo *grid = user->grid;
587 PetscScalar dx = grid->dx, dz = grid->dz;
588 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, jlid = grid->jlid;
589 PetscScalar TE, TN, TS, TW, residual;
590 PetscScalar uE, uW, wN, wS;
591 PetscScalar fN, fS, fE, fW, dTdxW, dTdxE, dTdzN, dTdzS;
592
593 dTdzN = (x[j + 1][i].T - x[j][i].T) / dz;
594 dTdzS = (x[j][i].T - x[j - 1][i].T) / dz;
595 dTdxE = (x[j][i + 1].T - x[j][i].T) / dx;
596 dTdxW = (x[j][i].T - x[j][i - 1].T) / dx;
597
598 residual = ((dTdzN - dTdzS) / dz + /* diffusion term */
599 (dTdxE - dTdxW) / dx) *
600 dx * dz / param->peclet;
601
602 if (j <= jlid && i >= j) {
603 /* don't advect in the lid */
604 return residual;
605 } else if (i < j) {
606 /* beneath the slab sfc */
607 uW = uE = param->cb;
608 wS = wN = param->sb;
609 } else {
610 /* advect in the slab and wedge */
611 uW = x[j][i - 1].u;
612 uE = x[j][i].u;
613 wS = x[j - 1][i].w;
614 wN = x[j][i].w;
615 }
616
617 if (param->adv_scheme == ADVECT_FV || i == ilim - 1 || j == jlim - 1 || i == 1 || j == 1) {
618 /* finite volume advection */
619 TS = (x[j][i].T + x[j - 1][i].T) / 2.0;
620 TN = (x[j][i].T + x[j + 1][i].T) / 2.0;
621 TE = (x[j][i].T + x[j][i + 1].T) / 2.0;
622 TW = (x[j][i].T + x[j][i - 1].T) / 2.0;
623 fN = wN * TN * dx;
624 fS = wS * TS * dx;
625 fE = uE * TE * dz;
626 fW = uW * TW * dz;
627
628 } else {
629 /* Fromm advection scheme */
630 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;
631 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;
632 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;
633 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;
634 }
635
636 residual -= (fE - fW + fN - fS);
637
638 return residual;
639 }
640
641 /* computes the shear stress---used on the boundaries */
ShearStress(Field ** x,PetscInt i,PetscInt j,PetscInt ipos,AppCtx * user)642 static inline PetscScalar ShearStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
643 {
644 Parameter *param = user->param;
645 GridInfo *grid = user->grid;
646 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1;
647 PetscScalar uN, uS, wE, wW;
648
649 if (j <= grid->jlid || i < j || i == ilim || j == jlim) return EPS_ZERO;
650
651 if (ipos == CELL_CENTER) { /* on cell center */
652
653 wE = WInterp(x, i, j - 1);
654 if (i == j) {
655 wW = param->sb;
656 uN = param->cb;
657 } else {
658 wW = WInterp(x, i - 1, j - 1);
659 uN = UInterp(x, i - 1, j);
660 }
661 if (j == grid->jlid + 1) uS = 0.0;
662 else uS = UInterp(x, i - 1, j - 1);
663
664 } else { /* on cell corner */
665
666 uN = x[j + 1][i].u;
667 uS = x[j][i].u;
668 wW = x[j][i].w;
669 wE = x[j][i + 1].w;
670 }
671
672 return (uN - uS) / grid->dz + (wE - wW) / grid->dx;
673 }
674
675 /* computes the normal stress---used on the boundaries */
XNormalStress(Field ** x,PetscInt i,PetscInt j,PetscInt ipos,AppCtx * user)676 static inline PetscScalar XNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
677 {
678 Parameter *param = user->param;
679 GridInfo *grid = user->grid;
680 PetscScalar dx = grid->dx, dz = grid->dz;
681 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
682 PetscScalar epsC = 0.0, etaC, TC, uE, uW, pC, z_scale;
683 if (i < j || j <= grid->jlid) return EPS_ZERO;
684
685 ivisc = param->ivisc;
686 z_scale = param->z_scale;
687
688 if (ipos == CELL_CENTER) { /* on cell center */
689
690 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
691 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
692 etaC = Viscosity(TC, epsC, dz * j, param);
693
694 uW = x[j][i - 1].u;
695 uE = x[j][i].u;
696 pC = x[j][i].p;
697
698 } else { /* on cell corner */
699 if (i == ilim || j == jlim) return EPS_ZERO;
700
701 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
702 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
703 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);
704
705 if (i == j) uW = param->sb;
706 else uW = UInterp(x, i - 1, j);
707 uE = UInterp(x, i, j);
708 pC = PInterp(x, i, j);
709 }
710
711 return 2.0 * etaC * (uE - uW) / dx - pC;
712 }
713
714 /* computes the normal stress---used on the boundaries */
ZNormalStress(Field ** x,PetscInt i,PetscInt j,PetscInt ipos,AppCtx * user)715 static inline PetscScalar ZNormalStress(Field **x, PetscInt i, PetscInt j, PetscInt ipos, AppCtx *user)
716 {
717 Parameter *param = user->param;
718 GridInfo *grid = user->grid;
719 PetscScalar dz = grid->dz;
720 PetscInt ilim = grid->ni - 1, jlim = grid->nj - 1, ivisc;
721 PetscScalar epsC = 0.0, etaC, TC;
722 PetscScalar pC, wN, wS, z_scale;
723 if (i < j || j <= grid->jlid) return EPS_ZERO;
724
725 ivisc = param->ivisc;
726 z_scale = param->z_scale;
727
728 if (ipos == CELL_CENTER) { /* on cell center */
729
730 TC = param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * z_scale);
731 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CENTER, user);
732 etaC = Viscosity(TC, epsC, dz * j, param);
733 wN = x[j][i].w;
734 wS = x[j - 1][i].w;
735 pC = x[j][i].p;
736
737 } else { /* on cell corner */
738 if ((i == ilim) || (j == jlim)) return EPS_ZERO;
739
740 TC = param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * z_scale);
741 if (ivisc >= VISC_DISL) epsC = CalcSecInv(x, i, j, CELL_CORNER, user);
742 etaC = Viscosity(TC, epsC, dz * (j + 0.5), param);
743 if (i == j) wN = param->sb;
744 else wN = WInterp(x, i, j);
745 wS = WInterp(x, i, j - 1);
746 pC = PInterp(x, i, j);
747 }
748
749 return 2.0 * etaC * (wN - wS) / dz - pC;
750 }
751
752 /*=====================================================================
753 INITIALIZATION, POST-PROCESSING AND OUTPUT FUNCTIONS
754 =====================================================================*/
755
756 /* initializes the problem parameters and checks for
757 command line changes */
SetParams(Parameter * param,GridInfo * grid)758 PetscErrorCode SetParams(Parameter *param, GridInfo *grid)
759 {
760 PetscReal SEC_PER_YR = 3600.00 * 24.00 * 365.2500;
761 PetscReal alpha_g_on_cp_units_inverse_km = 4.0e-5 * 9.8;
762
763 PetscFunctionBeginUser;
764 /* domain geometry */
765 param->slab_dip = 45.0;
766 param->width = 320.0; /* km */
767 param->depth = 300.0; /* km */
768 param->lid_depth = 35.0; /* km */
769 param->fault_depth = 35.0; /* km */
770
771 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_dip", ¶m->slab_dip, NULL));
772 PetscCall(PetscOptionsGetReal(NULL, NULL, "-width", ¶m->width, NULL));
773 PetscCall(PetscOptionsGetReal(NULL, NULL, "-depth", ¶m->depth, NULL));
774 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_depth", ¶m->lid_depth, NULL));
775 PetscCall(PetscOptionsGetReal(NULL, NULL, "-fault_depth", ¶m->fault_depth, NULL));
776
777 param->slab_dip = param->slab_dip * PETSC_PI / 180.0; /* radians */
778
779 /* grid information */
780 PetscCall(PetscOptionsGetInt(NULL, NULL, "-jfault", &grid->jfault, NULL));
781 grid->ni = 82;
782 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ni", &grid->ni, NULL));
783
784 grid->dx = param->width / ((PetscReal)(grid->ni - 2)); /* km */
785 grid->dz = grid->dx * PetscTanReal(param->slab_dip); /* km */
786 grid->nj = (PetscInt)(param->depth / grid->dz + 3.0); /* gridpoints*/
787 param->depth = grid->dz * (grid->nj - 2); /* km */
788 grid->inose = 0; /* gridpoints*/
789 PetscCall(PetscOptionsGetInt(NULL, NULL, "-inose", &grid->inose, NULL));
790 grid->bx = DM_BOUNDARY_NONE;
791 grid->by = DM_BOUNDARY_NONE;
792 grid->stencil = DMDA_STENCIL_BOX;
793 grid->dof = 4;
794 grid->stencil_width = 2;
795 grid->mglevels = 1;
796
797 /* boundary conditions */
798 param->pv_analytic = PETSC_FALSE;
799 param->ibound = BC_NOSTRESS;
800 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ibound", ¶m->ibound, NULL));
801
802 /* physical constants */
803 param->slab_velocity = 5.0; /* cm/yr */
804 param->slab_age = 50.0; /* Ma */
805 param->lid_age = 50.0; /* Ma */
806 param->kappa = 0.7272e-6; /* m^2/sec */
807 param->potentialT = 1300.0; /* degrees C */
808
809 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_velocity", ¶m->slab_velocity, NULL));
810 PetscCall(PetscOptionsGetReal(NULL, NULL, "-slab_age", ¶m->slab_age, NULL));
811 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lid_age", ¶m->lid_age, NULL));
812 PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", ¶m->kappa, NULL));
813 PetscCall(PetscOptionsGetReal(NULL, NULL, "-potentialT", ¶m->potentialT, NULL));
814
815 /* viscosity */
816 param->ivisc = 3; /* 0=isovisc, 1=difn creep, 2=disl creep, 3=full */
817 param->eta0 = 1e24; /* Pa-s */
818 param->visc_cutoff = 0.0; /* factor of eta_0 */
819 param->continuation = 1.0;
820
821 /* constants for diffusion creep */
822 param->diffusion.A = 1.8e7; /* Pa-s */
823 param->diffusion.n = 1.0; /* dim'less */
824 param->diffusion.Estar = 375e3; /* J/mol */
825 param->diffusion.Vstar = 5e-6; /* m^3/mol */
826
827 /* constants for param->dislocationocation creep */
828 param->dislocation.A = 2.8969e4; /* Pa-s */
829 param->dislocation.n = 3.5; /* dim'less */
830 param->dislocation.Estar = 530e3; /* J/mol */
831 param->dislocation.Vstar = 14e-6; /* m^3/mol */
832
833 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ivisc", ¶m->ivisc, NULL));
834 PetscCall(PetscOptionsGetReal(NULL, NULL, "-visc_cutoff", ¶m->visc_cutoff, NULL));
835
836 param->output_ivisc = param->ivisc;
837
838 PetscCall(PetscOptionsGetInt(NULL, NULL, "-output_ivisc", ¶m->output_ivisc, NULL));
839 PetscCall(PetscOptionsGetReal(NULL, NULL, "-vstar", ¶m->dislocation.Vstar, NULL));
840
841 /* output options */
842 param->quiet = PETSC_FALSE;
843 param->param_test = PETSC_FALSE;
844
845 PetscCall(PetscOptionsHasName(NULL, NULL, "-quiet", ¶m->quiet));
846 PetscCall(PetscOptionsHasName(NULL, NULL, "-test", ¶m->param_test));
847 PetscCall(PetscOptionsGetString(NULL, NULL, "-file", param->filename, sizeof(param->filename), ¶m->output_to_file));
848
849 /* advection */
850 param->adv_scheme = ADVECT_FROMM; /* advection scheme: 0=finite vol, 1=Fromm */
851
852 PetscCall(PetscOptionsGetInt(NULL, NULL, "-adv_scheme", ¶m->adv_scheme, NULL));
853
854 /* misc. flags */
855 param->stop_solve = PETSC_FALSE;
856 param->interrupted = PETSC_FALSE;
857 param->kspmon = PETSC_FALSE;
858 param->toggle_kspmon = PETSC_FALSE;
859
860 /* derived parameters for slab angle */
861 param->sb = PetscSinReal(param->slab_dip);
862 param->cb = PetscCosReal(param->slab_dip);
863 param->c = param->slab_dip * param->sb / (param->slab_dip * param->slab_dip - param->sb * param->sb);
864 param->d = (param->slab_dip * param->cb - param->sb) / (param->slab_dip * param->slab_dip - param->sb * param->sb);
865
866 /* length, velocity and time scale for non-dimensionalization */
867 param->L = PetscMin(param->width, param->depth); /* km */
868 param->V = param->slab_velocity / 100.0 / SEC_PER_YR; /* m/sec */
869
870 /* other unit conversions and derived parameters */
871 param->scaled_width = param->width / param->L; /* dim'less */
872 param->scaled_depth = param->depth / param->L; /* dim'less */
873 param->lid_depth = param->lid_depth / param->L; /* dim'less */
874 param->fault_depth = param->fault_depth / param->L; /* dim'less */
875 grid->dx = grid->dx / param->L; /* dim'less */
876 grid->dz = grid->dz / param->L; /* dim'less */
877 grid->jlid = (PetscInt)(param->lid_depth / grid->dz); /* gridcells */
878 grid->jfault = (PetscInt)(param->fault_depth / grid->dz); /* gridcells */
879 param->lid_depth = grid->jlid * grid->dz; /* dim'less */
880 param->fault_depth = grid->jfault * grid->dz; /* dim'less */
881 grid->corner = grid->jlid + 1; /* gridcells */
882 param->peclet = param->V /* m/sec */
883 * param->L * 1000.0 /* m */
884 / param->kappa; /* m^2/sec */
885 param->z_scale = param->L * alpha_g_on_cp_units_inverse_km;
886 param->skt = PetscSqrtReal(param->kappa * param->slab_age * SEC_PER_YR);
887 PetscCall(PetscOptionsGetReal(NULL, NULL, "-peclet", ¶m->peclet, NULL));
888 PetscFunctionReturn(PETSC_SUCCESS);
889 }
890
891 /* prints a report of the problem parameters to stdout */
ReportParams(Parameter * param,GridInfo * grid)892 PetscErrorCode ReportParams(Parameter *param, GridInfo *grid)
893 {
894 char date[30];
895
896 PetscFunctionBeginUser;
897 PetscCall(PetscGetDate(date, 30));
898
899 if (!param->quiet) {
900 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------BEGIN ex30 PARAM REPORT-------------------\n"));
901 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Domain: \n"));
902 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Width = %g km, Depth = %g km\n", (double)param->width, (double)param->depth));
903 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));
904 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)));
905
906 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nGrid: \n"));
907 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)));
908 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " jlid = %3" PetscInt_FMT " jfault = %3" PetscInt_FMT " \n", grid->jlid, grid->jfault));
909 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pe = %g\n", (double)param->peclet));
910
911 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nRheology:"));
912 if (param->ivisc == VISC_CONST) {
913 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Isoviscous \n"));
914 if (param->pv_analytic) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Pressure and Velocity prescribed! \n"));
915 } else if (param->ivisc == VISC_DIFN) {
916 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Diffusion Creep (T-Dependent Newtonian) \n"));
917 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
918 } else if (param->ivisc == VISC_DISL) {
919 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Dislocation Creep (T-Dependent Non-Newtonian) \n"));
920 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
921 } else if (param->ivisc == VISC_FULL) {
922 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Full Rheology \n"));
923 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Viscosity range: %g--%g Pa-sec \n", (double)param->eta0, (double)(param->visc_cutoff * param->eta0)));
924 } else {
925 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Invalid! \n"));
926 PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
927 }
928
929 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Boundary condition:"));
930 if (param->ibound == BC_ANALYTIC) {
931 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Isoviscous Analytic Dirichlet \n"));
932 } else if (param->ibound == BC_NOSTRESS) {
933 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Stress-Free (normal & shear stress)\n"));
934 } else if (param->ibound == BC_EXPERMNT) {
935 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Experimental boundary condition \n"));
936 } else {
937 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Invalid! \n"));
938 PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
939 }
940
941 if (param->output_to_file) {
942 #if defined(PETSC_HAVE_MATLAB)
943 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: Mat file \"%s\"\n", param->filename));
944 #else
945 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Output Destination: PETSc binary file \"%s\"\n", param->filename));
946 #endif
947 }
948 if (param->output_ivisc != param->ivisc) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Output viscosity: -ivisc %" PetscInt_FMT "\n", param->output_ivisc));
949
950 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "---------------------END ex30 PARAM REPORT---------------------\n"));
951 }
952 if (param->param_test) PetscCall(PetscEnd());
953 PetscFunctionReturn(PETSC_SUCCESS);
954 }
955
956 /* ------------------------------------------------------------------- */
957 /* generates an initial guess using the analytic solution for isoviscous
958 corner flow */
Initialize(DM da)959 PetscErrorCode Initialize(DM da)
960 /* ------------------------------------------------------------------- */
961 {
962 AppCtx *user;
963 Parameter *param;
964 GridInfo *grid;
965 PetscInt i, j, is, js, im, jm;
966 Field **x;
967 Vec Xguess;
968
969 PetscFunctionBeginUser;
970 /* Get the fine grid */
971 PetscCall(DMGetApplicationContext(da, &user));
972 Xguess = user->Xguess;
973 param = user->param;
974 grid = user->grid;
975 PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
976 PetscCall(DMDAVecGetArray(da, Xguess, (void **)&x));
977
978 /* Compute initial guess */
979 for (j = js; j < js + jm; j++) {
980 for (i = is; i < is + im; i++) {
981 if (i < j) x[j][i].u = param->cb;
982 else if (j <= grid->jlid) x[j][i].u = 0.0;
983 else x[j][i].u = HorizVelocity(i, j, user);
984
985 if (i <= j) x[j][i].w = param->sb;
986 else if (j <= grid->jlid) x[j][i].w = 0.0;
987 else x[j][i].w = VertVelocity(i, j, user);
988
989 if (i < j || j <= grid->jlid) x[j][i].p = 0.0;
990 else x[j][i].p = Pressure(i, j, user);
991
992 x[j][i].T = PetscMin(grid->dz * (j - 0.5), 1.0);
993 }
994 }
995
996 /* Restore x to Xguess */
997 PetscCall(DMDAVecRestoreArray(da, Xguess, (void **)&x));
998 PetscFunctionReturn(PETSC_SUCCESS);
999 }
1000
1001 /* controls output to a file */
DoOutput(SNES snes,PetscInt its)1002 PetscErrorCode DoOutput(SNES snes, PetscInt its)
1003 {
1004 AppCtx *user;
1005 Parameter *param;
1006 GridInfo *grid;
1007 PetscInt ivt;
1008 PetscMPIInt rank;
1009 PetscViewer viewer;
1010 Vec res, pars;
1011 MPI_Comm comm;
1012 DM da;
1013
1014 PetscFunctionBeginUser;
1015 PetscCall(SNESGetDM(snes, &da));
1016 PetscCall(DMGetApplicationContext(da, &user));
1017 param = user->param;
1018 grid = user->grid;
1019 ivt = param->ivisc;
1020
1021 param->ivisc = param->output_ivisc;
1022
1023 /* compute final residual and final viscosity/strain rate fields */
1024 PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
1025 PetscCall(ViscosityField(da, user->x, user->Xguess));
1026
1027 /* get the communicator and the rank of the processor */
1028 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
1029 PetscCallMPI(MPI_Comm_rank(comm, &rank));
1030
1031 if (param->output_to_file) { /* send output to binary file */
1032 PetscCall(VecCreate(comm, &pars));
1033 if (rank == 0) { /* on processor 0 */
1034 PetscCall(VecSetSizes(pars, 20, PETSC_DETERMINE));
1035 PetscCall(VecSetFromOptions(pars));
1036 PetscCall(VecSetValue(pars, 0, (PetscScalar)grid->ni, INSERT_VALUES));
1037 PetscCall(VecSetValue(pars, 1, (PetscScalar)grid->nj, INSERT_VALUES));
1038 PetscCall(VecSetValue(pars, 2, (PetscScalar)grid->dx, INSERT_VALUES));
1039 PetscCall(VecSetValue(pars, 3, (PetscScalar)grid->dz, INSERT_VALUES));
1040 PetscCall(VecSetValue(pars, 4, (PetscScalar)param->L, INSERT_VALUES));
1041 PetscCall(VecSetValue(pars, 5, (PetscScalar)param->V, INSERT_VALUES));
1042 /* skipped 6 intentionally */
1043 PetscCall(VecSetValue(pars, 7, (PetscScalar)param->slab_dip, INSERT_VALUES));
1044 PetscCall(VecSetValue(pars, 8, (PetscScalar)grid->jlid, INSERT_VALUES));
1045 PetscCall(VecSetValue(pars, 9, (PetscScalar)param->lid_depth, INSERT_VALUES));
1046 PetscCall(VecSetValue(pars, 10, (PetscScalar)grid->jfault, INSERT_VALUES));
1047 PetscCall(VecSetValue(pars, 11, (PetscScalar)param->fault_depth, INSERT_VALUES));
1048 PetscCall(VecSetValue(pars, 12, (PetscScalar)param->potentialT, INSERT_VALUES));
1049 PetscCall(VecSetValue(pars, 13, (PetscScalar)param->ivisc, INSERT_VALUES));
1050 PetscCall(VecSetValue(pars, 14, (PetscScalar)param->visc_cutoff, INSERT_VALUES));
1051 PetscCall(VecSetValue(pars, 15, (PetscScalar)param->ibound, INSERT_VALUES));
1052 PetscCall(VecSetValue(pars, 16, (PetscScalar)its, INSERT_VALUES));
1053 } else { /* on some other processor */
1054 PetscCall(VecSetSizes(pars, 0, PETSC_DETERMINE));
1055 PetscCall(VecSetFromOptions(pars));
1056 }
1057 PetscCall(VecAssemblyBegin(pars));
1058 PetscCall(VecAssemblyEnd(pars));
1059
1060 /* create viewer */
1061 #if defined(PETSC_HAVE_MATLAB)
1062 PetscCall(PetscViewerMatlabOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
1063 #else
1064 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, param->filename, FILE_MODE_WRITE, &viewer));
1065 #endif
1066
1067 /* send vectors to viewer */
1068 PetscCall(PetscObjectSetName((PetscObject)res, "res"));
1069 PetscCall(VecView(res, viewer));
1070 PetscCall(PetscObjectSetName((PetscObject)user->x, "out"));
1071 PetscCall(VecView(user->x, viewer));
1072 PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "aux"));
1073 PetscCall(VecView(user->Xguess, viewer));
1074 PetscCall(StressField(da)); /* compute stress fields */
1075 PetscCall(PetscObjectSetName((PetscObject)user->Xguess, "str"));
1076 PetscCall(VecView(user->Xguess, viewer));
1077 PetscCall(PetscObjectSetName((PetscObject)pars, "par"));
1078 PetscCall(VecView(pars, viewer));
1079
1080 /* destroy viewer and vector */
1081 PetscCall(PetscViewerDestroy(&viewer));
1082 PetscCall(VecDestroy(&pars));
1083 }
1084
1085 param->ivisc = ivt;
1086 PetscFunctionReturn(PETSC_SUCCESS);
1087 }
1088
1089 /* ------------------------------------------------------------------- */
1090 /* Compute both the second invariant of the strain rate tensor and the viscosity, at both cell centers and cell corners */
ViscosityField(DM da,Vec X,Vec V)1091 PetscErrorCode ViscosityField(DM da, Vec X, Vec V)
1092 /* ------------------------------------------------------------------- */
1093 {
1094 AppCtx *user;
1095 Parameter *param;
1096 GridInfo *grid;
1097 Vec localX;
1098 Field **v, **x;
1099 PetscReal eps, /* dx,*/ dz, T, epsC, TC;
1100 PetscInt i, j, is, js, im, jm, ilim, jlim, ivt;
1101
1102 PetscFunctionBeginUser;
1103 PetscCall(DMGetApplicationContext(da, &user));
1104 param = user->param;
1105 grid = user->grid;
1106 ivt = param->ivisc;
1107 param->ivisc = param->output_ivisc;
1108
1109 PetscCall(DMGetLocalVector(da, &localX));
1110 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
1111 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
1112 PetscCall(DMDAVecGetArray(da, localX, (void **)&x));
1113 PetscCall(DMDAVecGetArray(da, V, (void **)&v));
1114
1115 /* Parameters */
1116 /* dx = grid->dx; */ dz = grid->dz;
1117
1118 ilim = grid->ni - 1;
1119 jlim = grid->nj - 1;
1120
1121 /* Compute real temperature, strain rate and viscosity */
1122 PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
1123 for (j = js; j < js + jm; j++) {
1124 for (i = is; i < is + im; i++) {
1125 T = PetscRealPart(param->potentialT * x[j][i].T * PetscExpScalar((j - 0.5) * dz * param->z_scale));
1126 if (i < ilim && j < jlim) {
1127 TC = PetscRealPart(param->potentialT * TInterp(x, i, j) * PetscExpScalar(j * dz * param->z_scale));
1128 } else {
1129 TC = T;
1130 }
1131 eps = PetscRealPart(CalcSecInv(x, i, j, CELL_CENTER, user));
1132 epsC = PetscRealPart(CalcSecInv(x, i, j, CELL_CORNER, user));
1133
1134 v[j][i].u = eps;
1135 v[j][i].w = epsC;
1136 v[j][i].p = Viscosity(T, eps, dz * (j - 0.5), param);
1137 v[j][i].T = Viscosity(TC, epsC, dz * j, param);
1138 }
1139 }
1140 PetscCall(DMDAVecRestoreArray(da, V, (void **)&v));
1141 PetscCall(DMDAVecRestoreArray(da, localX, (void **)&x));
1142 PetscCall(DMRestoreLocalVector(da, &localX));
1143
1144 param->ivisc = ivt;
1145 PetscFunctionReturn(PETSC_SUCCESS);
1146 }
1147
1148 /* ------------------------------------------------------------------- */
1149 /* post-processing: compute stress everywhere */
StressField(DM da)1150 PetscErrorCode StressField(DM da)
1151 /* ------------------------------------------------------------------- */
1152 {
1153 AppCtx *user;
1154 PetscInt i, j, is, js, im, jm;
1155 Vec locVec;
1156 Field **x, **y;
1157
1158 PetscFunctionBeginUser;
1159 PetscCall(DMGetApplicationContext(da, &user));
1160
1161 /* Get the fine grid of Xguess and X */
1162 PetscCall(DMDAGetCorners(da, &is, &js, NULL, &im, &jm, NULL));
1163 PetscCall(DMDAVecGetArray(da, user->Xguess, (void **)&x));
1164
1165 PetscCall(DMGetLocalVector(da, &locVec));
1166 PetscCall(DMGlobalToLocalBegin(da, user->x, INSERT_VALUES, locVec));
1167 PetscCall(DMGlobalToLocalEnd(da, user->x, INSERT_VALUES, locVec));
1168 PetscCall(DMDAVecGetArray(da, locVec, (void **)&y));
1169
1170 /* Compute stress on the corner points */
1171 for (j = js; j < js + jm; j++) {
1172 for (i = is; i < is + im; i++) {
1173 x[j][i].u = ShearStress(y, i, j, CELL_CENTER, user);
1174 x[j][i].w = ShearStress(y, i, j, CELL_CORNER, user);
1175 x[j][i].p = XNormalStress(y, i, j, CELL_CENTER, user);
1176 x[j][i].T = ZNormalStress(y, i, j, CELL_CENTER, user);
1177 }
1178 }
1179
1180 /* Restore the fine grid of Xguess and X */
1181 PetscCall(DMDAVecRestoreArray(da, user->Xguess, (void **)&x));
1182 PetscCall(DMDAVecRestoreArray(da, locVec, (void **)&y));
1183 PetscCall(DMRestoreLocalVector(da, &locVec));
1184 PetscFunctionReturn(PETSC_SUCCESS);
1185 }
1186
1187 /*=====================================================================
1188 UTILITY FUNCTIONS
1189 =====================================================================*/
1190
1191 /* returns the velocity of the subducting slab and handles fault nodes for BC */
SlabVel(char c,PetscInt i,PetscInt j,AppCtx * user)1192 static inline PetscScalar SlabVel(char c, PetscInt i, PetscInt j, AppCtx *user)
1193 {
1194 Parameter *param = user->param;
1195 GridInfo *grid = user->grid;
1196
1197 if (c == 'U' || c == 'u') {
1198 if (i < j - 1) return param->cb;
1199 else if (j <= grid->jfault) return 0.0;
1200 else return param->cb;
1201
1202 } else {
1203 if (i < j) return param->sb;
1204 else if (j <= grid->jfault) return 0.0;
1205 else return param->sb;
1206 }
1207 }
1208
1209 /* solution to diffusive half-space cooling model for BC */
PlateModel(PetscInt j,PetscInt plate,AppCtx * user)1210 static inline PetscScalar PlateModel(PetscInt j, PetscInt plate, AppCtx *user)
1211 {
1212 Parameter *param = user->param;
1213 PetscScalar z;
1214 if (plate == PLATE_LID) z = (j - 0.5) * user->grid->dz;
1215 else z = (j - 0.5) * user->grid->dz * param->cb; /* PLATE_SLAB */
1216 #if defined(PETSC_HAVE_ERF)
1217 return (PetscReal)(erf((double)PetscRealPart(z * param->L / 2.0 / param->skt)));
1218 #else
1219 (*PetscErrorPrintf)("erf() not available on this machine\n");
1220 MPI_Abort(PETSC_COMM_SELF, 1);
1221 #endif
1222 }
1223
1224 /*=====================================================================
1225 INTERACTIVE SIGNAL HANDLING
1226 =====================================================================*/
1227
1228 /* ------------------------------------------------------------------- */
SNESConverged_Interactive(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason * reason,PetscCtx ctx)1229 PetscErrorCode SNESConverged_Interactive(SNES snes, PetscInt it, PetscReal xnorm, PetscReal snorm, PetscReal fnorm, SNESConvergedReason *reason, PetscCtx ctx)
1230 /* ------------------------------------------------------------------- */
1231 {
1232 AppCtx *user = (AppCtx *)ctx;
1233 Parameter *param = user->param;
1234 KSP ksp;
1235
1236 PetscFunctionBeginUser;
1237 if (param->interrupted) {
1238 param->interrupted = PETSC_FALSE;
1239 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: exiting SNES solve. \n"));
1240 *reason = SNES_CONVERGED_FNORM_ABS;
1241 PetscFunctionReturn(PETSC_SUCCESS);
1242 } else if (param->toggle_kspmon) {
1243 param->toggle_kspmon = PETSC_FALSE;
1244
1245 PetscCall(SNESGetKSP(snes, &ksp));
1246
1247 if (param->kspmon) {
1248 PetscCall(KSPMonitorCancel(ksp));
1249
1250 param->kspmon = PETSC_FALSE;
1251 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: deactivating ksp singular value monitor. \n"));
1252 } else {
1253 PetscViewerAndFormat *vf;
1254 PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
1255 PetscCall(KSPMonitorSet(ksp, (PetscErrorCode (*)(KSP, PetscInt, PetscReal, void *))KSPMonitorSingularValue, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
1256
1257 param->kspmon = PETSC_TRUE;
1258 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "USER SIGNAL: activating ksp singular value monitor. \n"));
1259 }
1260 }
1261 PetscCall(SNESConvergedDefault(snes, it, xnorm, snorm, fnorm, reason, ctx));
1262 PetscFunctionReturn(PETSC_SUCCESS);
1263 }
1264
1265 /* ------------------------------------------------------------------- */
1266 #include <signal.h>
InteractiveHandler(int signum,PetscCtx ctx)1267 PetscErrorCode InteractiveHandler(int signum, PetscCtx ctx)
1268 /* ------------------------------------------------------------------- */
1269 {
1270 AppCtx *user = (AppCtx *)ctx;
1271 Parameter *param = user->param;
1272
1273 if (signum == SIGILL) {
1274 param->toggle_kspmon = PETSC_TRUE;
1275 #if !defined(PETSC_MISSING_SIGCONT)
1276 } else if (signum == SIGCONT) {
1277 param->interrupted = PETSC_TRUE;
1278 #endif
1279 #if !defined(PETSC_MISSING_SIGURG)
1280 } else if (signum == SIGURG) {
1281 param->stop_solve = PETSC_TRUE;
1282 #endif
1283 }
1284 return PETSC_SUCCESS;
1285 }
1286
1287 /* main call-back function that computes the processor-local piece of the residual */
FormFunctionLocal(DMDALocalInfo * info,Field ** x,Field ** f,void * ptr)1288 PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, Field **x, Field **f, void *ptr)
1289 {
1290 AppCtx *user = (AppCtx *)ptr;
1291 Parameter *param = user->param;
1292 GridInfo *grid = user->grid;
1293 PetscScalar mag_w, mag_u;
1294 PetscInt i, j, mx, mz, ilim, jlim;
1295 PetscInt is, ie, js, je, ibound; /* ,ivisc */
1296
1297 PetscFunctionBeginUser;
1298 /* Define global and local grid parameters */
1299 mx = info->mx;
1300 mz = info->my;
1301 ilim = mx - 1;
1302 jlim = mz - 1;
1303 is = info->xs;
1304 ie = info->xs + info->xm;
1305 js = info->ys;
1306 je = info->ys + info->ym;
1307
1308 /* Define geometric and numeric parameters */
1309 /* ivisc = param->ivisc; */ ibound = param->ibound;
1310
1311 for (j = js; j < je; j++) {
1312 for (i = is; i < ie; i++) {
1313 /************* X-MOMENTUM/VELOCITY *************/
1314 if (i < j) f[j][i].u = x[j][i].u - SlabVel('U', i, j, user);
1315 else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1316 /* in the lithospheric lid */
1317 f[j][i].u = x[j][i].u - 0.0;
1318 } else if (i == ilim) {
1319 /* on the right side boundary */
1320 if (ibound == BC_ANALYTIC) {
1321 f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
1322 } else {
1323 f[j][i].u = XNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
1324 }
1325
1326 } else if (j == jlim) {
1327 /* on the bottom boundary */
1328 if (ibound == BC_ANALYTIC) {
1329 f[j][i].u = x[j][i].u - HorizVelocity(i, j, user);
1330 } else if (ibound == BC_NOSTRESS) {
1331 f[j][i].u = XMomentumResidual(x, i, j, user);
1332 } else {
1333 /* experimental boundary condition */
1334 }
1335
1336 } else {
1337 /* in the mantle wedge */
1338 f[j][i].u = XMomentumResidual(x, i, j, user);
1339 }
1340
1341 /************* Z-MOMENTUM/VELOCITY *************/
1342 if (i <= j) {
1343 f[j][i].w = x[j][i].w - SlabVel('W', i, j, user);
1344
1345 } else if (j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1346 /* in the lithospheric lid */
1347 f[j][i].w = x[j][i].w - 0.0;
1348
1349 } else if (j == jlim) {
1350 /* on the bottom boundary */
1351 if (ibound == BC_ANALYTIC) {
1352 f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
1353 } else {
1354 f[j][i].w = ZNormalStress(x, i, j, CELL_CENTER, user) - EPS_ZERO;
1355 }
1356
1357 } else if (i == ilim) {
1358 /* on the right side boundary */
1359 if (ibound == BC_ANALYTIC) {
1360 f[j][i].w = x[j][i].w - VertVelocity(i, j, user);
1361 } else if (ibound == BC_NOSTRESS) {
1362 f[j][i].w = ZMomentumResidual(x, i, j, user);
1363 } else {
1364 /* experimental boundary condition */
1365 }
1366
1367 } else {
1368 /* in the mantle wedge */
1369 f[j][i].w = ZMomentumResidual(x, i, j, user);
1370 }
1371
1372 /************* CONTINUITY/PRESSURE *************/
1373 if (i < j || j <= grid->jlid || (j < grid->corner + grid->inose && i < grid->corner + grid->inose)) {
1374 /* in the lid or slab */
1375 f[j][i].p = x[j][i].p;
1376
1377 } else if ((i == ilim || j == jlim) && ibound == BC_ANALYTIC) {
1378 /* on an analytic boundary */
1379 f[j][i].p = x[j][i].p - Pressure(i, j, user);
1380
1381 } else {
1382 /* in the mantle wedge */
1383 f[j][i].p = ContinuityResidual(x, i, j, user);
1384 }
1385
1386 /************* TEMPERATURE *************/
1387 if (j == 0) {
1388 /* on the surface */
1389 f[j][i].T = x[j][i].T + x[j + 1][i].T + PetscMax(PetscRealPart(x[j][i].T), 0.0);
1390
1391 } else if (i == 0) {
1392 /* slab inflow boundary */
1393 f[j][i].T = x[j][i].T - PlateModel(j, PLATE_SLAB, user);
1394
1395 } else if (i == ilim) {
1396 /* right side boundary */
1397 mag_u = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j][i - 1].u) / param->cb, 1.0), 0.0), 5);
1398 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);
1399
1400 } else if (j == jlim) {
1401 /* bottom boundary */
1402 mag_w = 1.0 - PetscPowRealInt(1.0 - PetscMax(PetscMin(PetscRealPart(x[j - 1][i].w) / param->sb, 1.0), 0.0), 5);
1403 f[j][i].T = x[j][i].T - mag_w * x[j - 1][i - 1].T - (1.0 - mag_w);
1404
1405 } else {
1406 /* in the mantle wedge */
1407 f[j][i].T = EnergyResidual(x, i, j, user);
1408 }
1409 }
1410 }
1411 PetscFunctionReturn(PETSC_SUCCESS);
1412 }
1413
1414 /*TEST
1415
1416 build:
1417 requires: !complex erf
1418
1419 test:
1420 args: -ni 18 -fp_trap 0
1421 filter: grep -v Destination
1422 requires: !single
1423
1424 TEST*/
1425