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 <petscdmmg.h>
48: #include <petscdmcomposite.h>
49: #include <ctype.h> /* toupper() */
51: #if defined __SSE2__
52: # include <emmintrin.h>
53: #endif
55: /* The SSE2 kernels are only for PetscScalar=double on architectures that support it */
56: #define USE_SSE2_KERNELS (!defined NO_SSE2 \
57: && !defined PETSC_USE_COMPLEX \
58: && !defined PETSC_USE_REAL_SINGLE \
59: && !defined PETSC_USE_REAL_LONG_DOUBLE \
60: && defined __SSE2__)
62: #if !defined __STDC_VERSION__ || __STDC_VERSION__ < 199901L
63: # if defined __cplusplus /* C++ restrict is nonstandard and compilers have inconsistent rules about where it can be used */
64: # define restrict
65: # else
66: # define restrict PETSC_RESTRICT
67: # endif
68: #endif
70: static PetscClassId THI_CLASSID;
72: typedef enum {QUAD_GAUSS,QUAD_LOBATTO} QuadratureType;
73: static const char *QuadratureTypes[] = {"gauss","lobatto","QuadratureType","QUAD_",0};
74: static const PetscReal HexQWeights[8] = {1,1,1,1,1,1,1,1};
75: static const PetscReal HexQNodes[] = {-0.57735026918962573, 0.57735026918962573};
76: #define G 0.57735026918962573
77: #define H (0.5*(1.+G))
78: #define L (0.5*(1.-G))
79: #define M (-0.5)
80: #define P (0.5)
81: /* Special quadrature: Lobatto in horizontal, Gauss in vertical */
82: static const PetscReal HexQInterp_Lobatto[8][8] = {{H,0,0,0,L,0,0,0},
83: {0,H,0,0,0,L,0,0},
84: {0,0,H,0,0,0,L,0},
85: {0,0,0,H,0,0,0,L},
86: {L,0,0,0,H,0,0,0},
87: {0,L,0,0,0,H,0,0},
88: {0,0,L,0,0,0,H,0},
89: {0,0,0,L,0,0,0,H}};
90: static const PetscReal HexQDeriv_Lobatto[8][8][3] = {
91: {{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} },
92: {{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} },
93: {{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} },
94: {{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}},
95: {{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} },
96: {{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} },
97: {{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} },
98: {{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}}};
99: /* Stanndard Gauss */
100: static const PetscReal HexQInterp_Gauss[8][8] = {{H*H*H,L*H*H,L*L*H,H*L*H, H*H*L,L*H*L,L*L*L,H*L*L},
101: {L*H*H,H*H*H,H*L*H,L*L*H, L*H*L,H*H*L,H*L*L,L*L*L},
102: {L*L*H,H*L*H,H*H*H,L*H*H, L*L*L,H*L*L,H*H*L,L*H*L},
103: {H*L*H,L*L*H,L*H*H,H*H*H, H*L*L,L*L*L,L*H*L,H*H*L},
104: {H*H*L,L*H*L,L*L*L,H*L*L, H*H*H,L*H*H,L*L*H,H*L*H},
105: {L*H*L,H*H*L,H*L*L,L*L*L, L*H*H,H*H*H,H*L*H,L*L*H},
106: {L*L*L,H*L*L,H*H*L,L*H*L, L*L*H,H*L*H,H*H*H,L*H*H},
107: {H*L*L,L*L*L,L*H*L,H*H*L, H*L*H,L*L*H,L*H*H,H*H*H}};
108: static const PetscReal HexQDeriv_Gauss[8][8][3] = {
109: {{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}},
110: {{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}},
111: {{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}},
112: {{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}},
113: {{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}},
114: {{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}},
115: {{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}},
116: {{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}}};
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] = {{H*H,L*H,L*L,H*L},
120: {L*H,H*H,H*L,L*L},
121: {L*L,H*L,H*H,L*H},
122: {H*L,L*L,L*H,H*H}};
123: static const PetscReal QuadQDeriv[4][4][2] = {
124: {{M*H,M*H},{P*H,M*L},{P*L,P*L},{M*L,P*H}},
125: {{M*H,M*L},{P*H,M*H},{P*L,P*H},{M*L,P*L}},
126: {{M*L,M*L},{P*L,M*H},{P*H,P*H},{M*H,P*L}},
127: {{M*L,M*H},{P*L,M*L},{P*H,P*L},{M*H,P*H}}};
128: #undef G
129: #undef H
130: #undef L
131: #undef M
132: #undef P
134: #define HexExtract(x,i,j,k,n) do { \
135: (n)[0] = (x)[i][j][k]; \
136: (n)[1] = (x)[i+1][j][k]; \
137: (n)[2] = (x)[i+1][j+1][k]; \
138: (n)[3] = (x)[i][j+1][k]; \
139: (n)[4] = (x)[i][j][k+1]; \
140: (n)[5] = (x)[i+1][j][k+1]; \
141: (n)[6] = (x)[i+1][j+1][k+1]; \
142: (n)[7] = (x)[i][j+1][k+1]; \
143: } while (0)
145: #define HexExtractRef(x,i,j,k,n) 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 QuadExtract(x,i,j,n) do { \
157: (n)[0] = (x)[i][j]; \
158: (n)[1] = (x)[i+1][j]; \
159: (n)[2] = (x)[i+1][j+1]; \
160: (n)[3] = (x)[i][j+1]; \
161: } while (0)
163: static PetscScalar Sqr(PetscScalar a) {return a*a;}
165: static void HexGrad(const PetscReal dphi[][3],const PetscReal zn[],PetscReal dz[])
166: {
167: PetscInt i;
168: dz[0] = dz[1] = dz[2] = 0;
169: for (i=0; i<8; i++) {
170: dz[0] += dphi[i][0] * zn[i];
171: dz[1] += dphi[i][1] * zn[i];
172: dz[2] += dphi[i][2] * zn[i];
173: }
174: }
176: static void HexComputeGeometry(PetscInt q,PetscReal hx,PetscReal hy,const PetscReal dz[restrict],PetscReal phi[restrict],PetscReal dphi[restrict][3],PetscReal *restrict jw)
177: {
178: const PetscReal
179: jac[3][3] = {{hx/2,0,0}, {0,hy/2,0}, {dz[0],dz[1],dz[2]}}
180: ,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]}}
181: ,jdet = jac[0][0]*jac[1][1]*jac[2][2];
182: PetscInt i;
184: for (i=0; i<8; i++) {
185: const PetscReal *dphir = HexQDeriv[q][i];
186: phi[i] = HexQInterp[q][i];
187: dphi[i][0] = dphir[0]*ijac[0][0] + dphir[1]*ijac[1][0] + dphir[2]*ijac[2][0];
188: dphi[i][1] = dphir[0]*ijac[0][1] + dphir[1]*ijac[1][1] + dphir[2]*ijac[2][1];
189: dphi[i][2] = dphir[0]*ijac[0][2] + dphir[1]*ijac[1][2] + dphir[2]*ijac[2][2];
190: }
191: *jw = 1.0 * jdet;
192: }
194: typedef struct _p_THI *THI;
195: typedef struct _n_Units *Units;
197: typedef struct {
198: PetscScalar u,v;
199: } Node;
201: typedef struct {
202: PetscScalar b; /* bed */
203: PetscScalar h; /* thickness */
204: PetscScalar beta2; /* friction */
205: } PrmNode;
207: #define FieldSize(ntype) ((PetscInt)(sizeof(ntype)/sizeof(PetscScalar)))
208: #define FieldOffset(ntype,member) ((PetscInt)(offsetof(ntype,member)/sizeof(PetscScalar)))
209: #define FieldIndex(ntype,i,member) ((PetscInt)((i)*FieldSize(ntype) + FieldOffset(ntype,member)))
210: #define NODE_SIZE FieldSize(Node)
211: #define PRMNODE_SIZE FieldSize(PrmNode)
213: typedef struct {
214: PetscReal min,max,cmin,cmax;
215: } PRange;
217: struct _p_THI {
218: PETSCHEADER(int);
219: void (*initialize)(THI,PetscReal x,PetscReal y,PrmNode *p);
220: PetscInt nlevels;
221: PetscInt zlevels;
222: PetscReal Lx,Ly,Lz; /* Model domain */
223: PetscReal alpha; /* Bed angle */
224: Units units;
225: PetscReal dirichlet_scale;
226: PetscReal ssa_friction_scale;
227: PetscReal inertia;
228: PRange eta;
229: PRange beta2;
230: struct {
231: PetscReal Bd2,eps,exponent,glen_n;
232: } viscosity;
233: struct {
234: PetscReal irefgam,eps2,exponent;
235: } friction;
236: struct {
237: PetscReal rate,exponent,refvel;
238: } erosion;
239: PetscReal rhog;
240: PetscBool no_slip;
241: PetscBool verbose;
242: MatType mattype;
243: char *monitor_basename;
244: PetscInt monitor_interval;
245: };
247: struct _n_Units {
248: /* fundamental */
249: PetscReal meter;
250: PetscReal kilogram;
251: PetscReal second;
252: /* derived */
253: PetscReal Pascal;
254: PetscReal year;
255: };
257: static void PrmHexGetZ(const PrmNode pn[],PetscInt k,PetscInt zm,PetscReal zn[])
258: {
259: const PetscScalar zm1 = zm-1,
260: znl[8] = {pn[0].b + pn[0].h*(PetscScalar)k/zm1,
261: pn[1].b + pn[1].h*(PetscScalar)k/zm1,
262: pn[2].b + pn[2].h*(PetscScalar)k/zm1,
263: pn[3].b + pn[3].h*(PetscScalar)k/zm1,
264: pn[0].b + pn[0].h*(PetscScalar)(k+1)/zm1,
265: pn[1].b + pn[1].h*(PetscScalar)(k+1)/zm1,
266: pn[2].b + pn[2].h*(PetscScalar)(k+1)/zm1,
267: pn[3].b + pn[3].h*(PetscScalar)(k+1)/zm1};
268: PetscInt i;
269: for (i=0; i<8; i++) zn[i] = PetscRealPart(znl[i]);
270: }
274: /* Compute a gradient of all the 2D fields at four quadrature points. Output for [quadrature_point][direction].field_name */
275: static PetscErrorCode QuadComputeGrad4(const PetscReal dphi[][4][2],PetscReal hx,PetscReal hy,const PrmNode pn[4],PrmNode dp[4][2])
276: {
278: PetscInt q,i,f;
279: const PetscScalar (*restrict pg)[PRMNODE_SIZE] = (const PetscScalar(*)[PRMNODE_SIZE])pn; /* Get generic array pointers to the node */
280: PetscScalar (*restrict dpg)[2][PRMNODE_SIZE] = (PetscScalar(*)[2][PRMNODE_SIZE])dp;
283: PetscMemzero(dpg,4*sizeof(dpg[0]));
284: for (q=0; q<4; q++) {
285: for (i=0; i<4; i++) {
286: for (f=0; f<PRMNODE_SIZE; f++) {
287: dpg[q][0][f] += dphi[q][i][0]/hx * pg[i][f];
288: dpg[q][1][f] += dphi[q][i][1]/hy * pg[i][f];
289: }
290: }
291: }
292: return(0);
293: }
295: static inline PetscReal StaggeredMidpoint2D(PetscScalar a,PetscScalar b,PetscScalar c,PetscScalar d)
296: {return 0.5*PetscRealPart(0.75*a + 0.75*b + 0.25*c + 0.25*d);}
297: static inline PetscReal UpwindFlux1D(PetscReal u,PetscReal hL,PetscReal hR)
298: {return (u > 0) ? hL*u : hR*u;}
300: #define UpwindFluxXW(x3,x2,h,i,j,k,dj) 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), \
301: PetscRealPart(0.75*x2[i-1][j ].h+0.25*x2[i-1][j+dj].h), PetscRealPart(0.75*x2[i ][j ].h+0.25*x2[i ][j+dj].h))
302: #define UpwindFluxXE(x3,x2,h,i,j,k,dj) 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), \
303: PetscRealPart(0.75*x2[i ][j ].h+0.25*x2[i ][j+dj].h), PetscRealPart(0.75*x2[i+1][j ].h+0.25*x2[i+1][j+dj].h))
304: #define UpwindFluxYS(x3,x2,h,i,j,k,di) 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), \
305: PetscRealPart(0.75*x2[i ][j-1].h+0.25*x2[i+di][j-1].h), PetscRealPart(0.75*x2[i ][j ].h+0.25*x2[i+di][j ].h))
306: #define UpwindFluxYN(x3,x2,h,i,j,k,di) 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), \
307: PetscRealPart(0.75*x2[i ][j ].h+0.25*x2[i+di][j ].h), PetscRealPart(0.75*x2[i ][j+1].h+0.25*x2[i+di][j+1].h))
309: static void PrmNodeGetFaceMeasure(const PrmNode **p,PetscInt i,PetscInt j,PetscScalar h[])
310: {
311: /* West */
312: h[0] = StaggeredMidpoint2D(p[i][j].h,p[i-1][j].h,p[i-1][j-1].h,p[i][j-1].h);
313: h[1] = StaggeredMidpoint2D(p[i][j].h,p[i-1][j].h,p[i-1][j+1].h,p[i][j+1].h);
314: /* East */
315: h[2] = StaggeredMidpoint2D(p[i][j].h,p[i+1][j].h,p[i+1][j+1].h,p[i][j+1].h);
316: h[3] = StaggeredMidpoint2D(p[i][j].h,p[i+1][j].h,p[i+1][j-1].h,p[i][j-1].h);
317: /* South */
318: h[4] = StaggeredMidpoint2D(p[i][j].h,p[i][j-1].h,p[i+1][j-1].h,p[i+1][j].h);
319: h[5] = StaggeredMidpoint2D(p[i][j].h,p[i][j-1].h,p[i-1][j-1].h,p[i-1][j].h);
320: /* North */
321: h[6] = StaggeredMidpoint2D(p[i][j].h,p[i][j+1].h,p[i-1][j+1].h,p[i-1][j].h);
322: h[7] = StaggeredMidpoint2D(p[i][j].h,p[i][j+1].h,p[i+1][j+1].h,p[i+1][j].h);
323: }
325: /* Tests A and C are from the ISMIP-HOM paper (Pattyn et al. 2008) */
326: static void THIInitialize_HOM_A(THI thi,PetscReal x,PetscReal y,PrmNode *p)
327: {
328: Units units = thi->units;
329: PetscReal s = -x*sin(thi->alpha);
330: p->b = s - 1000*units->meter + 500*units->meter * sin(x*2*PETSC_PI/thi->Lx) * sin(y*2*PETSC_PI/thi->Ly);
331: p->h = s - p->b;
332: p->beta2 = -1e-10; /* This value is not used, but it should not be huge because that would change the finite difference step size */
333: }
335: static void THIInitialize_HOM_C(THI thi,PetscReal x,PetscReal y,PrmNode *p)
336: {
337: Units units = thi->units;
338: PetscReal s = -x*sin(thi->alpha);
339: p->b = s - 1000*units->meter;
340: p->h = s - p->b;
341: /* tau_b = beta2 v is a stress (Pa).
342: * 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. */
343: p->beta2 = 1000 * (1 + sin(x*2*PETSC_PI/thi->Lx)*sin(y*2*PETSC_PI/thi->Ly)) * units->Pascal * units->year / units->meter / thi->rhog;
344: }
346: /* These are just toys */
348: /* From Fred Herman */
349: static void THIInitialize_HOM_F(THI thi,PetscReal x,PetscReal y,PrmNode *p)
350: {
351: Units units = thi->units;
352: PetscReal s = -x*sin(thi->alpha);
353: p->b = s - 1000*units->meter + 100*units->meter * sin(x*2*PETSC_PI/thi->Lx);// * sin(y*2*PETSC_PI/thi->Ly);
354: p->h = s - p->b;
355: p->h = (1-(atan((x-thi->Lx/2)/1.)+PETSC_PI/2.)/PETSC_PI)*500*units->meter+1*units->meter;
356: s = PetscRealPart(p->b + p->h);
357: p->beta2 = -1e-10;
358: // p->beta2 = 1000 * units->Pascal * units->year / units->meter;
359: }
361: /* Same bed as test A, free slip everywhere except for a discontinuous jump to a circular sticky region in the middle. */
362: static void THIInitialize_HOM_X(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
363: {
364: Units units = thi->units;
365: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
366: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
367: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
368: p->h = s - p->b;
369: p->beta2 = 1000 * (r < 1 ? 2 : 0) * units->Pascal * units->year / units->meter / thi->rhog;
370: }
372: /* Like Z, but with 200 meter cliffs */
373: static void THIInitialize_HOM_Y(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
374: {
375: Units units = thi->units;
376: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
377: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
378: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
379: if (PetscRealPart(p->b) > -700*units->meter) p->b += 200*units->meter;
380: p->h = s - p->b;
381: p->beta2 = 1000 * (1. + sin(sqrt(16*r))/sqrt(1e-2 + 16*r)*cos(x*3/2)*cos(y*3/2)) * units->Pascal * units->year / units->meter / thi->rhog;
382: }
384: /* Same bed as A, smoothly varying slipperiness, similar to MATLAB's "sombrero" (uncorrelated with bathymetry) */
385: static void THIInitialize_HOM_Z(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
386: {
387: Units units = thi->units;
388: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
389: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
390: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
391: p->h = s - p->b;
392: p->beta2 = 1000 * (1. + sin(sqrt(16*r))/sqrt(1e-2 + 16*r)*cos(x*3/2)*cos(y*3/2)) * units->Pascal * units->year / units->meter / thi->rhog;
393: }
395: static void THIFriction(THI thi,PetscReal rbeta2,PetscReal gam,PetscReal *beta2,PetscReal *dbeta2)
396: {
397: if (thi->friction.irefgam == 0) {
398: Units units = thi->units;
399: thi->friction.irefgam = 1./(0.5*PetscSqr(100 * units->meter / units->year));
400: thi->friction.eps2 = 0.5*PetscSqr(1.e-4 / thi->friction.irefgam);
401: }
402: if (thi->friction.exponent == 0) {
403: *beta2 = rbeta2;
404: *dbeta2 = 0;
405: } else {
406: *beta2 = rbeta2 * pow(thi->friction.eps2 + gam*thi->friction.irefgam,thi->friction.exponent);
407: *dbeta2 = thi->friction.exponent * *beta2 / (thi->friction.eps2 + gam*thi->friction.irefgam) * thi->friction.irefgam;
408: }
409: }
411: static void THIViscosity(THI thi,PetscReal gam,PetscReal *eta,PetscReal *deta)
412: {
413: PetscReal Bd2,eps,exponent;
414: if (thi->viscosity.Bd2 == 0) {
415: Units units = thi->units;
416: const PetscReal
417: n = thi->viscosity.glen_n, /* Glen exponent */
418: p = 1. + 1./n, /* for Stokes */
419: A = 1.e-16 * pow(units->Pascal,-n) / units->year, /* softness parameter (Pa^{-n}/s) */
420: B = pow(A,-1./n); /* hardness parameter */
421: thi->viscosity.Bd2 = B/2;
422: thi->viscosity.exponent = (p-2)/2;
423: thi->viscosity.eps = 0.5*PetscSqr(1e-5 / units->year);
424: }
425: Bd2 = thi->viscosity.Bd2;
426: exponent = thi->viscosity.exponent;
427: eps = thi->viscosity.eps;
428: *eta = Bd2 * pow(eps + gam,exponent);
429: *deta = exponent * (*eta) / (eps + gam);
430: }
432: static void THIErosion(THI thi,const Node *vel,PetscScalar *erate,Node *derate)
433: {
434: const PetscScalar magref2 = 1.e-10 + (PetscSqr(vel->u) + PetscSqr(vel->v)) / PetscSqr(thi->erosion.refvel),
435: rate = - thi->erosion.rate * PetscPowScalar(magref2, 0.5*thi->erosion.exponent);
436: if (erate) *erate = rate;
437: if (derate) {
438: if (thi->erosion.exponent == 1) {
439: derate->u = 0;
440: derate->v = 0;
441: } else {
442: derate->u = 0.5*thi->erosion.exponent * rate / magref2 * 2. * vel->u / PetscSqr(thi->erosion.refvel);
443: derate->v = 0.5*thi->erosion.exponent * rate / magref2 * 2. * vel->v / PetscSqr(thi->erosion.refvel);
444: }
445: }
446: }
448: static void RangeUpdate(PetscReal *min,PetscReal *max,PetscReal x)
449: {
450: if (x < *min) *min = x;
451: if (x > *max) *max = x;
452: }
454: static void PRangeClear(PRange *p)
455: {
456: p->cmin = p->min = 1e100;
457: p->cmax = p->max = -1e100;
458: }
462: static PetscErrorCode PRangeMinMax(PRange *p,PetscReal min,PetscReal max)
463: {
466: p->cmin = min;
467: p->cmax = max;
468: if (min < p->min) p->min = min;
469: if (max > p->max) p->max = max;
470: return(0);
471: }
475: static PetscErrorCode THIDestroy(THI *thi)
476: {
480: if (--((PetscObject)(*thi))->refct > 0) return(0);
481: PetscFree((*thi)->units);
482: PetscFree((*thi)->mattype);
483: PetscFree((*thi)->monitor_basename);
484: PetscHeaderDestroy(thi);
485: return(0);
486: }
490: static PetscErrorCode THICreate(MPI_Comm comm,THI *inthi)
491: {
492: static PetscBool registered = PETSC_FALSE;
493: THI thi;
494: Units units;
495: char monitor_basename[PETSC_MAX_PATH_LEN] = "thi-";
496: PetscErrorCode ierr;
499: *inthi = 0;
500: if (!registered) {
501: PetscClassIdRegister("Toy Hydrostatic Ice",&THI_CLASSID);
502: registered = PETSC_TRUE;
503: }
504: PetscHeaderCreate(thi,_p_THI,0,THI_CLASSID,-1,"THI","Toy Hydrostatic Ice","THI",comm,THIDestroy,0);
506: PetscNew(struct _n_Units,&thi->units);
507: units = thi->units;
508: units->meter = 1e-2;
509: units->second = 1e-7;
510: units->kilogram = 1e-12;
511: PetscOptionsBegin(comm,NULL,"Scaled units options","");
512: {
513: PetscOptionsReal("-units_meter","1 meter in scaled length units","",units->meter,&units->meter,NULL);
514: PetscOptionsReal("-units_second","1 second in scaled time units","",units->second,&units->second,NULL);
515: 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: PetscOptionsReal("-thi_L","Domain size (m)","",L,&L,&flg);
542: if (flg) thi->Lx = thi->Ly = L;
543: PetscOptionsReal("-thi_Lx","X Domain size (m)","",thi->Lx,&thi->Lx,NULL);
544: PetscOptionsReal("-thi_Ly","Y Domain size (m)","",thi->Ly,&thi->Ly,NULL);
545: PetscOptionsReal("-thi_Lz","Z Domain size (m)","",thi->Lz,&thi->Lz,NULL);
546: 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: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"HOM experiment '%c' not implemented",homexp[0]);
580: }
581: 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: PetscOptionsReal("-thi_alpha","Bed angle (degrees)","",thi->alpha,&thi->alpha,NULL);
593: 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: PetscOptionsReal("-thi_friction_m","Friction exponent, 0=Coulomb, 1=Navier","",m,&m,NULL);
595: thi->friction.exponent = (m-1)/2;
596: 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: PetscOptionsReal("-thi_erosion_exponent","Power of sliding velocity appearing in erosion relation",NULL,thi->erosion.exponent,&thi->erosion.exponent,NULL);
598: 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: PetscOptionsReal("-thi_dirichlet_scale","Scale Dirichlet boundary conditions by this factor","",thi->dirichlet_scale,&thi->dirichlet_scale,NULL);
602: 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: PetscOptionsReal("-thi_inertia","Coefficient of accelaration term in velocity system, physical is almost zero",NULL,thi->inertia,&thi->inertia,NULL);
604: PetscOptionsInt("-thi_nlevels","Number of levels of refinement","",thi->nlevels,&thi->nlevels,NULL);
605: PetscOptionsList("-thi_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)mtype,sizeof(mtype),NULL);
606: PetscStrallocpy(mtype,&thi->mattype);
607: PetscOptionsBool("-thi_verbose","Enable verbose output (like matrix sizes and statistics)","",thi->verbose,&thi->verbose,NULL);
608: PetscOptionsString("-thi_monitor","Basename to write state files to",NULL,monitor_basename,monitor_basename,sizeof monitor_basename,&flg);
609: if (flg) {
610: PetscStrallocpy(monitor_basename,&thi->monitor_basename);
611: thi->monitor_interval = 1;
612: 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),
628: gradu = u / (100*units->meter),eta,deta,
629: rho = 910 * units->kilogram/pow(units->meter,3),
630: grav = 9.81 * units->meter/PetscSqr(units->second),
631: driving = rho * grav * sin(thi->alpha) * 1000*units->meter;
632: THIViscosity(thi,0.5*gradu*gradu,&eta,&deta);
633: thi->rhog = rho * grav;
634: if (thi->verbose) {
635: PetscPrintf(((PetscObject)thi)->comm,"Units: meter %8.2g second %8.2g kg %8.2g Pa %8.2g\n",units->meter,units->second,units->kilogram,units->Pascal);
636: PetscPrintf(((PetscObject)thi)->comm,"Domain (%6.2g,%6.2g,%6.2g), pressure %8.2g, driving stress %8.2g\n",thi->Lx,thi->Ly,thi->Lz,rho*grav*1e3*units->meter,driving);
637: PetscPrintf(((PetscObject)thi)->comm,"Large velocity 1km/a %8.2g, velocity gradient %8.2g, eta %8.2g, stress %8.2g, ratio %8.2g\n",u,gradu,eta,2*eta*gradu,2*eta*gradu/driving);
638: THIViscosity(thi,0.5*PetscSqr(1e-3*gradu),&eta,&deta);
639: PetscPrintf(((PetscObject)thi)->comm,"Small velocity 1m/a %8.2g, velocity gradient %8.2g, eta %8.2g, stress %8.2g, ratio %8.2g\n",1e-3*u,1e-3*gradu,eta,2*eta*1e-3*gradu,2*eta*1e-3*gradu/driving);
640: }
641: }
643: *inthi = thi;
644: return(0);
645: }
649: /* Our problem is periodic, but the domain has a mean slope of alpha so the bed does not line up between the upstream
650: * and downstream ends of the domain. This function fixes the ghost values so that the domain appears truly periodic in
651: * the horizontal. */
652: static PetscErrorCode THIFixGhosts(THI thi,DM da3,DM da2,Vec X3,Vec X2)
653: {
655: DMDALocalInfo info;
656: PrmNode **x2;
657: PetscInt i,j;
660: DMDAGetLocalInfo(da3,&info);
661: //VecView(X2,PETSC_VIEWER_STDOUT_WORLD);
662: DMDAVecGetArray(da2,X2,&x2);
663: for (i=info.gzs; i<info.gzs+info.gzm; i++) {
664: if (i > -1 && i < info.mz) continue;
665: for (j=info.gys; j<info.gys+info.gym; j++) {
666: x2[i][j].b += sin(thi->alpha) * thi->Lx * (i<0 ? 1.0 : -1.0);
667: }
668: }
669: DMDAVecRestoreArray(da2,X2,&x2);
670: //VecView(X2,PETSC_VIEWER_STDOUT_WORLD);
671: return(0);
672: }
676: static PetscErrorCode THIInitializePrm(THI thi,DM da2prm,PrmNode **p)
677: {
678: PetscInt i,j,xs,xm,ys,ym,mx,my;
682: DMDAGetGhostCorners(da2prm,&ys,&xs,0,&ym,&xm,0);
683: DMDAGetInfo(da2prm,0, &my,&mx,0, 0,0,0, 0,0,0,0,0,0);
684: for (i=xs; i<xs+xm; i++) {
685: for (j=ys; j<ys+ym; j++) {
686: PetscReal xx = thi->Lx*i/mx,yy = thi->Ly*j/my;
687: thi->initialize(thi,xx,yy,&p[i][j]);
688: }
689: }
690: return(0);
691: }
695: static PetscErrorCode THIInitial(THI thi,DM pack,Vec X)
696: {
697: DM da3,da2;
698: PetscInt i,j,k,xs,xm,ys,ym,zs,zm,mx,my;
699: PetscReal hx,hy;
700: PrmNode **prm;
701: Node ***x;
702: Vec X3g,X2g,X2;
706: DMCompositeGetEntries(pack,&da3,&da2);
707: DMCompositeGetAccess(pack,X,&X3g,&X2g);
708: DMGetLocalVector(da2,&X2);
710: DMDAGetInfo(da3,0, 0,&my,&mx, 0,0,0, 0,0,0,0,0,0);
711: DMDAGetCorners(da3,&zs,&ys,&xs,&zm,&ym,&xm);
712: DMDAVecGetArray(da3,X3g,&x);
713: DMDAVecGetArray(da2,X2,&prm);
715: THIInitializePrm(thi,da2,prm);
717: hx = thi->Lx / mx;
718: hy = thi->Ly / my;
719: for (i=xs; i<xs+xm; i++) {
720: for (j=ys; j<ys+ym; j++) {
721: for (k=zs; k<zs+zm; k++) {
722: const PetscScalar zm1 = zm-1,
723: 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),
724: drivingy = thi->rhog * (prm[i][j+1].b+prm[i][j+1].h - prm[i][j-1].b-prm[i][j-1].h) / (2*hx);
725: x[i][j][k].u = 0. * drivingx * prm[i][j].h*(PetscScalar)k/zm1;
726: x[i][j][k].v = 0. * drivingy * prm[i][j].h*(PetscScalar)k/zm1;
727: }
728: }
729: }
731: DMDAVecRestoreArray(da3,X3g,&x);
732: DMDAVecRestoreArray(da2,X2,&prm);
734: DMLocalToGlobalBegin(da2,X2,INSERT_VALUES,X2g);
735: DMLocalToGlobalEnd (da2,X2,INSERT_VALUES,X2g);
736: DMRestoreLocalVector(da2,&X2);
738: DMCompositeRestoreAccess(pack,X,&X3g,&X2g);
739: return(0);
740: }
744: static PetscErrorCode THIInitialDMMG(DMMG dmmg,Vec X)
745: {
749: THIInitial((THI)dmmg->user,dmmg->dm,X);
750: return(0);
751: }
753: 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)
754: {
755: PetscInt l,ll;
756: PetscScalar gam;
758: du[0] = du[1] = du[2] = 0;
759: dv[0] = dv[1] = dv[2] = 0;
760: *u = 0;
761: *v = 0;
762: for (l=0; l<8; l++) {
763: *u += phi[l] * n[l].u;
764: *v += phi[l] * n[l].v;
765: for (ll=0; ll<3; ll++) {
766: du[ll] += dphi[l][ll] * n[l].u;
767: dv[ll] += dphi[l][ll] * n[l].v;
768: }
769: }
770: 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]);
771: THIViscosity(thi,PetscRealPart(gam),eta,deta);
772: }
776: static PetscErrorCode THIFunctionLocal_3D(DMDALocalInfo *info,const Node ***x,const PrmNode **prm,const Node ***xdot,Node ***f,THI thi)
777: {
778: PetscInt xs,ys,xm,ym,zm,i,j,k,q,l;
779: PetscReal hx,hy,etamin,etamax,beta2min,beta2max;
783: xs = info->zs;
784: ys = info->ys;
785: xm = info->zm;
786: ym = info->ym;
787: zm = info->xm;
788: hx = thi->Lx / info->mz;
789: hy = thi->Ly / info->my;
791: etamin = 1e100;
792: etamax = 0;
793: beta2min = 1e100;
794: beta2max = 0;
796: for (i=xs; i<xs+xm; i++) {
797: for (j=ys; j<ys+ym; j++) {
798: PrmNode pn[4],dpn[4][2];
799: QuadExtract(prm,i,j,pn);
800: QuadComputeGrad4(QuadQDeriv,hx,hy,pn,dpn);
801: for (k=0; k<zm-1; k++) {
802: PetscInt ls = 0;
803: Node n[8],ndot[8],*fn[8];
804: PetscReal zn[8],etabase = 0;
805: PrmHexGetZ(pn,k,zm,zn);
806: HexExtract(x,i,j,k,n);
807: HexExtract(xdot,i,j,k,ndot);
808: HexExtractRef(f,i,j,k,fn);
809: if (thi->no_slip && k == 0) {
810: for (l=0; l<4; l++) n[l].u = n[l].v = 0;
811: /* The first 4 basis functions lie on the bottom layer, so their contribution is exactly 0, hence we can skip them */
812: ls = 4;
813: }
814: for (q=0; q<8; q++) {
815: PetscReal dz[3],phi[8],dphi[8][3],jw,eta,deta;
816: PetscScalar du[3],dv[3],u,v,udot=0,vdot=0;
817: for (l=ls; l<8; l++) {
818: udot += HexQInterp[q][l]*ndot[l].u;
819: vdot += HexQInterp[q][l]*ndot[l].v;
820: }
821: HexGrad(HexQDeriv[q],zn,dz);
822: HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
823: PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
824: jw /= thi->rhog; /* scales residuals to be O(1) */
825: if (q == 0) etabase = eta;
826: RangeUpdate(&etamin,&etamax,eta);
827: for (l=ls; l<8; l++) { /* test functions */
828: const PetscScalar ds[2] = {dpn[q%4][0].h+dpn[q%4][0].b, dpn[q%4][1].h+dpn[q%4][1].b};
829: const PetscReal pp=phi[l],*dp = dphi[l];
830: 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];
831: 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];
832: fn[l]->u += pp*jw*udot*thi->inertia*pp;
833: fn[l]->v += pp*jw*vdot*thi->inertia*pp;
834: }
835: }
836: if (k == 0) { /* we are on a bottom face */
837: if (thi->no_slip) {
838: /* Note: Non-Galerkin coarse grid operators are very sensitive to the scaling of Dirichlet boundary
839: * conditions. After shenanigans above, etabase contains the effective viscosity at the closest quadrature
840: * point to the bed. We want the diagonal entry in the Dirichlet condition to have similar magnitude to the
841: * diagonal entry corresponding to the adjacent node. The fundamental scaling of the viscous part is in
842: * diagu, diagv below. This scaling is easy to recognize by considering the finite difference operator after
843: * scaling by element size. The no-slip Dirichlet condition is scaled by this factor, and also in the
844: * assembled matrix (see the similar block in THIJacobianLocal).
845: *
846: * Note that the residual at this Dirichlet node is linear in the state at this node, but also depends
847: * (nonlinearly in general) on the neighboring interior nodes through the local viscosity. This will make
848: * a matrix-free Jacobian have extra entries in the corresponding row. We assemble only the diagonal part,
849: * so the solution will exactly satisfy the boundary condition after the first linear iteration.
850: */
851: const PetscReal hz = PetscRealPart(pn[0].h)/(zm-1.);
852: 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);
853: fn[0]->u = thi->dirichlet_scale*diagu*x[i][j][k].u;
854: fn[0]->v = thi->dirichlet_scale*diagv*x[i][j][k].v;
855: } else { /* Integrate over bottom face to apply boundary condition */
856: 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) */
857: const PetscReal jw = 0.25*hx*hy,*phi = QuadQInterp[q];
858: PetscScalar u=0,v=0,rbeta2=0;
859: PetscReal beta2,dbeta2;
860: for (l=0; l<4; l++) {
861: u += phi[l]*n[l].u;
862: v += phi[l]*n[l].v;
863: rbeta2 += phi[l]*pn[l].beta2;
864: }
865: THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
866: RangeUpdate(&beta2min,&beta2max,beta2);
867: for (l=0; l<4; l++) {
868: const PetscReal pp = phi[l];
869: fn[ls+l]->u += pp*jw*beta2*u;
870: fn[ls+l]->v += pp*jw*beta2*v;
871: }
872: }
873: }
874: }
875: }
876: }
877: }
879: PRangeMinMax(&thi->eta,etamin,etamax);
880: PRangeMinMax(&thi->beta2,beta2min,beta2max);
881: return(0);
882: }
886: static PetscErrorCode THIFunctionLocal_2D(DMDALocalInfo *info,const Node ***x,const PrmNode **prm,const PrmNode **prmdot,PrmNode **f,THI thi)
887: {
888: PetscInt xs,ys,xm,ym,zm,i,j,k;
889: PetscReal hx,hy;
892: xs = info->zs;
893: ys = info->ys;
894: xm = info->zm;
895: ym = info->ym;
896: zm = info->xm;
897: hx = thi->Lx / info->mz;
898: hy = thi->Ly / info->my;
900: for (i=xs; i<xs+xm; i++) {
901: for (j=ys; j<ys+ym; j++) {
902: PetscScalar div = 0,erate,h[8];
903: PrmNodeGetFaceMeasure(prm,i,j,h);
904: for (k=0; k<zm; k++) {
905: PetscScalar weight = (k==0 || k == zm-1) ? 0.5/(zm-1) : 1.0/(zm-1);
906: if (0) { /* centered flux */
907: 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)
908: - 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)
909: + 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)
910: + 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)
911: - 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)
912: - 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)
913: + 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)
914: + 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));
915: } else { /* Upwind flux */
916: div += weight*(-UpwindFluxXW(x,prm,h,i,j,k, 1)
917: -UpwindFluxXW(x,prm,h,i,j,k,-1)
918: +UpwindFluxXE(x,prm,h,i,j,k, 1)
919: +UpwindFluxXE(x,prm,h,i,j,k,-1)
920: -UpwindFluxYS(x,prm,h,i,j,k, 1)
921: -UpwindFluxYS(x,prm,h,i,j,k,-1)
922: +UpwindFluxYN(x,prm,h,i,j,k, 1)
923: +UpwindFluxYN(x,prm,h,i,j,k,-1));
924: }
925: }
926: /* printf("div[%d][%d] %g\n",i,j,div); */
927: THIErosion(thi,&x[i][j][0],&erate,NULL);
928: f[i][j].b = prmdot[i][j].b - erate;
929: f[i][j].h = prmdot[i][j].h + div;
930: f[i][j].beta2 = prmdot[i][j].beta2;
931: }
932: }
933: return(0);
934: }
938: static PetscErrorCode THIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
939: {
941: THI thi = (THI)ctx;
942: DM pack,da3,da2;
943: Vec X3,X2,Xdot3,Xdot2,F3,F2,F3g,F2g;
944: const Node ***x3,***xdot3;
945: const PrmNode **x2,**xdot2;
946: Node ***f3;
947: PrmNode **f2;
948: DMDALocalInfo info3;
951: TSGetDM(ts,&pack);
952: DMCompositeGetEntries(pack,&da3,&da2);
953: DMDAGetLocalInfo(da3,&info3);
954: DMCompositeGetLocalVectors(pack,&X3,&X2);
955: DMCompositeGetLocalVectors(pack,&Xdot3,&Xdot2);
956: DMCompositeScatter(pack,X,X3,X2);
957: THIFixGhosts(thi,da3,da2,X3,X2);
958: DMCompositeScatter(pack,Xdot,Xdot3,Xdot2);
960: DMGetLocalVector(da3,&F3);
961: DMGetLocalVector(da2,&F2);
962: VecZeroEntries(F3);
964: DMDAVecGetArray(da3,X3,&x3);
965: DMDAVecGetArray(da2,X2,&x2);
966: DMDAVecGetArray(da3,Xdot3,&xdot3);
967: DMDAVecGetArray(da2,Xdot2,&xdot2);
968: DMDAVecGetArray(da3,F3,&f3);
969: DMDAVecGetArray(da2,F2,&f2);
971: THIFunctionLocal_3D(&info3,x3,x2,xdot3,f3,thi);
972: THIFunctionLocal_2D(&info3,x3,x2,xdot2,f2,thi);
974: DMDAVecRestoreArray(da3,X3,&x3);
975: DMDAVecRestoreArray(da2,X2,&x2);
976: DMDAVecRestoreArray(da3,Xdot3,&xdot3);
977: DMDAVecRestoreArray(da2,Xdot2,&xdot2);
978: DMDAVecRestoreArray(da3,F3,&f3);
979: DMDAVecRestoreArray(da2,F2,&f2);
981: DMCompositeRestoreLocalVectors(pack,&X3,&X2);
982: DMCompositeRestoreLocalVectors(pack,&Xdot3,&Xdot2);
984: VecZeroEntries(F);
985: DMCompositeGetAccess(pack,F,&F3g,&F2g);
986: DMLocalToGlobalBegin(da3,F3,ADD_VALUES,F3g);
987: DMLocalToGlobalEnd (da3,F3,ADD_VALUES,F3g);
988: DMLocalToGlobalBegin(da2,F2,INSERT_VALUES,F2g);
989: DMLocalToGlobalEnd (da2,F2,INSERT_VALUES,F2g);
991: if (thi->verbose) {
992: PetscViewer viewer;
993: PetscViewerASCIIGetStdout(((PetscObject)thi)->comm,&viewer);
994: PetscViewerASCIIPrintf(viewer,"3D_Velocity residual (bs=2):\n");
995: PetscViewerASCIIPushTab(viewer);
996: VecView(F3,viewer);
997: PetscViewerASCIIPopTab(viewer);
998: PetscViewerASCIIPrintf(viewer,"2D_Fields residual (bs=3):\n");
999: PetscViewerASCIIPushTab(viewer);
1000: VecView(F2,viewer);
1001: PetscViewerASCIIPopTab(viewer);
1002: }
1004: DMCompositeRestoreAccess(pack,F,&F3g,&F2g);
1006: DMRestoreLocalVector(da3,&F3);
1007: DMRestoreLocalVector(da2,&F2);
1009: return(0);
1010: }
1014: static PetscErrorCode THIMatrixStatistics(THI thi,Mat B,PetscViewer viewer)
1015: {
1017: PetscReal nrm;
1018: PetscInt m;
1019: PetscMPIInt rank;
1022: MatNorm(B,NORM_FROBENIUS,&nrm);
1023: MatGetSize(B,&m,0);
1024: MPI_Comm_rank(((PetscObject)B)->comm,&rank);
1025: if (!rank) {
1026: PetscScalar val0,val2;
1027: MatGetValue(B,0,0,&val0);
1028: MatGetValue(B,2,2,&val2);
1029: PetscViewerASCIIPrintf(viewer,"Matrix dim %8d norm %8.2e, (0,0) %8.2e (2,2) %8.2e, eta [%8.2e,%8.2e] beta2 [%8.2e,%8.2e]\n",m,nrm,PetscRealPart(val0),PetscRealPart(val2),thi->eta.cmin,thi->eta.cmax,thi->beta2.cmin,thi->beta2.cmax);
1030: }
1031: return(0);
1032: }
1036: static PetscErrorCode THISurfaceStatistics(DM pack,Vec X,PetscReal *min,PetscReal *max,PetscReal *mean)
1037: {
1039: DM da3,da2;
1040: Vec X3,X2;
1041: Node ***x;
1042: PetscInt i,j,xs,ys,zs,xm,ym,zm,mx,my,mz;
1043: PetscReal umin = 1e100,umax=-1e100;
1044: PetscScalar usum=0.0,gusum;
1047: DMCompositeGetEntries(pack,&da3,&da2);
1048: DMCompositeGetAccess(pack,X,&X3,&X2);
1049: *min = *max = *mean = 0;
1050: DMDAGetInfo(da3,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
1051: DMDAGetCorners(da3,&zs,&ys,&xs,&zm,&ym,&xm);
1052: if (zs != 0 || zm != mz) SETERRQ(PETSC_COMM_SELF,1,"Unexpected decomposition");
1053: DMDAVecGetArray(da3,X3,&x);
1054: for (i=xs; i<xs+xm; i++) {
1055: for (j=ys; j<ys+ym; j++) {
1056: PetscReal u = PetscRealPart(x[i][j][zm-1].u);
1057: RangeUpdate(&umin,&umax,u);
1058: usum += u;
1059: }
1060: }
1061: DMDAVecRestoreArray(da3,X3,&x);
1062: DMCompositeRestoreAccess(pack,X,&X3,&X2);
1064: MPI_Allreduce(&umin,min,1,MPIU_REAL,MPIU_MIN,((PetscObject)da3)->comm);
1065: MPI_Allreduce(&umax,max,1,MPIU_REAL,MPIU_MAX,((PetscObject)da3)->comm);
1066: MPI_Allreduce(&usum,&gusum,1,MPIU_SCALAR,MPIU_SUM,((PetscObject)da3)->comm);
1067: *mean = PetscRealPart(gusum) / (mx*my);
1068: return(0);
1069: }
1073: static PetscErrorCode THISolveStatistics(THI thi,DMMG *dmmg,PetscInt coarsened,const char name[])
1074: {
1075: MPI_Comm comm = ((PetscObject)thi)->comm;
1076: PetscInt nlevels = DMMGGetLevels(dmmg),level = nlevels-1-coarsened;
1077: SNES snes = dmmg[level]->snes;
1078: DM pack = dmmg[level]->dm;
1079: Vec X = dmmg[level]->x,X3,X2;
1083: DMCompositeGetAccess(pack,X,&X3,&X2);
1084: PetscPrintf(comm,"Solution statistics after solve: %s\n",name);
1085: {
1086: PetscInt its,lits;
1087: SNESConvergedReason reason;
1088: SNESGetIterationNumber(snes,&its);
1089: SNESGetConvergedReason(snes,&reason);
1090: SNESGetLinearSolveIterations(snes,&lits);
1091: PetscPrintf(comm,"%s: Number of Newton iterations = %d, total linear iterations = %d\n",SNESConvergedReasons[reason],its,lits);
1092: }
1093: {
1094: PetscReal nrm2,tmin[3]={1e100,1e100,1e100},tmax[3]={-1e100,-1e100,-1e100},min[3],max[3];
1095: PetscInt i,j,m;
1096: PetscScalar *x;
1097: VecNorm(X3,NORM_2,&nrm2);
1098: VecGetLocalSize(X3,&m);
1099: VecGetArray(X3,&x);
1100: for (i=0; i<m; i+=2) {
1101: PetscReal u = PetscRealPart(x[i]),v = PetscRealPart(x[i+1]),c = sqrt(u*u+v*v);
1102: tmin[0] = PetscMin(u,tmin[0]);
1103: tmin[1] = PetscMin(v,tmin[1]);
1104: tmin[2] = PetscMin(c,tmin[2]);
1105: tmax[0] = PetscMax(u,tmax[0]);
1106: tmax[1] = PetscMax(v,tmax[1]);
1107: tmax[2] = PetscMax(c,tmax[2]);
1108: }
1109: VecRestoreArray(X,&x);
1110: MPI_Allreduce(tmin,min,3,MPIU_REAL,MPIU_MIN,((PetscObject)thi)->comm);
1111: MPI_Allreduce(tmax,max,3,MPIU_REAL,MPIU_MAX,((PetscObject)thi)->comm);
1112: /* Dimensionalize to meters/year */
1113: nrm2 *= thi->units->year / thi->units->meter;
1114: for (j=0; j<3; j++) {
1115: min[j] *= thi->units->year / thi->units->meter;
1116: max[j] *= thi->units->year / thi->units->meter;
1117: }
1118: PetscPrintf(comm,"|X|_2 %g u in [%g, %g] v in [%g, %g] c in [%g, %g] \n",nrm2,min[0],max[0],min[1],max[1],min[2],max[2]);
1119: {
1120: PetscReal umin,umax,umean;
1121: THISurfaceStatistics(dmmg[level]->dm,X,&umin,&umax,&umean);
1122: umin *= thi->units->year / thi->units->meter;
1123: umax *= thi->units->year / thi->units->meter;
1124: umean *= thi->units->year / thi->units->meter;
1125: PetscPrintf(comm,"Surface statistics: u in [%12.6e, %12.6e] mean %12.6e\n",umin,umax,umean);
1126: }
1127: /* These values stay nondimensional */
1128: PetscPrintf(comm,"Global eta range [%g, %g], converged range [%g, %g]\n",thi->eta.min,thi->eta.max,thi->eta.cmin,thi->eta.cmax);
1129: PetscPrintf(comm,"Global beta2 range [%g, %g], converged range [%g, %g]\n",thi->beta2.min,thi->beta2.max,thi->beta2.cmin,thi->beta2.cmax);
1130: }
1131: PetscPrintf(comm,"\n");
1132: DMCompositeRestoreAccess(pack,X,&X3,&X2);
1133: return(0);
1134: }
1136: static inline PetscInt DMDALocalIndex3D(DMDALocalInfo *info,PetscInt i,PetscInt j,PetscInt k)
1137: {return ((i-info->gzs)*info->gym + (j-info->gys))*info->gxm + (k-info->gxs);}
1138: static inline PetscInt DMDALocalIndex2D(DMDALocalInfo *info,PetscInt i,PetscInt j)
1139: {return (i-info->gzs)*info->gym + (j-info->gys);}
1143: static PetscErrorCode THIJacobianLocal_Momentum(DMDALocalInfo *info,const Node ***x,const PrmNode **prm,Mat B,Mat Bcpl,THI thi)
1144: {
1145: PetscInt xs,ys,xm,ym,zm,i,j,k,q,l,ll;
1146: PetscReal hx,hy;
1150: xs = info->zs;
1151: ys = info->ys;
1152: xm = info->zm;
1153: ym = info->ym;
1154: zm = info->xm;
1155: hx = thi->Lx / info->mz;
1156: hy = thi->Ly / info->my;
1158: for (i=xs; i<xs+xm; i++) {
1159: for (j=ys; j<ys+ym; j++) {
1160: PrmNode pn[4],dpn[4][2];
1161: QuadExtract(prm,i,j,pn);
1162: QuadComputeGrad4(QuadQDeriv,hx,hy,pn,dpn);
1163: for (k=0; k<zm-1; k++) {
1164: Node n[8];
1165: PetscReal zn[8],etabase = 0;
1166: PetscScalar Ke[8*NODE_SIZE][8*NODE_SIZE],Kcpl[8*NODE_SIZE][4*PRMNODE_SIZE];
1167: PetscInt ls = 0;
1169: PrmHexGetZ(pn,k,zm,zn);
1170: HexExtract(x,i,j,k,n);
1171: PetscMemzero(Ke,sizeof(Ke));
1172: PetscMemzero(Kcpl,sizeof(Kcpl));
1173: if (thi->no_slip && k == 0) {
1174: for (l=0; l<4; l++) n[l].u = n[l].v = 0;
1175: ls = 4;
1176: }
1177: for (q=0; q<8; q++) {
1178: PetscReal dz[3],phi[8],dphi[8][3],jw,eta,deta;
1179: PetscScalar du[3],dv[3],u,v;
1180: HexGrad(HexQDeriv[q],zn,dz);
1181: HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
1182: PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
1183: jw /= thi->rhog; /* residuals are scaled by this factor */
1184: if (q == 0) etabase = eta;
1185: for (l=ls; l<8; l++) { /* test functions */
1186: const PetscReal pp=phi[l],*restrict dp = dphi[l];
1187: for (ll=ls; ll<8; ll++) { /* trial functions */
1188: const PetscReal *restrict dpl = dphi[ll];
1189: PetscScalar dgdu,dgdv;
1190: dgdu = 2.*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1]+dv[0])*dpl[1] + 0.5*du[2]*dpl[2];
1191: dgdv = 2.*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1]+dv[0])*dpl[0] + 0.5*dv[2]*dpl[2];
1192: /* Picard part */
1193: 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];
1194: Ke[l*2+0][ll*2+1] += dp[0]*jw*eta*2.*dpl[1] + dp[1]*jw*eta*dpl[0];
1195: Ke[l*2+1][ll*2+0] += dp[1]*jw*eta*2.*dpl[0] + dp[0]*jw*eta*dpl[1];
1196: 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];
1197: /* extra Newton terms */
1198: 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];
1199: 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];
1200: 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];
1201: 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];
1202: /* inertial part */
1203: Ke[l*2+0][ll*2+0] += pp*jw*thi->inertia*pp;
1204: Ke[l*2+1][ll*2+1] += pp*jw*thi->inertia*pp;
1205: }
1206: for (ll=0; ll<4; ll++) { /* Trial functions for surface/bed */
1207: const PetscReal dpl[] = {QuadQDeriv[q%4][ll][0]/hx, QuadQDeriv[q%4][ll][1]/hy}; /* surface = h + b */
1208: Kcpl[FieldIndex(Node,l,u)][FieldIndex(PrmNode,ll,h)] += pp*jw*thi->rhog*dpl[0];
1209: Kcpl[FieldIndex(Node,l,u)][FieldIndex(PrmNode,ll,b)] += pp*jw*thi->rhog*dpl[0];
1210: Kcpl[FieldIndex(Node,l,v)][FieldIndex(PrmNode,ll,h)] += pp*jw*thi->rhog*dpl[1];
1211: Kcpl[FieldIndex(Node,l,v)][FieldIndex(PrmNode,ll,b)] += pp*jw*thi->rhog*dpl[1];
1212: }
1213: }
1214: }
1215: if (k == 0) { /* on a bottom face */
1216: if (thi->no_slip) {
1217: const PetscReal hz = PetscRealPart(pn[0].h)/(zm-1);
1218: 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);
1219: Ke[0][0] = thi->dirichlet_scale*diagu;
1220: Ke[0][1] = 0;
1221: Ke[1][0] = 0;
1222: Ke[1][1] = thi->dirichlet_scale*diagv;
1223: } else {
1224: for (q=0; q<4; q++) { /* We remove the explicit scaling by 1/rhog because beta2 already has that scaling to be O(1) */
1225: const PetscReal jw = 0.25*hx*hy,*phi = QuadQInterp[q];
1226: PetscScalar u=0,v=0,rbeta2=0;
1227: PetscReal beta2,dbeta2;
1228: for (l=0; l<4; l++) {
1229: u += phi[l]*n[l].u;
1230: v += phi[l]*n[l].v;
1231: rbeta2 += phi[l]*pn[l].beta2;
1232: }
1233: THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
1234: for (l=0; l<4; l++) {
1235: const PetscReal pp = phi[l];
1236: for (ll=0; ll<4; ll++) {
1237: const PetscReal ppl = phi[ll];
1238: Ke[l*2+0][ll*2+0] += pp*jw*beta2*ppl + pp*jw*dbeta2*u*u*ppl;
1239: Ke[l*2+0][ll*2+1] += pp*jw*dbeta2*u*v*ppl;
1240: Ke[l*2+1][ll*2+0] += pp*jw*dbeta2*v*u*ppl;
1241: Ke[l*2+1][ll*2+1] += pp*jw*beta2*ppl + pp*jw*dbeta2*v*v*ppl;
1242: }
1243: }
1244: }
1245: }
1246: }
1247: {
1248: const PetscInt rc3blocked[8] = {
1249: DMDALocalIndex3D(info,i+0,j+0,k+0),
1250: DMDALocalIndex3D(info,i+1,j+0,k+0),
1251: DMDALocalIndex3D(info,i+1,j+1,k+0),
1252: DMDALocalIndex3D(info,i+0,j+1,k+0),
1253: DMDALocalIndex3D(info,i+0,j+0,k+1),
1254: DMDALocalIndex3D(info,i+1,j+0,k+1),
1255: DMDALocalIndex3D(info,i+1,j+1,k+1),
1256: DMDALocalIndex3D(info,i+0,j+1,k+1)
1257: },col2blocked[PRMNODE_SIZE*4] = {
1258: DMDALocalIndex2D(info,i+0,j+0),
1259: DMDALocalIndex2D(info,i+1,j+0),
1260: DMDALocalIndex2D(info,i+1,j+1),
1261: DMDALocalIndex2D(info,i+0,j+1)
1262: };
1263: #if !defined COMPUTE_LOWER_TRIANGULAR /* fill in lower-triangular part, this is really cheap compared to computing the entries */
1264: for (l=0; l<8; l++) {
1265: for (ll=l+1; ll<8; ll++) {
1266: Ke[ll*2+0][l*2+0] = Ke[l*2+0][ll*2+0];
1267: Ke[ll*2+1][l*2+0] = Ke[l*2+0][ll*2+1];
1268: Ke[ll*2+0][l*2+1] = Ke[l*2+1][ll*2+0];
1269: Ke[ll*2+1][l*2+1] = Ke[l*2+1][ll*2+1];
1270: }
1271: }
1272: #endif
1273: MatSetValuesBlockedLocal(B,8,rc3blocked,8,rc3blocked,&Ke[0][0],ADD_VALUES); /* velocity-velocity coupling can use blocked insertion */
1274: { /* The off-diagonal part cannot (yet) */
1275: PetscInt row3scalar[NODE_SIZE*8],col2scalar[PRMNODE_SIZE*4];
1276: for (l=0; l<8; l++) for (ll=0; ll<NODE_SIZE; ll++)
1277: row3scalar[l*NODE_SIZE+ll] = rc3blocked[l]*NODE_SIZE+ll;
1278: for (l=0; l<4; l++) for (ll=0; ll<PRMNODE_SIZE; ll++)
1279: col2scalar[l*PRMNODE_SIZE+ll] = col2blocked[l]*PRMNODE_SIZE+ll;
1280: MatSetValuesLocal(Bcpl,8*NODE_SIZE,row3scalar,4*PRMNODE_SIZE,col2scalar,&Kcpl[0][0],ADD_VALUES);
1281: }
1282: }
1283: }
1284: }
1285: }
1286: return(0);
1287: }
1291: static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo *info,const Node ***x3,const PrmNode **x2,const PrmNode **xdot2,PetscReal a,Mat B22,Mat B21,THI thi)
1292: {
1294: PetscInt xs,ys,xm,ym,zm,i,j,k;
1295: PetscReal hx,hy;
1298: xs = info->zs;
1299: ys = info->ys;
1300: xm = info->zm;
1301: ym = info->ym;
1302: zm = info->xm;
1303: hx = thi->Lx / info->mz;
1304: hy = thi->Ly / info->my;
1306: if (zm > 1024) SETERRQ(((PetscObject)info->da)->comm,PETSC_ERR_SUP,"Need to allocate more space");
1307: for (i=xs; i<xs+xm; i++) {
1308: for (j=ys; j<ys+ym; j++) {
1309: { /* Self-coupling */
1310: const PetscInt row[] = {DMDALocalIndex2D(info,i,j)};
1311: const PetscInt col[] = {DMDALocalIndex2D(info,i,j)};
1312: const PetscScalar vals[] = {
1313: a,0,0,
1314: 0,a,0,
1315: 0,0,a
1316: };
1317: MatSetValuesBlockedLocal(B22,1,row,1,col,vals,INSERT_VALUES);
1318: }
1319: for (k=0; k<zm; k++) { /* Coupling to velocity problem */
1320: /* Use a cheaper quadrature than for residual evaluation, because it is much sparser */
1321: const PetscInt row[] = {FieldIndex(PrmNode,DMDALocalIndex2D(info,i,j),h)};
1322: const PetscInt cols[] = {
1323: FieldIndex(Node,DMDALocalIndex3D(info,i-1,j,k),u),
1324: FieldIndex(Node,DMDALocalIndex3D(info,i ,j,k),u),
1325: FieldIndex(Node,DMDALocalIndex3D(info,i+1,j,k),u),
1326: FieldIndex(Node,DMDALocalIndex3D(info,i,j-1,k),v),
1327: FieldIndex(Node,DMDALocalIndex3D(info,i,j ,k),v),
1328: FieldIndex(Node,DMDALocalIndex3D(info,i,j+1,k),v)
1329: };
1330: const PetscScalar
1331: w = (k && k<zm-1) ? 0.5 : 0.25,
1332: hW = w*(x2[i-1][j ].h+x2[i ][j ].h)/(zm-1.),
1333: hE = w*(x2[i ][j ].h+x2[i+1][j ].h)/(zm-1.),
1334: hS = w*(x2[i ][j-1].h+x2[i ][j ].h)/(zm-1.),
1335: hN = w*(x2[i ][j ].h+x2[i ][j+1].h)/(zm-1.);
1336: PetscScalar *vals,
1337: vals_upwind[] = {((PetscRealPart(x3[i][j][k].u) > 0) ? -hW : 0),
1338: ((PetscRealPart(x3[i][j][k].u) > 0) ? +hE : -hW),
1339: ((PetscRealPart(x3[i][j][k].u) > 0) ? 0 : +hE),
1340: ((PetscRealPart(x3[i][j][k].v) > 0) ? -hS : 0),
1341: ((PetscRealPart(x3[i][j][k].v) > 0) ? +hN : -hS),
1342: ((PetscRealPart(x3[i][j][k].v) > 0) ? 0 : +hN)},
1343: vals_centered[] = {-0.5*hW, 0.5*(-hW+hE), 0.5*hE,
1344: -0.5*hS, 0.5*(-hS+hN), 0.5*hN};
1345: vals = 1 ? vals_upwind : vals_centered;
1346: if (k == 0) {
1347: Node derate;
1348: THIErosion(thi,&x3[i][j][0],NULL,&derate);
1349: vals[1] -= derate.u;
1350: vals[4] -= derate.v;
1351: }
1352: MatSetValuesLocal(B21,1,row,6,cols,vals,INSERT_VALUES);
1353: }
1354: }
1355: }
1356: return(0);
1357: }
1361: static PetscErrorCode THIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *A,Mat *B,MatStructure *mstr,void *ctx)
1362: {
1364: THI thi = (THI)ctx;
1365: DM pack,da3,da2;
1366: Vec X3,X2,Xdot2;
1367: Mat B11,B12,B21,B22;
1368: DMDALocalInfo info3;
1369: IS *isloc;
1370: const Node ***x3;
1371: const PrmNode **x2,**xdot2;
1374: TSGetDM(ts,&pack);
1375: DMCompositeGetEntries(pack,&da3,&da2);
1376: DMDAGetLocalInfo(da3,&info3);
1377: DMCompositeGetLocalVectors(pack,&X3,&X2);
1378: DMCompositeGetLocalVectors(pack,0,&Xdot2);
1379: DMCompositeScatter(pack,X,X3,X2);
1380: THIFixGhosts(thi,da3,da2,X3,X2);
1381: DMCompositeScatter(pack,Xdot,0,Xdot2);
1383: MatZeroEntries(*B);
1385: DMCompositeGetLocalISs(pack,&isloc);
1386: MatGetLocalSubMatrix(*B,isloc[0],isloc[0],&B11);
1387: MatGetLocalSubMatrix(*B,isloc[0],isloc[1],&B12);
1388: MatGetLocalSubMatrix(*B,isloc[1],isloc[0],&B21);
1389: MatGetLocalSubMatrix(*B,isloc[1],isloc[1],&B22);
1391: DMDAVecGetArray(da3,X3,&x3);
1392: DMDAVecGetArray(da2,X2,&x2);
1393: DMDAVecGetArray(da2,Xdot2,&xdot2);
1395: THIJacobianLocal_Momentum(&info3,x3,x2,B11,B12,thi);
1397: /* Need to switch from ADD_VALUES to INSERT_VALUES */
1398: MatAssemblyBegin(*B,MAT_FLUSH_ASSEMBLY);
1399: MatAssemblyEnd(*B,MAT_FLUSH_ASSEMBLY);
1401: THIJacobianLocal_2D(&info3,x3,x2,xdot2,a,B22,B21,thi);
1403: DMDAVecRestoreArray(da3,X3,&x3);
1404: DMDAVecRestoreArray(da2,X2,&x2);
1405: DMDAVecRestoreArray(da2,Xdot2,&xdot2);
1407: MatRestoreLocalSubMatrix(*B,isloc[0],isloc[0],&B11);
1408: MatRestoreLocalSubMatrix(*B,isloc[0],isloc[1],&B12);
1409: MatRestoreLocalSubMatrix(*B,isloc[1],isloc[0],&B21);
1410: MatRestoreLocalSubMatrix(*B,isloc[1],isloc[1],&B22);
1411: ISDestroy(&isloc[0]);
1412: ISDestroy(&isloc[1]);
1413: PetscFree(isloc);
1415: DMCompositeRestoreLocalVectors(pack,&X3,&X2);
1416: DMCompositeRestoreLocalVectors(pack,0,&Xdot2);
1418: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1419: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1420: if (*A != *B) {
1421: MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1422: MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1423: }
1424: *mstr = SAME_NONZERO_PATTERN;
1425: if (thi->verbose) {THIMatrixStatistics(thi,*B,PETSC_VIEWER_STDOUT_WORLD);}
1426: return(0);
1427: }
1431: /* VTK's XML formats are so brain-dead that they can't handle multiple grids in the same file. Since the communication
1432: * can be shared between the two grids, we write two files at once, one for velocity (living on a 3D grid defined by
1433: * h=thickness and b=bed) and another for all properties living on the 2D grid.
1434: */
1435: static PetscErrorCode THIDAVecView_VTK_XML(THI thi,DM pack,Vec X,const char filename[],const char filename2[])
1436: {
1437: const PetscInt dof = NODE_SIZE,dof2 = PRMNODE_SIZE;
1438: Units units = thi->units;
1439: MPI_Comm comm;
1441: PetscViewer viewer3,viewer2;
1442: PetscMPIInt rank,size,tag,nn,nmax,nn2,nmax2;
1443: PetscInt mx,my,mz,r,range[6];
1444: PetscScalar *x,*x2;
1445: DM da3,da2;
1446: Vec X3,X2;
1449: comm = ((PetscObject)thi)->comm;
1450: DMCompositeGetEntries(pack,&da3,&da2);
1451: DMCompositeGetAccess(pack,X,&X3,&X2);
1452: DMDAGetInfo(da3,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
1453: MPI_Comm_size(comm,&size);
1454: MPI_Comm_rank(comm,&rank);
1455: PetscViewerASCIIOpen(comm,filename,&viewer3);
1456: PetscViewerASCIIOpen(comm,filename2,&viewer2);
1457: PetscViewerASCIIPrintf(viewer3,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1458: PetscViewerASCIIPrintf(viewer2,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1459: PetscViewerASCIIPrintf(viewer3," <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n",0,mz-1,0,my-1,0,mx-1);
1460: PetscViewerASCIIPrintf(viewer2," <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n",0,0,0,my-1,0,mx-1);
1462: DMDAGetCorners(da3,range,range+1,range+2,range+3,range+4,range+5);
1463: nn = PetscMPIIntCast(range[3]*range[4]*range[5]*dof);
1464: MPI_Reduce(&nn,&nmax,1,MPI_INT,MPI_MAX,0,comm);
1465: nn2 = PetscMPIIntCast(range[4]*range[5]*dof2);
1466: MPI_Reduce(&nn2,&nmax2,1,MPI_INT,MPI_MAX,0,comm);
1467: tag = ((PetscObject)viewer3)->tag;
1468: VecGetArray(X3,&x);
1469: VecGetArray(X2,&x2);
1470: if (!rank) {
1471: PetscScalar *array,*array2;
1472: PetscMalloc2(nmax,PetscScalar,&array,nmax2,PetscScalar,&array2);
1473: for (r=0; r<size; r++) {
1474: PetscInt i,j,k,f,xs,xm,ys,ym,zs,zm;
1475: Node *y3;
1476: PetscScalar (*y2)[PRMNODE_SIZE];
1477: MPI_Status status;
1478: if (r) {
1479: MPI_Recv(range,6,MPIU_INT,r,tag,comm,MPI_STATUS_IGNORE);
1480: }
1481: zs = range[0];ys = range[1];xs = range[2];zm = range[3];ym = range[4];xm = range[5];
1482: if (xm*ym*zm*dof > nmax) SETERRQ(PETSC_COMM_SELF,1,"should not happen");
1483: if (r) {
1484: MPI_Recv(array,nmax,MPIU_SCALAR,r,tag,comm,&status);
1485: MPI_Get_count(&status,MPIU_SCALAR,&nn);
1486: if (nn != xm*ym*zm*dof) SETERRQ(PETSC_COMM_SELF,1,"corrupt da3 send");
1487: y3 = (Node*)array;
1488: MPI_Recv(array2,nmax2,MPIU_SCALAR,r,tag,comm,&status);
1489: MPI_Get_count(&status,MPIU_SCALAR,&nn2);
1490: if (nn2 != xm*ym*dof2) SETERRQ(PETSC_COMM_SELF,1,"corrupt da2 send");
1491: y2 = (PetscScalar(*)[PRMNODE_SIZE])array2;
1492: } else {
1493: y3 = (Node*)x;
1494: y2 = (PetscScalar(*)[PRMNODE_SIZE])x2;
1495: }
1496: PetscViewerASCIIPrintf(viewer3," <Piece Extent=\"%d %d %d %d %d %d\">\n",zs,zs+zm-1,ys,ys+ym-1,xs,xs+xm-1);
1497: PetscViewerASCIIPrintf(viewer2," <Piece Extent=\"%d %d %d %d %d %d\">\n",0,0,ys,ys+ym-1,xs,xs+xm-1);
1499: PetscViewerASCIIPrintf(viewer3," <Points>\n");
1500: PetscViewerASCIIPrintf(viewer2," <Points>\n");
1501: PetscViewerASCIIPrintf(viewer3," <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n");
1502: PetscViewerASCIIPrintf(viewer2," <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n");
1503: for (i=xs; i<xs+xm; i++) {
1504: for (j=ys; j<ys+ym; j++) {
1505: PetscReal
1506: xx = thi->Lx*i/mx,
1507: yy = thi->Ly*j/my,
1508: b = PetscRealPart(y2[i*ym+j][FieldOffset(PrmNode,b)]),
1509: h = PetscRealPart(y2[i*ym+j][FieldOffset(PrmNode,h)]);
1510: for (k=zs; k<zs+zm; k++) {
1511: PetscReal zz = b + h*k/(mz-1.);
1512: PetscViewerASCIIPrintf(viewer3,"%f %f %f\n",xx,yy,zz);
1513: }
1514: PetscViewerASCIIPrintf(viewer2,"%f %f %f\n",xx,yy,(double)0.0);
1515: }
1516: }
1517: PetscViewerASCIIPrintf(viewer3," </DataArray>\n");
1518: PetscViewerASCIIPrintf(viewer2," </DataArray>\n");
1519: PetscViewerASCIIPrintf(viewer3," </Points>\n");
1520: PetscViewerASCIIPrintf(viewer2," </Points>\n");
1522: { /* Velocity and rank (3D) */
1523: PetscViewerASCIIPrintf(viewer3," <PointData>\n");
1524: PetscViewerASCIIPrintf(viewer3," <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n");
1525: for (i=0; i<nn/dof; i++) {
1526: PetscViewerASCIIPrintf(viewer3,"%f %f %f\n",PetscRealPart(y3[i].u)*units->year/units->meter,PetscRealPart(y3[i].v)*units->year/units->meter,0.0);
1527: }
1528: PetscViewerASCIIPrintf(viewer3," </DataArray>\n");
1530: PetscViewerASCIIPrintf(viewer3," <DataArray type=\"Int32\" Name=\"rank\" NumberOfComponents=\"1\" format=\"ascii\">\n");
1531: for (i=0; i<nn; i+=dof) {
1532: PetscViewerASCIIPrintf(viewer3,"%d\n",r);
1533: }
1534: PetscViewerASCIIPrintf(viewer3," </DataArray>\n");
1535: PetscViewerASCIIPrintf(viewer3," </PointData>\n");
1536: }
1538: { /* 2D */
1539: PetscViewerASCIIPrintf(viewer2," <PointData>\n");
1540: for (f=0; f<PRMNODE_SIZE; f++) {
1541: const char *fieldname;
1542: DMDAGetFieldName(da2,f,&fieldname);
1543: PetscViewerASCIIPrintf(viewer2," <DataArray type=\"Float32\" Name=\"%s\" format=\"ascii\">\n",fieldname);
1544: for (i=0; i<nn2/PRMNODE_SIZE; i++) {
1545: PetscViewerASCIIPrintf(viewer2,"%g\n",y2[i][f]);
1546: }
1547: PetscViewerASCIIPrintf(viewer2," </DataArray>\n");
1548: }
1549: PetscViewerASCIIPrintf(viewer2," </PointData>\n");
1550: }
1552: PetscViewerASCIIPrintf(viewer3," </Piece>\n");
1553: PetscViewerASCIIPrintf(viewer2," </Piece>\n");
1554: }
1555: PetscFree2(array,array2);
1556: } else {
1557: MPI_Send(range,6,MPIU_INT,0,tag,comm);
1558: MPI_Send(x,nn,MPIU_SCALAR,0,tag,comm);
1559: MPI_Send(x2,nn2,MPIU_SCALAR,0,tag,comm);
1560: }
1561: VecRestoreArray(X3,&x);
1562: VecRestoreArray(X2,&x2);
1563: PetscViewerASCIIPrintf(viewer3," </StructuredGrid>\n");
1564: PetscViewerASCIIPrintf(viewer2," </StructuredGrid>\n");
1566: DMCompositeRestoreAccess(pack,X,&X3,&X2);
1567: PetscViewerASCIIPrintf(viewer3,"</VTKFile>\n");
1568: PetscViewerASCIIPrintf(viewer2,"</VTKFile>\n");
1569: PetscViewerDestroy(&viewer3);
1570: PetscViewerDestroy(&viewer2);
1571: return(0);
1572: }
1576: static PetscErrorCode THITSMonitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
1577: {
1579: THI thi = (THI)ctx;
1580: DM pack;
1581: char filename3[PETSC_MAX_PATH_LEN],filename2[PETSC_MAX_PATH_LEN];
1584: PetscPrintf(((PetscObject)ts)->comm,"%3D: t=%G\n",step,t);
1585: if (thi->monitor_interval && step % thi->monitor_interval) return(0);
1586: TSGetDM(ts,&pack);
1587: PetscSNPrintf(filename3,sizeof filename3,"%s-3d-%03d.vts",thi->monitor_basename,step);
1588: PetscSNPrintf(filename2,sizeof filename2,"%s-2d-%03d.vts",thi->monitor_basename,step);
1589: THIDAVecView_VTK_XML(thi,pack,X,filename3,filename2);
1590: return(0);
1591: }
1596: static PetscErrorCode THICreateDM3d(THI thi,DM *dm3d)
1597: {
1598: MPI_Comm comm = ((PetscObject)thi)->comm;
1599: PetscInt M = 3,N = 3,P = 2;
1600: DM da;
1604: PetscOptionsBegin(comm,NULL,"Grid resolution options","");
1605: {
1606: PetscOptionsInt("-M","Number of elements in x-direction on coarse level","",M,&M,NULL);
1607: N = M;
1608: PetscOptionsInt("-N","Number of elements in y-direction on coarse level (if different from M)","",N,&N,NULL);
1609: PetscOptionsInt("-P","Number of elements in z-direction on coarse level","",P,&P,NULL);
1610: }
1611: PetscOptionsEnd();
1612: DMDACreate3d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,P,N,M,1,PETSC_DETERMINE,PETSC_DETERMINE,sizeof(Node)/sizeof(PetscScalar),1,0,0,0,&da);
1613: DMDASetFieldName(da,0,"x-velocity");
1614: DMDASetFieldName(da,1,"y-velocity");
1615: *dm3d = da;
1616: return(0);
1617: }
1619: #if 0
1622: static PetscErrorCode SNESSolve_THI_DMMG(SNES snes)
1623: {
1624: PetscErrorCode ierr,(*realsolve)(SNES);
1625: DMMG *dmmg;
1628: PetscObjectQuery((PetscObject)snes,"DMMG",(PetscObject*)&dmmg);
1629: if (!dmmg) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONG,"No DMMG attached to this SNES");
1630: if (snes != DMMGGetSNES(dmmg)) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"DMMG does not have attached SNES");
1631: if (!DMMGGetx(dmmg)) {VecDuplicate(X,&dmmg->x);}
1632: VecCopy(X,DMMGGetx(dmmg));
1634: PetscObjectQueryFunction((PetscObject)dmmg,"SNESSolve_Real",(void(**)(void))&realsolve);
1635: if (!realsolve) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_INCOMP,"No cache solve function pointer");
1636: snes->ops->solve = realsolve;
1638: DMMGSolve(dmmg);
1639: snes->ops->solve = &SNESSolve_THI_DMMG;
1640: return(0);
1641: }
1642: #endif
1646: int main(int argc,char *argv[])
1647: {
1648: MPI_Comm comm;
1649: DM pack,da3,da2;
1650: DMMG *dmmg;
1651: TS ts;
1652: THI thi;
1653: Vec X;
1654: Mat B;
1655: PetscInt i,steps;
1656: PetscReal ftime;
1659: PetscInitialize(&argc,&argv,0,help);
1660: comm = PETSC_COMM_WORLD;
1662: THICreate(comm,&thi);
1663: THICreateDM3d(thi,&da3);
1664: {
1665: PetscInt Mx,My,mx,my,s;
1666: DMDAStencilType st;
1667: DMDAGetInfo(da3,0, 0,&My,&Mx, 0,&my,&mx, 0,&s,0,0,0,&st);
1668: DMDACreate2d(((PetscObject)thi)->comm,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,st,My,Mx,my,mx,sizeof(PrmNode)/sizeof(PetscScalar),s,0,0,&da2);
1669: }
1671: PetscObjectSetName((PetscObject)da3,"3D_Velocity");
1672: DMSetOptionsPrefix(da3,"f3d_");
1673: DMDASetFieldName(da3,0,"u");
1674: DMDASetFieldName(da3,1,"v");
1675: PetscObjectSetName((PetscObject)da2,"2D_Fields");
1676: DMSetOptionsPrefix(da2,"f2d_");
1677: DMDASetFieldName(da2,0,"b");
1678: DMDASetFieldName(da2,1,"h");
1679: DMDASetFieldName(da2,2,"beta2");
1680: DMCompositeCreate(comm,&pack);
1681: DMCompositeAddDM(pack,da3);
1682: DMCompositeAddDM(pack,da2);
1683: DMDestroy(&da3);
1684: DMDestroy(&da2);
1685: DMSetUp(pack);
1686: DMGetMatrix(pack,PETSC_NULL,&B);
1687: MatSetOptionsPrefix(B,"thi_");
1689: DMMGCreate(comm,thi->nlevels,thi,&dmmg);
1690: DMMGSetDM(dmmg,pack);
1692: for (i=0; i<thi->nlevels; i++) {
1693: PetscReal Lx = thi->Lx / thi->units->meter,Ly = thi->Ly / thi->units->meter,Lz = thi->Lz / thi->units->meter;
1694: PetscInt Mx,My,Mz;
1695: DMCompositeGetEntries(pack,&da3,&da2);
1696: DMDAGetInfo(da3,0, &Mz,&My,&Mx, 0,0,0, 0,0,0,0,0,0);
1697: PetscPrintf(((PetscObject)thi)->comm,"Level %d 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));
1698: }
1699: DMMGSetInitialGuess(dmmg,THIInitialDMMG);
1701: {
1702: /* Use the user-defined matrix type on all but the coarse level */
1703: DMMGSetMatType(dmmg,thi->mattype);
1704: /* PCREDUNDANT only works with AIJ, and so do the third-party direct solvers. So when running in parallel, we can't
1705: * use the faster (S)BAIJ formats on the coarse level. */
1706: PetscFree(dmmg[0]->mtype);
1707: PetscStrallocpy(MATAIJ,&dmmg[0]->mtype);
1708: }
1709: //DMMGSetSNES(dmmg,THIFunction,PETSC_NULL);
1710: //MatSetOptionsPrefix(DMMGGetB(dmmg),"thi_");
1711: //DMMGSetFromOptions(dmmg);
1713: DMCreateGlobalVector(pack,&X);
1714: THIInitial(thi,pack,X);
1716: TSCreate(comm,&ts);
1717: TSSetDM(ts,pack);
1718: TSSetProblemType(ts,TS_NONLINEAR);
1719: TSMonitorSet(ts,THITSMonitor,thi,NULL);
1720: TSSetType(ts,TSTHETA);
1721: TSSetIFunction(ts,PETSC_NULL,THIFunction,thi);
1722: TSSetIJacobian(ts,B,B,THIJacobian,thi);
1723: TSSetDuration(ts,100,10.0);
1724: TSSetSolution(ts,X);
1725: TSSetInitialTimeStep(ts,0.,1e-3);
1726: TSSetFromOptions(ts);
1728: TSSolve(ts,X,&ftime);
1729: TSGetTimeStepNumber(ts,&steps);
1730: PetscPrintf(PETSC_COMM_WORLD,"Steps %D final time %G\n",steps,ftime);
1732: if (0) {THISolveStatistics(thi,dmmg,0,"Full");}
1734: {
1735: PetscBool flg;
1736: char filename[PETSC_MAX_PATH_LEN] = "";
1737: PetscOptionsGetString(PETSC_NULL,"-o",filename,sizeof(filename),&flg);
1738: if (flg) {
1739: THIDAVecView_VTK_XML(thi,pack,X,filename,NULL);
1740: }
1741: }
1743: DMMGDestroy(dmmg);
1744: DMDestroy(&pack);
1745: TSDestroy(&ts);
1746: THIDestroy(&thi);
1747: PetscFinalize();
1748: return 0;
1749: }