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. (5) compute the residual x^+ - dt/2Ny^{k,i} - y^{k,i}, (6) solve the linear system for dy^{k,i+1}, (7) update y^{k,i+1} = y^{k,i} + dy^{k,i+1}. (8) Increment i. If i<i_max, go to 5. (8) increment k. If k<k_max, go to 3.

Code implementation

A timestepper class which has a run() method,

def run(self):
   t = self.t0
   dt = self.dt
   tmax = self.tmax

   self.xn = self.x_init

   while(t<tmax - 0.5*dt):
       t += dt 
       self.apply_nonlinear(self.xn,self.xstar)
       self.xnp1.assign(self.xn)

       for(k in range(self.maxk)):
           self.set_ubar()  #computes self.ubar from self.xn and self.xnp1
           self.rho_advection.step(self.xstar,self.xp) #advects rho from xstar and puts result in rho component of xp
           self.u_advection.step(self.xstar,self.xp) #same thing for u
           self.theta_advection.step(self.xstar,self.xp) #same thing for theta
           for(i in range(self.maxi)):
               self.xrhs.assign(0.) #xrhs is the residual which goes in the linear solve
               self.apply_nonlinear(self.xp,self.xrhs)
               self.xrhs -= self.xnp1
               self.linear_system.solve() # solves linear system and places result in self.dy
               self.xnp1 += self.dy

It's useful to have the advection methods as class objects as then we can test them independently of the rest of the solver. So, a usercode execution script needs to instantiate these class objects and then put them into the timestepper class object.

Q: is it better to do this with an init method or just to put them in the object from the execution script?