06 - Time integration

ReMKiT1D offers full control over time integration schemes used with models, as well as a wide range of detail for users of varying levels.

In this tutorial we cover:

1. Integrators, and in particular the implemented default Backwards Euler integrator options

2. IntegrationRules and IntegrationSteps

3. Timestep control

4. IntegrationSchemes and inline composition of integration steps

from RMK_support import Variable,Grid,Model,DiagonalStencil,IntegrationRule,IntegrationStep,IntegrationScheme,BDEIntegrator,ModelCollection, Timestep

import numpy as np

Let’s first build some Variables/Models (see previous tutorial for details)

grid = Grid(xGrid = 0.1 * np.ones(16),
            interpretXGridAsWidths = True,
            vGrid = 0.1 * np.ones(8),
            interpretVGridAsWidths = True,
            lMax = 3,
            )

a,b,c = (Variable(name,grid,data=(i+1)*np.ones(16)) for i,name in enumerate(["a","b","c"]))

model1 = Model("model1")

diag = DiagonalStencil()

model1.ddt[a] += 2*diag(b).rename("term1") - 0.1*diag(b*c).rename("term2")
model1.ddt[b] += 2*diag(c).rename("term3") - diag(a*b).rename("term4")

model2 = model1.rename("model2")
for term in model2.ddt[a].terms:
    term = term.regroup([2])

collection = ModelCollection(model1,model2)

In ReMKiT1D, a timestepping scheme consists of IntegrationSteps, each of which is associated with an Integrator, a set of Rules for each Model it evolves, as well as the fraction of the global Timestep.

While there are shortcuts one can use when building the most basic of integration schemes, here we go over all of the components and features, starting with the built-in Backwards Euler integrator with fixed-point iterations. For other implemented integrators see examples and documentation.

For the basics of how the integrator works see documentation and Appendix B in the code paper.

The integrator uses PETSc, and will, given a timestep length, attempt to take the step. If the solve fails, which can be because of linear solver failure (associated with a KSP failure code in PETSc), or because the maximum number of allowed fixed-point iterations is read, the solver can be instructed to retry the solve with half the requested step length. After a number of requested length steps the solver will try to consolidate back to a single substep. The properties of this internal step control can be set at solver construction - see BDEIntegrator docstring.

integrator1 = BDEIntegrator("BDE1", # Name of the integrator
                           convergenceVars = [a, b], # Variables used to test convergence
                           nonlinTol = 1e-12, # Relative convergence tolerance for fixed-point iterations
                           maxNonlinIters = 100, # Allowed number of failed fixed-point iterations before the solver declares non-linear convergence failure
                           associatedPETScGroup = 1, # If using multiple different BDEIntegrators to solve different model/term contributions they should be associated with different PETSc objects - Fortran 1-indexing
                           internalStepControl = True, # Allow the solver to substep in case of failed solves (see other doctrings for control customisation options)
                           relaxationWeight = 0.75 # Fixed-point relaxation weight <1 - under-relaxation, >1 over-relaxation
                        )

integrator2 = BDEIntegrator("BDE2", associatedPETScGroup = 2) # We will use this for demonstrating operator splitting

Integration steps and rules

Integrators are associated with IntegrationSteps, which can be chained (see last sectionbelow).

Let’s create two steps with our integrators

step1 = IntegrationStep("step1",
                        integrator1)

step2 = IntegrationStep("step2",
                        integrator2)

# By default each step in a sequence starts from the output state of the previous step and evolves the time variable - this can be changed in the constructor, or by using the following

step1 = step1.disableTimeEvo()
step1 = step1.enableTimeEvo()

step1 = step1.startFromZero() # Starts from the initial state of the system (at the start of the step sequence)
step1 = step1.startFromLast() # Starts from the output of the previous step in the sequence (this is the default)

Each step has a fraction of the global timestep (the timestep associated with the full step sequence) associated with it. This is by default 1.0, and can be automatically set when combining steps (see further down below)

print(step1.stepFraction)

1.0

In order to properly make use of IntegrationSteps they should be supplied with IntegrationRules for each Model the step evolves.

The default rule is to update and evaluate all term groups in the Model and to update any ModelboundData. If the default behaviour is acceptable, Models and ModelCollections can be added directly to the IntegrationStep, otherwise IntegrationRules should be specified.

step1.add(model1) # Will add a default rule for model1

customRule = IntegrationRule(model2,
                             updatedGroups=[1], # Updates only group 1, so that group 2 is only updated at the start of the step sequence
                             evaluatedGroups=[1,2], # Evaluate both term groups
                             updateModelData=True # Update modelbound data if present
                             )

step2.add(customRule)

print(step1.rules[model1].modelName) # rules can be directly indexed by evolved model

print(step2.rules.evolvedModels)
print(step2.rules[model2].evaluatedGroups)

model1
['model2']
[1, 2]

Timestep control

The global Timestep associated with the full IntegrationStepSequence can be set with either a fixed value, or a MultiplicativeArgument/Variable, using the min or max value of that argument within the spatial domain (should all be fluid variables).

fixedTimestep = Timestep(0.1) # Fixed timestep at 0.1 (in normalised time units)

variableTimestep = Timestep(0.1*a**2/b) # Timestep set to the minimum value of 0.1*a**2/b in the domain

variableTimestep = variableTimestep.max() # Now using the max value
print(variableTimestep.usingMaxVal)

variableTimestep = variableTimestep.min() # Now using the min value (default)
print(variableTimestep.usingMaxVal)

True
False

IntegrationScheme and setting a step sequence

The highest level time integration object is the IntegrationScheme, which contains an IntegrationStepSequence and Timestep, fully defining the integration scheme.

If only a single step is required (no splitting scheme), the scheme can be initialised directly with a Timestep and IntegrationStep

simpleScheme = IntegrationScheme(0.1,step1) # Uses a fixed Timestep and single step

We can also compose steps explicitly by calling them with the global step fraction we want to associate with them in the following way

compositeScheme = IntegrationScheme(variableTimestep) # Now using a variable step length

compositeScheme.steps = step1(0.8)*step2(1.0)*step1(0.2) # Steps are applied right to left, so the first one is step1 at 0.2 of the global timestep

for step in compositeScheme.steps:
    print(step.name) # Steps in the scheme get renamed to ensure no name conflicts
    print(step.stepFraction)

step12
0.2
step20
1.0
step11
0.8

Finally, the timestepping mode can be selected for the used scheme.

compositeScheme.setFixedNumTimesteps(numTimesteps=1000, # Run for 1000 (global) steps
                                     outputInterval=100 # Output every 100 steps
                                     )

# or

compositeScheme.setOutputPoints(outputPoints=[1.0,5.0,100.0]) # Output at given (normalised) time values - complete run when all outputs performed.