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