BISMG:StephenC

Friday 13th November
I've been investigating Chombo, attempting to decide how we might use it to build an ice sheet model with adaptive mesh refinement. Ideally, we would like to integrate it with the GLIMMER-CISM model, but at the moment I don't know enough about either to say how that might be done. I think that the best way to start would be to look at how a Chombo application needs to be written by trying to do just that: implementing a simple ice sheet model, essentially a 2D analogue of the model Rupert used in his 1D grounding line work.

The (nonlinear) PDEs are (1), a vertically integrated mass balance equation $$\frac{\partial H}{\partial t} + (Hu_j)_{,j} = a$$ and (2) a vertically integrated momentum balance equation $$\left [ \nu \left (u_{i,j} + u_{j,i} + 2 \delta _{ij} u_{k,k} \right ) \right]_{,j} = \beta ^2 u_i + \rho g H s_{,i}$$. $$H$$ is the ice height, and $$u$$ the vertically averaged velocity vector. $$ \nu = b \left [ (u_{i,j} + u_{j,i})(u_{j,i} + u_{i,j}) \right ]^\frac{n-1}{2n} $$ is a viscosity given by Glen's flow law for ice, $$\beta^2(x,y)$$ is a basal drag coefficient, $$s(x,y)$$ is the surface elevation of ice.

(1) is a first-order hyperbolic conservation equation and (2) is second-order elliptic. On top of this there will be second-order parabolic equations (heat equation), which can be broken into hyperbolic and elliptic parts, so I think that building a model that solves these should demonstrate that Chombo can do the basic job, and help me decide what kind of bridge between the two code bases might be needed. I expect that the most effort will be due to elliptic operators, which need to be inverted on the whole mesh. Chombo provides a fairly involved multigrid procedure (AMRVCycleMG(l)) which does the job, and only seems to need definitions of the operators on uniform grids, which looks promising. There is a C++ class, ViscousTensorOp which can be used to solve a Stokes problem, with $$A u_i + B \left [ C \left ( u_{i,j} + u_{j,i} \right ) + D \delta _{ij} u_{k,k} \right ] _{,j} = E $$ That is, a linear, constant coefficient elliptic PDE that is, basically, (2) linearised. Chombo also has an example application (SelfGravity) where a system of hyperbolic equations are to be advanced in time alongside solutions to one elliptic equation (Poisson's equation). So, it seems that I have a template to work from.

If anyone would like see the sort of thing Chombo can do, let me know. There is a nice advection-diffusion (heat equation) example included but it has a minor flaw which renders the results insensible - it uses $$|\phi _{,j} / \phi| > \epsilon$$ to decide where to refine the mesh. Since $$\phi$$ is close to zero across much of the domain, the mesh is refined everywhere. In my copy, I have changed the criterion to $$|\phi _{,j} / (1 + \phi)| > \epsilon$$, and the resulting adaptive refinement works nicely