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