Home

Welcome to the dcore wiki!

Here is a sketch of the proposed dcore code layout.

Mathematical Description

a cartoon of the compressible equations is

Dx/Dt + Nx = 0,

where x is the model prognostic variables (u,rho,theta), Dx/Dt = x_t + u.grad x, and N is a nonlinear operator with Nx = N(x_0) + L(x-x_0) + o(x^2), where x_0 is a reference profile with u=0, and given rho, theta. We apply further approximations to L by neglecting horizontal derivatives

Then, the timestepping scheme is an iterative method for solving the implicit system

x^* = x^n - (1-\alpha) dt Nx^n, x^+ = \Phi_{\bar{u},dt}x^*, x^{n+1} = x^+ - \alpha dt Nx^{n+1},

where \Phi_{\bar{u},dt}x^* is the numerical solution of the pure advection equation

Dx/Dt = 0,

solved over one timestep, with constant advection velocity \bar{u} = u^n + \alpha(u^{n+1}-u^n), with initial condition x^*.

To solve this iteratively, we calculate iterative approximations y^k to x^{n+1}, with y^k = x^n, using the incremental formulation

A dy^{k+1} = x^+ - dt/2Ny^k - y^k, y^{k+1} = y^k + dy^{k+1},

where x^+ is obtained by

x^+ = \Phi_{\bar{u}}x^*,

where \bar{u} = u^n + \alpha(v^k - u^n), and v^k is the velocity from y^k.

The operator A is I + alphadt(L - B),

where B is the linearisation of \Phi_{\bar{u},dt}x^* about the reference profile, i.e. \Phi_{\bar{u},x^} \approx x^ + B(x^n + \alphadt(y^k - x^n)).

Hence, the steps are:

1. Compute x^*.
2. set k=0, y^k = x^n.
3. compute \bar{u} from y^k.
4. calculate x^+ using \bar{u} by calling advection schemes independently for (u,\rho,\theta) [or possibly sequentially for conservative theta advection],
5. compute the residual x^+ - dt/2Ny^k - y^k,
6. solve the linear system for dy^{k+1},
7. update y^{k+1} = y^k + dy^{k+1}.
8. increment k and return to (3) if the max number of iterations has not been reached.

If the advection schemes are expensive to evaluate compared to evaluating N, then there is some benefit to returning to (5) instead of (3). Hence, we get the modified scheme:

Hence, the steps are:

1. Compute x^*.
2. set "outer iteration" counter k=0, y^0 = x^n.
3. compute \bar{u} from y^k.
4. calculate x^+ using \bar{u} by calling advection schemes independently for (u,\rho,\theta) [or possibly sequentially for conservative theta advection],
5. set "inner iteration" counter i=0, y^{k,0} = y^k.
6. compute the residual x^+ - dt/2Ny^{k,i} - y^{k,i},
7. solve the linear system for dy^{k,i+1},
8. update y^{k,i+1} = y^{k,i} + dy^{k,i+1}.
9. Increment i. If i<i_max, go to 5.
10. increment k. If k<k_max, go to 3.

Code implementation

Timestepper class

Timestepping is described on a separate page

Model state class

Given the mesh, and some function space parameters (horizontal/vertical degree) through state.__init__(), this class builds the function spaces V2,V3,Vt,... followed by the mixed function space M = V2*V3*Vt*... for the state variable. It also holds any other parameters such as dt etc.