Actual source code: ex48.c
petsc-3.4.2 2013-07-02

      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 <petscsnes.h>
     59: #include <ctype.h>              /* toupper() */
     60: #include <petsc-private/dmdaimpl.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 PETSC_USE_REAL___FLOAT128         \
     78:                           && defined __SSE2__)

     80: static PetscClassId THI_CLASSID;

     82: typedef enum {QUAD_GAUSS,QUAD_LOBATTO} QuadratureType;
     83: static const char      *QuadratureTypes[] = {"gauss","lobatto","QuadratureType","QUAD_",0};
     84: static const PetscReal HexQWeights[8]     = {1,1,1,1,1,1,1,1};
     85: static const PetscReal HexQNodes[]        = {-0.57735026918962573, 0.57735026918962573};
     86: #define G 0.57735026918962573
     87: #define H (0.5*(1.+G))
     88: #define L (0.5*(1.-G))
     89: #define M (-0.5)
     90: #define P (0.5)
     91: /* Special quadrature: Lobatto in horizontal, Gauss in vertical */
     92: static const PetscReal HexQInterp_Lobatto[8][8] = {{H,0,0,0,L,0,0,0},
     93:                                                    {0,H,0,0,0,L,0,0},
     94:                                                    {0,0,H,0,0,0,L,0},
     95:                                                    {0,0,0,H,0,0,0,L},
     96:                                                    {L,0,0,0,H,0,0,0},
     97:                                                    {0,L,0,0,0,H,0,0},
     98:                                                    {0,0,L,0,0,0,H,0},
     99:                                                    {0,0,0,L,0,0,0,H}};
    100: static const PetscReal HexQDeriv_Lobatto[8][8][3] = {
    101:   {{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}  },
    102:   {{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}    },
    103:   {{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}  },
    104:   {{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}},
    105:   {{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}  },
    106:   {{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}    },
    107:   {{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}  },
    108:   {{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}}};
    109: /* Stanndard Gauss */
    110: 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},
    111:                                                  {L*H*H,H*H*H,H*L*H,L*L*H, L*H*L,H*H*L,H*L*L,L*L*L},
    112:                                                  {L*L*H,H*L*H,H*H*H,L*H*H, L*L*L,H*L*L,H*H*L,L*H*L},
    113:                                                  {H*L*H,L*L*H,L*H*H,H*H*H, H*L*L,L*L*L,L*H*L,H*H*L},
    114:                                                  {H*H*L,L*H*L,L*L*L,H*L*L, H*H*H,L*H*H,L*L*H,H*L*H},
    115:                                                  {L*H*L,H*H*L,H*L*L,L*L*L, L*H*H,H*H*H,H*L*H,L*L*H},
    116:                                                  {L*L*L,H*L*L,H*H*L,L*H*L, L*L*H,H*L*H,H*H*H,L*H*H},
    117:                                                  {H*L*L,L*L*L,L*H*L,H*H*L, H*L*H,L*L*H,L*H*H,H*H*H}};
    118: static const PetscReal HexQDeriv_Gauss[8][8][3] = {
    119:   {{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}},
    120:   {{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}},
    121:   {{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}},
    122:   {{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}},
    123:   {{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}},
    124:   {{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}},
    125:   {{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}},
    126:   {{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}}};
    127: static const PetscReal (*HexQInterp)[8],(*HexQDeriv)[8][3];
    128: /* Standard 2x2 Gauss quadrature for the bottom layer. */
    129: static const PetscReal QuadQInterp[4][4] = {{H*H,L*H,L*L,H*L},
    130:                                             {L*H,H*H,H*L,L*L},
    131:                                             {L*L,H*L,H*H,L*H},
    132:                                             {H*L,L*L,L*H,H*H}};
    133: static const PetscReal QuadQDeriv[4][4][2] = {
    134:   {{M*H,M*H},{P*H,M*L},{P*L,P*L},{M*L,P*H}},
    135:   {{M*H,M*L},{P*H,M*H},{P*L,P*H},{M*L,P*L}},
    136:   {{M*L,M*L},{P*L,M*H},{P*H,P*H},{M*H,P*L}},
    137:   {{M*L,M*H},{P*L,M*L},{P*H,P*L},{M*H,P*H}}};
    138: #undef G
    139: #undef H
    140: #undef L
    141: #undef M
    142: #undef P

    144: #define HexExtract(x,i,j,k,n) do {              \
    145:     (n)[0] = (x)[i][j][k];                      \
    146:     (n)[1] = (x)[i+1][j][k];                    \
    147:     (n)[2] = (x)[i+1][j+1][k];                  \
    148:     (n)[3] = (x)[i][j+1][k];                    \
    149:     (n)[4] = (x)[i][j][k+1];                    \
    150:     (n)[5] = (x)[i+1][j][k+1];                  \
    151:     (n)[6] = (x)[i+1][j+1][k+1];                \
    152:     (n)[7] = (x)[i][j+1][k+1];                  \
    153:   } while (0)

    155: #define HexExtractRef(x,i,j,k,n) do {           \
    156:     (n)[0] = &(x)[i][j][k];                     \
    157:     (n)[1] = &(x)[i+1][j][k];                   \
    158:     (n)[2] = &(x)[i+1][j+1][k];                 \
    159:     (n)[3] = &(x)[i][j+1][k];                   \
    160:     (n)[4] = &(x)[i][j][k+1];                   \
    161:     (n)[5] = &(x)[i+1][j][k+1];                 \
    162:     (n)[6] = &(x)[i+1][j+1][k+1];               \
    163:     (n)[7] = &(x)[i][j+1][k+1];                 \
    164:   } while (0)

    166: #define QuadExtract(x,i,j,n) do {               \
    167:     (n)[0] = (x)[i][j];                         \
    168:     (n)[1] = (x)[i+1][j];                       \
    169:     (n)[2] = (x)[i+1][j+1];                     \
    170:     (n)[3] = (x)[i][j+1];                       \
    171:   } while (0)

    173: static PetscScalar Sqr(PetscScalar a) {return a*a;}

    175: static void HexGrad(const PetscReal dphi[][3],const PetscReal zn[],PetscReal dz[])
    176: {
    177:   PetscInt i;
    178:   dz[0] = dz[1] = dz[2] = 0;
    179:   for (i=0; i<8; i++) {
    180:     dz[0] += dphi[i][0] * zn[i];
    181:     dz[1] += dphi[i][1] * zn[i];
    182:     dz[2] += dphi[i][2] * zn[i];
    183:   }
    184: }

    186: static void HexComputeGeometry(PetscInt q,PetscReal hx,PetscReal hy,const PetscReal dz[restrict],PetscReal phi[restrict],PetscReal dphi[restrict][3],PetscReal *restrict jw)
    187: {
    188:   const PetscReal jac[3][3]  = {{hx/2,0,0}, {0,hy/2,0}, {dz[0],dz[1],dz[2]}};
    189:   const PetscReal 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:   const PetscReal 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  zlevels;
    226:   PetscReal Lx,Ly,Lz;           /* Model domain */
    227:   PetscReal alpha;              /* Bed angle */
    228:   Units     units;
    229:   PetscReal dirichlet_scale;
    230:   PetscReal ssa_friction_scale;
    231:   PRange    eta;
    232:   PRange    beta2;
    233:   struct {
    234:     PetscReal Bd2,eps,exponent;
    235:   } viscosity;
    236:   struct {
    237:     PetscReal irefgam,eps2,exponent,refvel,epsvel;
    238:   } friction;
    239:   PetscReal rhog;
    240:   PetscBool no_slip;
    241:   PetscBool tridiagonal;
    242:   PetscBool coarse2d;
    243:   PetscBool verbose;
    244:   MatType   mattype;
    245: };

    247: struct _n_Units {
    248:   /* fundamental */
    249:   PetscReal meter;
    250:   PetscReal kilogram;
    251:   PetscReal second;
    252:   /* derived */
    253:   PetscReal Pascal;
    254:   PetscReal year;
    255: };

    257: typedef PetscErrorCode (*DMDASNESJacobianLocal)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*);
    258: typedef PetscErrorCode (*DMDASNESFunctionLocal)(DMDALocalInfo*,void*,void*,void*);
    259: static PetscErrorCode THIJacobianLocal_3D_Full(DMDALocalInfo*,Node***,Mat,Mat,MatStructure*,THI);
    260: static PetscErrorCode THIJacobianLocal_3D_Tridiagonal(DMDALocalInfo*,Node***,Mat,THI);
    261: static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo*,Node**,Mat,THI);

    263: static void PrmHexGetZ(const PrmNode pn[],PetscInt k,PetscInt zm,PetscReal zn[])
    264: {
    265:   const PetscScalar zm1    = zm-1,
    266:                     znl[8] = {pn[0].b + pn[0].h*(PetscScalar)k/zm1,
    267:                               pn[1].b + pn[1].h*(PetscScalar)k/zm1,
    268:                               pn[2].b + pn[2].h*(PetscScalar)k/zm1,
    269:                               pn[3].b + pn[3].h*(PetscScalar)k/zm1,
    270:                               pn[0].b + pn[0].h*(PetscScalar)(k+1)/zm1,
    271:                               pn[1].b + pn[1].h*(PetscScalar)(k+1)/zm1,
    272:                               pn[2].b + pn[2].h*(PetscScalar)(k+1)/zm1,
    273:                               pn[3].b + pn[3].h*(PetscScalar)(k+1)/zm1};
    274:   PetscInt i;
    275:   for (i=0; i<8; i++) zn[i] = PetscRealPart(znl[i]);
    276: }

    278: /* Tests A and C are from the ISMIP-HOM paper (Pattyn et al. 2008) */
    279: static void THIInitialize_HOM_A(THI thi,PetscReal x,PetscReal y,PrmNode *p)
    280: {
    281:   Units     units = thi->units;
    282:   PetscReal s     = -x*sin(thi->alpha);

    284:   p->b     = s - 1000*units->meter + 500*units->meter * sin(x*2*PETSC_PI/thi->Lx) * sin(y*2*PETSC_PI/thi->Ly);
    285:   p->h     = s - p->b;
    286:   p->beta2 = 1e30;
    287: }

    289: static void THIInitialize_HOM_C(THI thi,PetscReal x,PetscReal y,PrmNode *p)
    290: {
    291:   Units     units = thi->units;
    292:   PetscReal s     = -x*sin(thi->alpha);

    294:   p->b = s - 1000*units->meter;
    295:   p->h = s - p->b;
    296:   /* tau_b = beta2 v   is a stress (Pa) */
    297:   p->beta2 = 1000 * (1 + sin(x*2*PETSC_PI/thi->Lx)*sin(y*2*PETSC_PI/thi->Ly)) * units->Pascal * units->year / units->meter;
    298: }

    300: /* These are just toys */

    302: /* Same bed as test A, free slip everywhere except for a discontinuous jump to a circular sticky region in the middle. */
    303: static void THIInitialize_HOM_X(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
    304: {
    305:   Units     units = thi->units;
    306:   PetscReal x     = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
    307:   PetscReal r     = PetscSqrtReal(x*x + y*y),s = -x*PetscSinReal(thi->alpha);
    308:   p->b     = s - 1000*units->meter + 500*units->meter*sin(x + PETSC_PI) * sin(y + PETSC_PI);
    309:   p->h     = s - p->b;
    310:   p->beta2 = 1000 * (r < 1 ? 2 : 0) * units->Pascal * units->year / units->meter;
    311: }

    313: /* Like Z, but with 200 meter cliffs */
    314: static void THIInitialize_HOM_Y(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
    315: {
    316:   Units     units = thi->units;
    317:   PetscReal x     = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
    318:   PetscReal r     = PetscSqrtReal(x*x + y*y),s = -x*PetscSinReal(thi->alpha);

    320:   p->b = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
    321:   if (PetscRealPart(p->b) > -700*units->meter) p->b += 200*units->meter;
    322:   p->h     = s - p->b;
    323:   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;
    324: }

    326: /* Same bed as A, smoothly varying slipperiness, similar to MATLAB's "sombrero" (uncorrelated with bathymetry) */
    327: static void THIInitialize_HOM_Z(THI thi,PetscReal xx,PetscReal yy,PrmNode *p)
    328: {
    329:   Units     units = thi->units;
    330:   PetscReal x     = xx*2*PETSC_PI/thi->Lx - PETSC_PI,y = yy*2*PETSC_PI/thi->Ly - PETSC_PI; /* [-pi,pi] */
    331:   PetscReal r     = PetscSqrtReal(x*x + y*y),s = -x*PetscSinReal(thi->alpha);

    333:   p->b     = s - 1000*units->meter + 500*units->meter * sin(x + PETSC_PI) * sin(y + PETSC_PI);
    334:   p->h     = s - p->b;
    335:   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;
    336: }

    338: static void THIFriction(THI thi,PetscReal rbeta2,PetscReal gam,PetscReal *beta2,PetscReal *dbeta2)
    339: {
    340:   if (thi->friction.irefgam == 0) {
    341:     Units units = thi->units;
    342:     thi->friction.irefgam = 1./(0.5*PetscSqr(thi->friction.refvel * units->meter / units->year));
    343:     thi->friction.eps2    = 0.5*PetscSqr(thi->friction.epsvel * units->meter / units->year) * thi->friction.irefgam;
    344:   }
    345:   if (thi->friction.exponent == 0) {
    346:     *beta2  = rbeta2;
    347:     *dbeta2 = 0;
    348:   } else {
    349:     *beta2  = rbeta2 * pow(thi->friction.eps2 + gam*thi->friction.irefgam,thi->friction.exponent);
    350:     *dbeta2 = thi->friction.exponent * *beta2 / (thi->friction.eps2 + gam*thi->friction.irefgam) * thi->friction.irefgam;
    351:   }
    352: }

    354: static void THIViscosity(THI thi,PetscReal gam,PetscReal *eta,PetscReal *deta)
    355: {
    356:   PetscReal Bd2,eps,exponent;
    357:   if (thi->viscosity.Bd2 == 0) {
    358:     Units units = thi->units;
    359:     const PetscReal
    360:       n = 3.,                                           /* Glen exponent */
    361:       p = 1. + 1./n,                                    /* for Stokes */
    362:       A = 1.e-16 * pow(units->Pascal,-n) / units->year, /* softness parameter (Pa^{-n}/s) */
    363:       B = pow(A,-1./n);                                 /* hardness parameter */
    364:     thi->viscosity.Bd2      = B/2;
    365:     thi->viscosity.exponent = (p-2)/2;
    366:     thi->viscosity.eps      = 0.5*PetscSqr(1e-5 / units->year);
    367:   }
    368:   Bd2      = thi->viscosity.Bd2;
    369:   exponent = thi->viscosity.exponent;
    370:   eps      = thi->viscosity.eps;
    371:   *eta     = Bd2 * pow(eps + gam,exponent);
    372:   *deta    = exponent * (*eta) / (eps + gam);
    373: }

    375: static void RangeUpdate(PetscReal *min,PetscReal *max,PetscReal x)
    376: {
    377:   if (x < *min) *min = x;
    378:   if (x > *max) *max = x;
    379: }

    381: static void PRangeClear(PRange *p)
    382: {
    383:   p->cmin = p->min = 1e100;
    384:   p->cmax = p->max = -1e100;
    385: }

    389: static PetscErrorCode PRangeMinMax(PRange *p,PetscReal min,PetscReal max)
    390: {

    393:   p->cmin = min;
    394:   p->cmax = max;
    395:   if (min < p->min) p->min = min;
    396:   if (max > p->max) p->max = max;
    397:   return(0);
    398: }

    402: static PetscErrorCode THIDestroy(THI *thi)
    403: {

    407:   if (!*thi) return(0);
    408:   if (--((PetscObject)(*thi))->refct > 0) {*thi = 0; return(0);}
    409:   PetscFree((*thi)->units);
    410:   PetscFree((*thi)->mattype);
    411:   PetscHeaderDestroy(thi);
    412:   return(0);
    413: }

    417: static PetscErrorCode THICreate(MPI_Comm comm,THI *inthi)
    418: {
    419:   static PetscBool registered = PETSC_FALSE;
    420:   THI              thi;
    421:   Units            units;
    422:   PetscErrorCode   ierr;

    425:   *inthi = 0;
    426:   if (!registered) {
    427:     PetscClassIdRegister("Toy Hydrostatic Ice",&THI_CLASSID);
    428:     registered = PETSC_TRUE;
    429:   }
    430:   PetscHeaderCreate(thi,_p_THI,0,THI_CLASSID,"THI","Toy Hydrostatic Ice","",comm,THIDestroy,0);

    432:   PetscNew(struct _n_Units,&thi->units);
    433:   units           = thi->units;
    434:   units->meter    = 1e-2;
    435:   units->second   = 1e-7;
    436:   units->kilogram = 1e-12;

    438:   PetscOptionsBegin(comm,NULL,"Scaled units options","");
    439:   {
    440:     PetscOptionsReal("-units_meter","1 meter in scaled length units","",units->meter,&units->meter,NULL);
    441:     PetscOptionsReal("-units_second","1 second in scaled time units","",units->second,&units->second,NULL);
    442:     PetscOptionsReal("-units_kilogram","1 kilogram in scaled mass units","",units->kilogram,&units->kilogram,NULL);
    443:   }
    444:   PetscOptionsEnd();
    445:   units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second));
    446:   units->year   = 31556926. * units->second, /* seconds per year */

    448:   thi->Lx              = 10.e3;
    449:   thi->Ly              = 10.e3;
    450:   thi->Lz              = 1000;
    451:   thi->dirichlet_scale = 1;
    452:   thi->verbose         = PETSC_FALSE;

    454:   PetscOptionsBegin(comm,NULL,"Toy Hydrostatic Ice options","");
    455:   {
    456:     QuadratureType quad       = QUAD_GAUSS;
    457:     char           homexp[]   = "A";
    458:     char           mtype[256] = MATSBAIJ;
    459:     PetscReal      L,m = 1.0;
    460:     PetscBool      flg;
    461:     L    = thi->Lx;
    462:     PetscOptionsReal("-thi_L","Domain size (m)","",L,&L,&flg);
    463:     if (flg) thi->Lx = thi->Ly = L;
    464:     PetscOptionsReal("-thi_Lx","X Domain size (m)","",thi->Lx,&thi->Lx,NULL);
    465:     PetscOptionsReal("-thi_Ly","Y Domain size (m)","",thi->Ly,&thi->Ly,NULL);
    466:     PetscOptionsReal("-thi_Lz","Z Domain size (m)","",thi->Lz,&thi->Lz,NULL);
    467:     PetscOptionsString("-thi_hom","ISMIP-HOM experiment (A or C)","",homexp,homexp,sizeof(homexp),NULL);
    468:     switch (homexp[0] = toupper(homexp[0])) {
    469:     case 'A':
    470:       thi->initialize = THIInitialize_HOM_A;
    471:       thi->no_slip    = PETSC_TRUE;
    472:       thi->alpha      = 0.5;
    473:       break;
    474:     case 'C':
    475:       thi->initialize = THIInitialize_HOM_C;
    476:       thi->no_slip    = PETSC_FALSE;
    477:       thi->alpha      = 0.1;
    478:       break;
    479:     case 'X':
    480:       thi->initialize = THIInitialize_HOM_X;
    481:       thi->no_slip    = PETSC_FALSE;
    482:       thi->alpha      = 0.3;
    483:       break;
    484:     case 'Y':
    485:       thi->initialize = THIInitialize_HOM_Y;
    486:       thi->no_slip    = PETSC_FALSE;
    487:       thi->alpha      = 0.5;
    488:       break;
    489:     case 'Z':
    490:       thi->initialize = THIInitialize_HOM_Z;
    491:       thi->no_slip    = PETSC_FALSE;
    492:       thi->alpha      = 0.5;
    493:       break;
    494:     default:
    495:       SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"HOM experiment '%c' not implemented",homexp[0]);
    496:     }
    497:     PetscOptionsEnum("-thi_quadrature","Quadrature to use for 3D elements","",QuadratureTypes,(PetscEnum)quad,(PetscEnum*)&quad,NULL);
    498:     switch (quad) {
    499:     case QUAD_GAUSS:
    500:       HexQInterp = HexQInterp_Gauss;
    501:       HexQDeriv  = HexQDeriv_Gauss;
    502:       break;
    503:     case QUAD_LOBATTO:
    504:       HexQInterp = HexQInterp_Lobatto;
    505:       HexQDeriv  = HexQDeriv_Lobatto;
    506:       break;
    507:     }
    508:     PetscOptionsReal("-thi_alpha","Bed angle (degrees)","",thi->alpha,&thi->alpha,NULL);

    510:     thi->friction.refvel = 100.;
    511:     thi->friction.epsvel = 1.;

    513:     PetscOptionsReal("-thi_friction_refvel","Reference velocity for sliding","",thi->friction.refvel,&thi->friction.refvel,NULL);
    514:     PetscOptionsReal("-thi_friction_epsvel","Regularization velocity for sliding","",thi->friction.epsvel,&thi->friction.epsvel,NULL);
    515:     PetscOptionsReal("-thi_friction_m","Friction exponent, 0=Coulomb, 1=Navier","",m,&m,NULL);

    517:     thi->friction.exponent = (m-1)/2;

    519:     PetscOptionsReal("-thi_dirichlet_scale","Scale Dirichlet boundary conditions by this factor","",thi->dirichlet_scale,&thi->dirichlet_scale,NULL);
    520:     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);
    521:     PetscOptionsBool("-thi_coarse2d","Use a 2D coarse space corresponding to SSA","",thi->coarse2d,&thi->coarse2d,NULL);
    522:     PetscOptionsBool("-thi_tridiagonal","Assemble a tridiagonal system (column coupling only) on the finest level","",thi->tridiagonal,&thi->tridiagonal,NULL);
    523:     PetscOptionsList("-thi_mat_type","Matrix type","MatSetType",MatList,mtype,(char*)mtype,sizeof(mtype),NULL);
    524:     PetscStrallocpy(mtype,(char**)&thi->mattype);
    525:     PetscOptionsBool("-thi_verbose","Enable verbose output (like matrix sizes and statistics)","",thi->verbose,&thi->verbose,NULL);
    526:   }
    527:   PetscOptionsEnd();

    529:   /* dimensionalize */
    530:   thi->Lx    *= units->meter;
    531:   thi->Ly    *= units->meter;
    532:   thi->Lz    *= units->meter;
    533:   thi->alpha *= PETSC_PI / 180;

    535:   PRangeClear(&thi->eta);
    536:   PRangeClear(&thi->beta2);

    538:   {
    539:     PetscReal u       = 1000*units->meter/(3e7*units->second),
    540:               gradu   = u / (100*units->meter),eta,deta,
    541:               rho     = 910 * units->kilogram/pow(units->meter,3),
    542:               grav    = 9.81 * units->meter/PetscSqr(units->second),
    543:               driving = rho * grav * sin(thi->alpha) * 1000*units->meter;
    544:     THIViscosity(thi,0.5*gradu*gradu,&eta,&deta);
    545:     thi->rhog = rho * grav;
    546:     if (thi->verbose) {
    547:       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);
    548:       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);
    549:       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);
    550:       THIViscosity(thi,0.5*PetscSqr(1e-3*gradu),&eta,&deta);
    551:       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);
    552:     }
    553:   }

    555:   *inthi = thi;
    556:   return(0);
    557: }

    561: static PetscErrorCode THIInitializePrm(THI thi,DM da2prm,Vec prm)
    562: {
    563:   PrmNode        **p;
    564:   PetscInt       i,j,xs,xm,ys,ym,mx,my;

    568:   DMDAGetGhostCorners(da2prm,&ys,&xs,0,&ym,&xm,0);
    569:   DMDAGetInfo(da2prm,0, &my,&mx,0, 0,0,0, 0,0,0,0,0,0);
    570:   DMDAVecGetArray(da2prm,prm,&p);
    571:   for (i=xs; i<xs+xm; i++) {
    572:     for (j=ys; j<ys+ym; j++) {
    573:       PetscReal xx = thi->Lx*i/mx,yy = thi->Ly*j/my;
    574:       thi->initialize(thi,xx,yy,&p[i][j]);
    575:     }
    576:   }
    577:   DMDAVecRestoreArray(da2prm,prm,&p);
    578:   return(0);
    579: }

    583: static PetscErrorCode THISetUpDM(THI thi,DM dm)
    584: {
    585:   PetscErrorCode  ierr;
    586:   PetscInt        refinelevel,coarsenlevel,level,dim,Mx,My,Mz,mx,my,s;
    587:   DMDAStencilType st;
    588:   DM              da2prm;
    589:   Vec             X;

    592:   DMDAGetInfo(dm,&dim, &Mz,&My,&Mx, 0,&my,&mx, 0,&s,0,0,0,&st);
    593:   if (dim == 2) {
    594:     DMDAGetInfo(dm,&dim, &My,&Mx,0, &my,&mx,0, 0,&s,0,0,0,&st);
    595:   }
    596:   DMGetRefineLevel(dm,&refinelevel);
    597:   DMGetCoarsenLevel(dm,&coarsenlevel);
    598:   level = refinelevel - coarsenlevel;
    599:   DMDACreate2d(PetscObjectComm((PetscObject)thi),DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,st,My,Mx,my,mx,sizeof(PrmNode)/sizeof(PetscScalar),s,0,0,&da2prm);
    600:   DMCreateLocalVector(da2prm,&X);
    601:   {
    602:     PetscReal Lx = thi->Lx / thi->units->meter,Ly = thi->Ly / thi->units->meter,Lz = thi->Lz / thi->units->meter;
    603:     if (dim == 2) {
    604:       PetscPrintf(PetscObjectComm((PetscObject)thi),"Level %D domain size (m) %8.2g x %8.2g, num elements %3d x %3d (%8d), size (m) %g x %g\n",level,Lx,Ly,Mx,My,Mx*My,Lx/Mx,Ly/My);
    605:     } else {
    606:       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",level,Lx,Ly,Lz,Mx,My,Mz,Mx*My*Mz,Lx/Mx,Ly/My,1000./(Mz-1));
    607:     }
    608:   }
    609:   THIInitializePrm(thi,da2prm,X);
    610:   if (thi->tridiagonal) {       /* Reset coarse Jacobian evaluation */
    611:     DMDASNESSetJacobianLocal(dm,(DMDASNESJacobianLocal)THIJacobianLocal_3D_Full,thi);
    612:   }
    613:   if (thi->coarse2d) {
    614:     DMDASNESSetJacobianLocal(dm,(DMDASNESJacobianLocal)THIJacobianLocal_2D,thi);
    615:   }
    616:   PetscObjectCompose((PetscObject)dm,"DMDA2Prm",(PetscObject)da2prm);
    617:   PetscObjectCompose((PetscObject)dm,"DMDA2Prm_Vec",(PetscObject)X);
    618:   DMDestroy(&da2prm);
    619:   VecDestroy(&X);
    620:   return(0);
    621: }

    625: static PetscErrorCode DMCoarsenHook_THI(DM dmf,DM dmc,void *ctx)
    626: {
    627:   THI            thi = (THI)ctx;
    629:   PetscInt       rlevel,clevel;

    632:   THISetUpDM(thi,dmc);
    633:   DMGetRefineLevel(dmc,&rlevel);
    634:   DMGetCoarsenLevel(dmc,&clevel);
    635:   if (rlevel-clevel == 0) {DMSetMatType(dmc,MATAIJ);}
    636:   DMCoarsenHookAdd(dmc,DMCoarsenHook_THI,NULL,thi);
    637:   return(0);
    638: }

    642: static PetscErrorCode DMRefineHook_THI(DM dmc,DM dmf,void *ctx)
    643: {
    644:   THI            thi = (THI)ctx;

    648:   THISetUpDM(thi,dmf);
    649:   DMSetMatType(dmf,thi->mattype);
    650:   DMRefineHookAdd(dmf,DMRefineHook_THI,NULL,thi);
    651:   /* With grid sequencing, a formerly-refined DM will later be coarsened by PCSetUp_MG */
    652:   DMCoarsenHookAdd(dmf,DMCoarsenHook_THI,NULL,thi);
    653:   return(0);
    654: }

    658: static PetscErrorCode THIDAGetPrm(DM da,PrmNode ***prm)
    659: {
    661:   DM             da2prm;
    662:   Vec            X;

    665:   PetscObjectQuery((PetscObject)da,"DMDA2Prm",(PetscObject*)&da2prm);
    666:   if (!da2prm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm composed with given DMDA");
    667:   PetscObjectQuery((PetscObject)da,"DMDA2Prm_Vec",(PetscObject*)&X);
    668:   if (!X) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm_Vec composed with given DMDA");
    669:   DMDAVecGetArray(da2prm,X,prm);
    670:   return(0);
    671: }

    675: static PetscErrorCode THIDARestorePrm(DM da,PrmNode ***prm)
    676: {
    678:   DM             da2prm;
    679:   Vec            X;

    682:   PetscObjectQuery((PetscObject)da,"DMDA2Prm",(PetscObject*)&da2prm);
    683:   if (!da2prm) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm composed with given DMDA");
    684:   PetscObjectQuery((PetscObject)da,"DMDA2Prm_Vec",(PetscObject*)&X);
    685:   if (!X) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No DMDA2Prm_Vec composed with given DMDA");
    686:   DMDAVecRestoreArray(da2prm,X,prm);
    687:   return(0);
    688: }

    692: static PetscErrorCode THIInitial(SNES snes,Vec X,void *ctx)
    693: {
    694:   THI            thi;
    695:   PetscInt       i,j,k,xs,xm,ys,ym,zs,zm,mx,my;
    696:   PetscReal      hx,hy;
    697:   PrmNode        **prm;
    698:   Node           ***x;
    700:   DM             da;

    703:   SNESGetDM(snes,&da);
    704:   DMGetApplicationContext(da,&thi);
    705:   DMDAGetInfo(da,0, 0,&my,&mx, 0,0,0, 0,0,0,0,0,0);
    706:   DMDAGetCorners(da,&zs,&ys,&xs,&zm,&ym,&xm);
    707:   DMDAVecGetArray(da,X,&x);
    708:   THIDAGetPrm(da,&prm);
    709:   hx   = thi->Lx / mx;
    710:   hy   = thi->Ly / my;
    711:   for (i=xs; i<xs+xm; i++) {
    712:     for (j=ys; j<ys+ym; j++) {
    713:       for (k=zs; k<zs+zm; k++) {
    714:         const PetscScalar zm1      = zm-1,
    715:                           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),
    716:                           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);
    717:         x[i][j][k].u = 0. * drivingx * prm[i][j].h*(PetscScalar)k/zm1;
    718:         x[i][j][k].v = 0. * drivingy * prm[i][j].h*(PetscScalar)k/zm1;
    719:       }
    720:     }
    721:   }
    722:   DMDAVecRestoreArray(da,X,&x);
    723:   THIDARestorePrm(da,&prm);
    724:   return(0);
    725: }

    727: 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)
    728: {
    729:   PetscInt    l,ll;
    730:   PetscScalar gam;

    732:   du[0] = du[1] = du[2] = 0;
    733:   dv[0] = dv[1] = dv[2] = 0;
    734:   *u    = 0;
    735:   *v    = 0;
    736:   for (l=0; l<8; l++) {
    737:     *u += phi[l] * n[l].u;
    738:     *v += phi[l] * n[l].v;
    739:     for (ll=0; ll<3; ll++) {
    740:       du[ll] += dphi[l][ll] * n[l].u;
    741:       dv[ll] += dphi[l][ll] * n[l].v;
    742:     }
    743:   }
    744:   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]);
    745:   THIViscosity(thi,PetscRealPart(gam),eta,deta);
    746: }

    748: 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)
    749: {
    750:   PetscInt    l,ll;
    751:   PetscScalar gam;

    753:   du[0] = du[1] = 0;
    754:   dv[0] = dv[1] = 0;
    755:   *u    = 0;
    756:   *v    = 0;
    757:   for (l=0; l<4; l++) {
    758:     *u += phi[l] * n[l].u;
    759:     *v += phi[l] * n[l].v;
    760:     for (ll=0; ll<2; ll++) {
    761:       du[ll] += dphi[l][ll] * n[l].u;
    762:       dv[ll] += dphi[l][ll] * n[l].v;
    763:     }
    764:   }
    765:   gam = Sqr(du[0]) + Sqr(dv[1]) + du[0]*dv[1] + 0.25*Sqr(du[1]+dv[0]);
    766:   THIViscosity(thi,PetscRealPart(gam),eta,deta);
    767: }

    771: static PetscErrorCode THIFunctionLocal(DMDALocalInfo *info,Node ***x,Node ***f,THI thi)
    772: {
    773:   PetscInt       xs,ys,xm,ym,zm,i,j,k,q,l;
    774:   PetscReal      hx,hy,etamin,etamax,beta2min,beta2max;
    775:   PrmNode        **prm;

    779:   xs = info->zs;
    780:   ys = info->ys;
    781:   xm = info->zm;
    782:   ym = info->ym;
    783:   zm = info->xm;
    784:   hx = thi->Lx / info->mz;
    785:   hy = thi->Ly / info->my;

    787:   etamin   = 1e100;
    788:   etamax   = 0;
    789:   beta2min = 1e100;
    790:   beta2max = 0;

    792:   THIDAGetPrm(info->da,&prm);

    794:   for (i=xs; i<xs+xm; i++) {
    795:     for (j=ys; j<ys+ym; j++) {
    796:       PrmNode pn[4];
    797:       QuadExtract(prm,i,j,pn);
    798:       for (k=0; k<zm-1; k++) {
    799:         PetscInt  ls = 0;
    800:         Node      n[8],*fn[8];
    801:         PetscReal zn[8],etabase = 0;
    802:         PrmHexGetZ(pn,k,zm,zn);
    803:         HexExtract(x,i,j,k,n);
    804:         HexExtractRef(f,i,j,k,fn);
    805:         if (thi->no_slip && k == 0) {
    806:           for (l=0; l<4; l++) n[l].u = n[l].v = 0;
    807:           /* The first 4 basis functions lie on the bottom layer, so their contribution is exactly 0, hence we can skip them */
    808:           ls = 4;
    809:         }
    810:         for (q=0; q<8; q++) {
    811:           PetscReal   dz[3],phi[8],dphi[8][3],jw,eta,deta;
    812:           PetscScalar du[3],dv[3],u,v;
    813:           HexGrad(HexQDeriv[q],zn,dz);
    814:           HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
    815:           PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
    816:           jw /= thi->rhog;      /* scales residuals to be O(1) */
    817:           if (q == 0) etabase = eta;
    818:           RangeUpdate(&etamin,&etamax,eta);
    819:           for (l=ls; l<8; l++) { /* test functions */
    820:             const PetscReal ds[2] = {-sin(thi->alpha),0};
    821:             const PetscReal pp    = phi[l],*dp = dphi[l];
    822:             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];
    823:             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];
    824:           }
    825:         }
    826:         if (k == 0) { /* we are on a bottom face */
    827:           if (thi->no_slip) {
    828:             /* Note: Non-Galerkin coarse grid operators are very sensitive to the scaling of Dirichlet boundary
    829:             * conditions.  After shenanigans above, etabase contains the effective viscosity at the closest quadrature
    830:             * point to the bed.  We want the diagonal entry in the Dirichlet condition to have similar magnitude to the
    831:             * diagonal entry corresponding to the adjacent node.  The fundamental scaling of the viscous part is in
    832:             * diagu, diagv below.  This scaling is easy to recognize by considering the finite difference operator after
    833:             * scaling by element size.  The no-slip Dirichlet condition is scaled by this factor, and also in the
    834:             * assembled matrix (see the similar block in THIJacobianLocal).
    835:             *
    836:             * Note that the residual at this Dirichlet node is linear in the state at this node, but also depends
    837:             * (nonlinearly in general) on the neighboring interior nodes through the local viscosity.  This will make
    838:             * a matrix-free Jacobian have extra entries in the corresponding row.  We assemble only the diagonal part,
    839:             * so the solution will exactly satisfy the boundary condition after the first linear iteration.
    840:             */
    841:             const PetscReal   hz    = PetscRealPart(pn[0].h)/(zm-1.);
    842:             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);
    843:             fn[0]->u = thi->dirichlet_scale*diagu*x[i][j][k].u;
    844:             fn[0]->v = thi->dirichlet_scale*diagv*x[i][j][k].v;
    845:           } else {              /* Integrate over bottom face to apply boundary condition */
    846:             for (q=0; q<4; q++) {
    847:               const PetscReal jw = 0.25*hx*hy/thi->rhog,*phi = QuadQInterp[q];
    848:               PetscScalar     u  =0,v=0,rbeta2=0;
    849:               PetscReal       beta2,dbeta2;
    850:               for (l=0; l<4; l++) {
    851:                 u      += phi[l]*n[l].u;
    852:                 v      += phi[l]*n[l].v;
    853:                 rbeta2 += phi[l]*pn[l].beta2;
    854:               }
    855:               THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
    856:               RangeUpdate(&beta2min,&beta2max,beta2);
    857:               for (l=0; l<4; l++) {
    858:                 const PetscReal pp = phi[l];
    859:                 fn[ls+l]->u += pp*jw*beta2*u;
    860:                 fn[ls+l]->v += pp*jw*beta2*v;
    861:               }
    862:             }
    863:           }
    864:         }
    865:       }
    866:     }
    867:   }

    869:   THIDARestorePrm(info->da,&prm);

    871:   PRangeMinMax(&thi->eta,etamin,etamax);
    872:   PRangeMinMax(&thi->beta2,beta2min,beta2max);
    873:   return(0);
    874: }

    878: static PetscErrorCode THIMatrixStatistics(THI thi,Mat B,PetscViewer viewer)
    879: {
    881:   PetscReal      nrm;
    882:   PetscInt       m;
    883:   PetscMPIInt    rank;

    886:   MatNorm(B,NORM_FROBENIUS,&nrm);
    887:   MatGetSize(B,&m,0);
    888:   MPI_Comm_rank(PetscObjectComm((PetscObject)B),&rank);
    889:   if (!rank) {
    890:     PetscScalar val0,val2;
    891:     MatGetValue(B,0,0,&val0);
    892:     MatGetValue(B,2,2,&val2);
    893:     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);
    894:   }
    895:   return(0);
    896: }

    900: static PetscErrorCode THISurfaceStatistics(DM da,Vec X,PetscReal *min,PetscReal *max,PetscReal *mean)
    901: {
    903:   Node           ***x;
    904:   PetscInt       i,j,xs,ys,zs,xm,ym,zm,mx,my,mz;
    905:   PetscReal      umin = 1e100,umax=-1e100;
    906:   PetscScalar    usum = 0.0,gusum;

    909:   *min = *max = *mean = 0;
    910:   DMDAGetInfo(da,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
    911:   DMDAGetCorners(da,&zs,&ys,&xs,&zm,&ym,&xm);
    912:   if (zs != 0 || zm != mz) SETERRQ(PETSC_COMM_SELF,1,"Unexpected decomposition");
    913:   DMDAVecGetArray(da,X,&x);
    914:   for (i=xs; i<xs+xm; i++) {
    915:     for (j=ys; j<ys+ym; j++) {
    916:       PetscReal u = PetscRealPart(x[i][j][zm-1].u);
    917:       RangeUpdate(&umin,&umax,u);
    918:       usum += u;
    919:     }
    920:   }
    921:   DMDAVecRestoreArray(da,X,&x);
    922:   MPI_Allreduce(&umin,min,1,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)da));
    923:   MPI_Allreduce(&umax,max,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));
    924:   MPI_Allreduce(&usum,&gusum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));
    925:   *mean = PetscRealPart(gusum) / (mx*my);
    926:   return(0);
    927: }

    931: static PetscErrorCode THISolveStatistics(THI thi,SNES snes,PetscInt coarsened,const char name[])
    932: {
    933:   MPI_Comm       comm;
    934:   Vec            X;
    935:   DM             dm;

    939:   PetscObjectGetComm((PetscObject)thi,&comm);
    940:   SNESGetSolution(snes,&X);
    941:   SNESGetDM(snes,&dm);
    942:   PetscPrintf(comm,"Solution statistics after solve: %s\n",name);
    943:   {
    944:     PetscInt            its,lits;
    945:     SNESConvergedReason reason;
    946:     SNESGetIterationNumber(snes,&its);
    947:     SNESGetConvergedReason(snes,&reason);
    948:     SNESGetLinearSolveIterations(snes,&lits);
    949:     PetscPrintf(comm,"%s: Number of SNES iterations = %d, total linear iterations = %d\n",SNESConvergedReasons[reason],its,lits);
    950:   }
    951:   {
    952:     PetscReal   nrm2,tmin[3]={1e100,1e100,1e100},tmax[3]={-1e100,-1e100,-1e100},min[3],max[3];
    953:     PetscInt    i,j,m;
    954:     PetscScalar *x;
    955:     VecNorm(X,NORM_2,&nrm2);
    956:     VecGetLocalSize(X,&m);
    957:     VecGetArray(X,&x);
    958:     for (i=0; i<m; i+=2) {
    959:       PetscReal u = PetscRealPart(x[i]),v = PetscRealPart(x[i+1]),c = PetscSqrtReal(u*u+v*v);
    960:       tmin[0] = PetscMin(u,tmin[0]);
    961:       tmin[1] = PetscMin(v,tmin[1]);
    962:       tmin[2] = PetscMin(c,tmin[2]);
    963:       tmax[0] = PetscMax(u,tmax[0]);
    964:       tmax[1] = PetscMax(v,tmax[1]);
    965:       tmax[2] = PetscMax(c,tmax[2]);
    966:     }
    967:     VecRestoreArray(X,&x);
    968:     MPI_Allreduce(tmin,min,3,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)thi));
    969:     MPI_Allreduce(tmax,max,3,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)thi));
    970:     /* Dimensionalize to meters/year */
    971:     nrm2 *= thi->units->year / thi->units->meter;
    972:     for (j=0; j<3; j++) {
    973:       min[j] *= thi->units->year / thi->units->meter;
    974:       max[j] *= thi->units->year / thi->units->meter;
    975:     }
    976:     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]);
    977:     {
    978:       PetscReal umin,umax,umean;
    979:       THISurfaceStatistics(dm,X,&umin,&umax,&umean);
    980:       umin  *= thi->units->year / thi->units->meter;
    981:       umax  *= thi->units->year / thi->units->meter;
    982:       umean *= thi->units->year / thi->units->meter;
    983:       PetscPrintf(comm,"Surface statistics: u in [%12.6e, %12.6e] mean %12.6e\n",umin,umax,umean);
    984:     }
    985:     /* These values stay nondimensional */
    986:     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);
    987:     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);
    988:   }
    989:   PetscPrintf(comm,"\n");
    990:   return(0);
    991: }

    995: static PetscErrorCode THIJacobianLocal_2D(DMDALocalInfo *info,Node **x,Mat B,THI thi)
    996: {
    997:   PetscInt       xs,ys,xm,ym,i,j,q,l,ll;
    998:   PetscReal      hx,hy;
    999:   PrmNode        **prm;

    1003:   xs = info->ys;
    1004:   ys = info->xs;
    1005:   xm = info->ym;
    1006:   ym = info->xm;
    1007:   hx = thi->Lx / info->my;
    1008:   hy = thi->Ly / info->mx;

    1010:   MatZeroEntries(B);
    1011:   THIDAGetPrm(info->da,&prm);

    1013:   for (i=xs; i<xs+xm; i++) {
    1014:     for (j=ys; j<ys+ym; j++) {
    1015:       Node        n[4];
    1016:       PrmNode     pn[4];
    1017:       PetscScalar Ke[4*2][4*2];
    1018:       QuadExtract(prm,i,j,pn);
    1019:       QuadExtract(x,i,j,n);
    1020:       PetscMemzero(Ke,sizeof(Ke));
    1021:       for (q=0; q<4; q++) {
    1022:         PetscReal   phi[4],dphi[4][2],jw,eta,deta,beta2,dbeta2;
    1023:         PetscScalar u,v,du[2],dv[2],h = 0,rbeta2 = 0;
    1024:         for (l=0; l<4; l++) {
    1025:           phi[l]     = QuadQInterp[q][l];
    1026:           dphi[l][0] = QuadQDeriv[q][l][0]*2./hx;
    1027:           dphi[l][1] = QuadQDeriv[q][l][1]*2./hy;
    1028:           h         += phi[l] * pn[l].h;
    1029:           rbeta2    += phi[l] * pn[l].beta2;
    1030:         }
    1031:         jw = 0.25*hx*hy / thi->rhog; /* rhog is only scaling */
    1032:         PointwiseNonlinearity2D(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
    1033:         THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
    1034:         for (l=0; l<4; l++) {
    1035:           const PetscReal pp = phi[l],*dp = dphi[l];
    1036:           for (ll=0; ll<4; ll++) {
    1037:             const PetscReal ppl = phi[ll],*dpl = dphi[ll];
    1038:             PetscScalar     dgdu,dgdv;
    1039:             dgdu = 2.*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1]+dv[0])*dpl[1];
    1040:             dgdv = 2.*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1]+dv[0])*dpl[0];
    1041:             /* Picard part */
    1042:             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;
    1043:             Ke[l*2+0][ll*2+1] += dp[0]*jw*eta*2.*dpl[1] + dp[1]*jw*eta*dpl[0];
    1044:             Ke[l*2+1][ll*2+0] += dp[1]*jw*eta*2.*dpl[0] + dp[0]*jw*eta*dpl[1];
    1045:             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;
    1046:             /* extra Newton terms */
    1047:             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;
    1048:             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;
    1049:             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;
    1050:             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;
    1051:           }
    1052:         }
    1053:       }
    1054:       {
    1055:         const MatStencil rc[4] = {{0,i,j,0},{0,i+1,j,0},{0,i+1,j+1,0},{0,i,j+1,0}};
    1056:         MatSetValuesBlockedStencil(B,4,rc,4,rc,&Ke[0][0],ADD_VALUES);
    1057:       }
    1058:     }
    1059:   }
    1060:   THIDARestorePrm(info->da,&prm);

    1062:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    1063:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    1064:   MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
    1065:   if (thi->verbose) {THIMatrixStatistics(thi,B,PETSC_VIEWER_STDOUT_WORLD);}
    1066:   return(0);
    1067: }

    1071: static PetscErrorCode THIJacobianLocal_3D(DMDALocalInfo *info,Node ***x,Mat B,THI thi,THIAssemblyMode amode)
    1072: {
    1073:   PetscInt       xs,ys,xm,ym,zm,i,j,k,q,l,ll;
    1074:   PetscReal      hx,hy;
    1075:   PrmNode        **prm;

    1079:   xs = info->zs;
    1080:   ys = info->ys;
    1081:   xm = info->zm;
    1082:   ym = info->ym;
    1083:   zm = info->xm;
    1084:   hx = thi->Lx / info->mz;
    1085:   hy = thi->Ly / info->my;

    1087:   MatZeroEntries(B);
    1088:   THIDAGetPrm(info->da,&prm);

    1090:   for (i=xs; i<xs+xm; i++) {
    1091:     for (j=ys; j<ys+ym; j++) {
    1092:       PrmNode pn[4];
    1093:       QuadExtract(prm,i,j,pn);
    1094:       for (k=0; k<zm-1; k++) {
    1095:         Node        n[8];
    1096:         PetscReal   zn[8],etabase = 0;
    1097:         PetscScalar Ke[8*2][8*2];
    1098:         PetscInt    ls = 0;

    1100:         PrmHexGetZ(pn,k,zm,zn);
    1101:         HexExtract(x,i,j,k,n);
    1102:         PetscMemzero(Ke,sizeof(Ke));
    1103:         if (thi->no_slip && k == 0) {
    1104:           for (l=0; l<4; l++) n[l].u = n[l].v = 0;
    1105:           ls = 4;
    1106:         }
    1107:         for (q=0; q<8; q++) {
    1108:           PetscReal   dz[3],phi[8],dphi[8][3],jw,eta,deta;
    1109:           PetscScalar du[3],dv[3],u,v;
    1110:           HexGrad(HexQDeriv[q],zn,dz);
    1111:           HexComputeGeometry(q,hx,hy,dz,phi,dphi,&jw);
    1112:           PointwiseNonlinearity(thi,n,phi,dphi,&u,&v,du,dv,&eta,&deta);
    1113:           jw /= thi->rhog;      /* residuals are scaled by this factor */
    1114:           if (q == 0) etabase = eta;
    1115:           for (l=ls; l<8; l++) { /* test functions */
    1116:             const PetscReal *restrict dp = dphi[l];
    1117: #if USE_SSE2_KERNELS
    1118:             /* gcc (up to my 4.5 snapshot) is really bad at hoisting intrinsics so we do it manually */
    1119:             __m128d
    1120:               p4         = _mm_set1_pd(4),p2 = _mm_set1_pd(2),p05 = _mm_set1_pd(0.5),
    1121:               p42        = _mm_setr_pd(4,2),p24 = _mm_shuffle_pd(p42,p42,_MM_SHUFFLE2(0,1)),
    1122:               du0        = _mm_set1_pd(du[0]),du1 = _mm_set1_pd(du[1]),du2 = _mm_set1_pd(du[2]),
    1123:               dv0        = _mm_set1_pd(dv[0]),dv1 = _mm_set1_pd(dv[1]),dv2 = _mm_set1_pd(dv[2]),
    1124:               jweta      = _mm_set1_pd(jw*eta),jwdeta = _mm_set1_pd(jw*deta),
    1125:               dp0        = _mm_set1_pd(dp[0]),dp1 = _mm_set1_pd(dp[1]),dp2 = _mm_set1_pd(dp[2]),
    1126:               dp0jweta   = _mm_mul_pd(dp0,jweta),dp1jweta = _mm_mul_pd(dp1,jweta),dp2jweta = _mm_mul_pd(dp2,jweta),
    1127:               p4du0p2dv1 = _mm_add_pd(_mm_mul_pd(p4,du0),_mm_mul_pd(p2,dv1)), /* 4 du0 + 2 dv1 */
    1128:               p4dv1p2du0 = _mm_add_pd(_mm_mul_pd(p4,dv1),_mm_mul_pd(p2,du0)), /* 4 dv1 + 2 du0 */
    1129:               pdu2dv2    = _mm_unpacklo_pd(du2,dv2),                          /* [du2, dv2] */
    1130:               du1pdv0    = _mm_add_pd(du1,dv0),                               /* du1 + dv0 */
    1131:               t1         = _mm_mul_pd(dp0,p4du0p2dv1),                        /* dp0 (4 du0 + 2 dv1) */
    1132:               t2         = _mm_mul_pd(dp1,p4dv1p2du0);                        /* dp1 (4 dv1 + 2 du0) */

    1134: #endif
    1135: #if defined COMPUTE_LOWER_TRIANGULAR  /* The element matrices are always symmetric so computing the lower-triangular part is not necessary */
    1136:             for (ll=ls; ll<8; ll++) { /* trial functions */
    1137: #else
    1138:             for (ll=l; ll<8; ll++) {
    1139: #endif
    1140:               const PetscReal *restrict dpl = dphi[ll];
    1141:               if (amode == THIASSEMBLY_TRIDIAGONAL && (l-ll)%4) continue; /* these entries would not be inserted */
    1142: #if !USE_SSE2_KERNELS
    1143:               /* The analytic Jacobian in nice, easy-to-read form */
    1144:               {
    1145:                 PetscScalar dgdu,dgdv;
    1146:                 dgdu = 2.*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1]+dv[0])*dpl[1] + 0.5*du[2]*dpl[2];
    1147:                 dgdv = 2.*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1]+dv[0])*dpl[0] + 0.5*dv[2]*dpl[2];
    1148:                 /* Picard part */
    1149:                 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];
    1150:                 Ke[l*2+0][ll*2+1] += dp[0]*jw*eta*2.*dpl[1] + dp[1]*jw*eta*dpl[0];
    1151:                 Ke[l*2+1][ll*2+0] += dp[1]*jw*eta*2.*dpl[0] + dp[0]*jw*eta*dpl[1];
    1152:                 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];
    1153:                 /* extra Newton terms */
    1154:                 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];
    1155:                 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];
    1156:                 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];
    1157:                 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];
    1158:               }
    1159: #else
    1160:               /* This SSE2 code is an exact replica of above, but uses explicit packed instructions for some speed
    1161:               * benefit.  On my hardware, these intrinsics are almost twice as fast as above, reducing total assembly cost
    1162:               * by 25 to 30 percent. */
    1163:               {
    1164:                 __m128d
    1165:                   keu   = _mm_loadu_pd(&Ke[l*2+0][ll*2+0]),
    1166:                   kev   = _mm_loadu_pd(&Ke[l*2+1][ll*2+0]),
    1167:                   dpl01 = _mm_loadu_pd(&dpl[0]),dpl10 = _mm_shuffle_pd(dpl01,dpl01,_MM_SHUFFLE2(0,1)),dpl2 = _mm_set_sd(dpl[2]),
    1168:                   t0,t3,pdgduv;
    1169:                 keu = _mm_add_pd(keu,_mm_add_pd(_mm_mul_pd(_mm_mul_pd(dp0jweta,p42),dpl01),
    1170:                                                 _mm_add_pd(_mm_mul_pd(dp1jweta,dpl10),
    1171:                                                            _mm_mul_pd(dp2jweta,dpl2))));
    1172:                 kev = _mm_add_pd(kev,_mm_add_pd(_mm_mul_pd(_mm_mul_pd(dp1jweta,p24),dpl01),
    1173:                                                 _mm_add_pd(_mm_mul_pd(dp0jweta,dpl10),
    1174:                                                            _mm_mul_pd(dp2jweta,_mm_shuffle_pd(dpl2,dpl2,_MM_SHUFFLE2(0,1))))));
    1175:                 pdgduv = _mm_mul_pd(p05,_mm_add_pd(_mm_add_pd(_mm_mul_pd(p42,_mm_mul_pd(du0,dpl01)),
    1176:                                                               _mm_mul_pd(p24,_mm_mul_pd(dv1,dpl01))),
    1177:                                                    _mm_add_pd(_mm_mul_pd(du1pdv0,dpl10),
    1178:                                                               _mm_mul_pd(pdu2dv2,_mm_set1_pd(dpl[2]))))); /* [dgdu, dgdv] */
    1179:                 t0 = _mm_mul_pd(jwdeta,pdgduv);  /* jw deta [dgdu, dgdv] */
    1180:                 t3 = _mm_mul_pd(t0,du1pdv0);     /* t0 (du1 + dv0) */
    1181:                 _mm_storeu_pd(&Ke[l*2+0][ll*2+0],_mm_add_pd(keu,_mm_add_pd(_mm_mul_pd(t1,t0),
    1182:                                                                            _mm_add_pd(_mm_mul_pd(dp1,t3),
    1183:                                                                                       _mm_mul_pd(t0,_mm_mul_pd(dp2,du2))))));
    1184:                 _mm_storeu_pd(&Ke[l*2+1][ll*2+0],_mm_add_pd(kev,_mm_add_pd(_mm_mul_pd(t2,t0),
    1185:                                                                            _mm_add_pd(_mm_mul_pd(dp0,t3),
    1186:                                                                                       _mm_mul_pd(t0,_mm_mul_pd(dp2,dv2))))));
    1187:               }
    1188: #endif
    1189:             }
    1190:           }
    1191:         }
    1192:         if (k == 0) { /* on a bottom face */
    1193:           if (thi->no_slip) {
    1194:             const PetscReal   hz    = PetscRealPart(pn[0].h)/(zm-1);
    1195:             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);
    1196:             Ke[0][0] = thi->dirichlet_scale*diagu;
    1197:             Ke[1][1] = thi->dirichlet_scale*diagv;
    1198:           } else {
    1199:             for (q=0; q<4; q++) {
    1200:               const PetscReal jw = 0.25*hx*hy/thi->rhog,*phi = QuadQInterp[q];
    1201:               PetscScalar     u  =0,v=0,rbeta2=0;
    1202:               PetscReal       beta2,dbeta2;
    1203:               for (l=0; l<4; l++) {
    1204:                 u      += phi[l]*n[l].u;
    1205:                 v      += phi[l]*n[l].v;
    1206:                 rbeta2 += phi[l]*pn[l].beta2;
    1207:               }
    1208:               THIFriction(thi,PetscRealPart(rbeta2),PetscRealPart(u*u+v*v)/2,&beta2,&dbeta2);
    1209:               for (l=0; l<4; l++) {
    1210:                 const PetscReal pp = phi[l];
    1211:                 for (ll=0; ll<4; ll++) {
    1212:                   const PetscReal ppl = phi[ll];
    1213:                   Ke[l*2+0][ll*2+0] += pp*jw*beta2*ppl + pp*jw*dbeta2*u*u*ppl;
    1214:                   Ke[l*2+0][ll*2+1] +=                   pp*jw*dbeta2*u*v*ppl;
    1215:                   Ke[l*2+1][ll*2+0] +=                   pp*jw*dbeta2*v*u*ppl;
    1216:                   Ke[l*2+1][ll*2+1] += pp*jw*beta2*ppl + pp*jw*dbeta2*v*v*ppl;
    1217:                 }
    1218:               }
    1219:             }
    1220:           }
    1221:         }
    1222:         {
    1223:           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}};
    1224:           if (amode == THIASSEMBLY_TRIDIAGONAL) {
    1225:             for (l=0; l<4; l++) { /* Copy out each of the blocks, discarding horizontal coupling */
    1226:               const PetscInt   l4     = l+4;
    1227:               const MatStencil rcl[2] = {{rc[l].k,rc[l].j,rc[l].i,0},{rc[l4].k,rc[l4].j,rc[l4].i,0}};
    1228: #if defined COMPUTE_LOWER_TRIANGULAR
    1229:               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]},
    1230:                                              {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]},
    1231:                                              {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]},
    1232:                                              {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]}};
    1233: #else
    1234:               /* Same as above except for the lower-left block */
    1235:               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]},
    1236:                                              {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]},
    1237:                                              {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]},
    1238:                                              {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]}};
    1239: #endif
    1240:               MatSetValuesBlockedStencil(B,2,rcl,2,rcl,&Kel[0][0],ADD_VALUES);
    1241:             }
    1242:           } else {
    1243: #if !defined COMPUTE_LOWER_TRIANGULAR /* fill in lower-triangular part, this is really cheap compared to computing the entries */
    1244:             for (l=0; l<8; l++) {
    1245:               for (ll=l+1; ll<8; ll++) {
    1246:                 Ke[ll*2+0][l*2+0] = Ke[l*2+0][ll*2+0];
    1247:                 Ke[ll*2+1][l*2+0] = Ke[l*2+0][ll*2+1];
    1248:                 Ke[ll*2+0][l*2+1] = Ke[l*2+1][ll*2+0];
    1249:                 Ke[ll*2+1][l*2+1] = Ke[l*2+1][ll*2+1];
    1250:               }
    1251:             }
    1252: #endif
    1253:             MatSetValuesBlockedStencil(B,8,rc,8,rc,&Ke[0][0],ADD_VALUES);
    1254:           }
    1255:         }
    1256:       }
    1257:     }
    1258:   }
    1259:   THIDARestorePrm(info->da,&prm);

    1261:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    1262:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    1263:   MatSetOption(B,MAT_SYMMETRIC,PETSC_TRUE);
    1264:   if (thi->verbose) {THIMatrixStatistics(thi,B,PETSC_VIEWER_STDOUT_WORLD);}
    1265:   return(0);
    1266: }

    1270: static PetscErrorCode THIJacobianLocal_3D_Full(DMDALocalInfo *info,Node ***x,Mat A,Mat B,MatStructure *mstr,THI thi)
    1271: {

    1275:   THIJacobianLocal_3D(info,x,B,thi,THIASSEMBLY_FULL);
    1276:   *mstr = SAME_NONZERO_PATTERN;
    1277:   return(0);
    1278: }

    1282: static PetscErrorCode THIJacobianLocal_3D_Tridiagonal(DMDALocalInfo *info,Node ***x,Mat B,THI thi)
    1283: {

    1287:   THIJacobianLocal_3D(info,x,B,thi,THIASSEMBLY_TRIDIAGONAL);
    1288:   return(0);
    1289: }

    1293: static PetscErrorCode DMRefineHierarchy_THI(DM dac0,PetscInt nlevels,DM hierarchy[])
    1294: {
    1295:   PetscErrorCode  ierr;
    1296:   THI             thi;
    1297:   PetscInt        dim,M,N,m,n,s,dof;
    1298:   DM              dac,daf;
    1299:   DMDAStencilType st;
    1300:   DM_DA           *ddf,*ddc;

    1303:   PetscObjectQuery((PetscObject)dac0,"THI",(PetscObject*)&thi);
    1304:   if (!thi) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot refine this DMDA, missing composed THI instance");
    1305:   if (nlevels > 1) {
    1306:     DMRefineHierarchy(dac0,nlevels-1,hierarchy);
    1307:     dac  = hierarchy[nlevels-2];
    1308:   } else {
    1309:     dac = dac0;
    1310:   }
    1311:   DMDAGetInfo(dac,&dim, &N,&M,0, &n,&m,0, &dof,&s,0,0,0,&st);
    1312:   if (dim != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"This function can only refine 2D DMDAs");

    1314:   /* Creates a 3D DMDA with the same map-plane layout as the 2D one, with contiguous columns */
    1315:   DMDACreate3d(PetscObjectComm((PetscObject)dac),DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,st,thi->zlevels,N,M,1,n,m,dof,s,NULL,NULL,NULL,&daf);

    1317:   daf->ops->creatematrix        = dac->ops->creatematrix;
    1318:   daf->ops->createinterpolation = dac->ops->createinterpolation;
    1319:   daf->ops->getcoloring         = dac->ops->getcoloring;
    1320:   ddf                           = (DM_DA*)daf->data;
    1321:   ddc                           = (DM_DA*)dac->data;
    1322:   ddf->interptype               = ddc->interptype;

    1324:   DMDASetFieldName(daf,0,"x-velocity");
    1325:   DMDASetFieldName(daf,1,"y-velocity");

    1327:   hierarchy[nlevels-1] = daf;
    1328:   return(0);
    1329: }

    1333: static PetscErrorCode DMCreateInterpolation_DA_THI(DM dac,DM daf,Mat *A,Vec *scale)
    1334: {
    1336:   PetscInt       dim;

    1343:   DMDAGetInfo(daf,&dim,0,0,0,0,0,0,0,0,0,0,0,0);
    1344:   if (dim  == 2) {
    1345:     /* We are in the 2D problem and use normal DMDA interpolation */
    1346:     DMCreateInterpolation(dac,daf,A,scale);
    1347:   } else {
    1348:     PetscInt i,j,k,xs,ys,zs,xm,ym,zm,mx,my,mz,rstart,cstart;
    1349:     Mat      B;

    1351:     DMDAGetInfo(daf,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
    1352:     DMDAGetCorners(daf,&zs,&ys,&xs,&zm,&ym,&xm);
    1353:     if (zs != 0) SETERRQ(PETSC_COMM_SELF,1,"unexpected");
    1354:     MatCreate(PetscObjectComm((PetscObject)daf),&B);
    1355:     MatSetSizes(B,xm*ym*zm,xm*ym,mx*my*mz,mx*my);

    1357:     MatSetType(B,MATAIJ);
    1358:     MatSeqAIJSetPreallocation(B,1,NULL);
    1359:     MatMPIAIJSetPreallocation(B,1,NULL,0,NULL);
    1360:     MatGetOwnershipRange(B,&rstart,NULL);
    1361:     MatGetOwnershipRangeColumn(B,&cstart,NULL);
    1362:     for (i=xs; i<xs+xm; i++) {
    1363:       for (j=ys; j<ys+ym; j++) {
    1364:         for (k=zs; k<zs+zm; k++) {
    1365:           PetscInt    i2  = i*ym+j,i3 = i2*zm+k;
    1366:           PetscScalar val = ((k == 0 || k == mz-1) ? 0.5 : 1.) / (mz-1.); /* Integration using trapezoid rule */
    1367:           MatSetValue(B,cstart+i3,rstart+i2,val,INSERT_VALUES);
    1368:         }
    1369:       }
    1370:     }
    1371:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
    1372:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
    1373:     MatCreateMAIJ(B,sizeof(Node)/sizeof(PetscScalar),A);
    1374:     MatDestroy(&B);
    1375:   }
    1376:   return(0);
    1377: }

    1381: static PetscErrorCode DMCreateMatrix_THI_Tridiagonal(DM da,MatType mtype,Mat *J)
    1382: {
    1383:   PetscErrorCode         ierr;
    1384:   Mat                    A;
    1385:   PetscInt               xm,ym,zm,dim,dof = 2,starts[3],dims[3];
    1386:   ISLocalToGlobalMapping ltog,ltogb;

    1389:   DMDAGetInfo(da,&dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
    1390:   if (dim != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected DMDA to be 3D");
    1391:   DMDAGetCorners(da,0,0,0,&zm,&ym,&xm);
    1392:   DMGetLocalToGlobalMapping(da,&ltog);
    1393:   DMGetLocalToGlobalMappingBlock(da,&ltogb);
    1394:   MatCreate(PetscObjectComm((PetscObject)da),&A);
    1395:   MatSetSizes(A,dof*xm*ym*zm,dof*xm*ym*zm,PETSC_DETERMINE,PETSC_DETERMINE);
    1396:   MatSetType(A,mtype);
    1397:   MatSetFromOptions(A);
    1398:   MatSeqAIJSetPreallocation(A,3*2,NULL);
    1399:   MatMPIAIJSetPreallocation(A,3*2,NULL,0,NULL);
    1400:   MatSeqBAIJSetPreallocation(A,2,3,NULL);
    1401:   MatMPIBAIJSetPreallocation(A,2,3,NULL,0,NULL);
    1402:   MatSeqSBAIJSetPreallocation(A,2,2,NULL);
    1403:   MatMPISBAIJSetPreallocation(A,2,2,NULL,0,NULL);
    1404:   MatSetLocalToGlobalMapping(A,ltog,ltog);
    1405:   MatSetLocalToGlobalMappingBlock(A,ltogb,ltogb);
    1406:   DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
    1407:   MatSetStencil(A,dim,dims,starts,dof);
    1408:   *J   = A;
    1409:   return(0);
    1410: }

    1414: static PetscErrorCode THIDAVecView_VTK_XML(THI thi,DM da,Vec X,const char filename[])
    1415: {
    1416:   const PetscInt dof   = 2;
    1417:   Units          units = thi->units;
    1418:   MPI_Comm       comm;
    1420:   PetscViewer    viewer;
    1421:   PetscMPIInt    rank,size,tag,nn,nmax;
    1422:   PetscInt       mx,my,mz,r,range[6];
    1423:   PetscScalar    *x;

    1426:   PetscObjectGetComm((PetscObject)thi,&comm);
    1427:   DMDAGetInfo(da,0, &mz,&my,&mx, 0,0,0, 0,0,0,0,0,0);
    1428:   MPI_Comm_size(comm,&size);
    1429:   MPI_Comm_rank(comm,&rank);
    1430:   PetscViewerASCIIOpen(comm,filename,&viewer);
    1431:   PetscViewerASCIIPrintf(viewer,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
    1432:   PetscViewerASCIIPrintf(viewer,"  <StructuredGrid WholeExtent=\"%d %d %d %d %d %d\">\n",0,mz-1,0,my-1,0,mx-1);

    1434:   DMDAGetCorners(da,range,range+1,range+2,range+3,range+4,range+5);
    1435:   PetscMPIIntCast(range[3]*range[4]*range[5]*dof,&nn);
    1436:   MPI_Reduce(&nn,&nmax,1,MPI_INT,MPI_MAX,0,comm);
    1437:   tag  = ((PetscObject) viewer)->tag;
    1438:   VecGetArray(X,&x);
    1439:   if (!rank) {
    1440:     PetscScalar *array;
    1441:     PetscMalloc(nmax*sizeof(PetscScalar),&array);
    1442:     for (r=0; r<size; r++) {
    1443:       PetscInt    i,j,k,xs,xm,ys,ym,zs,zm;
    1444:       PetscScalar *ptr;
    1445:       MPI_Status  status;
    1446:       if (r) {
    1447:         MPI_Recv(range,6,MPIU_INT,r,tag,comm,MPI_STATUS_IGNORE);
    1448:       }
    1449:       zs = range[0];ys = range[1];xs = range[2];zm = range[3];ym = range[4];xm = range[5];
    1450:       if (xm*ym*zm*dof > nmax) SETERRQ(PETSC_COMM_SELF,1,"should not happen");
    1451:       if (r) {
    1452:         MPI_Recv(array,nmax,MPIU_SCALAR,r,tag,comm,&status);
    1453:         MPI_Get_count(&status,MPIU_SCALAR,&nn);
    1454:         if (nn != xm*ym*zm*dof) SETERRQ(PETSC_COMM_SELF,1,"should not happen");
    1455:         ptr = array;
    1456:       } else ptr = x;
    1457:       PetscViewerASCIIPrintf(viewer,"    <Piece Extent=\"%d %d %d %d %d %d\">\n",zs,zs+zm-1,ys,ys+ym-1,xs,xs+xm-1);

    1459:       PetscViewerASCIIPrintf(viewer,"      <Points>\n");
    1460:       PetscViewerASCIIPrintf(viewer,"        <DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n");
    1461:       for (i=xs; i<xs+xm; i++) {
    1462:         for (j=ys; j<ys+ym; j++) {
    1463:           for (k=zs; k<zs+zm; k++) {
    1464:             PrmNode   p;
    1465:             PetscReal xx = thi->Lx*i/mx,yy = thi->Ly*j/my,zz;
    1466:             thi->initialize(thi,xx,yy,&p);
    1467:             zz   = PetscRealPart(p.b) + PetscRealPart(p.h)*k/(mz-1);
    1468:             PetscViewerASCIIPrintf(viewer,"%f %f %f\n",xx,yy,zz);
    1469:           }
    1470:         }
    1471:       }
    1472:       PetscViewerASCIIPrintf(viewer,"        </DataArray>\n");
    1473:       PetscViewerASCIIPrintf(viewer,"      </Points>\n");

    1475:       PetscViewerASCIIPrintf(viewer,"      <PointData>\n");
    1476:       PetscViewerASCIIPrintf(viewer,"        <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n");
    1477:       for (i=0; i<nn; i+=dof) {
    1478:         PetscViewerASCIIPrintf(viewer,"%f %f %f\n",PetscRealPart(ptr[i])*units->year/units->meter,PetscRealPart(ptr[i+1])*units->year/units->meter,0.0);
    1479:       }
    1480:       PetscViewerASCIIPrintf(viewer,"        </DataArray>\n");

    1482:       PetscViewerASCIIPrintf(viewer,"        <DataArray type=\"Int32\" Name=\"rank\" NumberOfComponents=\"1\" format=\"ascii\">\n");
    1483:       for (i=0; i<nn; i+=dof) {
    1484:         PetscViewerASCIIPrintf(viewer,"%d\n",r);
    1485:       }
    1486:       PetscViewerASCIIPrintf(viewer,"        </DataArray>\n");
    1487:       PetscViewerASCIIPrintf(viewer,"      </PointData>\n");

    1489:       PetscViewerASCIIPrintf(viewer,"    </Piece>\n");
    1490:     }
    1491:     PetscFree(array);
    1492:   } else {
    1493:     MPI_Send(range,6,MPIU_INT,0,tag,comm);
    1494:     MPI_Send(x,nn,MPIU_SCALAR,0,tag,comm);
    1495:   }
    1496:   VecRestoreArray(X,&x);
    1497:   PetscViewerASCIIPrintf(viewer,"  </StructuredGrid>\n");
    1498:   PetscViewerASCIIPrintf(viewer,"</VTKFile>\n");
    1499:   PetscViewerDestroy(&viewer);
    1500:   return(0);
    1501: }

    1505: int main(int argc,char *argv[])
    1506: {
    1507:   MPI_Comm       comm;
    1508:   THI            thi;
    1510:   DM             da;
    1511:   SNES           snes;

    1513:   PetscInitialize(&argc,&argv,0,help);
    1514:   comm = PETSC_COMM_WORLD;

    1516:   THICreate(comm,&thi);
    1517:   {
    1518:     PetscInt M = 3,N = 3,P = 2;
    1519:     PetscOptionsBegin(comm,NULL,"Grid resolution options","");
    1520:     {
    1521:       PetscOptionsInt("-M","Number of elements in x-direction on coarse level","",M,&M,NULL);
    1522:       N    = M;
    1523:       PetscOptionsInt("-N","Number of elements in y-direction on coarse level (if different from M)","",N,&N,NULL);
    1524:       if (thi->coarse2d) {
    1525:         PetscOptionsInt("-zlevels","Number of elements in z-direction on fine level","",thi->zlevels,&thi->zlevels,NULL);
    1526:       } else {
    1527:         PetscOptionsInt("-P","Number of elements in z-direction on coarse level","",P,&P,NULL);
    1528:       }
    1529:     }
    1530:     PetscOptionsEnd();
    1531:     if (thi->coarse2d) {
    1532:       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);

    1534:       da->ops->refinehierarchy     = DMRefineHierarchy_THI;
    1535:       da->ops->createinterpolation = DMCreateInterpolation_DA_THI;

    1537:       PetscObjectCompose((PetscObject)da,"THI",(PetscObject)thi);
    1538:     } else {
    1539:       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);
    1540:     }
    1541:     DMDASetFieldName(da,0,"x-velocity");
    1542:     DMDASetFieldName(da,1,"y-velocity");
    1543:   }
    1544:   THISetUpDM(thi,da);
    1545:   if (thi->tridiagonal) da->ops->creatematrix = DMCreateMatrix_THI_Tridiagonal;

    1547:   {                             /* Set the fine level matrix type if -da_refine */
    1548:     PetscInt rlevel,clevel;
    1549:     DMGetRefineLevel(da,&rlevel);
    1550:     DMGetCoarsenLevel(da,&clevel);
    1551:     if (rlevel - clevel > 0) {DMSetMatType(da,thi->mattype);}
    1552:   }

    1554:   DMDASNESSetFunctionLocal(da,ADD_VALUES,(DMDASNESFunctionLocal)THIFunctionLocal,thi);
    1555:   if (thi->tridiagonal) {
    1556:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobianLocal)THIJacobianLocal_3D_Tridiagonal,thi);
    1557:   } else {
    1558:     DMDASNESSetJacobianLocal(da,(DMDASNESJacobianLocal)THIJacobianLocal_3D_Full,thi);
    1559:   }
    1560:   DMCoarsenHookAdd(da,DMCoarsenHook_THI,NULL,thi);
    1561:   DMRefineHookAdd(da,DMRefineHook_THI,NULL,thi);

    1563:   DMSetApplicationContext(da,thi);

    1565:   SNESCreate(comm,&snes);
    1566:   SNESSetDM(snes,da);
    1567:   DMDestroy(&da);
    1568:   SNESSetComputeInitialGuess(snes,THIInitial,NULL);
    1569:   SNESSetFromOptions(snes);

    1571:   SNESSolve(snes,NULL,NULL);

    1573:   THISolveStatistics(thi,snes,0,"Full");

    1575:   {
    1576:     PetscBool flg;
    1577:     char      filename[PETSC_MAX_PATH_LEN] = "";
    1578:     PetscOptionsGetString(NULL,"-o",filename,sizeof(filename),&flg);
    1579:     if (flg) {
    1580:       Vec X;
    1581:       DM  dm;
    1582:       SNESGetSolution(snes,&X);
    1583:       SNESGetDM(snes,&dm);
    1584:       THIDAVecView_VTK_XML(thi,dm,X,filename);
    1585:     }
    1586:   }

    1588:   DMDestroy(&da);
    1589:   SNESDestroy(&snes);
    1590:   THIDestroy(&thi);
    1591:   PetscFinalize();
    1592:   return 0;
    1593: }