Actual source code: ex14.c
1: static const char help[] = "Toy hydrostatic ice flow with multigrid in 3D.\n\
2: \n\
3: Solves the hydrostatic (aka Blatter/Pattyn/First Order) equations for ice sheet flow\n\
4: using multigrid. The ice uses a power-law rheology with \"Glen\" exponent 3 (corresponds\n\
5: to p=4/3 in a p-Laplacian). The focus is on ISMIP-HOM experiments which assume periodic\n\
6: boundary conditions in the x- and y-directions.\n\
7: \n\
8: Equations are rescaled so that the domain size and solution are O(1), details of this scaling\n\
9: can be controlled by the options -units_meter, -units_second, and -units_kilogram.\n\
10: \n\
11: A VTK StructuredGrid output file can be written using the option -o filename.vts\n\
12: \n\n";
14: /*
15: The equations for horizontal velocity (u,v) are
17: - [eta (4 u_x + 2 v_y)]_x - [eta (u_y + v_x)]_y - [eta u_z]_z + rho g s_x = 0
18: - [eta (4 v_y + 2 u_x)]_y - [eta (u_y + v_x)]_x - [eta v_z]_z + rho g s_y = 0
20: where
22: eta = B/2 (epsilon + gamma)^((p-2)/2)
24: is the nonlinear effective viscosity with regularization epsilon and hardness parameter B,
25: written in terms of the second invariant
27: gamma = u_x^2 + v_y^2 + u_x v_y + (1/4) (u_y + v_x)^2 + (1/4) u_z^2 + (1/4) v_z^2
29: The surface boundary conditions are the natural conditions. The basal boundary conditions
30: are either no-slip, or Navier (linear) slip with spatially variant friction coefficient beta^2.
32: In the code, the equations for (u,v) are multiplied through by 1/(rho g) so that residuals are O(1).
34: The discretization is Q1 finite elements, managed by a DMDA. The grid is never distorted in the
35: map (x,y) plane, but the bed and surface may be bumpy. This is handled as usual in FEM, through
36: the Jacobian of the coordinate transformation from a reference element to the physical element.
38: Since ice-flow is tightly coupled in the z-direction (within columns), the DMDA is managed
39: specially so that columns are never distributed, and are always contiguous in memory.
40: This amounts to reversing the meaning of X,Y,Z compared to the DMDA's internal interpretation,
41: and then indexing as vec[i][j][k]. The exotic coarse spaces require 2D DMDAs which are made to
42: use compatible domain decomposition relative to the 3D DMDAs.
44: */
46: #include <petscts.h>
47: #include <petscdm.h>
48: #include <petscdmda.h>
49: #include <petscdmcomposite.h>
50: #include <ctype.h> /* toupper() */
51: #include <petsc/private/petscimpl.h>
54: #if defined(__cplusplus) /* C++ restrict is nonstandard and compilers have inconsistent rules about where it can be used */
55: #define restrict
56: #else
57: #define restrict PETSC_RESTRICT
58: #endif
59: #endif
61: static PetscClassId THI_CLASSID;
63: typedef enum {
64: QUAD_GAUSS,
65: QUAD_LOBATTO
66: } QuadratureType;
67: static const char *QuadratureTypes[] = {"gauss", "lobatto", "QuadratureType", "QUAD_", 0};
68: static const PetscReal HexQWeights[8] = {1, 1, 1, 1, 1, 1, 1, 1};
69: static const PetscReal HexQNodes[] = {-0.57735026918962573, 0.57735026918962573};
70: #define G 0.57735026918962573
71: #define H (0.5 * (1. + G))
72: #define L (0.5 * (1. - G))
73: #define M (-0.5)
74: #define P (0.5)
75: /* Special quadrature: Lobatto in horizontal, Gauss in vertical */
76: static const PetscReal HexQInterp_Lobatto[8][8] = {
77: {H, 0, 0, 0, L, 0, 0, 0},
78: {0, H, 0, 0, 0, L, 0, 0},
79: {0, 0, H, 0, 0, 0, L, 0},
80: {0, 0, 0, H, 0, 0, 0, L},
81: {L, 0, 0, 0, H, 0, 0, 0},
82: {0, L, 0, 0, 0, H, 0, 0},
83: {0, 0, L, 0, 0, 0, H, 0},
84: {0, 0, 0, L, 0, 0, 0, H}
85: };
86: static const PetscReal HexQDeriv_Lobatto[8][8][3] = {
87: {{M * H, M * H, M}, {P * H, 0, 0}, {0, 0, 0}, {0, P * H, 0}, {M * L, M * L, P}, {P * L, 0, 0}, {0, 0, 0}, {0, P * L, 0} },
88: {{M * H, 0, 0}, {P * H, M * H, M}, {0, P * H, 0}, {0, 0, 0}, {M * L, 0, 0}, {P * L, M * L, P}, {0, P * L, 0}, {0, 0, 0} },
89: {{0, 0, 0}, {0, M * H, 0}, {P * H, P * H, M}, {M * H, 0, 0}, {0, 0, 0}, {0, M * L, 0}, {P * L, P * L, P}, {M * L, 0, 0} },
90: {{0, M * H, 0}, {0, 0, 0}, {P * H, 0, 0}, {M * H, P * H, M}, {0, M * L, 0}, {0, 0, 0}, {P * L, 0, 0}, {M * L, P * L, P}},
91: {{M * L, M * L, M}, {P * L, 0, 0}, {0, 0, 0}, {0, P * L, 0}, {M * H, M * H, P}, {P * H, 0, 0}, {0, 0, 0}, {0, P * H, 0} },
92: {{M * L, 0, 0}, {P * L, M * L, M}, {0, P * L, 0}, {0, 0, 0}, {M * H, 0, 0}, {P * H, M * H, P}, {0, P * H, 0}, {0, 0, 0} },
93: {{0, 0, 0}, {0, M * L, 0}, {P * L, P * L, M}, {M * L, 0, 0}, {0, 0, 0}, {0, M * H, 0}, {P * H, P * H, P}, {M * H, 0, 0} },
94: {{0, M * L, 0}, {0, 0, 0}, {P * L, 0, 0}, {M * L, P * L, M}, {0, M * H, 0}, {0, 0, 0}, {P * H, 0, 0}, {M * H, P * H, P}}
95: };
96: /* Standard Gauss */
97: static const PetscReal HexQInterp_Gauss[8][8] = {
98: {H * H * H, L * H * H, L * L * H, H * L * H, H * H * L, L * H * L, L * L * L, H * L * L},
99: {L * H * H, H * H * H, H * L * H, L * L * H, L * H * L, H * H * L, H * L * L, L * L * L},
100: {L * L * H, H * L * H, H * H * H, L * H * H, L * L * L, H * L * L, H * H * L, L * H * L},
101: {H * L * H, L * L * H, L * H * H, H * H * H, H * L * L, L * L * L, L * H * L, H * H * L},
102: {H * H * L, L * H * L, L * L * L, H * L * L, H * H * H, L * H * H, L * L * H, H * L * H},
103: {L * H * L, H * H * L, H * L * L, L * L * L, L * H * H, H * H * H, H * L * H, L * L * H},
104: {L * L * L, H * L * L, H * H * L, L * H * L, L * L * H, H * L * H, H * H * H, L * H * H},
105: {H * L * L, L * L * L, L * H * L, H * H * L, H * L * H, L * L * H, L * H * H, H * H * H}
106: };
107: static const PetscReal HexQDeriv_Gauss[8][8][3] = {
108: {{M * H * H, H * M * H, H * H * M}, {P * H * H, L * M * H, L * H * M}, {P * L * H, L * P * H, L * L * M}, {M * L * H, H * P * H, H * L * M}, {M * H * L, H * M * L, H * H * P}, {P * H * L, L * M * L, L * H * P}, {P * L * L, L * P * L, L * L * P}, {M * L * L, H * P * L, H * L * P}},
109: {{M * H * H, L * M * H, L * H * M}, {P * H * H, H * M * H, H * H * M}, {P * L * H, H * P * H, H * L * M}, {M * L * H, L * P * H, L * L * M}, {M * H * L, L * M * L, L * H * P}, {P * H * L, H * M * L, H * H * P}, {P * L * L, H * P * L, H * L * P}, {M * L * L, L * P * L, L * L * P}},
110: {{M * L * H, L * M * H, L * L * M}, {P * L * H, H * M * H, H * L * M}, {P * H * H, H * P * H, H * H * M}, {M * H * H, L * P * H, L * H * M}, {M * L * L, L * M * L, L * L * P}, {P * L * L, H * M * L, H * L * P}, {P * H * L, H * P * L, H * H * P}, {M * H * L, L * P * L, L * H * P}},
111: {{M * L * H, H * M * H, H * L * M}, {P * L * H, L * M * H, L * L * M}, {P * H * H, L * P * H, L * H * M}, {M * H * H, H * P * H, H * H * M}, {M * L * L, H * M * L, H * L * P}, {P * L * L, L * M * L, L * L * P}, {P * H * L, L * P * L, L * H * P}, {M * H * L, H * P * L, H * H * P}},
112: {{M * H * L, H * M * L, H * H * M}, {P * H * L, L * M * L, L * H * M}, {P * L * L, L * P * L, L * L * M}, {M * L * L, H * P * L, H * L * M}, {M * H * H, H * M * H, H * H * P}, {P * H * H, L * M * H, L * H * P}, {P * L * H, L * P * H, L * L * P}, {M * L * H, H * P * H, H * L * P}},
113: {{M * H * L, L * M * L, L * H * M}, {P * H * L, H * M * L, H * H * M}, {P * L * L, H * P * L, H * L * M}, {M * L * L, L * P * L, L * L * M}, {M * H * H, L * M * H, L * H * P}, {P * H * H, H * M * H, H * H * P}, {P * L * H, H * P * H, H * L * P}, {M * L * H, L * P * H, L * L * P}},
114: {{M * L * L, L * M * L, L * L * M}, {P * L * L, H * M * L, H * L * M}, {P * H * L, H * P * L, H * H * M}, {M * H * L, L * P * L, L * H * M}, {M * L * H, L * M * H, L * L * P}, {P * L * H, H * M * H, H * L * P}, {P * H * H, H * P * H, H * H * P}, {M * H * H, L * P * H, L * H * P}},
115: {{M * L * L, H * M * L, H * L * M}, {P * L * L, L * M * L, L * L * M}, {P * H * L, L * P * L, L * H * M}, {M * H * L, H * P * L, H * H * M}, {M * L * H, H * M * H, H * L * P}, {P * L * H, L * M * H, L * L * P}, {P * H * H, L * P * H, L * H * P}, {M * H * H, H * P * H, H * H * P}}
116: };
117: static const PetscReal (*HexQInterp)[8], (*HexQDeriv)[8][3];
118: /* Standard 2x2 Gauss quadrature for the bottom layer. */
119: static const PetscReal QuadQInterp[4][4] = {
120: {H * H, L * H, L * L, H * L},
121: {L * H, H * H, H * L, L * L},
122: {L * L, H * L, H * H, L * H},
123: {H * L, L * L, L * H, H * H}
124: };
125: static const PetscReal QuadQDeriv[4][4][2] = {
126: {{M * H, M * H}, {P * H, M * L}, {P * L, P * L}, {M * L, P * H}},
127: {{M * H, M * L}, {P * H, M * H}, {P * L, P * H}, {M * L, P * L}},
128: {{M * L, M * L}, {P * L, M * H}, {P * H, P * H}, {M * H, P * L}},
129: {{M * L, M * H}, {P * L, M * L}, {P * H, P * L}, {M * H, P * H}}
130: };
131: #undef G
132: #undef H
133: #undef L
134: #undef M
135: #undef P
137: #define HexExtract(x, i, j, k, n) \
138: do { \
139: (n)[0] = (x)[i][j][k]; \
140: (n)[1] = (x)[(i) + 1][j][k]; \
141: (n)[2] = (x)[(i) + 1][(j) + 1][k]; \
142: (n)[3] = (x)[i][(j) + 1][k]; \
143: (n)[4] = (x)[i][j][(k) + 1]; \
144: (n)[5] = (x)[(i) + 1][j][(k) + 1]; \
145: (n)[6] = (x)[(i) + 1][(j) + 1][(k) + 1]; \
146: (n)[7] = (x)[i][(j) + 1][(k) + 1]; \
147: } while (0)
149: #define HexExtractRef(x, i, j, k, n) \
150: do { \
151: (n)[0] = &(x)[i][j][k]; \
152: (n)[1] = &(x)[(i) + 1][j][k]; \
153: (n)[2] = &(x)[(i) + 1][(j) + 1][k]; \
154: (n)[3] = &(x)[i][(j) + 1][k]; \
155: (n)[4] = &(x)[i][j][(k) + 1]; \
156: (n)[5] = &(x)[(i) + 1][j][(k) + 1]; \
157: (n)[6] = &(x)[(i) + 1][(j) + 1][(k) + 1]; \
158: (n)[7] = &(x)[i][(j) + 1][(k) + 1]; \
159: } while (0)
161: #define QuadExtract(x, i, j, n) \
162: do { \
163: (n)[0] = (x)[i][j]; \
164: (n)[1] = (x)[(i) + 1][j]; \
165: (n)[2] = (x)[(i) + 1][(j) + 1]; \
166: (n)[3] = (x)[i][(j) + 1]; \
167: } while (0)
169: static PetscScalar Sqr(PetscScalar a)
170: {
171: return a * a;
172: }
174: static void HexGrad(const PetscReal dphi[][3], const PetscReal zn[], PetscReal dz[])
175: {
176: PetscInt i;
177: dz[0] = dz[1] = dz[2] = 0;
178: for (i = 0; i < 8; i++) {
179: dz[0] += dphi[i][0] * zn[i];
180: dz[1] += dphi[i][1] * zn[i];
181: dz[2] += dphi[i][2] * zn[i];
182: }
183: }
185: static void HexComputeGeometry(PetscInt q, PetscReal hx, PetscReal hy, const PetscReal dz[restrict], PetscReal phi[restrict], PetscReal dphi[restrict][3], PetscReal *restrict jw)
186: {
187: const PetscReal jac[3][3] =
188: {
189: {hx / 2, 0, 0 },
190: {0, hy / 2, 0 },
191: {dz[0], dz[1], dz[2]}
192: },
193: ijac[3][3] = {{1 / jac[0][0], 0, 0}, {0, 1 / jac[1][1], 0}, {-jac[2][0] / (jac[0][0] * jac[2][2]), -jac[2][1] / (jac[1][1] * jac[2][2]), 1 / jac[2][2]}}, jdet = jac[0][0] * jac[1][1] * jac[2][2];
194: PetscInt i;
196: for (i = 0; i < 8; i++) {
197: const PetscReal *dphir = HexQDeriv[q][i];
198: phi[i] = HexQInterp[q][i];
199: dphi[i][0] = dphir[0] * ijac[0][0] + dphir[1] * ijac[1][0] + dphir[2] * ijac[2][0];
200: dphi[i][1] = dphir[0] * ijac[0][1] + dphir[1] * ijac[1][1] + dphir[2] * ijac[2][1];
201: dphi[i][2] = dphir[0] * ijac[0][2] + dphir[1] * ijac[1][2] + dphir[2] * ijac[2][2];
202: }
203: *jw = 1.0 * jdet;
204: }
206: typedef struct _p_THI *THI;
207: typedef struct _n_Units *Units;
209: typedef struct {
210: PetscScalar u, v;
211: } Node;
213: typedef struct {
214: PetscScalar b; /* bed */
215: PetscScalar h; /* thickness */
216: PetscScalar beta2; /* friction */
217: } PrmNode;
219: #define FieldSize(ntype) ((PetscInt)(sizeof(ntype) / sizeof(PetscScalar)))
220: #define FieldOffset(ntype, member) ((PetscInt)(offsetof(ntype, member) / sizeof(PetscScalar)))
221: #define FieldIndex(ntype, i, member) ((PetscInt)((i) * FieldSize(ntype) + FieldOffset(ntype, member)))
222: #define NODE_SIZE FieldSize(Node)
223: #define PRMNODE_SIZE FieldSize(PrmNode)
225: typedef struct {
226: PetscReal min, max, cmin, cmax;
227: } PRange;
229: struct _p_THI {
230: PETSCHEADER(int);
231: void (*initialize)(THI, PetscReal x, PetscReal y, PrmNode *p);
232: PetscInt nlevels;
233: PetscInt zlevels;
234: PetscReal Lx, Ly, Lz; /* Model domain */
235: PetscReal alpha; /* Bed angle */
236: Units units;
237: PetscReal dirichlet_scale;
238: PetscReal ssa_friction_scale;
239: PetscReal inertia;
240: PRange eta;
241: PRange beta2;
242: struct {
243: PetscReal Bd2, eps, exponent, glen_n;
244: } viscosity;
245: struct {
246: PetscReal irefgam, eps2, exponent;
247: } friction;
248: struct {
249: PetscReal rate, exponent, refvel;
250: } erosion;
251: PetscReal rhog;
252: PetscBool no_slip;
253: PetscBool verbose;
254: char *mattype;
255: char *monitor_basename;
256: PetscInt monitor_interval;
257: };
259: struct _n_Units {
260: /* fundamental */
261: PetscReal meter;
262: PetscReal kilogram;
263: PetscReal second;
264: /* derived */
265: PetscReal Pascal;
266: PetscReal year;
267: };
269: static void PrmHexGetZ(const PrmNode pn[], PetscInt k, PetscInt zm, PetscReal zn[])
270: {
271: const PetscScalar zm1 = zm - 1, znl[8] = {pn[0].b + pn[0].h * (PetscScalar)k / zm1, pn[1].b + pn[1].h * (PetscScalar)k / zm1, pn[2].b + pn[2].h * (PetscScalar)k / zm1, pn[3].b + pn[3].h * (PetscScalar)k / zm1,
272: pn[0].b + pn[0].h * (PetscScalar)(k + 1) / zm1, pn[1].b + pn[1].h * (PetscScalar)(k + 1) / zm1, pn[2].b + pn[2].h * (PetscScalar)(k + 1) / zm1, pn[3].b + pn[3].h * (PetscScalar)(k + 1) / zm1};
273: PetscInt i;
274: for (i = 0; i < 8; i++) zn[i] = PetscRealPart(znl[i]);
275: }
277: /* Compute a gradient of all the 2D fields at four quadrature points. Output for [quadrature_point][direction].field_name */
278: static PetscErrorCode QuadComputeGrad4(const PetscReal dphi[][4][2], PetscReal hx, PetscReal hy, const PrmNode pn[4], PrmNode dp[4][2])
279: {
280: PetscInt q, i, f;
281: const PetscScalar (*restrict pg)[PRMNODE_SIZE] = (const PetscScalar (*)[PRMNODE_SIZE])pn; /* Get generic array pointers to the node */
282: PetscScalar (*restrict dpg)[2][PRMNODE_SIZE] = (PetscScalar (*)[2][PRMNODE_SIZE])dp;
284: PetscFunctionBeginUser;
285: PetscCall(PetscArrayzero(dpg, 4));
286: for (q = 0; q < 4; q++) {
287: for (i = 0; i < 4; i++) {
288: for (f = 0; f < PRMNODE_SIZE; f++) {
289: dpg[q][0][f] += dphi[q][i][0] / hx * pg[i][f];
290: dpg[q][1][f] += dphi[q][i][1] / hy * pg[i][f];
291: }
292: }
293: }
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static inline PetscReal StaggeredMidpoint2D(PetscScalar a, PetscScalar b, PetscScalar c, PetscScalar d)
298: {
299: return 0.5 * PetscRealPart(0.75 * a + 0.75 * b + 0.25 * c + 0.25 * d);
300: }
301: static inline PetscReal UpwindFlux1D(PetscReal u, PetscReal hL, PetscReal hR)
302: {
303: return (u > 0) ? hL * u : hR * u;
304: }
306: #define UpwindFluxXW(x3, x2, h, i, j, k, dj) \
307: UpwindFlux1D(StaggeredMidpoint2D((x3)[i][j][k].u, (x3)[(i) - 1][j][k].u, (x3)[(i) - 1][(j) + (dj)][k].u, (x3)[i][(k) + (dj)][k].u), PetscRealPart(0.75 * (x2)[(i) - 1][j].h + 0.25 * (x2)[(i) - 1][(j) + (dj)].h), \
308: PetscRealPart(0.75 * (x2)[i][j].h + 0.25 * (x2)[i][(j) + (dj)].h))
309: #define UpwindFluxXE(x3, x2, h, i, j, k, dj) \
310: UpwindFlux1D(StaggeredMidpoint2D((x3)[i][j][k].u, (x3)[(i) + 1][j][k].u, (x3)[(i) + 1][(j) + (dj)][k].u, (x3)[i][(k) + (dj)][k].u), PetscRealPart(0.75 * (x2)[i][j].h + 0.25 * (x2)[i][(j) + (dj)].h), \
311: PetscRealPart(0.75 * (x2)[(i) + 1][j].h + 0.25 * (x2)[(i) + 1][(j) + (dj)].h))
312: #define UpwindFluxYS(x3, x2, h, i, j, k, di) \
313: UpwindFlux1D(StaggeredMidpoint2D((x3)[i][j][k].v, (x3)[i][(j) - 1][k].v, (x3)[(i) + (di)][(j) - 1][k].v, (x3)[(i) + (di)][j][k].v), PetscRealPart(0.75 * (x2)[i][(j) - 1].h + 0.25 * (x2)[(i) + (di)][(j) - 1].h), \
314: PetscRealPart(0.75 * (x2)[i][j].h + 0.25 * (x2)[(i) + (di)][j].h))
315: #define UpwindFluxYN(x3, x2, h, i, j, k, di) \
316: UpwindFlux1D(StaggeredMidpoint2D((x3)[i][j][k].v, (x3)[i][(j) + 1][k].v, (x3)[(i) + (di)][(j) + 1][k].v, (x3)[(i) + (di)][j][k].v), PetscRealPart(0.75 * (x2)[i][j].h + 0.25 * (x2)[(i) + (di)][j].h), \
317: PetscRealPart(0.75 * (x2)[i][(j) + 1].h + 0.25 * (x2)[(i) + (di)][(j) + 1].h))
319: static void PrmNodeGetFaceMeasure(const PrmNode **p, PetscInt i, PetscInt j, PetscScalar h[])
320: {
321: /* West */
322: h[0] = StaggeredMidpoint2D(p[i][j].h, p[i - 1][j].h, p[i - 1][j - 1].h, p[i][j - 1].h);
323: h[1] = StaggeredMidpoint2D(p[i][j].h, p[i - 1][j].h, p[i - 1][j + 1].h, p[i][j + 1].h);
324: /* East */
325: h[2] = StaggeredMidpoint2D(p[i][j].h, p[i + 1][j].h, p[i + 1][j + 1].h, p[i][j + 1].h);
326: h[3] = StaggeredMidpoint2D(p[i][j].h, p[i + 1][j].h, p[i + 1][j - 1].h, p[i][j - 1].h);
327: /* South */
328: h[4] = StaggeredMidpoint2D(p[i][j].h, p[i][j - 1].h, p[i + 1][j - 1].h, p[i + 1][j].h);
329: h[5] = StaggeredMidpoint2D(p[i][j].h, p[i][j - 1].h, p[i - 1][j - 1].h, p[i - 1][j].h);
330: /* North */
331: h[6] = StaggeredMidpoint2D(p[i][j].h, p[i][j + 1].h, p[i - 1][j + 1].h, p[i - 1][j].h);
332: h[7] = StaggeredMidpoint2D(p[i][j].h, p[i][j + 1].h, p[i + 1][j + 1].h, p[i + 1][j].h);
333: }
335: /* Tests A and C are from the ISMIP-HOM paper (Pattyn et al. 2008) */
336: static void THIInitialize_HOM_A(THI thi, PetscReal x, PetscReal y, PrmNode *p)
337: {
338: Units units = thi->units;
339: PetscReal s = -x * PetscSinReal(thi->alpha);
340: p->b = s - 1000 * units->meter + 500 * units->meter * PetscSinReal(x * 2 * PETSC_PI / thi->Lx) * PetscSinReal(y * 2 * PETSC_PI / thi->Ly);
341: p->h = s - p->b;
342: p->beta2 = -1e-10; /* This value is not used, but it should not be huge because that would change the finite difference step size */
343: }
345: static void THIInitialize_HOM_C(THI thi, PetscReal x, PetscReal y, PrmNode *p)
346: {
347: Units units = thi->units;
348: PetscReal s = -x * PetscSinReal(thi->alpha);
349: p->b = s - 1000 * units->meter;
350: p->h = s - p->b;
351: /* tau_b = beta2 v is a stress (Pa).
352: * This is a big number in our units (it needs to balance the driving force from the surface), so we scale it by 1/rhog, just like the residual. */
353: p->beta2 = 1000 * (1 + PetscSinReal(x * 2 * PETSC_PI / thi->Lx) * PetscSinReal(y * 2 * PETSC_PI / thi->Ly)) * units->Pascal * units->year / units->meter / thi->rhog;
354: }
356: /* These are just toys */
358: /* From Fred Herman */
359: static void THIInitialize_HOM_F(THI thi, PetscReal x, PetscReal y, PrmNode *p)
360: {
361: Units units = thi->units;
362: PetscReal s = -x * PetscSinReal(thi->alpha);
363: p->b = s - 1000 * units->meter + 100 * units->meter * PetscSinReal(x * 2 * PETSC_PI / thi->Lx); /* * sin(y*2*PETSC_PI/thi->Ly); */
364: p->h = s - p->b;
365: p->h = (1 - (atan((x - thi->Lx / 2) / 1.) + PETSC_PI / 2.) / PETSC_PI) * 500 * units->meter + 1 * units->meter;
366: s = PetscRealPart(p->b + p->h);
367: p->beta2 = -1e-10;
368: /* p->beta2 = 1000 * units->Pascal * units->year / units->meter; */
369: }
371: /* Same bed as test A, free slip everywhere except for a discontinuous jump to a circular sticky region in the middle. */
372: static void THIInitialize_HOM_X(THI thi, PetscReal xx, PetscReal yy, PrmNode *p)
373: {
374: Units units = thi->units;
375: PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; /* [-pi,pi] */
376: PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha);
377: p->b = s - 1000 * units->meter + 500 * units->meter * PetscSinReal(x + PETSC_PI) * PetscSinReal(y + PETSC_PI);
378: p->h = s - p->b;
379: p->beta2 = 1000 * (r < 1 ? 2 : 0) * units->Pascal * units->year / units->meter / thi->rhog;
380: }
382: /* Like Z, but with 200 meter cliffs */
383: static void THIInitialize_HOM_Y(THI thi, PetscReal xx, PetscReal yy, PrmNode *p)
384: {
385: Units units = thi->units;
386: PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; /* [-pi,pi] */
387: PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha);
388: p->b = s - 1000 * units->meter + 500 * units->meter * PetscSinReal(x + PETSC_PI) * PetscSinReal(y + PETSC_PI);
389: if (PetscRealPart(p->b) > -700 * units->meter) p->b += 200 * units->meter;
390: p->h = s - p->b;
391: p->beta2 = 1000 * (1. + PetscSinReal(PetscSqrtReal(16 * r)) / PetscSqrtReal(1e-2 + 16 * r) * PetscCosReal(x * 3 / 2) * PetscCosReal(y * 3 / 2)) * units->Pascal * units->year / units->meter / thi->rhog;
392: }
394: /* Same bed as A, smoothly varying slipperiness, similar to MATLAB's "sombrero" (uncorrelated with bathymetry) */
395: static void THIInitialize_HOM_Z(THI thi, PetscReal xx, PetscReal yy, PrmNode *p)
396: {
397: Units units = thi->units;
398: PetscReal x = xx * 2 * PETSC_PI / thi->Lx - PETSC_PI, y = yy * 2 * PETSC_PI / thi->Ly - PETSC_PI; /* [-pi,pi] */
399: PetscReal r = PetscSqrtReal(x * x + y * y), s = -x * PetscSinReal(thi->alpha);
400: p->b = s - 1000 * units->meter + 500 * units->meter * PetscSinReal(x + PETSC_PI) * PetscSinReal(y + PETSC_PI);
401: p->h = s - p->b;
402: p->beta2 = 1000 * (1. + PetscSinReal(PetscSqrtReal(16 * r)) / PetscSqrtReal(1e-2 + 16 * r) * PetscCosReal(x * 3 / 2) * PetscCosReal(y * 3 / 2)) * units->Pascal * units->year / units->meter / thi->rhog;
403: }
405: static void THIFriction(THI thi, PetscReal rbeta2, PetscReal gam, PetscReal *beta2, PetscReal *dbeta2)
406: {
407: if (thi->friction.irefgam == 0) {
408: Units units = thi->units;
409: thi->friction.irefgam = 1. / (0.5 * PetscSqr(100 * units->meter / units->year));
410: thi->friction.eps2 = 0.5 * PetscSqr(1.e-4 / thi->friction.irefgam);
411: }
412: if (thi->friction.exponent == 0) {
413: *beta2 = rbeta2;
414: *dbeta2 = 0;
415: } else {
416: *beta2 = rbeta2 * PetscPowReal(thi->friction.eps2 + gam * thi->friction.irefgam, thi->friction.exponent);
417: *dbeta2 = thi->friction.exponent * *beta2 / (thi->friction.eps2 + gam * thi->friction.irefgam) * thi->friction.irefgam;
418: }
419: }
421: static void THIViscosity(THI thi, PetscReal gam, PetscReal *eta, PetscReal *deta)
422: {
423: PetscReal Bd2, eps, exponent;
424: if (thi->viscosity.Bd2 == 0) {
425: Units units = thi->units;
426: const PetscReal n = thi->viscosity.glen_n, /* Glen exponent */
427: p = 1. + 1. / n, /* for Stokes */
428: A = 1.e-16 * PetscPowReal(units->Pascal, -n) / units->year, /* softness parameter (Pa^{-n}/s) */
429: B = PetscPowReal(A, -1. / n); /* hardness parameter */
430: thi->viscosity.Bd2 = B / 2;
431: thi->viscosity.exponent = (p - 2) / 2;
432: thi->viscosity.eps = 0.5 * PetscSqr(1e-5 / units->year);
433: }
434: Bd2 = thi->viscosity.Bd2;
435: exponent = thi->viscosity.exponent;
436: eps = thi->viscosity.eps;
437: *eta = Bd2 * PetscPowReal(eps + gam, exponent);
438: *deta = exponent * (*eta) / (eps + gam);
439: }
441: static void THIErosion(THI thi, const Node *vel, PetscScalar *erate, Node *derate)
442: {
443: const PetscScalar magref2 = 1.e-10 + (PetscSqr(vel->u) + PetscSqr(vel->v)) / PetscSqr(thi->erosion.refvel), rate = -thi->erosion.rate * PetscPowScalar(magref2, 0.5 * thi->erosion.exponent);
444: if (erate) *erate = rate;
445: if (derate) {
446: if (thi->erosion.exponent == 1) {
447: derate->u = 0;
448: derate->v = 0;
449: } else {
450: derate->u = 0.5 * thi->erosion.exponent * rate / magref2 * 2. * vel->u / PetscSqr(thi->erosion.refvel);
451: derate->v = 0.5 * thi->erosion.exponent * rate / magref2 * 2. * vel->v / PetscSqr(thi->erosion.refvel);
452: }
453: }
454: }
456: static void RangeUpdate(PetscReal *min, PetscReal *max, PetscReal x)
457: {
458: if (x < *min) *min = x;
459: if (x > *max) *max = x;
460: }
462: static void PRangeClear(PRange *p)
463: {
464: p->cmin = p->min = 1e100;
465: p->cmax = p->max = -1e100;
466: }
468: static PetscErrorCode PRangeMinMax(PRange *p, PetscReal min, PetscReal max)
469: {
470: PetscFunctionBeginUser;
471: p->cmin = min;
472: p->cmax = max;
473: if (min < p->min) p->min = min;
474: if (max > p->max) p->max = max;
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode THIDestroy(THI *thi)
479: {
480: PetscFunctionBeginUser;
481: if (--((PetscObject)*thi)->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
482: PetscCall(PetscFree((*thi)->units));
483: PetscCall(PetscFree((*thi)->mattype));
484: PetscCall(PetscFree((*thi)->monitor_basename));
485: PetscCall(PetscHeaderDestroy(thi));
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode THICreate(MPI_Comm comm, THI *inthi)
490: {
491: static PetscBool registered = PETSC_FALSE;
492: THI thi;
493: Units units;
494: char monitor_basename[PETSC_MAX_PATH_LEN] = "thi-";
496: PetscFunctionBeginUser;
497: *inthi = 0;
498: if (!registered) {
499: PetscCall(PetscClassIdRegister("Toy Hydrostatic Ice", &THI_CLASSID));
500: registered = PETSC_TRUE;
501: }
502: PetscCall(PetscHeaderCreate(thi, THI_CLASSID, "THI", "Toy Hydrostatic Ice", "THI", comm, THIDestroy, 0));
504: PetscCall(PetscNew(&thi->units));
506: units = thi->units;
507: units->meter = 1e-2;
508: units->second = 1e-7;
509: units->kilogram = 1e-12;
511: PetscOptionsBegin(comm, NULL, "Scaled units options", "");
512: {
513: PetscCall(PetscOptionsReal("-units_meter", "1 meter in scaled length units", "", units->meter, &units->meter, NULL));
514: PetscCall(PetscOptionsReal("-units_second", "1 second in scaled time units", "", units->second, &units->second, NULL));
515: PetscCall(PetscOptionsReal("-units_kilogram", "1 kilogram in scaled mass units", "", units->kilogram, &units->kilogram, NULL));
516: }
517: PetscOptionsEnd();
518: units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second));
519: units->year = 31556926. * units->second, /* seconds per year */
521: thi->Lx = 10.e3;
522: thi->Ly = 10.e3;
523: thi->Lz = 1000;
524: thi->nlevels = 1;
525: thi->dirichlet_scale = 1;
526: thi->verbose = PETSC_FALSE;
528: thi->viscosity.glen_n = 3.;
529: thi->erosion.rate = 1e-3; /* m/a */
530: thi->erosion.exponent = 1.;
531: thi->erosion.refvel = 1.; /* m/a */
533: PetscOptionsBegin(comm, NULL, "Toy Hydrostatic Ice options", "");
534: {
535: QuadratureType quad = QUAD_GAUSS;
536: char homexp[] = "A";
537: char mtype[256] = MATSBAIJ;
538: PetscReal L, m = 1.0;
539: PetscBool flg;
540: L = thi->Lx;
541: PetscCall(PetscOptionsReal("-thi_L", "Domain size (m)", "", L, &L, &flg));
542: if (flg) thi->Lx = thi->Ly = L;
543: PetscCall(PetscOptionsReal("-thi_Lx", "X Domain size (m)", "", thi->Lx, &thi->Lx, NULL));
544: PetscCall(PetscOptionsReal("-thi_Ly", "Y Domain size (m)", "", thi->Ly, &thi->Ly, NULL));
545: PetscCall(PetscOptionsReal("-thi_Lz", "Z Domain size (m)", "", thi->Lz, &thi->Lz, NULL));
546: PetscCall(PetscOptionsString("-thi_hom", "ISMIP-HOM experiment (A or C)", "", homexp, homexp, sizeof(homexp), NULL));
547: switch (homexp[0] = toupper(homexp[0])) {
548: case 'A':
549: thi->initialize = THIInitialize_HOM_A;
550: thi->no_slip = PETSC_TRUE;
551: thi->alpha = 0.5;
552: break;
553: case 'C':
554: thi->initialize = THIInitialize_HOM_C;
555: thi->no_slip = PETSC_FALSE;
556: thi->alpha = 0.1;
557: break;
558: case 'F':
559: thi->initialize = THIInitialize_HOM_F;
560: thi->no_slip = PETSC_FALSE;
561: thi->alpha = 0.5;
562: break;
563: case 'X':
564: thi->initialize = THIInitialize_HOM_X;
565: thi->no_slip = PETSC_FALSE;
566: thi->alpha = 0.3;
567: break;
568: case 'Y':
569: thi->initialize = THIInitialize_HOM_Y;
570: thi->no_slip = PETSC_FALSE;
571: thi->alpha = 0.5;
572: break;
573: case 'Z':
574: thi->initialize = THIInitialize_HOM_Z;
575: thi->no_slip = PETSC_FALSE;
576: thi->alpha = 0.5;
577: break;
578: default:
579: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "HOM experiment '%c' not implemented", homexp[0]);
580: }
581: PetscCall(PetscOptionsEnum("-thi_quadrature", "Quadrature to use for 3D elements", "", QuadratureTypes, (PetscEnum)quad, (PetscEnum *)&quad, NULL));
582: switch (quad) {
583: case QUAD_GAUSS:
584: HexQInterp = HexQInterp_Gauss;
585: HexQDeriv = HexQDeriv_Gauss;
586: break;
587: case QUAD_LOBATTO:
588: HexQInterp = HexQInterp_Lobatto;
589: HexQDeriv = HexQDeriv_Lobatto;
590: break;
591: }
592: PetscCall(PetscOptionsReal("-thi_alpha", "Bed angle (degrees)", "", thi->alpha, &thi->alpha, NULL));
593: PetscCall(PetscOptionsReal("-thi_viscosity_glen_n", "Exponent in Glen flow law, 1=linear, infty=ideal plastic", NULL, thi->viscosity.glen_n, &thi->viscosity.glen_n, NULL));
594: PetscCall(PetscOptionsReal("-thi_friction_m", "Friction exponent, 0=Coulomb, 1=Navier", "", m, &m, NULL));
595: thi->friction.exponent = (m - 1) / 2;
596: PetscCall(PetscOptionsReal("-thi_erosion_rate", "Rate of erosion relative to sliding velocity at reference velocity (m/a)", NULL, thi->erosion.rate, &thi->erosion.rate, NULL));
597: PetscCall(PetscOptionsReal("-thi_erosion_exponent", "Power of sliding velocity appearing in erosion relation", NULL, thi->erosion.exponent, &thi->erosion.exponent, NULL));
598: PetscCall(PetscOptionsReal("-thi_erosion_refvel", "Reference sliding velocity for erosion (m/a)", NULL, thi->erosion.refvel, &thi->erosion.refvel, NULL));
599: thi->erosion.rate *= units->meter / units->year;
600: thi->erosion.refvel *= units->meter / units->year;
601: PetscCall(PetscOptionsReal("-thi_dirichlet_scale", "Scale Dirichlet boundary conditions by this factor", "", thi->dirichlet_scale, &thi->dirichlet_scale, NULL));
602: PetscCall(PetscOptionsReal("-thi_ssa_friction_scale", "Scale slip boundary conditions by this factor in SSA (2D) assembly", "", thi->ssa_friction_scale, &thi->ssa_friction_scale, NULL));
603: PetscCall(PetscOptionsReal("-thi_inertia", "Coefficient of acceleration term in velocity system, physical is almost zero", NULL, thi->inertia, &thi->inertia, NULL));
604: PetscCall(PetscOptionsInt("-thi_nlevels", "Number of levels of refinement", "", thi->nlevels, &thi->nlevels, NULL));
605: PetscCall(PetscOptionsFList("-thi_mat_type", "Matrix type", "MatSetType", MatList, mtype, (char *)mtype, sizeof(mtype), NULL));
606: PetscCall(PetscStrallocpy(mtype, &thi->mattype));
607: PetscCall(PetscOptionsBool("-thi_verbose", "Enable verbose output (like matrix sizes and statistics)", "", thi->verbose, &thi->verbose, NULL));
608: PetscCall(PetscOptionsString("-thi_monitor", "Basename to write state files to", NULL, monitor_basename, monitor_basename, sizeof(monitor_basename), &flg));
609: if (flg) {
610: PetscCall(PetscStrallocpy(monitor_basename, &thi->monitor_basename));
611: thi->monitor_interval = 1;
612: PetscCall(PetscOptionsInt("-thi_monitor_interval", "Frequency at which to write state files", NULL, thi->monitor_interval, &thi->monitor_interval, NULL));
613: }
614: }
615: PetscOptionsEnd();
617: /* dimensionalize */
618: thi->Lx *= units->meter;
619: thi->Ly *= units->meter;
620: thi->Lz *= units->meter;
621: thi->alpha *= PETSC_PI / 180;
623: PRangeClear(&thi->eta);
624: PRangeClear(&thi->beta2);
626: {
627: PetscReal u = 1000 * units->meter / (3e7 * units->second), gradu = u / (100 * units->meter), eta, deta, rho = 910 * units->kilogram / PetscPowRealInt(units->meter, 3), grav = 9.81 * units->meter / PetscSqr(units->second),
628: driving = rho * grav * PetscSinReal(thi->alpha) * 1000 * units->meter;
629: THIViscosity(thi, 0.5 * gradu * gradu, &eta, &deta);
630: thi->rhog = rho * grav;
631: if (thi->verbose) {
632: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Units: meter %8.2g second %8.2g kg %8.2g Pa %8.2g\n", (double)units->meter, (double)units->second, (double)units->kilogram, (double)units->Pascal));
633: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Domain (%6.2g,%6.2g,%6.2g), pressure %8.2g, driving stress %8.2g\n", (double)thi->Lx, (double)thi->Ly, (double)thi->Lz, (double)(rho * grav * 1e3 * units->meter), (double)driving));
634: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Large velocity 1km/a %8.2g, velocity gradient %8.2g, eta %8.2g, stress %8.2g, ratio %8.2g\n", (double)u, (double)gradu, (double)eta, (double)(2 * eta * gradu, 2 * eta * gradu / driving)));
635: THIViscosity(thi, 0.5 * PetscSqr(1e-3 * gradu), &eta, &deta);
636: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Small velocity 1m/a %8.2g, velocity gradient %8.2g, eta %8.2g, stress %8.2g, ratio %8.2g\n", (double)(1e-3 * u), (double)(1e-3 * gradu), (double)eta, (double)(2 * eta * 1e-3 * gradu, 2 * eta * 1e-3 * gradu / driving)));
637: }
638: }
640: *inthi = thi;
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: /* Our problem is periodic, but the domain has a mean slope of alpha so the bed does not line up between the upstream
645: * and downstream ends of the domain. This function fixes the ghost values so that the domain appears truly periodic in
646: * the horizontal. */
647: static PetscErrorCode THIFixGhosts(THI thi, DM da3, DM da2, Vec X3, Vec X2)
648: {
649: DMDALocalInfo info;
650: PrmNode **x2;
651: PetscInt i, j;
653: PetscFunctionBeginUser;
654: PetscCall(DMDAGetLocalInfo(da3, &info));
655: /* PetscCall(VecView(X2,PETSC_VIEWER_STDOUT_WORLD)); */
656: PetscCall(DMDAVecGetArray(da2, X2, &x2));
657: for (i = info.gzs; i < info.gzs + info.gzm; i++) {
658: if (i > -1 && i < info.mz) continue;
659: for (j = info.gys; j < info.gys + info.gym; j++) x2[i][j].b += PetscSinReal(thi->alpha) * thi->Lx * (i < 0 ? 1.0 : -1.0);
660: }
661: PetscCall(DMDAVecRestoreArray(da2, X2, &x2));
662: /* PetscCall(VecView(X2,PETSC_VIEWER_STDOUT_WORLD)); */
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode THIInitializePrm(THI thi, DM da2prm, PrmNode **p)
667: {
668: PetscInt i, j, xs, xm, ys, ym, mx, my;
670: PetscFunctionBeginUser;
671: PetscCall(DMDAGetGhostCorners(da2prm, &ys, &xs, 0, &ym, &xm, 0));
672: PetscCall(DMDAGetInfo(da2prm, 0, &my, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
673: for (i = xs; i < xs + xm; i++) {
674: for (j = ys; j < ys + ym; j++) {
675: PetscReal xx = thi->Lx * i / mx, yy = thi->Ly * j / my;
676: thi->initialize(thi, xx, yy, &p[i][j]);
677: }
678: }
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode THIInitial(THI thi, DM pack, Vec X)
683: {
684: DM da3, da2;
685: PetscInt i, j, k, xs, xm, ys, ym, zs, zm, mx, my;
686: PetscReal hx, hy;
687: PrmNode **prm;
688: Node ***x;
689: Vec X3g, X2g, X2;
691: PetscFunctionBeginUser;
692: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
693: PetscCall(DMCompositeGetAccess(pack, X, &X3g, &X2g));
694: PetscCall(DMGetLocalVector(da2, &X2));
696: PetscCall(DMDAGetInfo(da3, 0, 0, &my, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0));
697: PetscCall(DMDAGetCorners(da3, &zs, &ys, &xs, &zm, &ym, &xm));
698: PetscCall(DMDAVecGetArray(da3, X3g, &x));
699: PetscCall(DMDAVecGetArray(da2, X2, &prm));
701: PetscCall(THIInitializePrm(thi, da2, prm));
703: hx = thi->Lx / mx;
704: hy = thi->Ly / my;
705: for (i = xs; i < xs + xm; i++) {
706: for (j = ys; j < ys + ym; j++) {
707: for (k = zs; k < zs + zm; k++) {
708: const PetscScalar zm1 = zm - 1, drivingx = thi->rhog * (prm[i + 1][j].b + prm[i + 1][j].h - prm[i - 1][j].b - prm[i - 1][j].h) / (2 * hx), drivingy = thi->rhog * (prm[i][j + 1].b + prm[i][j + 1].h - prm[i][j - 1].b - prm[i][j - 1].h) / (2 * hy);
709: x[i][j][k].u = 0. * drivingx * prm[i][j].h * (PetscScalar)k / zm1;
710: x[i][j][k].v = 0. * drivingy * prm[i][j].h * (PetscScalar)k / zm1;
711: }
712: }
713: }
715: PetscCall(DMDAVecRestoreArray(da3, X3g, &x));
716: PetscCall(DMDAVecRestoreArray(da2, X2, &prm));
718: PetscCall(DMLocalToGlobalBegin(da2, X2, INSERT_VALUES, X2g));
719: PetscCall(DMLocalToGlobalEnd(da2, X2, INSERT_VALUES, X2g));
720: PetscCall(DMRestoreLocalVector(da2, &X2));
722: PetscCall(DMCompositeRestoreAccess(pack, X, &X3g, &X2g));
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: static void PointwiseNonlinearity(THI thi, const Node n[restrict 8], const PetscReal phi[restrict 3], PetscReal dphi[restrict 8][3], PetscScalar *restrict u, PetscScalar *restrict v, PetscScalar du[restrict 3], PetscScalar dv[restrict 3], PetscReal *eta, PetscReal *deta)
727: {
728: PetscInt l;
729: PetscScalar gam;
731: du[0] = du[1] = du[2] = 0;
732: dv[0] = dv[1] = dv[2] = 0;
733: *u = 0;
734: *v = 0;
735: for (l = 0; l < 8; l++) {
736: *u += phi[l] * n[l].u;
737: *v += phi[l] * n[l].v;
738: for (PetscInt ll = 0; ll < 3; ll++) {
739: du[ll] += dphi[l][ll] * n[l].u;
740: dv[ll] += dphi[l][ll] * n[l].v;
741: }
742: }
743: gam = Sqr(du[0]) + Sqr(dv[1]) + du[0] * dv[1] + 0.25 * Sqr(du[1] + dv[0]) + 0.25 * Sqr(du[2]) + 0.25 * Sqr(dv[2]);
744: THIViscosity(thi, PetscRealPart(gam), eta, deta);
745: }
747: static PetscErrorCode THIFunctionLocal_3D(DMDALocalInfo *info, const Node ***x, const PrmNode **prm, const Node ***xdot, Node ***f, THI thi)
748: {
749: PetscInt xs, ys, xm, ym, zm, i, j, k, q, l;
750: PetscReal hx, hy, etamin, etamax, beta2min, beta2max;
752: PetscFunctionBeginUser;
753: xs = info->zs;
754: ys = info->ys;
755: xm = info->zm;
756: ym = info->ym;
757: zm = info->xm;
758: hx = thi->Lx / info->mz;
759: hy = thi->Ly / info->my;
761: etamin = 1e100;
762: etamax = 0;
763: beta2min = 1e100;
764: beta2max = 0;
766: for (i = xs; i < xs + xm; i++) {
767: for (j = ys; j < ys + ym; j++) {
768: PrmNode pn[4], dpn[4][2];
769: QuadExtract(prm, i, j, pn);
770: PetscCall(QuadComputeGrad4(QuadQDeriv, hx, hy, pn, dpn));
771: for (k = 0; k < zm - 1; k++) {
772: PetscInt ls = 0;
773: Node n[8], ndot[8], *fn[8];
774: PetscReal zn[8], etabase = 0;
776: PrmHexGetZ(pn, k, zm, zn);
777: HexExtract(x, i, j, k, n);
778: HexExtract(xdot, i, j, k, ndot);
779: HexExtractRef(f, i, j, k, fn);
780: if (thi->no_slip && k == 0) {
781: for (l = 0; l < 4; l++) n[l].u = n[l].v = 0;
782: /* The first 4 basis functions lie on the bottom layer, so their contribution is exactly 0, hence we can skip them */
783: ls = 4;
784: }
785: for (q = 0; q < 8; q++) {
786: PetscReal dz[3], phi[8], dphi[8][3], jw, eta, deta;
787: PetscScalar du[3], dv[3], u, v, udot = 0, vdot = 0;
788: for (l = ls; l < 8; l++) {
789: udot += HexQInterp[q][l] * ndot[l].u;
790: vdot += HexQInterp[q][l] * ndot[l].v;
791: }
792: HexGrad(HexQDeriv[q], zn, dz);
793: HexComputeGeometry(q, hx, hy, dz, phi, dphi, &jw);
794: PointwiseNonlinearity(thi, n, phi, dphi, &u, &v, du, dv, &eta, &deta);
795: jw /= thi->rhog; /* scales residuals to be O(1) */
796: if (q == 0) etabase = eta;
797: RangeUpdate(&etamin, &etamax, eta);
798: for (l = ls; l < 8; l++) { /* test functions */
799: const PetscScalar ds[2] = {dpn[q % 4][0].h + dpn[q % 4][0].b, dpn[q % 4][1].h + dpn[q % 4][1].b};
800: const PetscReal pp = phi[l], *dp = dphi[l];
801: fn[l]->u += dp[0] * jw * eta * (4. * du[0] + 2. * dv[1]) + dp[1] * jw * eta * (du[1] + dv[0]) + dp[2] * jw * eta * du[2] + pp * jw * thi->rhog * ds[0];
802: fn[l]->v += dp[1] * jw * eta * (2. * du[0] + 4. * dv[1]) + dp[0] * jw * eta * (du[1] + dv[0]) + dp[2] * jw * eta * dv[2] + pp * jw * thi->rhog * ds[1];
803: fn[l]->u += pp * jw * udot * thi->inertia * pp;
804: fn[l]->v += pp * jw * vdot * thi->inertia * pp;
805: }
806: }
807: if (k == 0) { /* we are on a bottom face */
808: if (thi->no_slip) {
809: /* Note: Non-Galerkin coarse grid operators are very sensitive to the scaling of Dirichlet boundary
810: * conditions. After shenanigans above, etabase contains the effective viscosity at the closest quadrature
811: * point to the bed. We want the diagonal entry in the Dirichlet condition to have similar magnitude to the
812: * diagonal entry corresponding to the adjacent node. The fundamental scaling of the viscous part is in
813: * diagu, diagv below. This scaling is easy to recognize by considering the finite difference operator after
814: * scaling by element size. The no-slip Dirichlet condition is scaled by this factor, and also in the
815: * assembled matrix (see the similar block in THIJacobianLocal).
816: *
817: * Note that the residual at this Dirichlet node is linear in the state at this node, but also depends
818: * (nonlinearly in general) on the neighboring interior nodes through the local viscosity. This will make
819: * a matrix-free Jacobian have extra entries in the corresponding row. We assemble only the diagonal part,
820: * so the solution will exactly satisfy the boundary condition after the first linear iteration.
821: */
822: const PetscReal hz = PetscRealPart(pn[0].h) / (zm - 1.);
823: const PetscScalar diagu = 2 * etabase / thi->rhog * (hx * hy / hz + hx * hz / hy + 4 * hy * hz / hx), diagv = 2 * etabase / thi->rhog * (hx * hy / hz + 4 * hx * hz / hy + hy * hz / hx);
824: fn[0]->u = thi->dirichlet_scale * diagu * x[i][j][k].u;
825: fn[0]->v = thi->dirichlet_scale * diagv * x[i][j][k].v;
826: } else { /* Integrate over bottom face to apply boundary condition */
827: for (q = 0; q < 4; q++) { /* We remove the explicit scaling of the residual by 1/rhog because beta2 already has that scaling to be O(1) */
828: const PetscReal jw = 0.25 * hx * hy, *phi = QuadQInterp[q];
829: PetscScalar u = 0, v = 0, rbeta2 = 0;
830: PetscReal beta2, dbeta2;
831: for (l = 0; l < 4; l++) {
832: u += phi[l] * n[l].u;
833: v += phi[l] * n[l].v;
834: rbeta2 += phi[l] * pn[l].beta2;
835: }
836: THIFriction(thi, PetscRealPart(rbeta2), PetscRealPart(u * u + v * v) / 2, &beta2, &dbeta2);
837: RangeUpdate(&beta2min, &beta2max, beta2);
838: for (l = 0; l < 4; l++) {
839: const PetscReal pp = phi[l];
840: fn[ls + l]->u += pp * jw * beta2 * u;
841: fn[ls + l]->v += pp * jw * beta2 * v;
842: }
843: }
844: }
845: }
846: }
847: }
848: }
850: PetscCall(PRangeMinMax(&thi->eta, etamin, etamax));
851: PetscCall(PRangeMinMax(&thi->beta2, beta2min, beta2max));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode THIFunctionLocal_2D(DMDALocalInfo *info, const Node ***x, const PrmNode **prm, const PrmNode **prmdot, PrmNode **f, THI thi)
856: {
857: PetscInt xs, ys, xm, ym, zm, i, j, k;
859: PetscFunctionBeginUser;
860: xs = info->zs;
861: ys = info->ys;
862: xm = info->zm;
863: ym = info->ym;
864: zm = info->xm;
866: for (i = xs; i < xs + xm; i++) {
867: for (j = ys; j < ys + ym; j++) {
868: PetscScalar div = 0, erate, h[8];
869: PrmNodeGetFaceMeasure(prm, i, j, h);
870: for (k = 0; k < zm; k++) {
871: PetscScalar weight = (k == 0 || k == zm - 1) ? 0.5 / (zm - 1) : 1.0 / (zm - 1);
872: if (0) { /* centered flux */
873: div += (-weight * h[0] * StaggeredMidpoint2D(x[i][j][k].u, x[i - 1][j][k].u, x[i - 1][j - 1][k].u, x[i][j - 1][k].u) - weight * h[1] * StaggeredMidpoint2D(x[i][j][k].u, x[i - 1][j][k].u, x[i - 1][j + 1][k].u, x[i][j + 1][k].u) +
874: weight * h[2] * StaggeredMidpoint2D(x[i][j][k].u, x[i + 1][j][k].u, x[i + 1][j + 1][k].u, x[i][j + 1][k].u) + weight * h[3] * StaggeredMidpoint2D(x[i][j][k].u, x[i + 1][j][k].u, x[i + 1][j - 1][k].u, x[i][j - 1][k].u) -
875: weight * h[4] * StaggeredMidpoint2D(x[i][j][k].v, x[i][j - 1][k].v, x[i + 1][j - 1][k].v, x[i + 1][j][k].v) - weight * h[5] * StaggeredMidpoint2D(x[i][j][k].v, x[i][j - 1][k].v, x[i - 1][j - 1][k].v, x[i - 1][j][k].v) +
876: weight * h[6] * StaggeredMidpoint2D(x[i][j][k].v, x[i][j + 1][k].v, x[i - 1][j + 1][k].v, x[i - 1][j][k].v) + weight * h[7] * StaggeredMidpoint2D(x[i][j][k].v, x[i][j + 1][k].v, x[i + 1][j + 1][k].v, x[i + 1][j][k].v));
877: } else { /* Upwind flux */
878: div += weight * (-UpwindFluxXW(x, prm, h, i, j, k, 1) - UpwindFluxXW(x, prm, h, i, j, k, -1) + UpwindFluxXE(x, prm, h, i, j, k, 1) + UpwindFluxXE(x, prm, h, i, j, k, -1) - UpwindFluxYS(x, prm, h, i, j, k, 1) - UpwindFluxYS(x, prm, h, i, j, k, -1) + UpwindFluxYN(x, prm, h, i, j, k, 1) + UpwindFluxYN(x, prm, h, i, j, k, -1));
879: }
880: }
881: /* printf("div[%d][%d] %g\n",i,j,div); */
882: THIErosion(thi, &x[i][j][0], &erate, NULL);
883: f[i][j].b = prmdot[i][j].b - erate;
884: f[i][j].h = prmdot[i][j].h + div;
885: f[i][j].beta2 = prmdot[i][j].beta2;
886: }
887: }
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: static PetscErrorCode THIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
892: {
893: THI thi = (THI)ctx;
894: DM pack, da3, da2;
895: Vec X3, X2, Xdot3, Xdot2, F3, F2, F3g, F2g;
896: const Node ***x3, ***xdot3;
897: const PrmNode **x2, **xdot2;
898: Node ***f3;
899: PrmNode **f2;
900: DMDALocalInfo info3;
902: PetscFunctionBeginUser;
903: PetscCall(TSGetDM(ts, &pack));
904: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
905: PetscCall(DMDAGetLocalInfo(da3, &info3));
906: PetscCall(DMCompositeGetLocalVectors(pack, &X3, &X2));
907: PetscCall(DMCompositeGetLocalVectors(pack, &Xdot3, &Xdot2));
908: PetscCall(DMCompositeScatter(pack, X, X3, X2));
909: PetscCall(THIFixGhosts(thi, da3, da2, X3, X2));
910: PetscCall(DMCompositeScatter(pack, Xdot, Xdot3, Xdot2));
912: PetscCall(DMGetLocalVector(da3, &F3));
913: PetscCall(DMGetLocalVector(da2, &F2));
914: PetscCall(VecZeroEntries(F3));
916: PetscCall(DMDAVecGetArray(da3, X3, &x3));
917: PetscCall(DMDAVecGetArray(da2, X2, &x2));
918: PetscCall(DMDAVecGetArray(da3, Xdot3, &xdot3));
919: PetscCall(DMDAVecGetArray(da2, Xdot2, &xdot2));
920: PetscCall(DMDAVecGetArray(da3, F3, &f3));
921: PetscCall(DMDAVecGetArray(da2, F2, &f2));
923: PetscCall(THIFunctionLocal_3D(&info3, x3, x2, xdot3, f3, thi));
924: PetscCall(THIFunctionLocal_2D(&info3, x3, x2, xdot2, f2, thi));
926: PetscCall(DMDAVecRestoreArray(da3, X3, &x3));
927: PetscCall(DMDAVecRestoreArray(da2, X2, &x2));
928: PetscCall(DMDAVecRestoreArray(da3, Xdot3, &xdot3));
929: PetscCall(DMDAVecRestoreArray(da2, Xdot2, &xdot2));
930: PetscCall(DMDAVecRestoreArray(da3, F3, &f3));
931: PetscCall(DMDAVecRestoreArray(da2, F2, &f2));
933: PetscCall(DMCompositeRestoreLocalVectors(pack, &X3, &X2));
934: PetscCall(DMCompositeRestoreLocalVectors(pack, &Xdot3, &Xdot2));
936: PetscCall(VecZeroEntries(F));
937: PetscCall(DMCompositeGetAccess(pack, F, &F3g, &F2g));
938: PetscCall(DMLocalToGlobalBegin(da3, F3, ADD_VALUES, F3g));
939: PetscCall(DMLocalToGlobalEnd(da3, F3, ADD_VALUES, F3g));
940: PetscCall(DMLocalToGlobalBegin(da2, F2, INSERT_VALUES, F2g));
941: PetscCall(DMLocalToGlobalEnd(da2, F2, INSERT_VALUES, F2g));
943: if (thi->verbose) {
944: PetscViewer viewer;
945: PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)thi), &viewer));
946: PetscCall(PetscViewerASCIIPrintf(viewer, "3D_Velocity residual (bs=2):\n"));
947: PetscCall(PetscViewerASCIIPushTab(viewer));
948: PetscCall(VecView(F3, viewer));
949: PetscCall(PetscViewerASCIIPopTab(viewer));
950: PetscCall(PetscViewerASCIIPrintf(viewer, "2D_Fields residual (bs=3):\n"));
951: PetscCall(PetscViewerASCIIPushTab(viewer));
952: PetscCall(VecView(F2, viewer));
953: PetscCall(PetscViewerASCIIPopTab(viewer));
954: }
956: PetscCall(DMCompositeRestoreAccess(pack, F, &F3g, &F2g));
958: PetscCall(DMRestoreLocalVector(da3, &F3));
959: PetscCall(DMRestoreLocalVector(da2, &F2));
960: PetscFunctionReturn(PETSC_SUCCESS);
961: }
963: static PetscErrorCode THIMatrixStatistics(THI thi, Mat B, PetscViewer viewer)
964: {
965: PetscReal nrm;
966: PetscInt m;
967: PetscMPIInt rank;
969: PetscFunctionBeginUser;
970: PetscCall(MatNorm(B, NORM_FROBENIUS, &nrm));
971: PetscCall(MatGetSize(B, &m, 0));
972: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)B), &rank));
973: if (rank == 0) {
974: PetscScalar val0, val2;
975: PetscCall(MatGetValue(B, 0, 0, &val0));
976: PetscCall(MatGetValue(B, 2, 2, &val2));
977: PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix dim %8" PetscInt_FMT " norm %8.2e, (0,0) %8.2e (2,2) %8.2e, eta [%8.2e,%8.2e] beta2 [%8.2e,%8.2e]\n", m, (double)nrm, (double)PetscRealPart(val0), (double)PetscRealPart(val2), (double)thi->eta.cmin,
978: (double)thi->eta.cmax, (double)thi->beta2.cmin, (double)thi->beta2.cmax));
979: }
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: static PetscErrorCode THISurfaceStatistics(DM pack, Vec X, PetscReal *min, PetscReal *max, PetscReal *mean)
984: {
985: DM da3, da2;
986: Vec X3, X2;
987: Node ***x;
988: PetscInt i, j, xs, ys, zs, xm, ym, zm, mx, my, mz;
989: PetscScalar gusum = 0.0;
991: PetscFunctionBeginUser;
992: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
993: PetscCall(DMCompositeGetAccess(pack, X, &X3, &X2));
994: *min = 1e100;
995: *max = -1e100;
996: *mean = 0;
997: PetscCall(DMDAGetInfo(da3, 0, &mz, &my, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0));
998: PetscCall(DMDAGetCorners(da3, &zs, &ys, &xs, &zm, &ym, &xm));
999: PetscCheck(zs == 0 && zm == mz, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected decomposition");
1000: PetscCall(DMDAVecGetArray(da3, X3, &x));
1001: for (i = xs; i < xs + xm; i++) {
1002: for (j = ys; j < ys + ym; j++) {
1003: PetscReal u = PetscRealPart(x[i][j][zm - 1].u);
1004: RangeUpdate(min, max, u);
1005: gusum += u;
1006: }
1007: }
1008: PetscCall(DMDAVecRestoreArray(da3, X3, &x));
1009: PetscCall(DMCompositeRestoreAccess(pack, X, &X3, &X2));
1011: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, min, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)da3)));
1012: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, max, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da3)));
1013: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gusum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da3)));
1014: *mean = PetscRealPart(gusum) / (mx * my);
1015: PetscFunctionReturn(PETSC_SUCCESS);
1016: }
1018: static PetscErrorCode THISolveStatistics(THI thi, TS ts, PetscInt coarsened, const char name[])
1019: {
1020: MPI_Comm comm;
1021: DM pack;
1022: Vec X, X3, X2;
1024: PetscFunctionBeginUser;
1025: PetscCall(PetscObjectGetComm((PetscObject)thi, &comm));
1026: PetscCall(TSGetDM(ts, &pack));
1027: PetscCall(TSGetSolution(ts, &X));
1028: PetscCall(DMCompositeGetAccess(pack, X, &X3, &X2));
1029: PetscCall(PetscPrintf(comm, "Solution statistics after solve: %s\n", name));
1030: {
1031: PetscInt its, lits;
1032: SNESConvergedReason reason;
1033: SNES snes;
1034: PetscCall(TSGetSNES(ts, &snes));
1035: PetscCall(SNESGetIterationNumber(snes, &its));
1036: PetscCall(SNESGetConvergedReason(snes, &reason));
1037: PetscCall(SNESGetLinearSolveIterations(snes, &lits));
1038: PetscCall(PetscPrintf(comm, "%s: Number of SNES iterations = %" PetscInt_FMT ", total linear iterations = %" PetscInt_FMT "\n", SNESConvergedReasons[reason], its, lits));
1039: }
1040: {
1041: PetscReal nrm2, min[3] = {1e100, 1e100, 1e100}, max[3] = {-1e100, -1e100, -1e100};
1042: PetscInt i, j, m;
1043: PetscScalar *x;
1044: PetscCall(VecNorm(X3, NORM_2, &nrm2));
1045: PetscCall(VecGetLocalSize(X3, &m));
1046: PetscCall(VecGetArray(X3, &x));
1047: for (i = 0; i < m; i += 2) {
1048: PetscReal u = PetscRealPart(x[i]), v = PetscRealPart(x[i + 1]), c = PetscSqrtReal(u * u + v * v);
1049: min[0] = PetscMin(u, min[0]);
1050: min[1] = PetscMin(v, min[1]);
1051: min[2] = PetscMin(c, min[2]);
1052: max[0] = PetscMax(u, max[0]);
1053: max[1] = PetscMax(v, max[1]);
1054: max[2] = PetscMax(c, max[2]);
1055: }
1056: PetscCall(VecRestoreArray(X, &x));
1057: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, min, 3, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)thi)));
1058: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, max, 3, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)thi)));
1059: /* Dimensionalize to meters/year */
1060: nrm2 *= thi->units->year / thi->units->meter;
1061: for (j = 0; j < 3; j++) {
1062: min[j] *= thi->units->year / thi->units->meter;
1063: max[j] *= thi->units->year / thi->units->meter;
1064: }
1065: PetscCall(PetscPrintf(comm, "|X|_2 %g u in [%g, %g] v in [%g, %g] c in [%g, %g] \n", (double)nrm2, (double)min[0], (double)max[0], (double)min[1], (double)max[1], (double)min[2], (double)max[2]));
1066: {
1067: PetscReal umin, umax, umean;
1068: PetscCall(THISurfaceStatistics(pack, X, &umin, &umax, &umean));
1069: umin *= thi->units->year / thi->units->meter;
1070: umax *= thi->units->year / thi->units->meter;
1071: umean *= thi->units->year / thi->units->meter;
1072: PetscCall(PetscPrintf(comm, "Surface statistics: u in [%12.6e, %12.6e] mean %12.6e\n", (double)umin, (double)umax, (double)umean));
1073: }
1074: /* These values stay nondimensional */
1075: PetscCall(PetscPrintf(comm, "Global eta range [%g, %g], converged range [%g, %g]\n", (double)thi->eta.min, (double)thi->eta.max, (double)thi->eta.cmin, (double)thi->eta.cmax));
1076: PetscCall(PetscPrintf(comm, "Global beta2 range [%g, %g], converged range [%g, %g]\n", (double)thi->beta2.min, (double)thi->beta2.max, (double)thi->beta2.cmin, (double)thi->beta2.cmax));
1077: }
1078: PetscCall(PetscPrintf(comm, "\n"));
1079: PetscCall(DMCompositeRestoreAccess(pack, X, &X3, &X2));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: static inline PetscInt DMDALocalIndex3D(DMDALocalInfo *info, PetscInt i, PetscInt j, PetscInt k)
1084: {
1085: return ((i - info->gzs) * info->gym + (j - info->gys)) * info->gxm + (k - info->gxs);
1086: }
1087: static inline PetscInt DMDALocalIndex2D(DMDALocalInfo *info, PetscInt i, PetscInt j)
1088: {
1089: return (i - info->gzs) * info->gym + (j - info->gys);
1090: }
1092: static PetscErrorCode THIJacobianLocal_Momentum(DMDALocalInfo *info, const Node ***x, const PrmNode **prm, Mat B, Mat Bcpl, THI thi)
1093: {
1094: PetscInt xs, ys, xm, ym, zm, i, j, k, q, l, ll;
1095: PetscReal hx, hy;
1097: PetscFunctionBeginUser;
1098: xs = info->zs;
1099: ys = info->ys;
1100: xm = info->zm;
1101: ym = info->ym;
1102: zm = info->xm;
1103: hx = thi->Lx / info->mz;
1104: hy = thi->Ly / info->my;
1106: for (i = xs; i < xs + xm; i++) {
1107: for (j = ys; j < ys + ym; j++) {
1108: PrmNode pn[4], dpn[4][2];
1109: QuadExtract(prm, i, j, pn);
1110: PetscCall(QuadComputeGrad4(QuadQDeriv, hx, hy, pn, dpn));
1111: for (k = 0; k < zm - 1; k++) {
1112: Node n[8];
1113: PetscReal zn[8], etabase = 0;
1114: PetscScalar Ke[8 * NODE_SIZE][8 * NODE_SIZE], Kcpl[8 * NODE_SIZE][4 * PRMNODE_SIZE];
1115: PetscInt ls = 0;
1117: PrmHexGetZ(pn, k, zm, zn);
1118: HexExtract(x, i, j, k, n);
1119: PetscCall(PetscMemzero(Ke, sizeof(Ke)));
1120: PetscCall(PetscMemzero(Kcpl, sizeof(Kcpl)));
1121: if (thi->no_slip && k == 0) {
1122: for (l = 0; l < 4; l++) n[l].u = n[l].v = 0;
1123: ls = 4;
1124: }
1125: for (q = 0; q < 8; q++) {
1126: PetscReal dz[3], phi[8], dphi[8][3], jw, eta, deta;
1127: PetscScalar du[3], dv[3], u, v;
1128: HexGrad(HexQDeriv[q], zn, dz);
1129: HexComputeGeometry(q, hx, hy, dz, phi, dphi, &jw);
1130: PointwiseNonlinearity(thi, n, phi, dphi, &u, &v, du, dv, &eta, &deta);
1131: jw /= thi->rhog; /* residuals are scaled by this factor */
1132: if (q == 0) etabase = eta;
1133: for (l = ls; l < 8; l++) { /* test functions */
1134: const PetscReal pp = phi[l], *restrict dp = dphi[l];
1135: for (ll = ls; ll < 8; ll++) { /* trial functions */
1136: const PetscReal *restrict dpl = dphi[ll];
1137: PetscScalar dgdu, dgdv;
1138: dgdu = 2. * du[0] * dpl[0] + dv[1] * dpl[0] + 0.5 * (du[1] + dv[0]) * dpl[1] + 0.5 * du[2] * dpl[2];
1139: dgdv = 2. * dv[1] * dpl[1] + du[0] * dpl[1] + 0.5 * (du[1] + dv[0]) * dpl[0] + 0.5 * dv[2] * dpl[2];
1140: /* Picard part */
1141: Ke[l * 2 + 0][ll * 2 + 0] += dp[0] * jw * eta * 4. * dpl[0] + dp[1] * jw * eta * dpl[1] + dp[2] * jw * eta * dpl[2];
1142: Ke[l * 2 + 0][ll * 2 + 1] += dp[0] * jw * eta * 2. * dpl[1] + dp[1] * jw * eta * dpl[0];
1143: Ke[l * 2 + 1][ll * 2 + 0] += dp[1] * jw * eta * 2. * dpl[0] + dp[0] * jw * eta * dpl[1];
1144: Ke[l * 2 + 1][ll * 2 + 1] += dp[1] * jw * eta * 4. * dpl[1] + dp[0] * jw * eta * dpl[0] + dp[2] * jw * eta * dpl[2];
1145: /* extra Newton terms */
1146: Ke[l * 2 + 0][ll * 2 + 0] += dp[0] * jw * deta * dgdu * (4. * du[0] + 2. * dv[1]) + dp[1] * jw * deta * dgdu * (du[1] + dv[0]) + dp[2] * jw * deta * dgdu * du[2];
1147: Ke[l * 2 + 0][ll * 2 + 1] += dp[0] * jw * deta * dgdv * (4. * du[0] + 2. * dv[1]) + dp[1] * jw * deta * dgdv * (du[1] + dv[0]) + dp[2] * jw * deta * dgdv * du[2];
1148: Ke[l * 2 + 1][ll * 2 + 0] += dp[1] * jw * deta * dgdu * (4. * dv[1] + 2. * du[0]) + dp[0] * jw * deta * dgdu * (du[1] + dv[0]) + dp[2] * jw * deta * dgdu * dv[2];
1149: Ke[l * 2 + 1][ll * 2 + 1] += dp[1] * jw * deta * dgdv * (4. * dv[1] + 2. * du[0]) + dp[0] * jw * deta * dgdv * (du[1] + dv[0]) + dp[2] * jw * deta * dgdv * dv[2];
1150: /* inertial part */
1151: Ke[l * 2 + 0][ll * 2 + 0] += pp * jw * thi->inertia * pp;
1152: Ke[l * 2 + 1][ll * 2 + 1] += pp * jw * thi->inertia * pp;
1153: }
1154: for (ll = 0; ll < 4; ll++) { /* Trial functions for surface/bed */
1155: const PetscReal dpl[] = {QuadQDeriv[q % 4][ll][0] / hx, QuadQDeriv[q % 4][ll][1] / hy}; /* surface = h + b */
1156: Kcpl[FieldIndex(Node, l, u)][FieldIndex(PrmNode, ll, h)] += pp * jw * thi->rhog * dpl[0];
1157: Kcpl[FieldIndex(Node, l, u)][FieldIndex(PrmNode, ll, b)] += pp * jw * thi->rhog * dpl[0];
1158: Kcpl[FieldIndex(Node, l, v)][FieldIndex(PrmNode, ll, h)] += pp * jw * thi->rhog * dpl[1];
1159: Kcpl[FieldIndex(Node, l, v)][FieldIndex(PrmNode, ll, b)] += pp * jw * thi->rhog * dpl[1];
1160: }
1161: }
1162: }
1163: if (k == 0) { /* on a bottom face */
1164: if (thi->no_slip) {
1165: const PetscReal hz = PetscRealPart(pn[0].h) / (zm - 1);
1166: const PetscScalar diagu = 2 * etabase / thi->rhog * (hx * hy / hz + hx * hz / hy + 4 * hy * hz / hx), diagv = 2 * etabase / thi->rhog * (hx * hy / hz + 4 * hx * hz / hy + hy * hz / hx);
1167: Ke[0][0] = thi->dirichlet_scale * diagu;
1168: Ke[0][1] = 0;
1169: Ke[1][0] = 0;
1170: Ke[1][1] = thi->dirichlet_scale * diagv;
1171: } else {
1172: for (q = 0; q < 4; q++) { /* We remove the explicit scaling by 1/rhog because beta2 already has that scaling to be O(1) */
1173: const PetscReal jw = 0.25 * hx * hy, *phi = QuadQInterp[q];
1174: PetscScalar u = 0, v = 0, rbeta2 = 0;
1175: PetscReal beta2, dbeta2;
1176: for (l = 0; l < 4; l++) {
1177: u += phi[l] * n[l].u;
1178: v += phi[l] * n[l].v;
1179: rbeta2 += phi[l] * pn[l].beta2;
1180: }
1181: THIFriction(thi, PetscRealPart(rbeta2), PetscRealPart(u * u + v * v) / 2, &beta2, &dbeta2);
1182: for (l = 0; l < 4; l++) {
1183: const PetscReal pp = phi[l];
1184: for (ll = 0; ll < 4; ll++) {
1185: const PetscReal ppl = phi[ll];
1186: Ke[l * 2 + 0][ll * 2 + 0] += pp * jw * beta2 * ppl + pp * jw * dbeta2 * u * u * ppl;
1187: Ke[l * 2 + 0][ll * 2 + 1] += pp * jw * dbeta2 * u * v * ppl;
1188: Ke[l * 2 + 1][ll * 2 + 0] += pp * jw * dbeta2 * v * u * ppl;
1189: Ke[l * 2 + 1][ll * 2 + 1] += pp * jw * beta2 * ppl + pp * jw * dbeta2 * v * v * ppl;
1190: }
1191: }
1192: }
1193: }
1194: }
1195: {
1196: const PetscInt rc3blocked[8] = {DMDALocalIndex3D(info, i + 0, j + 0, k + 0), DMDALocalIndex3D(info, i + 1, j + 0, k + 0), DMDALocalIndex3D(info, i + 1, j + 1, k + 0), DMDALocalIndex3D(info, i + 0, j + 1, k + 0),
1197: DMDALocalIndex3D(info, i + 0, j + 0, k + 1), DMDALocalIndex3D(info, i + 1, j + 0, k + 1), DMDALocalIndex3D(info, i + 1, j + 1, k + 1), DMDALocalIndex3D(info, i + 0, j + 1, k + 1)},
1198: col2blocked[PRMNODE_SIZE * 4] = {DMDALocalIndex2D(info, i + 0, j + 0), DMDALocalIndex2D(info, i + 1, j + 0), DMDALocalIndex2D(info, i + 1, j + 1), DMDALocalIndex2D(info, i + 0, j + 1)};
1199: #if !defined(COMPUTE_LOWER_TRIANGULAR) /* fill in lower-triangular part, this is really cheap compared to computing the entries */
1200: for (l = 0; l < 8; l++) {
1201: for (ll = l + 1; ll < 8; ll++) {
1202: Ke[ll * 2 + 0][l * 2 + 0] = Ke[l * 2 + 0][ll * 2 + 0];
1203: Ke[ll * 2 + 1][l * 2 + 0] = Ke[l * 2 + 0][ll * 2 + 1];
1204: Ke[ll * 2 + 0][l * 2 + 1] = Ke[l * 2 + 1][ll * 2 + 0];
1205: Ke[ll * 2 + 1][l * 2 + 1] = Ke[l * 2 + 1][ll * 2 + 1];
1206: }
1207: }
1208: #endif
1209: PetscCall(MatSetValuesBlockedLocal(B, 8, rc3blocked, 8, rc3blocked, &Ke[0][0], ADD_VALUES)); /* velocity-velocity coupling can use blocked insertion */
1210: { /* The off-diagonal part cannot (yet) */
1211: PetscInt row3scalar[NODE_SIZE * 8], col2scalar[PRMNODE_SIZE * 4];
1212: for (l = 0; l < 8; l++)
1213: for (ll = 0; ll < NODE_SIZE; ll++) row3scalar[l * NODE_SIZE + ll] = rc3blocked[l] * NODE_SIZE + ll;
1214: for (l = 0; l < 4; l++)
1215: for (ll = 0; ll < PRMNODE_SIZE; ll++) col2scalar[l * PRMNODE_SIZE + ll] = col2blocked[l] * PRMNODE_SIZE + ll;
1216: PetscCall(MatSetValuesLocal(Bcpl, 8 * NODE_SIZE, row3scalar, 4 * PRMNODE_SIZE, col2scalar, &Kcpl[0][0], ADD_VALUES));
1217: }
1218: }
1219: }
1220: }
1221: }
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo *info, const Node ***x3, const PrmNode **x2, const PrmNode **xdot2, PetscReal a, Mat B22, Mat B21, THI thi)
1226: {
1227: PetscInt xs, ys, xm, ym, zm, i, j, k;
1229: PetscFunctionBeginUser;
1230: xs = info->zs;
1231: ys = info->ys;
1232: xm = info->zm;
1233: ym = info->ym;
1234: zm = info->xm;
1236: PetscCheck(zm <= 1024, ((PetscObject)info->da)->comm, PETSC_ERR_SUP, "Need to allocate more space");
1237: for (i = xs; i < xs + xm; i++) {
1238: for (j = ys; j < ys + ym; j++) {
1239: { /* Self-coupling */
1240: const PetscInt row[] = {DMDALocalIndex2D(info, i, j)};
1241: const PetscInt col[] = {DMDALocalIndex2D(info, i, j)};
1242: const PetscScalar vals[] = {a, 0, 0, 0, a, 0, 0, 0, a};
1243: PetscCall(MatSetValuesBlockedLocal(B22, 1, row, 1, col, vals, INSERT_VALUES));
1244: }
1245: for (k = 0; k < zm; k++) { /* Coupling to velocity problem */
1246: /* Use a cheaper quadrature than for residual evaluation, because it is much sparser */
1247: const PetscInt row[] = {FieldIndex(PrmNode, DMDALocalIndex2D(info, i, j), h)};
1248: const PetscInt cols[] = {FieldIndex(Node, DMDALocalIndex3D(info, i - 1, j, k), u), FieldIndex(Node, DMDALocalIndex3D(info, i, j, k), u), FieldIndex(Node, DMDALocalIndex3D(info, i + 1, j, k), u),
1249: FieldIndex(Node, DMDALocalIndex3D(info, i, j - 1, k), v), FieldIndex(Node, DMDALocalIndex3D(info, i, j, k), v), FieldIndex(Node, DMDALocalIndex3D(info, i, j + 1, k), v)};
1250: const PetscScalar w = (k && k < zm - 1) ? 0.5 : 0.25, hW = w * (x2[i - 1][j].h + x2[i][j].h) / (zm - 1.), hE = w * (x2[i][j].h + x2[i + 1][j].h) / (zm - 1.), hS = w * (x2[i][j - 1].h + x2[i][j].h) / (zm - 1.),
1251: hN = w * (x2[i][j].h + x2[i][j + 1].h) / (zm - 1.);
1252: PetscScalar *vals, vals_upwind[] = {((PetscRealPart(x3[i][j][k].u) > 0) ? -hW : 0), (PetscRealPart(x3[i][j][k].u) > 0) ? +hE : -hW, (PetscRealPart(x3[i][j][k].u) > 0) ? 0 : +hE,
1253: (PetscRealPart(x3[i][j][k].v) > 0) ? -hS : 0, (PetscRealPart(x3[i][j][k].v) > 0) ? +hN : -hS, ((PetscRealPart(x3[i][j][k].v) > 0) ? 0 : +hN)},
1254: vals_centered[] = {-0.5 * hW, 0.5 * (-hW + hE), 0.5 * hE, -0.5 * hS, 0.5 * (-hS + hN), 0.5 * hN};
1255: vals = 1 ? vals_upwind : vals_centered;
1256: if (k == 0) {
1257: Node derate;
1258: THIErosion(thi, &x3[i][j][0], NULL, &derate);
1259: vals[1] -= derate.u;
1260: vals[4] -= derate.v;
1261: }
1262: PetscCall(MatSetValuesLocal(B21, 1, row, 6, cols, vals, INSERT_VALUES));
1263: }
1264: }
1265: }
1266: PetscFunctionReturn(PETSC_SUCCESS);
1267: }
1269: static PetscErrorCode THIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
1270: {
1271: THI thi = (THI)ctx;
1272: DM pack, da3, da2;
1273: Vec X3, X2, Xdot2;
1274: Mat B11, B12, B21, B22;
1275: DMDALocalInfo info3;
1276: IS *isloc;
1277: const Node ***x3;
1278: const PrmNode **x2, **xdot2;
1280: PetscFunctionBeginUser;
1281: PetscCall(TSGetDM(ts, &pack));
1282: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
1283: PetscCall(DMDAGetLocalInfo(da3, &info3));
1284: PetscCall(DMCompositeGetLocalVectors(pack, &X3, &X2));
1285: PetscCall(DMCompositeGetLocalVectors(pack, NULL, &Xdot2));
1286: PetscCall(DMCompositeScatter(pack, X, X3, X2));
1287: PetscCall(THIFixGhosts(thi, da3, da2, X3, X2));
1288: PetscCall(DMCompositeScatter(pack, Xdot, NULL, Xdot2));
1290: PetscCall(MatZeroEntries(B));
1292: PetscCall(DMCompositeGetLocalISs(pack, &isloc));
1293: PetscCall(MatGetLocalSubMatrix(B, isloc[0], isloc[0], &B11));
1294: PetscCall(MatGetLocalSubMatrix(B, isloc[0], isloc[1], &B12));
1295: PetscCall(MatGetLocalSubMatrix(B, isloc[1], isloc[0], &B21));
1296: PetscCall(MatGetLocalSubMatrix(B, isloc[1], isloc[1], &B22));
1298: PetscCall(DMDAVecGetArray(da3, X3, &x3));
1299: PetscCall(DMDAVecGetArray(da2, X2, &x2));
1300: PetscCall(DMDAVecGetArray(da2, Xdot2, &xdot2));
1302: PetscCall(THIJacobianLocal_Momentum(&info3, x3, x2, B11, B12, thi));
1304: /* Need to switch from ADD_VALUES to INSERT_VALUES */
1305: PetscCall(MatAssemblyBegin(B, MAT_FLUSH_ASSEMBLY));
1306: PetscCall(MatAssemblyEnd(B, MAT_FLUSH_ASSEMBLY));
1308: PetscCall(THIJacobianLocal_2D(&info3, x3, x2, xdot2, a, B22, B21, thi));
1310: PetscCall(DMDAVecRestoreArray(da3, X3, &x3));
1311: PetscCall(DMDAVecRestoreArray(da2, X2, &x2));
1312: PetscCall(DMDAVecRestoreArray(da2, Xdot2, &xdot2));
1314: PetscCall(MatRestoreLocalSubMatrix(B, isloc[0], isloc[0], &B11));
1315: PetscCall(MatRestoreLocalSubMatrix(B, isloc[0], isloc[1], &B12));
1316: PetscCall(MatRestoreLocalSubMatrix(B, isloc[1], isloc[0], &B21));
1317: PetscCall(MatRestoreLocalSubMatrix(B, isloc[1], isloc[1], &B22));
1318: PetscCall(ISDestroy(&isloc[0]));
1319: PetscCall(ISDestroy(&isloc[1]));
1320: PetscCall(PetscFree(isloc));
1322: PetscCall(DMCompositeRestoreLocalVectors(pack, &X3, &X2));
1323: PetscCall(DMCompositeRestoreLocalVectors(pack, 0, &Xdot2));
1325: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1326: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1327: if (A != B) {
1328: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1329: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1330: }
1331: if (thi->verbose) PetscCall(THIMatrixStatistics(thi, B, PETSC_VIEWER_STDOUT_WORLD));
1332: PetscFunctionReturn(PETSC_SUCCESS);
1333: }
1335: /* VTK's XML formats are so brain-dead that they can't handle multiple grids in the same file. Since the communication
1336: * can be shared between the two grids, we write two files at once, one for velocity (living on a 3D grid defined by
1337: * h=thickness and b=bed) and another for all properties living on the 2D grid.
1338: */
1339: static PetscErrorCode THIDAVecView_VTK_XML(THI thi, DM pack, Vec X, const char filename[], const char filename2[])
1340: {
1341: const PetscInt dof = NODE_SIZE, dof2 = PRMNODE_SIZE;
1342: Units units = thi->units;
1343: MPI_Comm comm;
1344: PetscViewer viewer3, viewer2;
1345: PetscMPIInt rank, size, tag, nn, nmax, nn2, nmax2;
1346: PetscInt mx, my, mz, r, range[6];
1347: PetscScalar *x, *x2;
1348: DM da3, da2;
1349: Vec X3, X2;
1351: PetscFunctionBeginUser;
1352: PetscCall(PetscObjectGetComm((PetscObject)thi, &comm));
1353: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
1354: PetscCall(DMCompositeGetAccess(pack, X, &X3, &X2));
1355: PetscCall(DMDAGetInfo(da3, 0, &mz, &my, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1356: PetscCallMPI(MPI_Comm_size(comm, &size));
1357: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1358: PetscCall(PetscViewerASCIIOpen(comm, filename, &viewer3));
1359: PetscCall(PetscViewerASCIIOpen(comm, filename2, &viewer2));
1360: PetscCall(PetscViewerASCIIPrintf(viewer3, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n"));
1361: PetscCall(PetscViewerASCIIPrintf(viewer2, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n"));
1362: PetscCall(PetscViewerASCIIPrintf(viewer3, " <StructuredGrid WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, mz - 1, 0, my - 1, 0, mx - 1));
1363: PetscCall(PetscViewerASCIIPrintf(viewer2, " <StructuredGrid WholeExtent=\"%d %d %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, 0, 0, my - 1, 0, mx - 1));
1365: PetscCall(DMDAGetCorners(da3, range, range + 1, range + 2, range + 3, range + 4, range + 5));
1366: PetscCall(PetscMPIIntCast(range[3] * range[4] * range[5] * dof, &nn));
1367: PetscCallMPI(MPI_Reduce(&nn, &nmax, 1, MPI_INT, MPI_MAX, 0, comm));
1368: PetscCall(PetscMPIIntCast(range[4] * range[5] * dof2, &nn2));
1369: PetscCallMPI(MPI_Reduce(&nn2, &nmax2, 1, MPI_INT, MPI_MAX, 0, comm));
1370: tag = ((PetscObject)viewer3)->tag;
1371: PetscCall(VecGetArrayRead(X3, (const PetscScalar **)&x));
1372: PetscCall(VecGetArrayRead(X2, (const PetscScalar **)&x2));
1373: if (rank == 0) {
1374: PetscScalar *array, *array2;
1375: PetscCall(PetscMalloc2(nmax, &array, nmax2, &array2));
1376: for (r = 0; r < size; r++) {
1377: PetscInt i, j, k, f, xs, xm, ys, ym, zs, zm;
1378: Node *y3;
1379: PetscScalar (*y2)[PRMNODE_SIZE];
1380: MPI_Status status;
1382: if (r) PetscCallMPI(MPI_Recv(range, 6, MPIU_INT, r, tag, comm, MPI_STATUS_IGNORE));
1383: zs = range[0];
1384: ys = range[1];
1385: xs = range[2];
1386: zm = range[3];
1387: ym = range[4];
1388: xm = range[5];
1389: PetscCheck(xm * ym * zm * dof <= nmax, PETSC_COMM_SELF, PETSC_ERR_PLIB, "should not happen");
1390: if (r) {
1391: PetscCallMPI(MPI_Recv(array, nmax, MPIU_SCALAR, r, tag, comm, &status));
1392: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn));
1393: PetscCheck(nn == xm * ym * zm * dof, PETSC_COMM_SELF, PETSC_ERR_PLIB, "corrupt da3 send");
1394: y3 = (Node *)array;
1395: PetscCallMPI(MPI_Recv(array2, nmax2, MPIU_SCALAR, r, tag, comm, &status));
1396: PetscCallMPI(MPI_Get_count(&status, MPIU_SCALAR, &nn2));
1397: PetscCheck(nn2 == xm * ym * dof2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "corrupt da2 send");
1398: y2 = (PetscScalar (*)[PRMNODE_SIZE])array2;
1399: } else {
1400: y3 = (Node *)x;
1401: y2 = (PetscScalar (*)[PRMNODE_SIZE])x2;
1402: }
1403: PetscCall(PetscViewerASCIIPrintf(viewer3, " <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", zs, zs + zm - 1, ys, ys + ym - 1, xs, xs + xm - 1));
1404: PetscCall(PetscViewerASCIIPrintf(viewer2, " <Piece Extent=\"%d %d %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", 0, 0, ys, ys + ym - 1, xs, xs + xm - 1));
1406: PetscCall(PetscViewerASCIIPrintf(viewer3, " <Points>\n"));
1407: PetscCall(PetscViewerASCIIPrintf(viewer2, " <Points>\n"));
1408: PetscCall(PetscViewerASCIIPrintf(viewer3, " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n"));
1409: PetscCall(PetscViewerASCIIPrintf(viewer2, " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n"));
1410: for (i = xs; i < xs + xm; i++) {
1411: for (j = ys; j < ys + ym; j++) {
1412: PetscReal xx = thi->Lx * i / mx, yy = thi->Ly * j / my, b = PetscRealPart(y2[i * ym + j][FieldOffset(PrmNode, b)]), h = PetscRealPart(y2[i * ym + j][FieldOffset(PrmNode, h)]);
1413: for (k = zs; k < zs + zm; k++) {
1414: PetscReal zz = b + h * k / (mz - 1.);
1415: PetscCall(PetscViewerASCIIPrintf(viewer3, "%f %f %f\n", (double)xx, (double)yy, (double)zz));
1416: }
1417: PetscCall(PetscViewerASCIIPrintf(viewer2, "%f %f %f\n", (double)xx, (double)yy, (double)0.0));
1418: }
1419: }
1420: PetscCall(PetscViewerASCIIPrintf(viewer3, " </DataArray>\n"));
1421: PetscCall(PetscViewerASCIIPrintf(viewer2, " </DataArray>\n"));
1422: PetscCall(PetscViewerASCIIPrintf(viewer3, " </Points>\n"));
1423: PetscCall(PetscViewerASCIIPrintf(viewer2, " </Points>\n"));
1425: { /* Velocity and rank (3D) */
1426: PetscCall(PetscViewerASCIIPrintf(viewer3, " <PointData>\n"));
1427: PetscCall(PetscViewerASCIIPrintf(viewer3, " <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n"));
1428: for (i = 0; i < nn / dof; i++) PetscCall(PetscViewerASCIIPrintf(viewer3, "%f %f %f\n", (double)(PetscRealPart(y3[i].u) * units->year / units->meter), (double)(PetscRealPart(y3[i].v) * units->year / units->meter), 0.0));
1429: PetscCall(PetscViewerASCIIPrintf(viewer3, " </DataArray>\n"));
1431: PetscCall(PetscViewerASCIIPrintf(viewer3, " <DataArray type=\"Int32\" Name=\"rank\" NumberOfComponents=\"1\" format=\"ascii\">\n"));
1432: for (i = 0; i < nn; i += dof) PetscCall(PetscViewerASCIIPrintf(viewer3, "%" PetscInt_FMT "\n", r));
1433: PetscCall(PetscViewerASCIIPrintf(viewer3, " </DataArray>\n"));
1434: PetscCall(PetscViewerASCIIPrintf(viewer3, " </PointData>\n"));
1435: }
1437: { /* 2D */
1438: PetscCall(PetscViewerASCIIPrintf(viewer2, " <PointData>\n"));
1439: for (f = 0; f < PRMNODE_SIZE; f++) {
1440: const char *fieldname;
1441: PetscCall(DMDAGetFieldName(da2, f, &fieldname));
1442: PetscCall(PetscViewerASCIIPrintf(viewer2, " <DataArray type=\"Float32\" Name=\"%s\" format=\"ascii\">\n", fieldname));
1443: for (i = 0; i < nn2 / PRMNODE_SIZE; i++) PetscCall(PetscViewerASCIIPrintf(viewer2, "%g\n", (double)y2[i][f]));
1444: PetscCall(PetscViewerASCIIPrintf(viewer2, " </DataArray>\n"));
1445: }
1446: PetscCall(PetscViewerASCIIPrintf(viewer2, " </PointData>\n"));
1447: }
1449: PetscCall(PetscViewerASCIIPrintf(viewer3, " </Piece>\n"));
1450: PetscCall(PetscViewerASCIIPrintf(viewer2, " </Piece>\n"));
1451: }
1452: PetscCall(PetscFree2(array, array2));
1453: } else {
1454: PetscCallMPI(MPI_Send(range, 6, MPIU_INT, 0, tag, comm));
1455: PetscCallMPI(MPI_Send(x, nn, MPIU_SCALAR, 0, tag, comm));
1456: PetscCallMPI(MPI_Send(x2, nn2, MPIU_SCALAR, 0, tag, comm));
1457: }
1458: PetscCall(VecRestoreArrayRead(X3, (const PetscScalar **)&x));
1459: PetscCall(VecRestoreArrayRead(X2, (const PetscScalar **)&x2));
1460: PetscCall(PetscViewerASCIIPrintf(viewer3, " </StructuredGrid>\n"));
1461: PetscCall(PetscViewerASCIIPrintf(viewer2, " </StructuredGrid>\n"));
1463: PetscCall(DMCompositeRestoreAccess(pack, X, &X3, &X2));
1464: PetscCall(PetscViewerASCIIPrintf(viewer3, "</VTKFile>\n"));
1465: PetscCall(PetscViewerASCIIPrintf(viewer2, "</VTKFile>\n"));
1466: PetscCall(PetscViewerDestroy(&viewer3));
1467: PetscCall(PetscViewerDestroy(&viewer2));
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: static PetscErrorCode THITSMonitor(TS ts, PetscInt step, PetscReal t, Vec X, PetscCtx ctx)
1472: {
1473: THI thi = (THI)ctx;
1474: DM pack;
1475: char filename3[PETSC_MAX_PATH_LEN], filename2[PETSC_MAX_PATH_LEN];
1477: PetscFunctionBeginUser;
1478: if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* negative one is used to indicate an interpolated solution */
1479: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%3" PetscInt_FMT ": t=%g\n", step, (double)t));
1480: if (thi->monitor_interval && step % thi->monitor_interval) PetscFunctionReturn(PETSC_SUCCESS);
1481: PetscCall(TSGetDM(ts, &pack));
1482: PetscCall(PetscSNPrintf(filename3, sizeof(filename3), "%s-3d-%03" PetscInt_FMT ".vts", thi->monitor_basename, step));
1483: PetscCall(PetscSNPrintf(filename2, sizeof(filename2), "%s-2d-%03" PetscInt_FMT ".vts", thi->monitor_basename, step));
1484: PetscCall(THIDAVecView_VTK_XML(thi, pack, X, filename3, filename2));
1485: PetscFunctionReturn(PETSC_SUCCESS);
1486: }
1488: static PetscErrorCode THICreateDM3d(THI thi, DM *dm3d)
1489: {
1490: MPI_Comm comm;
1491: PetscInt M = 3, N = 3, P = 2;
1492: DM da;
1494: PetscFunctionBeginUser;
1495: PetscCall(PetscObjectGetComm((PetscObject)thi, &comm));
1496: PetscOptionsBegin(comm, NULL, "Grid resolution options", "");
1497: {
1498: PetscCall(PetscOptionsInt("-M", "Number of elements in x-direction on coarse level", "", M, &M, NULL));
1499: N = M;
1500: PetscCall(PetscOptionsInt("-N", "Number of elements in y-direction on coarse level (if different from M)", "", N, &N, NULL));
1501: PetscCall(PetscOptionsInt("-P", "Number of elements in z-direction on coarse level", "", P, &P, NULL));
1502: }
1503: PetscOptionsEnd();
1504: PetscCall(DMDACreate3d(comm, DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_BOX, P, N, M, 1, PETSC_DETERMINE, PETSC_DETERMINE, sizeof(Node) / sizeof(PetscScalar), 1, 0, 0, 0, &da));
1505: PetscCall(DMSetFromOptions(da));
1506: PetscCall(DMSetUp(da));
1507: PetscCall(DMDASetFieldName(da, 0, "x-velocity"));
1508: PetscCall(DMDASetFieldName(da, 1, "y-velocity"));
1509: *dm3d = da;
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }
1513: int main(int argc, char *argv[])
1514: {
1515: MPI_Comm comm;
1516: DM pack, da3, da2;
1517: TS ts;
1518: THI thi;
1519: Vec X;
1520: Mat B;
1521: PetscInt steps;
1522: PetscReal ftime;
1524: PetscFunctionBeginUser;
1525: PetscCall(PetscInitialize(&argc, &argv, 0, help));
1526: comm = PETSC_COMM_WORLD;
1528: PetscCall(THICreate(comm, &thi));
1529: PetscCall(THICreateDM3d(thi, &da3));
1530: {
1531: PetscInt Mx, My, mx, my, s;
1532: DMDAStencilType st;
1533: PetscCall(DMDAGetInfo(da3, 0, 0, &My, &Mx, 0, &my, &mx, 0, &s, 0, 0, 0, &st));
1534: PetscCall(DMDACreate2d(PetscObjectComm((PetscObject)thi), DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, st, My, Mx, my, mx, sizeof(PrmNode) / sizeof(PetscScalar), s, 0, 0, &da2));
1535: PetscCall(DMSetUp(da2));
1536: }
1538: PetscCall(PetscObjectSetName((PetscObject)da3, "3D_Velocity"));
1539: PetscCall(DMSetOptionsPrefix(da3, "f3d_"));
1540: PetscCall(DMDASetFieldName(da3, 0, "u"));
1541: PetscCall(DMDASetFieldName(da3, 1, "v"));
1542: PetscCall(PetscObjectSetName((PetscObject)da2, "2D_Fields"));
1543: PetscCall(DMSetOptionsPrefix(da2, "f2d_"));
1544: PetscCall(DMDASetFieldName(da2, 0, "b"));
1545: PetscCall(DMDASetFieldName(da2, 1, "h"));
1546: PetscCall(DMDASetFieldName(da2, 2, "beta2"));
1547: PetscCall(DMCompositeCreate(comm, &pack));
1548: PetscCall(DMCompositeAddDM(pack, da3));
1549: PetscCall(DMCompositeAddDM(pack, da2));
1550: PetscCall(DMDestroy(&da3));
1551: PetscCall(DMDestroy(&da2));
1552: PetscCall(DMSetUp(pack));
1553: PetscCall(DMCreateMatrix(pack, &B));
1554: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
1555: PetscCall(MatSetOptionsPrefix(B, "thi_"));
1557: for (PetscInt i = 0; i < thi->nlevels; i++) {
1558: PetscReal Lx = thi->Lx / thi->units->meter, Ly = thi->Ly / thi->units->meter, Lz = thi->Lz / thi->units->meter;
1559: PetscInt Mx, My, Mz;
1560: PetscCall(DMCompositeGetEntries(pack, &da3, &da2));
1561: PetscCall(DMDAGetInfo(da3, 0, &Mz, &My, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1562: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)thi), "Level %" PetscInt_FMT " domain size (m) %8.2g x %8.2g x %8.2g, num elements %3d x %3d x %3d (%8d), size (m) %g x %g x %g\n", i, Lx, Ly, Lz, Mx, My, Mz, Mx * My * Mz, Lx / Mx, Ly / My, 1000. / (Mz - 1)));
1563: }
1565: PetscCall(DMCreateGlobalVector(pack, &X));
1566: PetscCall(THIInitial(thi, pack, X));
1568: PetscCall(TSCreate(comm, &ts));
1569: PetscCall(TSSetDM(ts, pack));
1570: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
1571: PetscCall(TSMonitorSet(ts, THITSMonitor, thi, NULL));
1572: PetscCall(TSSetType(ts, TSTHETA));
1573: PetscCall(TSSetIFunction(ts, NULL, THIFunction, thi));
1574: PetscCall(TSSetIJacobian(ts, B, B, THIJacobian, thi));
1575: PetscCall(TSSetMaxTime(ts, 10.0));
1576: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1577: PetscCall(TSSetSolution(ts, X));
1578: PetscCall(TSSetTimeStep(ts, 1e-3));
1579: PetscCall(TSSetFromOptions(ts));
1581: PetscCall(TSSolve(ts, X));
1582: PetscCall(TSGetSolveTime(ts, &ftime));
1583: PetscCall(TSGetStepNumber(ts, &steps));
1584: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Steps %" PetscInt_FMT " final time %g\n", steps, (double)ftime));
1586: if (0) PetscCall(THISolveStatistics(thi, ts, 0, "Full"));
1588: {
1589: PetscBool flg;
1590: char filename[PETSC_MAX_PATH_LEN] = "";
1591: PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &flg));
1592: if (flg) PetscCall(THIDAVecView_VTK_XML(thi, pack, X, filename, NULL));
1593: }
1595: PetscCall(VecDestroy(&X));
1596: PetscCall(MatDestroy(&B));
1597: PetscCall(DMDestroy(&pack));
1598: PetscCall(TSDestroy(&ts));
1599: PetscCall(THIDestroy(&thi));
1600: PetscCall(PetscFinalize());
1601: return 0;
1602: }