Time Stepping

An overview of the time step algorithms used in the dedicated solvers implemented in DORiE.

Model Run Algorithm

All solvers in DORiE are derived from the Model base class which implements the way a model “runs” a simulation. The run() method works as a simple loop over the step() method described below for each model separately.

All model run() algorithms proceed as follows:

  1. Write the output if enabled via the model’s output.policy.

  2. Check if the model has reached the simulation end time. If yes, return. Otherwise, continue.

  3. Perform a model time step.

  4. Adapt the grid if enabled via the adaptivity.policy.

  5. Repeat from 2.

Model Step Methods

Every model implements its own time step algorithm variant.

Richards Solver

Because the Richards equation is non-linear, a Newton solver is used to compute a time step. Since the equation is also very stiff, we only use an implicit time step scheme. The Newton solver needs to iteratively linearise the equation system in order to solve it. Depending on the “local” non-linearity of the problem, this operation might fail. However, for implicit time step schemes there is no a-priori estimate on a maximum time step.

We tackle these issues with a heuristic approach: If the equation system is highly non-linear, the linearisation in the Newton solver is ineffective and it requires more iteration steps. High non-linearities also hint at strongly varying soil water dynamics during that time step. In this case, we want the time steps to be small. If the dynamics are nearly constant, we want the time steps to be large and expect very few Newton solver iterations because the problem becomes nearly linear.

Our time step size heurisic is therefore as follows: Depending on the current time step size, we compute a maximum number of Newton iterations allowed. If the solver does not converge within this range, it is aborted and the time step is reduced. Whenever computing a solution succeeds, the time step is increased. The number of Newton solver iterations allowed is given by

\[N_\text{Newton} = \left\lfloor (N_\text{min} - N_\text{max}) \frac{\lg (T - T_\text{min} + 1)} {\lg (T_\text{max} - T_\text{min} + 1)} + N_\text{max} \right\rceil ,\]

where \(N_\text{min}\) and \(N_\text{max}\) are the allowed iterations at the minimum and maximum time step, respectively, \(T_\text{min}\) and \(T_\text{max}\) are the minimum and maximum time steps, respectively, and \(T\) is the current time step. These values may be adjusted in the Configuration File. Note that the result \(N_\text{Newton}\) is rounded.

The Richards solver step() algorithm proceeds as follows:

  1. Limit the current time step according to any boundary condition change.

  2. Compute the number of allowed Newton solver iterations with the equation above.

  3. Apply the non-linear solver.

    If this fails, reduce the time step by the user-supplied factor and repeat from 2.

  4. Increment the simulation time and increase the time step by the user-supplied factor.

  5. Write the current data depending on the output policy.

Transport Solver

The Transport solver usually does not work on its own because it requires information on soil water content and water flux which need to be computed by the Richards solver. However, it is implemented as an independent solver with its own time step algorithm.

The main difference to the Richards solver is that the equation is linear, which means that applying the solver only fails due to numerical errors. Therefore, we have devised an a-posteriori solution check which can be controlled by the user. This check may reject solutions and force a reduction of the time step. Additionally, the solver supports explicit time steps. In this case, the maximum time step is given by the CFL condition which can be computed based on local velocities before applying the solver.

The Transport solver step() algorithm proceeds as follows:

  1. If the time step scheme is explicit, compute the maximum time step with the CFL criterion and the user-supplied Courant number.

  2. Limit the time step according to any boundary condition change.

  3. Apply the linear solver.

    If this fails, reduce the time step by the user-supplied factor and repeat.

  4. Perform a check of the computed solution.

    Depending on the user settings, this may reject the solution, in which case the time step is reduced and step 3 is repeated.

  5. Increment the simulation time and increase the time step by the user-supplied factor.

  6. Depending on the output policy, write the current data.

Coupled Solver

If the Transport solver is enabled, the software actually runs a coupled solver which combines both Richards and Transport solver and manages their interaction. Whenever it calls one its subordinate solvers, they perform their respective time step algorithm as depicted above.

The coupled solver step() algorithm proceeds as follows:

  1. Retrieve water content and water flux states from the Richards solver at the start of the time step.

  2. Compute a time step with the Richards solver, calling its algorithm.

  3. Retrieve water content and water flux states from the Richards solver at the end of the time step.

  4. Pass water content and water flux states to the Transport solver.

  5. Suggest the difference between the current Transport solver time and the Richards solver time as maxiumum time step in the Transport solver.

  6. Compute a time step with the Transport solver, calling its algorithm.

    Repeat from 5 until both solvers are synchronized.

Note

Depending on the grid adaptivity policy, the task of marking and adapting the grid may be forwarded to one of the integrated models while both model solutions are projected onto the new grid.