BISMG:StephenC

Friday 20 Nov
I'm still getting to grips with Chombo. I've written a very simple application, approximately solving a first order equation

$$\frac{\partial u}{\partial t} = \sin(10(x+y) - t)$$

with initial values

$$u(x,y,t=0) = \cos(10(x+y))$$

To implement this in Chombo, I wrote a C++ class, AMRLevelSinWave which is a subclass of Chombo's AMRLevel base class. Any subclass of AMRLevel can be passed to the main driver for time dependent calculations, and needs to implement a few methods: advance, initialGrid(g), regrid(g), initialTagCells(a), tagCells(a).

Of these, it is advance that would interface with Glimmer-CISM. In this app, it updates the solution at each level like so

DataIterator dit(_uNew.disjointBoxLayout); for (dit.begin; dit.ok; ++dit){ const Box& b = _uNew.disjointBoxLayout[dit]; BoxIterator bit(b); for (bit.begin; bit.ok; ++bit){ const IntVect& iv = bit; Real x = (iv[0] + 0.5) * dx; Real y = (iv[1] + 0.5) * dx; _uNew[dit](iv,0) += dt * 0.5*(sin( 10.0 * (x+y) - time) + sin( 10.0 * (x+y) - time - dt)); } }

So, iterate through each of a level's set of rectangular regions, using the DataIterator dit. Then fish out _uNew[dit], which contains the data for just one rectangular box (conveniently including a block of data arranged as a FORTRAN multidimensional array would be), and updates it. It is easy enough to imagine passing this block (plus some boundary data, which could be derived from the solution on a coarser level) to some fixed mesh code to update (although something extra needs to be done to match both solution and fluxes at the boundaries)

By implementing the other methods, a series of refined meshes are generated as the solution is evolved, three timesteps worth shown below







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