Timestepping

Encapsulating the timestepping of the model

A timestepper class which has a run() method,

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

   self.xn.assign(self.x_init)

   xstar_fields = self.xstar.split()
   xp_fields = self.xp.split()

   while(t<tmax - 0.5*dt):
       t += dt 
       self.apply_forcing(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
           for advection, index in self.advection_list:
               advection.apply(xstar_fields[index],xp_fields[index]) #advects a field from xstar and puts result in xp
           for(i in range(self.maxi)):
               self.xrhs.assign(0.) #xrhs is the residual which goes in the linear solve
               self.apply_forcing(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.

The __init__ method of the timestepper has the following interface.

def __init__(self, state, advection_list):

where state is a model state class object, and advection_list is a list of tuples, the first being a class with an apply method, and the second being an index specifying which component of the mixed function space the field should be taken from.

Lawrence also wrote:

Rather than having a monolithic timestepper that encodes the scheme, we should think a little about how to encapsulate the set of operations required. The timestepping library the James Maddison wrote for dolfin-adjoint (http://www.sciencedirect.com/science/article/pii/S0045782514000966) may provide inspiration. Having loooked at the code in the past, I do not think we can lift it, however we may wish to consider some of the ideas.

The proposed timestepping run method is not extensible, or modifiable, without editing the main code base. This is a recipe for problems. James' solution was to have a timestepping object to which one adds a sequence of solves. He then picks these apart to find the dependencies and the updates (the assign statements above). Perhaps a halfway house if to have a sketch like this:

stepper = TimeStepper()

stepper.add_cycle('time', tmin, tmax)

t, cycle = stepper.get_cycle('time')

cycle.add_expression(lambda: apply_nonlinear(...))

cycle.add_cycle('k', 0, maxk)

k, kcycle = cycle.get_cycle('k')

kcycle.add_expression(...)

I have not thought hard about this