Actual source code: ex48.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: There are two compile-time options:
46: NO_SSE2:
47: If the host supports SSE2, we use integration code that has been vectorized with SSE2
48: intrinsics, unless this macro is defined. The intrinsics speed up integration by about
49: 30% on my architecture (P8700, gcc-4.5 snapshot).
51: COMPUTE_LOWER_TRIANGULAR:
52: The element matrices we assemble are lower-triangular so it is not necessary to compute
53: all entries explicitly. If this macro is defined, the lower-triangular entries are
54: computed explicitly.
56: */
58: #include <petscdmmg.h>
59: #include <ctype.h> /* toupper() */
60: #include <private/daimpl.h> /* There is not yet a public interface to manipulate dm->ops */
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
69: #if defined __SSE2__
70: # include <emmintrin.h>
71: #endif
73: /* The SSE2 kernels are only for PetscScalar=double on architectures that support it */
74: #define USE_SSE2_KERNELS (!defined NO_SSE2 \
75: && !defined PETSC_USE_COMPLEX \
76: && !defined PETSC_USE_REAL_SINGLE \
77: && defined __SSE2__)
79: static PetscClassId THI_CLASSID;
81: typedef enum {QUAD_GAUSS,QUAD_LOBATTO} QuadratureType;
82: static const char *QuadratureTypes[] = {"gauss","lobatto","QuadratureType","QUAD_",0};
83: static const PetscReal HexQWeights[8] = {1,1,1,1,1,1,1,1};
84: static const PetscReal HexQNodes[] = {-0.57735026918962573, 0.57735026918962573};
85: #define G 0.57735026918962573
86: #define H (0.5*(1.+G))
87: #define L (0.5*(1.-G))
88: #define M (-0.5)
89: #define P (0.5)
90: /* Special quadrature: Lobatto in horizontal, Gauss in vertical */
91: static const PetscReal HexQInterp_Lobatto[8][8] = {{H,0,0,0,L,0,0,0},
92: {0,H,0,0,0,L,0,0},
93: {0,0,H,0,0,0,L,0},
94: {0,0,0,H,0,0,0,L},
95: {L,0,0,0,H,0,0,0},
96: {0,L,0,0,0,H,0,0},
97: {0,0,L,0,0,0,H,0},
98: {0,0,0,L,0,0,0,H}};
99: static const PetscReal HexQDeriv_Lobatto[8][8][3] = {
100: {{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} },
101: {{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} },
102: {{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} },
103: {{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}},
104: {{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} },
105: {{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} },
106: {{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} },
107: {{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}}};
108: /* Stanndard Gauss */
109: 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},
110: {L*H*H,H*H*H,H*L*H,L*L*H, L*H*L,H*H*L,H*L*L,L*L*L},
111: {L*L*H,H*L*H,H*H*H,L*H*H, L*L*L,H*L*L,H*H*L,L*H*L},
112: {H*L*H,L*L*H,L*H*H,H*H*H, H*L*L,L*L*L,L*H*L,H*H*L},
113: {H*H*L,L*H*L,L*L*L,H*L*L, H*H*H,L*H*H,L*L*H,H*L*H},
114: {L*H*L,H*H*L,H*L*L,L*L*L, L*H*H,H*H*H,H*L*H,L*L*H},
115: {L*L*L,H*L*L,H*H*L,L*H*L, L*L*H,H*L*H,H*H*H,L*H*H},
116: {H*L*L,L*L*L,L*H*L,H*H*L, H*L*H,L*L*H,L*H*H,H*H*H}};
117: static const PetscReal HexQDeriv_Gauss[8][8][3] = {
118: {{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}},
119: {{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}},
120: {{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}},
121: {{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}},
122: {{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}},
123: {{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}},
124: {{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}},
125: {{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}}};
126: static const PetscReal (*HexQInterp)[8],(*HexQDeriv)[8][3];
127: /* Standard 2x2 Gauss quadrature for the bottom layer. */
128: static const PetscReal QuadQInterp[4][4] = {{H*H,L*H,L*L,H*L},
129: {L*H,H*H,H*L,L*L},
130: {L*L,H*L,H*H,L*H},
131: {H*L,L*L,L*H,H*H}};
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: #undef G
138: #undef H
139: #undef L
140: #undef M
141: #undef P
143: #define HexExtract(x,i,j,k,n) do { \
144: (n)[0] = (x)[i][j][k]; \
145: (n)[1] = (x)[i+1][j][k]; \
146: (n)[2] = (x)[i+1][j+1][k]; \
147: (n)[3] = (x)[i][j+1][k]; \
148: (n)[4] = (x)[i][j][k+1]; \
149: (n)[5] = (x)[i+1][j][k+1]; \
150: (n)[6] = (x)[i+1][j+1][k+1]; \
151: (n)[7] = (x)[i][j+1][k+1]; \
152: } while (0)
154: #define HexExtractRef(x,i,j,k,n) do { \
155: (n)[0] = &(x)[i][j][k]; \
156: (n)[1] = &(x)[i+1][j][k]; \
157: (n)[2] = &(x)[i+1][j+1][k]; \
158: (n)[3] = &(x)[i][j+1][k]; \
159: (n)[4] = &(x)[i][j][k+1]; \
160: (n)[5] = &(x)[i+1][j][k+1]; \
161: (n)[6] = &(x)[i+1][j+1][k+1]; \
162: (n)[7] = &(x)[i][j+1][k+1]; \
163: } while (0)
165: #define QuadExtract(x,i,j,n) do { \
166: (n)[0] = (x)[i][j]; \
167: (n)[1] = (x)[i+1][j]; \
168: (n)[2] = (x)[i+1][j+1]; \
169: (n)[3] = (x)[i][j+1]; \
170: } while (0)
172: static PetscScalar Sqr(PetscScalar a) {return a*a;}
174: static void HexGrad(const PetscReal dphi[][3],const PetscReal zn[],PetscReal dz[])
175: {
176: PetscInt i;
177: dz[0] = dz[1] = dz[2] = 0;
178: for (i=0; i<8; i++) {
179: dz[0] += dphi[i][0] * zn[i];
180: dz[1] += dphi[i][1] * zn[i];
181: dz[2] += dphi[i][2] * zn[i];
182: }
183: }
185: static void HexComputeGeometry(PetscInt q,PetscReal hx,PetscReal hy,const PetscReal dz[restrict],PetscReal phi[restrict],PetscReal dphi[restrict][3],PetscReal *restrict jw)
186: {
187: const PetscReal
188: jac[3][3] = {{hx/2,0,0}, {0,hy/2,0}, {dz[0],dz[1],dz[2]}}
189: ,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]}}
190: ,jdet = jac[0][0]*jac[1][1]*jac[2][2];
191: PetscInt i;
193: for (i=0; i<8; i++) {
194: const PetscReal *dphir = HexQDeriv[q][i];
195: phi[i] = HexQInterp[q][i];
196: dphi[i][0] = dphir[0]*ijac[0][0] + dphir[1]*ijac[1][0] + dphir[2]*ijac[2][0];
197: dphi[i][1] = dphir[0]*ijac[0][1] + dphir[1]*ijac[1][1] + dphir[2]*ijac[2][1];
198: dphi[i][2] = dphir[0]*ijac[0][2] + dphir[1]*ijac[1][2] + dphir[2]*ijac[2][2];
199: }
200: *jw = 1.0 * jdet;
201: }
203: typedef struct _p_THI *THI;
204: typedef struct _n_Units *Units;
206: typedef struct {
207: PetscScalar u,v;
208: } Node;
210: typedef struct {
211: PetscScalar b; /* bed */
212: PetscScalar h; /* thickness */
213: PetscScalar beta2; /* friction */
214: } PrmNode;
216: typedef struct {
217: PetscReal min,max,cmin,cmax;
218: } PRange;
220: typedef enum {THIASSEMBLY_TRIDIAGONAL,THIASSEMBLY_FULL} THIAssemblyMode;
222: struct _p_THI {
223: PETSCHEADER(int);
224: void (*initialize)(THI,PetscReal x,PetscReal y,PrmNode *p);
225: PetscInt nlevels;
226: PetscInt zlevels;
227: PetscReal Lx,Ly,Lz; /* Model domain */
228: PetscReal alpha; /* Bed angle */
229: Units units;
230: PetscReal dirichlet_scale;
231: PetscReal ssa_friction_scale;
232: PRange eta;
233: PRange beta2;
234: struct {
235: PetscReal Bd2,eps,exponent;
236: } viscosity;
237: struct {
238: PetscReal irefgam,eps2,exponent,refvel,epsvel;
239: } friction;
240: PetscReal rhog;
241: PetscBool no_slip;
242: PetscBool tridiagonal;
243: PetscBool coarse2d;
244: PetscBool verbose;
245: MatType mattype;
246: };
248: struct _n_Units {
249: /* fundamental */
250: PetscReal meter;
251: PetscReal kilogram;
252: PetscReal second;
253: /* derived */
254: PetscReal Pascal;
255: PetscReal year;
256: };
258: static void PrmHexGetZ(const PrmNode pn[],PetscInt k,PetscInt zm,PetscReal zn[])
259: {
260: const PetscScalar zm1 = zm-1,
261: znl[8] = {pn[0].b + pn[0].h*(PetscScalar)k/zm1,
262: pn[1].b + pn[1].h*(PetscScalar)k/zm1,
263: pn[2].b + pn[2].h*(PetscScalar)k/zm1,
264: pn[3].b + pn[3].h*(PetscScalar)k/zm1,
265: pn[0].b + pn[0].h*(PetscScalar)(k+1)/zm1,
266: pn[1].b + pn[1].h*(PetscScalar)(k+1)/zm1,
267: pn[2].b + pn[2].h*(PetscScalar)(k+1)/zm1,
268: pn[3].b + pn[3].h*(PetscScalar)(k+1)/zm1};
269: PetscInt i;
270: for (i=0; i<8; i++) zn[i] = PetscRealPart(znl[i]);
271: }
273: /* Tests A and C are from the ISMIP-HOM paper (Pattyn et al. 2008) */
274: static void THIInitialize_HOM_A(THI thi,PetscReal x,PetscReal y,PrmNode *p)
275: {
276: Units units = thi->units;
277: PetscReal s = -x*sin(thi->alpha);
278: p->b = s - 1000*units->meter + 500*units->meter * sin(x*2*PETSC_PI/thi->Lx) * sin(y*2*PETSC_PI/thi->Ly);
279: p->h = s - p->b;
280: p->beta2 = 1e30;
281: }
283: static void THIInitialize_HOM_C(THI thi,PetscReal x,PetscReal y,PrmNode *p)
284: {
285: Units units = thi->units;
286: PetscReal s = -x*sin(thi->alpha);
287: p->b = s - 1000*units->meter;
288: p->h = s - p->b;
289: /* tau_b = beta2 v is a stress (Pa) */
290: p->beta2 = 1000 * (1 + sin(x*2*PETSC_PI/thi->Lx)*sin(y*2*PETSC_PI/thi->Ly)) * units->Pascal * units->year / units->meter;
291: }
293: /* These are just toys */
295: /* Same bed as test A, free slip everywhere except for a discontinuous jump to a circular sticky region in the middle. */
296: static void THIInitialize_HOM_X(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
297: {
298: Units units = thi->units;
299: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
300: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
301: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
302: p->h = s - p->b;
303: p->beta2 = 1000 * (r < 1 ? 2 : 0) * units->Pascal * units->year / units->meter;
304: }
306: /* Like Z, but with 200 meter cliffs */
307: static void THIInitialize_HOM_Y(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
308: {
309: Units units = thi->units;
310: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
311: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
312: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
313: if (PetscRealPart(p->b) > -700*units->meter) p->b += 200*units->meter;
314: p->h = s - p->b;
315: 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;
316: }
318: /* Same bed as A, smoothly varying slipperiness, similar to MATLAB's "sombrero" (uncorrelated with bathymetry) */
319: static void THIInitialize_HOM_Z(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
320: {
321: Units units = thi->units;
322: PetscReal x = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
323: PetscReal r = sqrt(x*x + y*y),s = -x*sin(thi->alpha);
324: p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
325: p->h = s - p->b;
326: 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;
327: }
329: static void THIFriction(THI thi,PetscReal rbeta2,PetscReal gam,PetscReal *beta2,PetscReal *dbeta2)
330: {
331: if (thi->friction.irefgam == 0) {
332: Units units = thi->units;
333: thi->friction.irefgam = 1./(0.5*PetscSqr(thi->friction.refvel * units->meter / units->year));
334: thi->friction.eps2 = 0.5*PetscSqr(thi->friction.epsvel * units->meter / units->year) * thi->friction.irefgam;
335: }
336: if (thi->friction.exponent == 0) {
337: *beta2 = rbeta2;
338: *dbeta2 = 0;
339: } else {
340: *beta2 = rbeta2 * pow(thi->friction.eps2 + gam*thi->friction.irefgam,thi->friction.exponent);
341: *dbeta2 = thi->friction.exponent * *beta2 / (thi->friction.eps2 + gam*thi->friction.irefgam) * thi->friction.irefgam;
342: }
343: }
345: static void THIViscosity(THI thi,PetscReal gam,PetscReal *eta,PetscReal *deta)
346: {
347: PetscReal Bd2,eps,exponent;
348: if (thi->viscosity.Bd2 == 0) {
349: Units units = thi->units;
350: const PetscReal
351: n = 3., /* Glen exponent */
352: p = 1. + 1./n, /* for Stokes */
353: A = 1.e-16 * pow(units->Pascal,-n) / units->year, /* softness parameter (Pa^{-n}/s) */
354: B = pow(A,-1./n); /* hardness parameter */
355: thi->viscosity.Bd2 = B/2;
356: thi->viscosity.exponent = (p-2)/2;
357: thi->viscosity.eps = 0.5*PetscSqr(1e-5 / units->year);
358: }
359: Bd2 = thi->viscosity.Bd2;
360: exponent = thi->viscosity.exponent;
361: eps = thi->viscosity.eps;
362: *eta = Bd2 * pow(eps + gam,exponent);
363: *deta = exponent * (*eta) / (eps + gam);
364: }
366: static void RangeUpdate(PetscReal *min,PetscReal *max,PetscReal x)
367: {
368: if (x < *min) *min = x;
369: if (x > *max) *max = x;
370: }
372: static void PRangeClear(PRange *p)
373: {
374: p->cmin = p->min = 1e100;
375: p->cmax = p->max = -1e100;
376: }
380: static PetscErrorCode PRangeMinMax(PRange *p,PetscReal min,PetscReal max)
381: {
384: p->cmin = min;
385: p->cmax = max;
386: if (min < p->min) p->min = min;
387: if (max > p->max) p->max = max;
388: return(0);
389: }
393: static PetscErrorCode THIDestroy(THI *thi)
394: {
398: if (!*thi) return(0);
399: if (--((PetscObject)(*thi))->refct > 0) {*thi = 0; return(0);}
400: PetscFree((*thi)->units);
401: PetscFree((*thi)->mattype);
402: PetscHeaderDestroy(thi);
403: return(0);
404: }
408: static PetscErrorCode THICreate(MPI_Comm comm,THI *inthi)
409: {
410: static PetscBool registered = PETSC_FALSE;
411: THI thi;
412: Units units;
416: *inthi = 0;
417: if (!registered) {
418: PetscClassIdRegister("Toy Hydrostatic Ice",&THI_CLASSID);
419: registered = PETSC_TRUE;
420: }
421: PetscHeaderCreate(thi,_p_THI,0,THI_CLASSID,-1,"THI","Toy Hydrostatic Ice","",comm,THIDestroy,0);
423: PetscNew(struct _n_Units,&thi->units);
424: units = thi->units;
425: units->meter = 1e-2;
426: units->second = 1e-7;
427: units->kilogram = 1e-12;
428: PetscOptionsBegin(comm,NULL,"Scaled units options","");
429: {
430: PetscOptionsReal("-units_meter","1 meter in scaled length units","",units->meter,&units->meter,NULL);
431: PetscOptionsReal("-units_second","1 second in scaled time units","",units->second,&units->second,NULL);
432: PetscOptionsReal("-units_kilogram","1 kilogram in scaled mass units","",units->kilogram,&units->kilogram,NULL);
433: }
434: PetscOptionsEnd();
435: units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second));
436: units->year = 31556926. * units->second, /* seconds per year */
438: thi->Lx = 10.e3;
439: thi->Ly = 10.e3;
440: thi->Lz = 1000;
441: thi->nlevels = 1;
442: thi->dirichlet_scale = 1;
443: thi->verbose = PETSC_FALSE;
445: PetscOptionsBegin(comm,NULL,"Toy Hydrostatic Ice options","");
446: {
447: QuadratureType quad = QUAD_GAUSS;
448: char homexp[] = "A";
449: char mtype[256] = MATSBAIJ;
450: PetscReal L,m = 1.0;
451: PetscBool flg;
452: L = thi->Lx;
453: PetscOptionsReal("-thi_L","Domain size (m)","",L,&L,&flg);
454: if (flg) thi->Lx = thi->Ly = L;
455: PetscOptionsReal("-thi_Lx","X Domain size (m)","",thi->Lx,&thi->Lx,NULL);
456: PetscOptionsReal("-thi_Ly","Y Domain size (m)","",thi->Ly,&thi->Ly,NULL);
457: PetscOptionsReal("-thi_Lz","Z Domain size (m)","",thi->Lz,&thi->Lz,NULL);
458: PetscOptionsString("-thi_hom","ISMIP-HOM experiment (A or C)","",homexp,homexp,sizeof(homexp),NULL);
459: switch (homexp[0] = toupper(homexp[0])) {
460: case 'A':
461: thi->initialize = THIInitialize_HOM_A;
462: thi->no_slip = PETSC_TRUE;
463: thi->alpha = 0.5;
464: break;
465: case 'C':
466: thi->initialize = THIInitialize_HOM_C;
467: thi->no_slip = PETSC_FALSE;
468: thi->alpha = 0.1;
469: break;
470: case 'X':
471: thi->initialize = THIInitialize_HOM_X;
472: thi->no_slip = PETSC_FALSE;
473: thi->alpha = 0.3;
474: break;
475: case 'Y':
476: thi->initialize = THIInitialize_HOM_Y;
477: thi->no_slip = PETSC_FALSE;
478: thi->alpha = 0.5;
479: break;
480: case 'Z':
481: thi->initialize = THIInitialize_HOM_Z;
482: thi->no_slip = PETSC_FALSE;
483: thi->alpha = 0.5;
484: break;
485: default:
486: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"HOM experiment '%c' not implemented",homexp[0]);
487: }
488: PetscOptionsEnum("-thi_quadrature","Quadrature to use for 3D elements","",QuadratureTypes,(PetscEnum)quad,(PetscEnum*)&quad,NULL);
489: switch (quad) {
490: case QUAD_GAUSS:
491: HexQInterp = HexQInterp_Gauss;
492: HexQDeriv = HexQDeriv_Gauss;
493: break;
494: case QUAD_LOBATTO:
495: HexQInterp = HexQInterp_Lobatto;
496: HexQDeriv = HexQDeriv_Lobatto;
497: break;
498: }
499: PetscOptionsReal("-thi_alpha","Bed angle (degrees)","",thi->alpha,&thi->alpha,NULL);
500: thi->friction.refvel = 100.;
501: thi->friction.epsvel = 1.;
502: PetscOptionsReal("-thi_friction_refvel","Reference velocity for sliding","",thi->friction.refvel,&thi->friction.refvel,NULL);
503: PetscOptionsReal("-thi_friction_epsvel","Regularization velocity for sliding","",thi->friction.epsvel,&thi->friction.epsvel,NULL);
504: PetscOptionsReal("-thi_friction_m","Friction exponent, 0=Coulomb, 1=Navier","",m,&m,NULL);
505: thi->friction.exponent = (m-1)/2;
506: PetscOptionsReal("-thi_dirichlet_scale","Scale Dirichlet boundary conditions by this factor","",thi->dirichlet_scale,&thi->dirichlet_scale,NULL);
507: 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);
508: PetscOptionsInt("-thi_nlevels","Number of levels of refinement","",thi->nlevels,&thi->nlevels,NULL);
509: PetscOptionsBool("-thi_coarse2d","Use a 2D coarse space corresponding to SSA","",thi->coarse2d,&thi->coarse2d,NULL);
510: PetscOptionsBool("-thi_tridiagonal","Assemble a tridiagonal system (column coupling only) on the finest level","",thi->tridiagonal,&thi->tridiagonal,NULL);
511: PetscOptionsList("-thi_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)mtype,sizeof(mtype),NULL);
512: PetscStrallocpy(mtype,&thi->mattype);
513: PetscOptionsBool("-thi_verbose","Enable verbose output (like matrix sizes and statistics)","",thi->verbose,&thi->verbose,NULL);
514: }
515: PetscOptionsEnd();
517: /* dimensionalize */
518: thi->Lx *= units->meter;
519: thi->Ly *= units->meter;
520: thi->Lz *= units->meter;
521: thi->alpha *= PETSC_PI / 180;
523: PRangeClear(&thi->eta);
524: PRangeClear(&thi->beta2);
526: {
527: PetscReal u = 1000*units->meter/(3e7*units->second),
528: gradu = u / (100*units->meter),eta,deta,
529: rho = 910 * units->kilogram/pow(units->meter,3),
530: grav = 9.81 * units->meter/PetscSqr(units->second),
531: driving = rho * grav * sin(thi->alpha) * 1000*units->meter;
532: THIViscosity(thi,0.5*gradu*gradu,&eta,&deta);
533: thi->rhog = rho * grav;
534: if (thi->verbose) {
535: 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);
536: 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);
537: 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);
538: THIViscosity(thi,0.5*PetscSqr(1e-3*gradu),&eta,&deta);
539: 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);
540: }
541: }
543: *inthi = thi;
544: return(0);
545: }
549: static PetscErrorCode THIInitializePrm(THI thi,DM da2prm,Vec prm)
550: {
551: PrmNode **p;
552: PetscInt i,j,xs,xm,ys,ym,mx,my;
556: DMDAGetGhostCorners(da2prm,&ys,&xs,0,&ym,&xm,0);
557: DMDAGetInfo(da2prm,0, &my,&mx,0, 0,0,0, 0,0,0,0,0,0);
558: DMDAVecGetArray(da2prm,prm,&p);
559: for (i=xs; i<xs+xm; i++) {
560: for (j=ys; j<ys+ym; j++) {
561: PetscReal xx = thi->Lx*i/mx,yy = thi->Ly*j/my;
562: thi->initialize(thi,xx,yy,&p[i][j]);
563: }
564: }
565: DMDAVecRestoreArray(da2prm,prm,&p);
566: return(0);
567: }
571: static PetscErrorCode THISetDMMG(THI thi,DMMG *dmmg)
572: {
574: PetscInt i;
577: if (DMMGGetLevels(dmmg) != thi->nlevels) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMMG nlevels does not agree with THI");
578: for (i=0; i<thi->nlevels; i++) {
579: PetscInt Mx,My,Mz,mx,my,s,dim;
580: DMDAStencilType st;
581: DM da = dmmg[i]->dm,da2prm;
582: Vec X;
583: DMDAGetInfo(da,&dim, &Mz,&My,&Mx, 0,&my,&mx, 0,&s,0,0,0,&st);
584: if (dim == 2) {
585: DMDAGetInfo(da,&dim, &My,&Mx,0, &my,&mx,0, 0,&s,0,0,0,&st);
586: }
587: DMDACreate2d(((PetscObject)thi)->comm,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,st,My,Mx,my,mx,sizeof(PrmNode)/sizeof(PetscScalar),s,0,0,&da2prm);
588: DMCreateLocalVector(da2prm,&X);
589: {
590: PetscReal Lx = thi->Lx / thi->units->meter,Ly = thi->Ly / thi->units->meter,Lz = thi->Lz / thi->units->meter;
591: if (dim == 2) {
592: PetscPrintf(((PetscObject)thi)->comm,"Level %d domain size (m) %8.2g x %8.2g, num elements %3d x %3d (%8d), size (m) %g x %g\n",i,Lx,Ly,Mx,My,Mx*My,Lx/Mx,Ly/My);
593: } else {
594: 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));
595: }
596: }
597: THIInitializePrm(thi,da2prm,X);
598: PetscObjectCompose((PetscObject)da,"DMDA2Prm",(PetscObject)da2prm);
599: PetscObjectCompose((PetscObject)da,"DMDA2Prm_Vec",(PetscObject)X);
600: DMDestroy(&da2prm);
601: VecDestroy(&X);
602: }
603: return(0);
604: }
608: static PetscErrorCode THIDAGetPrm(DM da,PrmNode ***prm)
609: {
611: DM da2prm;
612: Vec X;
615: PetscObjectQuery((PetscObject)da,"DMDA2Prm",(PetscObject*)&da2prm);
616: if (!da2prm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm composed with given DMDA");
617: PetscObjectQuery((PetscObject)da,"DMDA2Prm_Vec",(PetscObject*)&X);
618: if (!X) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm_Vec composed with given DMDA");
619: DMDAVecGetArray(da2prm,X,prm);
620: return(0);
621: }
625: static PetscErrorCode THIDARestorePrm(DM da,PrmNode ***prm)
626: {
628: DM da2prm;
629: Vec X;
632: PetscObjectQuery((PetscObject)da,"DMDA2Prm",(PetscObject*)&da2prm);
633: if (!da2prm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm composed with given DMDA");
634: PetscObjectQuery((PetscObject)da,"DMDA2Prm_Vec",(PetscObject*)&X);
635: if (!X) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm_Vec composed with given DMDA");
636: DMDAVecRestoreArray(da2prm,X,prm);
637: return(0);
638: }
642: static PetscErrorCode THIInitial(DMMG dmmg,Vec X)
643: {
644: THI thi = (THI)dmmg->user;
645: DM da = dmmg->dm;
646: PetscInt i,j,k,xs,xm,ys,ym,zs,zm,mx,my;
647: PetscReal hx,hy;
648: PrmNode **prm;
649: Node ***x;
653: DMDAGetInfo(da,0, 0,&my,&mx, 0,0,0, 0,0,0,0,0,0);
654: DMDAGetCorners(da,&zs,&ys,&xs,&zm,&ym,&xm);
655: DMDAVecGetArray(da,X,&x);
656: THIDAGetPrm(da,&prm);
657: hx = thi->Lx / mx;
658: hy = thi->Ly / my;
659: for (i=xs; i<xs+xm; i++) {
660: for (j=ys; j<ys+ym; j++) {
661: for (k=zs; k<zs+zm; k++) {
662: const PetscScalar zm1 = zm-1,
663: 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),
664: 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);
665: x[i][j][k].u = 0. * drivingx * prm[i][j].h*(PetscScalar)k/zm1;
666: x[i][j][k].v = 0. * drivingy * prm[i][j].h*(PetscScalar)k/zm1;
667: }
668: }
669: }
670: DMDAVecRestoreArray(da,X,&x);
671: THIDARestorePrm(da,&prm);
672: return(0);
673: }
675: static void PointwiseNonlinearity(THI thi,const Node n[restrict],const PetscReal phi[restrict],PetscReal dphi[restrict][3],PetscScalar *restrict u,PetscScalar *restrict v,PetscScalar du[restrict],PetscScalar dv[restrict],PetscReal *eta,PetscReal *deta)
676: {
677: PetscInt l,ll;
678: PetscScalar gam;
680: du[0] = du[1] = du[2] = 0;
681: dv[0] = dv[1] = dv[2] = 0;
682: *u = 0;
683: *v = 0;
684: for (l=0; l<8; l++) {
685: *u += phi[l] * n[l].u;
686: *v += phi[l] * n[l].v;
687: for (ll=0; ll<3; ll++) {
688: du[ll] += dphi[l][ll] * n[l].u;
689: dv[ll] += dphi[l][ll] * n[l].v;
690: }
691: }
692: 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]);
693: THIViscosity(thi,PetscRealPart(gam),eta,deta);
694: }
696: static void PointwiseNonlinearity2D(THI thi,Node n[],PetscReal phi[],PetscReal dphi[4][2],PetscScalar *u,PetscScalar *v,PetscScalar du[],PetscScalar dv[],PetscReal *eta,PetscReal *deta)
697: {
698: PetscInt l,ll;
699: PetscScalar gam;
701: du[0] = du[1] = 0;
702: dv[0] = dv[1] = 0;
703: *u = 0;
704: *v = 0;
705: for (l=0; l<4; l++) {
706: *u += phi[l] * n[l].u;
707: *v += phi[l] * n[l].v;
708: for (ll=0; ll<2; ll++) {
709: du[ll] += dphi[l][ll] * n[l].u;
710: dv[ll] += dphi[l][ll] * n[l].v;
711: }
712: }
713: gam = Sqr(du[0]) + Sqr(dv[1]) + du[0]*dv[1] + 0.25*Sqr(du[1]+dv[0]);
714: THIViscosity(thi,PetscRealPart(gam),eta,deta);
715: }
719: static PetscErrorCode THIFunctionLocal(DMDALocalInfo *info,Node ***x,Node ***f,THI thi)
720: {
721: PetscInt xs,ys,xm,ym,zm,i,j,k,q,l;
722: PetscReal hx,hy,etamin,etamax,beta2min,beta2max;
723: PrmNode **prm;
727: xs = info->zs;
728: ys = info->ys;
729: xm = info->zm;
730: ym = info->ym;
731: zm = info->xm;
732: hx = thi->Lx / info->mz;
733: hy = thi->Ly / info->my;
735: etamin = 1e100;
736: etamax = 0;
737: beta2min = 1e100;
738: beta2max = 0;
740: THIDAGetPrm(info->da,&prm);
742: for (i=xs; i<xs+xm; i++) {
743: for (j=ys; j<ys+ym; j++) {
744: PrmNode pn[4];
745: QuadExtract(prm,i,j,pn);
746: for (k=0; k<zm-1; k++) {
747: PetscInt ls = 0;
748: Node n[8],*fn[8];
749: PetscReal zn[8],etabase = 0;
750: PrmHexGetZ(pn,k,zm,zn);
751: HexExtract(x,i,j,k,n);
752: HexExtractRef(f,i,j,k,fn);
753: if (thi->no_slip && k == 0) {
754: for (l=0; l<4; l++) n[l].u = n[l].v = 0;
755: /* The first 4 basis functions lie on the bottom layer, so their contribution is exactly 0, hence we can skip them */
756: ls = 4;
757: }
758: for (q=0; q<8; q++) {
759: PetscReal dz[3],phi[8],dphi[8][3],jw,eta,deta;
760: PetscScalar du[3],dv[3],u,v;
761: HexGrad(HexQDeriv[q],zn,dz);
762: HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
763: PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
764: jw /= thi->rhog; /* scales residuals to be O(1) */
765: if (q == 0) etabase = eta;
766: RangeUpdate(&etamin,&etamax,eta);
767: for (l=ls; l<8; l++) { /* test functions */
768: const PetscReal ds[2] = {-sin(thi->alpha),0};
769: const PetscReal pp=phi[l],*dp = dphi[l];
770: 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];
771: 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];
772: }
773: }
774: if (k == 0) { /* we are on a bottom face */
775: if (thi->no_slip) {
776: /* Note: Non-Galerkin coarse grid operators are very sensitive to the scaling of Dirichlet boundary
777: * conditions. After shenanigans above, etabase contains the effective viscosity at the closest quadrature
778: * point to the bed. We want the diagonal entry in the Dirichlet condition to have similar magnitude to the
779: * diagonal entry corresponding to the adjacent node. The fundamental scaling of the viscous part is in
780: * diagu, diagv below. This scaling is easy to recognize by considering the finite difference operator after
781: * scaling by element size. The no-slip Dirichlet condition is scaled by this factor, and also in the
782: * assembled matrix (see the similar block in THIJacobianLocal).
783: *
784: * Note that the residual at this Dirichlet node is linear in the state at this node, but also depends
785: * (nonlinearly in general) on the neighboring interior nodes through the local viscosity. This will make
786: * a matrix-free Jacobian have extra entries in the corresponding row. We assemble only the diagonal part,
787: * so the solution will exactly satisfy the boundary condition after the first linear iteration.
788: */
789: const PetscReal hz = PetscRealPart(pn[0].h)/(zm-1.);
790: 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);
791: fn[0]->u = thi->dirichlet_scale*diagu*x[i][j][k].u;
792: fn[0]->v = thi->dirichlet_scale*diagv*x[i][j][k].v;
793: } else { /* Integrate over bottom face to apply boundary condition */
794: for (q=0; q<4; q++) {
795: const PetscReal jw = 0.25*hx*hy/thi->rhog,*phi = QuadQInterp[q];
796: PetscScalar u=0,v=0,rbeta2=0;
797: PetscReal beta2,dbeta2;
798: for (l=0; l<4; l++) {
799: u += phi[l]*n[l].u;
800: v += phi[l]*n[l].v;
801: rbeta2 += phi[l]*pn[l].beta2;
802: }
803: THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
804: RangeUpdate(&beta2min,&beta2max,beta2);
805: for (l=0; l<4; l++) {
806: const PetscReal pp = phi[l];
807: fn[ls+l]->u += pp*jw*beta2*u;
808: fn[ls+l]->v += pp*jw*beta2*v;
809: }
810: }
811: }
812: }
813: }
814: }
815: }
817: THIDARestorePrm(info->da,&prm);
819: PRangeMinMax(&thi->eta,etamin,etamax);
820: PRangeMinMax(&thi->beta2,beta2min,beta2max);
821: return(0);
822: }
826: static PetscErrorCode THIMatrixStatistics(THI thi,Mat B,PetscViewer viewer)
827: {
829: PetscReal nrm;
830: PetscInt m;
831: PetscMPIInt rank;
834: MatNorm(B,NORM_FROBENIUS,&nrm);
835: MatGetSize(B,&m,0);
836: MPI_Comm_rank(((PetscObject)B)->comm,&rank);
837: if (!rank) {
838: PetscScalar val0,val2;
839: MatGetValue(B,0,0,&val0);
840: MatGetValue(B,2,2,&val2);
841: PetscViewerASCIIPrintf(viewer,"Matrix dim %8d norm %8.2e (0,0) %8.2e (2,2) %8.2e %8.2e <= eta <= %8.2e %8.2e <= beta2 <= %8.2e\n",m,nrm,PetscRealPart(val0),PetscRealPart(val2),thi->eta.cmin,thi->eta.cmax,thi->beta2.cmin,thi->beta2.cmax);
842: }
843: return(0);
844: }
848: static PetscErrorCode THISurfaceStatistics(DM da,Vec X,PetscReal *min,PetscReal *max,PetscReal *mean)
849: {
851: Node ***x;
852: PetscInt i,j,xs,ys,zs,xm,ym,zm,mx,my,mz;
853: PetscReal umin = 1e100,umax=-1e100;
854: PetscScalar usum=0.0,gusum;
857: *min = *max = *mean = 0;
858: DMDAGetInfo(da,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
859: DMDAGetCorners(da,&zs,&ys,&xs,&zm,&ym,&xm);
860: if (zs != 0 || zm != mz) SETERRQ(PETSC_COMM_SELF,1,"Unexpected decomposition");
861: DMDAVecGetArray(da,X,&x);
862: for (i=xs; i<xs+xm; i++) {
863: for (j=ys; j<ys+ym; j++) {
864: PetscReal u = PetscRealPart(x[i][j][zm-1].u);
865: RangeUpdate(&umin,&umax,u);
866: usum += u;
867: }
868: }
869: DMDAVecRestoreArray(da,X,&x);
870: MPI_Allreduce(&umin,min,1,MPIU_REAL,MPIU_MIN,((PetscObject)da)->comm);
871: MPI_Allreduce(&umax,max,1,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);
872: MPI_Allreduce(&usum,&gusum,1,MPIU_SCALAR,MPIU_SUM,((PetscObject)da)->comm);
873: *mean = PetscRealPart(gusum) / (mx*my);
874: return(0);
875: }
879: static PetscErrorCode THISolveStatistics(THI thi,DMMG *dmmg,PetscInt coarsened,const char name[])
880: {
881: MPI_Comm comm = ((PetscObject)thi)->comm;
882: PetscInt nlevels = DMMGGetLevels(dmmg),level = nlevels-1-coarsened;
883: SNES snes = dmmg[level]->snes;
884: Vec X = dmmg[level]->x;
888: PetscPrintf(comm,"Solution statistics after solve: %s\n",name);
889: {
890: PetscInt its,lits;
891: SNESConvergedReason reason;
892: SNESGetIterationNumber(snes,&its);
893: SNESGetConvergedReason(snes,&reason);
894: SNESGetLinearSolveIterations(snes,&lits);
895: PetscPrintf(comm,"%s: Number of SNES iterations = %d, total linear iterations = %d\n",SNESConvergedReasons[reason],its,lits);
896: }
897: {
898: PetscReal nrm2,tmin[3]={1e100,1e100,1e100},tmax[3]={-1e100,-1e100,-1e100},min[3],max[3];
899: PetscInt i,j,m;
900: PetscScalar *x;
901: VecNorm(X,NORM_2,&nrm2);
902: VecGetLocalSize(X,&m);
903: VecGetArray(X,&x);
904: for (i=0; i<m; i+=2) {
905: PetscReal u = PetscRealPart(x[i]),v = PetscRealPart(x[i+1]),c = sqrt(u*u+v*v);
906: tmin[0] = PetscMin(u,tmin[0]);
907: tmin[1] = PetscMin(v,tmin[1]);
908: tmin[2] = PetscMin(c,tmin[2]);
909: tmax[0] = PetscMax(u,tmax[0]);
910: tmax[1] = PetscMax(v,tmax[1]);
911: tmax[2] = PetscMax(c,tmax[2]);
912: }
913: VecRestoreArray(X,&x);
914: MPI_Allreduce(tmin,min,3,MPIU_REAL,MPIU_MIN,((PetscObject)thi)->comm);
915: MPI_Allreduce(tmax,max,3,MPIU_REAL,MPIU_MAX,((PetscObject)thi)->comm);
916: /* Dimensionalize to meters/year */
917: nrm2 *= thi->units->year / thi->units->meter;
918: for (j=0; j<3; j++) {
919: min[j] *= thi->units->year / thi->units->meter;
920: max[j] *= thi->units->year / thi->units->meter;
921: }
922: PetscPrintf(comm,"|X|_2 %g %g <= u <= %g %g <= v <= %g %g <= c <= %g \n",nrm2,min[0],max[0],min[1],max[1],min[2],max[2]);
923: {
924: PetscReal umin,umax,umean;
925: THISurfaceStatistics(dmmg[level]->dm,X,&umin,&umax,&umean);
926: umin *= thi->units->year / thi->units->meter;
927: umax *= thi->units->year / thi->units->meter;
928: umean *= thi->units->year / thi->units->meter;
929: PetscPrintf(comm,"Surface statistics: u in [%12.6e, %12.6e] mean %12.6e\n",umin,umax,umean);
930: }
931: /* These values stay nondimensional */
932: PetscPrintf(comm,"Global eta range %g to %g converged range %g to %g\n",thi->eta.min,thi->eta.max,thi->eta.cmin,thi->eta.cmax);
933: PetscPrintf(comm,"Global beta2 range %g to %g converged range %g to %g\n",thi->beta2.min,thi->beta2.max,thi->beta2.cmin,thi->beta2.cmax);
934: }
935: PetscPrintf(comm,"\n");
936: return(0);
937: }
941: static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo *info,Node **x,Mat B,THI thi)
942: {
943: PetscInt xs,ys,xm,ym,i,j,q,l,ll;
944: PetscReal hx,hy;
945: PrmNode **prm;
949: xs = info->ys;
950: ys = info->xs;
951: xm = info->ym;
952: ym = info->xm;
953: hx = thi->Lx / info->my;
954: hy = thi->Ly / info->mx;
956: MatZeroEntries(B);
957: THIDAGetPrm(info->da,&prm);
959: for (i=xs; i<xs+xm; i++) {
960: for (j=ys; j<ys+ym; j++) {
961: Node n[4];
962: PrmNode pn[4];
963: PetscScalar Ke[4*2][4*2];
964: QuadExtract(prm,i,j,pn);
965: QuadExtract(x,i,j,n);
966: PetscMemzero(Ke,sizeof(Ke));
967: for (q=0; q<4; q++) {
968: PetscReal phi[4],dphi[4][2],jw,eta,deta,beta2,dbeta2;
969: PetscScalar u,v,du[2],dv[2],h = 0,rbeta2 = 0;
970: for (l=0; l<4; l++) {
971: phi[l] = QuadQInterp[q][l];
972: dphi[l][0] = QuadQDeriv[q][l][0]*2./hx;
973: dphi[l][1] = QuadQDeriv[q][l][1]*2./hy;
974: h += phi[l] * pn[l].h;
975: rbeta2 += phi[l] * pn[l].beta2;
976: }
977: jw = 0.25*hx*hy / thi->rhog; /* rhog is only scaling */
978: PointwiseNonlinearity2D(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
979: THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
980: for (l=0; l<4; l++) {
981: const PetscReal pp = phi[l],*dp = dphi[l];
982: for (ll=0; ll<4; ll++) {
983: const PetscReal ppl = phi[ll],*dpl = dphi[ll];
984: PetscScalar dgdu,dgdv;
985: dgdu = 2.*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1]+dv[0])*dpl[1];
986: dgdv = 2.*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1]+dv[0])*dpl[0];
987: /* Picard part */
988: Ke[l*2+0][ll*2+0] += dp[0]*jw*eta*4.*dpl[0] + dp[1]*jw*eta*dpl[1] + pp*jw*(beta2/h)*ppl*thi->ssa_friction_scale;
989: Ke[l*2+0][ll*2+1] += dp[0]*jw*eta*2.*dpl[1] + dp[1]*jw*eta*dpl[0];
990: Ke[l*2+1][ll*2+0] += dp[1]*jw*eta*2.*dpl[0] + dp[0]*jw*eta*dpl[1];
991: Ke[l*2+1][ll*2+1] += dp[1]*jw*eta*4.*dpl[1] + dp[0]*jw*eta*dpl[0] + pp*jw*(beta2/h)*ppl*thi->ssa_friction_scale;
992: /* extra Newton terms */
993: 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]) + pp*jw*(dbeta2/h)*u*u*ppl*thi->ssa_friction_scale;
994: 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]) + pp*jw*(dbeta2/h)*u*v*ppl*thi->ssa_friction_scale;
995: 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]) + pp*jw*(dbeta2/h)*v*u*ppl*thi->ssa_friction_scale;
996: 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]) + pp*jw*(dbeta2/h)*v*v*ppl*thi->ssa_friction_scale;
997: }
998: }
999: }
1000: {
1001: const MatStencil rc[4] = {{0,i,j,0},{0,i+1,j,0},{0,i+1,j+1,0},{0,i,j+1,0}};
1002: MatSetValuesBlockedStencil(B,4,rc,4,rc,&Ke[0][0],ADD_VALUES);
1003: }
1004: }
1005: }
1006: THIDARestorePrm(info->da,&prm);
1008: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1009: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1010: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
1011: if (thi->verbose) {THIMatrixStatistics(thi,B,PETSC_VIEWER_STDOUT_WORLD);}
1012: return(0);
1013: }
1017: static PetscErrorCode THIJacobianLocal_3D(DMDALocalInfo *info,Node ***x,Mat B,THI thi,THIAssemblyMode amode)
1018: {
1019: PetscInt xs,ys,xm,ym,zm,i,j,k,q,l,ll;
1020: PetscReal hx,hy;
1021: PrmNode **prm;
1025: xs = info->zs;
1026: ys = info->ys;
1027: xm = info->zm;
1028: ym = info->ym;
1029: zm = info->xm;
1030: hx = thi->Lx / info->mz;
1031: hy = thi->Ly / info->my;
1033: MatZeroEntries(B);
1034: THIDAGetPrm(info->da,&prm);
1036: for (i=xs; i<xs+xm; i++) {
1037: for (j=ys; j<ys+ym; j++) {
1038: PrmNode pn[4];
1039: QuadExtract(prm,i,j,pn);
1040: for (k=0; k<zm-1; k++) {
1041: Node n[8];
1042: PetscReal zn[8],etabase = 0;
1043: PetscScalar Ke[8*2][8*2];
1044: PetscInt ls = 0;
1046: PrmHexGetZ(pn,k,zm,zn);
1047: HexExtract(x,i,j,k,n);
1048: PetscMemzero(Ke,sizeof(Ke));
1049: if (thi->no_slip && k == 0) {
1050: for (l=0; l<4; l++) n[l].u = n[l].v = 0;
1051: ls = 4;
1052: }
1053: for (q=0; q<8; q++) {
1054: PetscReal dz[3],phi[8],dphi[8][3],jw,eta,deta;
1055: PetscScalar du[3],dv[3],u,v;
1056: HexGrad(HexQDeriv[q],zn,dz);
1057: HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
1058: PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
1059: jw /= thi->rhog; /* residuals are scaled by this factor */
1060: if (q == 0) etabase = eta;
1061: for (l=ls; l<8; l++) { /* test functions */
1062: const PetscReal *restrict dp = dphi[l];
1063: #if USE_SSE2_KERNELS
1064: /* gcc (up to my 4.5 snapshot) is really bad at hoisting intrinsics so we do it manually */
1065: __m128d
1066: p4 = _mm_set1_pd(4),p2 = _mm_set1_pd(2),p05 = _mm_set1_pd(0.5),
1067: p42 = _mm_setr_pd(4,2),p24 = _mm_shuffle_pd(p42,p42,_MM_SHUFFLE2(0,1)),
1068: du0 = _mm_set1_pd(du[0]),du1 = _mm_set1_pd(du[1]),du2 = _mm_set1_pd(du[2]),
1069: dv0 = _mm_set1_pd(dv[0]),dv1 = _mm_set1_pd(dv[1]),dv2 = _mm_set1_pd(dv[2]),
1070: jweta = _mm_set1_pd(jw*eta),jwdeta = _mm_set1_pd(jw*deta),
1071: dp0 = _mm_set1_pd(dp[0]),dp1 = _mm_set1_pd(dp[1]),dp2 = _mm_set1_pd(dp[2]),
1072: dp0jweta = _mm_mul_pd(dp0,jweta),dp1jweta = _mm_mul_pd(dp1,jweta),dp2jweta = _mm_mul_pd(dp2,jweta),
1073: p4du0p2dv1 = _mm_add_pd(_mm_mul_pd(p4,du0),_mm_mul_pd(p2,dv1)), /* 4 du0 + 2 dv1 */
1074: p4dv1p2du0 = _mm_add_pd(_mm_mul_pd(p4,dv1),_mm_mul_pd(p2,du0)), /* 4 dv1 + 2 du0 */
1075: pdu2dv2 = _mm_unpacklo_pd(du2,dv2), /* [du2, dv2] */
1076: du1pdv0 = _mm_add_pd(du1,dv0), /* du1 + dv0 */
1077: t1 = _mm_mul_pd(dp0,p4du0p2dv1), /* dp0 (4 du0 + 2 dv1) */
1078: t2 = _mm_mul_pd(dp1,p4dv1p2du0); /* dp1 (4 dv1 + 2 du0) */
1080: #endif
1081: #if defined COMPUTE_LOWER_TRIANGULAR /* The element matrices are always symmetric so computing the lower-triangular part is not necessary */
1082: for (ll=ls; ll<8; ll++) { /* trial functions */
1083: #else
1084: for (ll=l; ll<8; ll++) {
1085: #endif
1086: const PetscReal *restrict dpl = dphi[ll];
1087: if (amode == THIASSEMBLY_TRIDIAGONAL && (l-ll)%4) continue; /* these entries would not be inserted */
1088: #if !USE_SSE2_KERNELS
1089: /* The analytic Jacobian in nice, easy-to-read form */
1090: {
1091: PetscScalar dgdu,dgdv;
1092: dgdu = 2.*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1]+dv[0])*dpl[1] + 0.5*du[2]*dpl[2];
1093: dgdv = 2.*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1]+dv[0])*dpl[0] + 0.5*dv[2]*dpl[2];
1094: /* Picard part */
1095: 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];
1096: Ke[l*2+0][ll*2+1] += dp[0]*jw*eta*2.*dpl[1] + dp[1]*jw*eta*dpl[0];
1097: Ke[l*2+1][ll*2+0] += dp[1]*jw*eta*2.*dpl[0] + dp[0]*jw*eta*dpl[1];
1098: 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];
1099: /* extra Newton terms */
1100: 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];
1101: 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];
1102: 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];
1103: 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];
1104: }
1105: #else
1106: /* This SSE2 code is an exact replica of above, but uses explicit packed instructions for some speed
1107: * benefit. On my hardware, these intrinsics are almost twice as fast as above, reducing total assembly cost
1108: * by 25 to 30 percent. */
1109: {
1110: __m128d
1111: keu = _mm_loadu_pd(&Ke[l*2+0][ll*2+0]),
1112: kev = _mm_loadu_pd(&Ke[l*2+1][ll*2+0]),
1113: dpl01 = _mm_loadu_pd(&dpl[0]),dpl10 = _mm_shuffle_pd(dpl01,dpl01,_MM_SHUFFLE2(0,1)),dpl2 = _mm_set_sd(dpl[2]),
1114: t0,t3,pdgduv;
1115: keu = _mm_add_pd(keu,_mm_add_pd(_mm_mul_pd(_mm_mul_pd(dp0jweta,p42),dpl01),
1116: _mm_add_pd(_mm_mul_pd(dp1jweta,dpl10),
1117: _mm_mul_pd(dp2jweta,dpl2))));
1118: kev = _mm_add_pd(kev,_mm_add_pd(_mm_mul_pd(_mm_mul_pd(dp1jweta,p24),dpl01),
1119: _mm_add_pd(_mm_mul_pd(dp0jweta,dpl10),
1120: _mm_mul_pd(dp2jweta,_mm_shuffle_pd(dpl2,dpl2,_MM_SHUFFLE2(0,1))))));
1121: pdgduv = _mm_mul_pd(p05,_mm_add_pd(_mm_add_pd(_mm_mul_pd(p42,_mm_mul_pd(du0,dpl01)),
1122: _mm_mul_pd(p24,_mm_mul_pd(dv1,dpl01))),
1123: _mm_add_pd(_mm_mul_pd(du1pdv0,dpl10),
1124: _mm_mul_pd(pdu2dv2,_mm_set1_pd(dpl[2]))))); /* [dgdu, dgdv] */
1125: t0 = _mm_mul_pd(jwdeta,pdgduv); /* jw deta [dgdu, dgdv] */
1126: t3 = _mm_mul_pd(t0,du1pdv0); /* t0 (du1 + dv0) */
1127: _mm_storeu_pd(&Ke[l*2+0][ll*2+0],_mm_add_pd(keu,_mm_add_pd(_mm_mul_pd(t1,t0),
1128: _mm_add_pd(_mm_mul_pd(dp1,t3),
1129: _mm_mul_pd(t0,_mm_mul_pd(dp2,du2))))));
1130: _mm_storeu_pd(&Ke[l*2+1][ll*2+0],_mm_add_pd(kev,_mm_add_pd(_mm_mul_pd(t2,t0),
1131: _mm_add_pd(_mm_mul_pd(dp0,t3),
1132: _mm_mul_pd(t0,_mm_mul_pd(dp2,dv2))))));
1133: }
1134: #endif
1135: }
1136: }
1137: }
1138: if (k == 0) { /* on a bottom face */
1139: if (thi->no_slip) {
1140: const PetscReal hz = PetscRealPart(pn[0].h)/(zm-1);
1141: 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);
1142: Ke[0][0] = thi->dirichlet_scale*diagu;
1143: Ke[1][1] = thi->dirichlet_scale*diagv;
1144: } else {
1145: for (q=0; q<4; q++) {
1146: const PetscReal jw = 0.25*hx*hy/thi->rhog,*phi = QuadQInterp[q];
1147: PetscScalar u=0,v=0,rbeta2=0;
1148: PetscReal beta2,dbeta2;
1149: for (l=0; l<4; l++) {
1150: u += phi[l]*n[l].u;
1151: v += phi[l]*n[l].v;
1152: rbeta2 += phi[l]*pn[l].beta2;
1153: }
1154: THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
1155: for (l=0; l<4; l++) {
1156: const PetscReal pp = phi[l];
1157: for (ll=0; ll<4; ll++) {
1158: const PetscReal ppl = phi[ll];
1159: Ke[l*2+0][ll*2+0] += pp*jw*beta2*ppl + pp*jw*dbeta2*u*u*ppl;
1160: Ke[l*2+0][ll*2+1] += pp*jw*dbeta2*u*v*ppl;
1161: Ke[l*2+1][ll*2+0] += pp*jw*dbeta2*v*u*ppl;
1162: Ke[l*2+1][ll*2+1] += pp*jw*beta2*ppl + pp*jw*dbeta2*v*v*ppl;
1163: }
1164: }
1165: }
1166: }
1167: }
1168: {
1169: const MatStencil rc[8] = {{i,j,k,0},{i+1,j,k,0},{i+1,j+1,k,0},{i,j+1,k,0},{i,j,k+1,0},{i+1,j,k+1,0},{i+1,j+1,k+1,0},{i,j+1,k+1,0}};
1170: if (amode == THIASSEMBLY_TRIDIAGONAL) {
1171: for (l=0; l<4; l++) { /* Copy out each of the blocks, discarding horizontal coupling */
1172: const PetscInt l4 = l+4;
1173: const MatStencil rcl[2] = {{rc[l].k,rc[l].j,rc[l].i,0},{rc[l4].k,rc[l4].j,rc[l4].i,0}};
1174: #if defined COMPUTE_LOWER_TRIANGULAR
1175: const PetscScalar Kel[4][4] = {{Ke[2*l+0][2*l+0] ,Ke[2*l+0][2*l+1] ,Ke[2*l+0][2*l4+0] ,Ke[2*l+0][2*l4+1]},
1176: {Ke[2*l+1][2*l+0] ,Ke[2*l+1][2*l+1] ,Ke[2*l+1][2*l4+0] ,Ke[2*l+1][2*l4+1]},
1177: {Ke[2*l4+0][2*l+0],Ke[2*l4+0][2*l+1],Ke[2*l4+0][2*l4+0],Ke[2*l4+0][2*l4+1]},
1178: {Ke[2*l4+1][2*l+0],Ke[2*l4+1][2*l+1],Ke[2*l4+1][2*l4+0],Ke[2*l4+1][2*l4+1]}};
1179: #else
1180: /* Same as above except for the lower-left block */
1181: const PetscScalar Kel[4][4] = {{Ke[2*l+0][2*l+0] ,Ke[2*l+0][2*l+1] ,Ke[2*l+0][2*l4+0] ,Ke[2*l+0][2*l4+1]},
1182: {Ke[2*l+1][2*l+0] ,Ke[2*l+1][2*l+1] ,Ke[2*l+1][2*l4+0] ,Ke[2*l+1][2*l4+1]},
1183: {Ke[2*l+0][2*l4+0],Ke[2*l+1][2*l4+0],Ke[2*l4+0][2*l4+0],Ke[2*l4+0][2*l4+1]},
1184: {Ke[2*l+0][2*l4+1],Ke[2*l+1][2*l4+1],Ke[2*l4+1][2*l4+0],Ke[2*l4+1][2*l4+1]}};
1185: #endif
1186: MatSetValuesBlockedStencil(B,2,rcl,2,rcl,&Kel[0][0],ADD_VALUES);
1187: }
1188: } else {
1189: #if !defined COMPUTE_LOWER_TRIANGULAR /* fill in lower-triangular part, this is really cheap compared to computing the entries */
1190: for (l=0; l<8; l++) {
1191: for (ll=l+1; ll<8; ll++) {
1192: Ke[ll*2+0][l*2+0] = Ke[l*2+0][ll*2+0];
1193: Ke[ll*2+1][l*2+0] = Ke[l*2+0][ll*2+1];
1194: Ke[ll*2+0][l*2+1] = Ke[l*2+1][ll*2+0];
1195: Ke[ll*2+1][l*2+1] = Ke[l*2+1][ll*2+1];
1196: }
1197: }
1198: #endif
1199: MatSetValuesBlockedStencil(B,8,rc,8,rc,&Ke[0][0],ADD_VALUES);
1200: }
1201: }
1202: }
1203: }
1204: }
1205: THIDARestorePrm(info->da,&prm);
1207: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1208: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1209: MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
1210: if (thi->verbose) {THIMatrixStatistics(thi,B,PETSC_VIEWER_STDOUT_WORLD);}
1211: return(0);
1212: }
1216: static PetscErrorCode THIJacobianLocal_3D_Full(DMDALocalInfo *info,Node ***x,Mat B,THI thi)
1217: {
1221: THIJacobianLocal_3D(info,x,B,thi,THIASSEMBLY_FULL);
1222: return(0);
1223: }
1227: static PetscErrorCode THIJacobianLocal_3D_Tridiagonal(DMDALocalInfo *info,Node ***x,Mat B,THI thi)
1228: {
1232: THIJacobianLocal_3D(info,x,B,thi,THIASSEMBLY_TRIDIAGONAL);
1233: return(0);
1234: }
1238: static PetscErrorCode DMRefineHierarchy_THI(DM dac0,PetscInt nlevels,DM hierarchy[])
1239: {
1241: THI thi;
1242: PetscInt dim,M,N,m,n,s,dof;
1243: DM dac,daf;
1244: DMDAStencilType st;
1245: DM_DA *ddf,*ddc;
1248: PetscObjectQuery((PetscObject)dac0,"THI",(PetscObject*)&thi);
1249: if (!thi) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot refine this DMDA, missing composed THI instance");
1250: if (nlevels > 1) {
1251: DMRefineHierarchy(dac0,nlevels-1,hierarchy);
1252: dac = hierarchy[nlevels-2];
1253: } else {
1254: dac = dac0;
1255: }
1256: DMDAGetInfo(dac,&dim, &N,&M,0, &n,&m,0, &dof,&s,0,0,0,&st);
1257: if (dim != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"This function can only refine 2D DMDAs");
1258: /* Creates a 3D DMDA with the same map-plane layout as the 2D one, with contiguous columns */
1259: DMDACreate3d(((PetscObject)dac)->comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,st,thi->zlevels,N,M,1,n,m,dof,s,PETSC_NULL,PETSC_NULL,PETSC_NULL,&daf);
1260: daf->ops->creatematrix = dac->ops->creatematrix;
1261: daf->ops->createinterpolation = dac->ops->createinterpolation;
1262: daf->ops->getcoloring = dac->ops->getcoloring;
1263: ddf = (DM_DA*)daf->data;
1264: ddc = (DM_DA*)dac->data;
1265: ddf->interptype = ddc->interptype;
1267: DMDASetFieldName(daf,0,"x-velocity");
1268: DMDASetFieldName(daf,1,"y-velocity");
1269: hierarchy[nlevels-1] = daf;
1270: return(0);
1271: }
1275: static PetscErrorCode DMCreateInterpolation_DA_THI(DM dac,DM daf,Mat *A,Vec *scale)
1276: {
1278: PetscInt dim;
1285: DMDAGetInfo(daf,&dim,0,0,0,0,0,0,0,0,0,0,0,0);
1286: if (dim == 2) {
1287: /* We are in the 2D problem and use normal DMDA interpolation */
1288: DMCreateInterpolation(dac,daf,A,scale);
1289: } else {
1290: PetscInt i,j,k,xs,ys,zs,xm,ym,zm,mx,my,mz,rstart,cstart;
1291: Mat B;
1293: DMDAGetInfo(daf,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
1294: DMDAGetCorners(daf,&zs,&ys,&xs,&zm,&ym,&xm);
1295: if (zs != 0) SETERRQ(PETSC_COMM_SELF,1,"unexpected");
1296: MatCreate(((PetscObject)daf)->comm,&B);
1297: MatSetSizes(B,xm*ym*zm,xm*ym,mx*my*mz,mx*my);
1298:
1299: MatSetType(B,MATAIJ);
1300: MatSeqAIJSetPreallocation(B,1,NULL);
1301: MatMPIAIJSetPreallocation(B,1,NULL,0,NULL);
1302: MatGetOwnershipRange(B,&rstart,NULL);
1303: MatGetOwnershipRangeColumn(B,&cstart,NULL);
1304: for (i=xs; i<xs+xm; i++) {
1305: for (j=ys; j<ys+ym; j++) {
1306: for (k=zs; k<zs+zm; k++) {
1307: PetscInt i2 = i*ym+j,i3 = i2*zm+k;
1308: PetscScalar val = ((k == 0 || k == mz-1) ? 0.5 : 1.) / (mz-1.); /* Integration using trapezoid rule */
1309: MatSetValue(B,cstart+i3,rstart+i2,val,INSERT_VALUES);
1310: }
1311: }
1312: }
1313: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1314: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1315: MatCreateMAIJ(B,sizeof(Node)/sizeof(PetscScalar),A);
1316: MatDestroy(&B);
1317: }
1318: return(0);
1319: }
1323: static PetscErrorCode DMCreateMatrix_THI_Tridiagonal(DM da,const MatType mtype,Mat *J)
1324: {
1326: Mat A;
1327: PetscInt xm,ym,zm,dim,dof = 2,starts[3],dims[3];
1328: ISLocalToGlobalMapping ltog,ltogb;
1331: DMDAGetInfo(da,&dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
1332: if (dim != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected DMDA to be 3D");
1333: DMDAGetCorners(da,0,0,0,&zm,&ym,&xm);
1334: DMGetLocalToGlobalMapping(da,<og);
1335: DMGetLocalToGlobalMappingBlock(da,<ogb);
1336: MatCreate(((PetscObject)da)->comm,&A);
1337: MatSetSizes(A,dof*xm*ym*zm,dof*xm*ym*zm,PETSC_DETERMINE,PETSC_DETERMINE);
1338: MatSetType(A,mtype);
1339: MatSetFromOptions(A);
1340: MatXAIJSetPreallocation(A,2,3,PETSC_NULL,0,PETSC_NULL,2,PETSC_NULL,0,PETSC_NULL);
1341: MatSetLocalToGlobalMapping(A,ltog,ltog);
1342: MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);
1343: DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
1344: MatSetStencil(A,dim,dims,starts,dof);
1345: *J = A;
1346: return(0);
1347: }
1351: static PetscErrorCode THIDAVecView_VTK_XML(THI thi,DM da,Vec X,const char filename[])
1352: {
1353: const PetscInt dof = 2;
1354: Units units = thi->units;
1355: MPI_Comm comm;
1357: PetscViewer viewer;
1358: PetscMPIInt rank,size,tag,nn,nmax;
1359: PetscInt mx,my,mz,r,range[6];
1360: PetscScalar *x;
1363: comm = ((PetscObject)thi)->comm;
1364: DMDAGetInfo(da,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
1365: MPI_Comm_size(comm,&size);
1366: MPI_Comm_rank(comm,&rank);
1367: PetscViewerASCIIOpen(comm,filename,&viewer);
1368: PetscViewerASCIIPrintf(viewer,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1369: PetscViewerASCIIPrintf(viewer," <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n",0,mz-1,0,my-1,0,mx-1);
1371: DMDAGetCorners(da,range,range+1,range+2,range+3,range+4,range+5);
1372: nn = PetscMPIIntCast(range[3]*range[4]*range[5]*dof);
1373: MPI_Reduce(&nn,&nmax,1,MPI_INT,MPI_MAX,0,comm);
1374: tag = ((PetscObject) viewer)->tag;
1375: VecGetArray(X,&x);
1376: if (!rank) {
1377: PetscScalar *array;
1378: PetscMalloc(nmax*sizeof(PetscScalar),&array);
1379: for (r=0; r<size; r++) {
1380: PetscInt i,j,k,xs,xm,ys,ym,zs,zm;
1381: PetscScalar *ptr;
1382: MPI_Status status;
1383: if (r) {
1384: MPI_Recv(range,6,MPIU_INT,r,tag,comm,MPI_STATUS_IGNORE);
1385: }
1386: zs = range[0];ys = range[1];xs = range[2];zm = range[3];ym = range[4];xm = range[5];
1387: if (xm*ym*zm*dof > nmax) SETERRQ(PETSC_COMM_SELF,1,"should not happen");
1388: if (r) {
1389: MPI_Recv(array,nmax,MPIU_SCALAR,r,tag,comm,&status);
1390: MPI_Get_count(&status,MPIU_SCALAR,&nn);
1391: if (nn != xm*ym*zm*dof) SETERRQ(PETSC_COMM_SELF,1,"should not happen");
1392: ptr = array;
1393: } else ptr = x;
1394: PetscViewerASCIIPrintf(viewer," <Piece Extent=\"%d %d %d %d %d %d\">\n",zs,zs+zm-1,ys,ys+ym-1,xs,xs+xm-1);
1396: PetscViewerASCIIPrintf(viewer," <Points>\n");
1397: PetscViewerASCIIPrintf(viewer," <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n");
1398: for (i=xs; i<xs+xm; i++) {
1399: for (j=ys; j<ys+ym; j++) {
1400: for (k=zs; k<zs+zm; k++) {
1401: PrmNode p;
1402: PetscReal xx = thi->Lx*i/mx,yy = thi->Ly*j/my,zz;
1403: thi->initialize(thi,xx,yy,&p);
1404: zz = PetscRealPart(p.b) + PetscRealPart(p.h)*k/(mz-1);
1405: PetscViewerASCIIPrintf(viewer,"%f %f %f\n",xx,yy,zz);
1406: }
1407: }
1408: }
1409: PetscViewerASCIIPrintf(viewer," </DataArray>\n");
1410: PetscViewerASCIIPrintf(viewer," </Points>\n");
1412: PetscViewerASCIIPrintf(viewer," <PointData>\n");
1413: PetscViewerASCIIPrintf(viewer," <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n");
1414: for (i=0; i<nn; i+=dof) {
1415: PetscViewerASCIIPrintf(viewer,"%f %f %f\n",PetscRealPart(ptr[i])*units->year/units->meter,PetscRealPart(ptr[i+1])*units->year/units->meter,0.0);
1416: }
1417: PetscViewerASCIIPrintf(viewer," </DataArray>\n");
1419: PetscViewerASCIIPrintf(viewer," <DataArray type=\"Int32\" Name=\"rank\" NumberOfComponents=\"1\" format=\"ascii\">\n");
1420: for (i=0; i<nn; i+=dof) {
1421: PetscViewerASCIIPrintf(viewer,"%d\n",r);
1422: }
1423: PetscViewerASCIIPrintf(viewer," </DataArray>\n");
1424: PetscViewerASCIIPrintf(viewer," </PointData>\n");
1426: PetscViewerASCIIPrintf(viewer," </Piece>\n");
1427: }
1428: PetscFree(array);
1429: } else {
1430: MPI_Send(range,6,MPIU_INT,0,tag,comm);
1431: MPI_Send(x,nn,MPIU_SCALAR,0,tag,comm);
1432: }
1433: VecRestoreArray(X,&x);
1434: PetscViewerASCIIPrintf(viewer," </StructuredGrid>\n");
1435: PetscViewerASCIIPrintf(viewer,"</VTKFile>\n");
1436: PetscViewerDestroy(&viewer);
1437: return(0);
1438: }
1442: int main(int argc,char *argv[])
1443: {
1444: MPI_Comm comm;
1445: DMMG *dmmg;
1446: THI thi;
1447: PetscInt i;
1449: PETSC_UNUSED PetscLogStage stages[3];
1450: PetscBool repeat_fine_solve = PETSC_FALSE;
1452: PetscInitialize(&argc,&argv,0,help);
1453: comm = PETSC_COMM_WORLD;
1455: /* We define two stages. The first includes all setup costs and solves from a naive initial guess. The second solve
1456: * is more indicative of what might occur during time-stepping. The initial guess is interpolated from the next
1457: * coarser (as in the last step of grid sequencing), and so requires fewer Newton steps. */
1458: PetscOptionsGetBool(NULL,"-repeat_fine_solve",&repeat_fine_solve,NULL);
1459: PetscLogStageRegister("Full solve",&stages[0]);
1460: if (repeat_fine_solve) {
1461: PetscLogStageRegister("Fine-1 solve",&stages[1]);
1462: PetscLogStageRegister("Fine-only solve",&stages[2]);
1463: }
1465: PetscLogStagePush(stages[0]);
1467: THICreate(comm,&thi);
1468: DMMGCreate(PETSC_COMM_WORLD,thi->nlevels,thi,&dmmg);
1469: {
1470: DM da;
1471: PetscInt M = 3,N = 3,P = 2;
1472: PetscOptionsBegin(comm,NULL,"Grid resolution options","");
1473: {
1474: PetscOptionsInt("-M","Number of elements in x-direction on coarse level","",M,&M,NULL);
1475: N = M;
1476: PetscOptionsInt("-N","Number of elements in y-direction on coarse level (if different from M)","",N,&N,NULL);
1477: if (thi->coarse2d) {
1478: PetscOptionsInt("-zlevels","Number of elements in z-direction on fine level","",thi->zlevels,&thi->zlevels,NULL);
1479: } else {
1480: PetscOptionsInt("-P","Number of elements in z-direction on coarse level","",P,&P,NULL);
1481: }
1482: }
1483: PetscOptionsEnd();
1484: if (thi->coarse2d) {
1485: DMDACreate2d(comm,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,N,M,PETSC_DETERMINE,PETSC_DETERMINE,sizeof(Node)/sizeof(PetscScalar),1,0,0,&da);
1486: da->ops->refinehierarchy = DMRefineHierarchy_THI;
1487: da->ops->createinterpolation = DMCreateInterpolation_DA_THI;
1488: PetscObjectCompose((PetscObject)da,"THI",(PetscObject)thi);
1489: } else {
1490: 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);
1491: }
1492: DMDASetFieldName(da,0,"x-velocity");
1493: DMDASetFieldName(da,1,"y-velocity");
1494: DMMGSetDM(dmmg,(DM)da);
1495: DMDestroy(&da);
1496: }
1497: if (thi->tridiagonal) {
1498: (DMMGGetDM(dmmg))->ops->creatematrix = DMCreateMatrix_THI_Tridiagonal;
1499: }
1500: {
1501: /* Use the user-defined matrix type on all but the coarse level */
1502: DMMGSetMatType(dmmg,thi->mattype);
1503: /* PCREDUNDANT only works with AIJ, and so do the third-party direct solvers. So when running in parallel, we can't
1504: * use the faster (S)BAIJ formats on the coarse level. */
1505: PetscFree(dmmg[0]->mtype);
1506: PetscStrallocpy(MATAIJ,&dmmg[0]->mtype);
1507: }
1508: PetscOptionsSetValue("-dmmg_form_function_ghost","1"); /* Spectacularly ugly API, our function evaluation provides ghost values */
1509: DMMGSetSNESLocal(dmmg,THIFunctionLocal,THIJacobianLocal_3D_Full,0,0);
1510: if (thi->tridiagonal) {
1511: DMDASetLocalJacobian(DMMGGetDM(dmmg),(DMDALocalFunction1)THIJacobianLocal_3D_Tridiagonal);
1512: }
1513: if (thi->coarse2d) {
1514: for (i=0; i<DMMGGetLevels(dmmg)-1; i++) {
1515: DMDASetLocalJacobian(dmmg[i]->dm,(DMDALocalFunction1)THIJacobianLocal_2D);
1516: }
1517: }
1518: for (i=0; i<DMMGGetLevels(dmmg); i++) {
1519: /* This option is only valid for the SBAIJ format. The matrices we assemble are symmetric, but the SBAIJ assembly
1520: * functions will complain if we provide lower-triangular entries without setting this option. */
1521: Mat B = dmmg[i]->B;
1522: PetscBool flg1,flg2;
1523: PetscTypeCompare((PetscObject)B,MATSEQSBAIJ,&flg1);
1524: PetscTypeCompare((PetscObject)B,MATMPISBAIJ,&flg2);
1525: if (flg1 || flg2) {
1526: MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
1527: }
1528: }
1529: MatSetOptionsPrefix(DMMGGetB(dmmg),"thi_");
1530: DMMGSetFromOptions(dmmg);
1531: THISetDMMG(thi,dmmg);
1533: DMMGSetInitialGuess(dmmg,THIInitial);
1534: DMMGSolve(dmmg);
1536: PetscLogStagePop();
1537: THISolveStatistics(thi,dmmg,0,"Full");
1538: /* The first solve is complete */
1540: if (repeat_fine_solve && DMMGGetLevels(dmmg) > 1) {
1541: PetscInt nlevels = DMMGGetLevels(dmmg);
1542: DMMG dmmgc = dmmg[nlevels-2],dmmgf = dmmg[nlevels-1];
1543: Vec Xc = dmmgc->x,Xf = dmmgf->x;
1544: MatRestrict(dmmgf->R,Xf,Xc);
1545: VecPointwiseMult(Xc,Xc,dmmgf->Rscale);
1547: /* Solve on the level with one coarsening, this is a more stringent test of latency */
1548: PetscLogStagePush(stages[1]);
1549: (*dmmgc->solve)(dmmg,nlevels-2);
1550: PetscLogStagePop();
1551: THISolveStatistics(thi,dmmg,1,"Fine-1");
1553: MatInterpolate(dmmgf->R,Xc,Xf);
1555: /* Solve again on the finest level, this is representative of what is needed in a time-stepping code */
1556: PetscLogStagePush(stages[2]);
1557: (*dmmgf->solve)(dmmg,nlevels-1);
1558: PetscLogStagePop();
1559: THISolveStatistics(thi,dmmg,0,"Fine");
1560: }
1562: {
1563: PetscBool flg;
1564: char filename[PETSC_MAX_PATH_LEN] = "";
1565: PetscOptionsGetString(PETSC_NULL,"-o",filename,sizeof(filename),&flg);
1566: if (flg) {
1567: THIDAVecView_VTK_XML(thi,DMMGGetDM(dmmg),DMMGGetx(dmmg),filename);
1568: }
1569: }
1571: DMMGDestroy(dmmg);
1572: THIDestroy(&thi);
1573: PetscFinalize();
1574: return 0;
1575: }