Adaptive time stepping

[SteinStoter]

Hi,

I would like to gradually increase the time-step size used by the thetamethod classes. In the docstring it says that the timestep provided to the constructor is the "Initial time step, will scale up as residual decreases". Does that actually happen though? If so, how do I keep track of and tweak that? If not, what is my best option for forcing that?

Best,
Stein

[joostvanzwieten]

The docstring of thetamethod is actually incorrect. This is fixed by #686. The steps yielded by thetamethod have a fixed timestep. However, if the newton solver fails to find a solution for a full timestep, thetamethod internally reduces the timestep by one half, repeatedly, until a solution is found.

With Nutils v7 you can use tensorial topologies to build a time integrator with variable step size yourself. The idea is to create a time topology with a single element, a time geometry that is scaled and shifted for each time step and space-time fields that are linear interpolations in time of spatial fields. As an example (for Nutils master, but adaptable to v7), here is an implementation of the time-dependent heat equation with doubling of the time step size at every iteration (disclaimer: I've not verified that results actually make sense):

import numpy
from nutils import export, function, mesh, solver
from nutils.expression_v2 import Namespace
import treelog

ns = Namespace()

# Space topology.
X, ns.x = mesh.rectilinear([numpy.linspace(-numpy.pi / 2, numpy.pi / 2, 8)]*2, space='X')
ns.define_for('x', gradient='∇', jacobians=('dV', 'dS'), normal='n')

# Time topology. The topology consists of a single element with linear, unit
# geometry `τ`. The geometry `t`, an affine transformation of `τ`, is the
# actual time for a timeslab, with argument `t0` defining the start of the
# timeslab and `Δt` the length.
T, ns.τ = mesh.line(1, space='T', bnames=('past', 'future'))
ns.add_field(('Δt', 't0'))
ns.t = 't0 + Δt τ'
ns.define_for('t', gradient='∂t', jacobians=('dt',), normal='nt')

# Fields. The field `u` is a linear interpolation of `u0` and `u1` in time.
ns.add_field(('u1', 'u0', 'v'), X.basis('spline', 2))
ns.u = 'u1 τ + u0 (1 - τ)'

# Samples of topology `T` at `τ ∈ {0, θ, 1}`.
θ = 1
Tat0 = T.boundary['past'].sample('gauss', 0)
Tatθ = T.locate(ns.τ, [θ], tol=1e-10, weights=numpy.array([1.]))
Tat1 = T.boundary['future'].sample('gauss', 0)

# Space-time weak form of the heat equation.
res = Tatθ.integral(X.integral('∇_i(v) ∇_i(u) dV dt' @ ns, degree=4))
res += (T.boundary * X).integral('v u dV nt' @ ns, degree=4)

# Dirichlet boundary condition.
cons = solver.optimize('u1,', X.boundary.integral('u1^2 dS' @ ns, degree=4), droptol=1e-10)

# Initial condition.
res0 = X.integral('(u0 - cos(x_0) cos(x_1))^2 dV' @ ns, degree=4)
args = solver.solve_linear('u0:u0,', res0)

# Sample for plotting.
smpl = X.sample('bezier', 5)
x = smpl.eval(ns.x)

# Time loop.
args.update(t0=0, Δt=0.01)
while args['t0'] < 1:
    with treelog.context(f't = {args["t0"]:.2f}'):
        # Solve for the next `u`.
        args = solver.solve_linear('u1:v,', res, constrain=cons, arguments=args)
        # Plot the old `u`.
        export.triplot('u.png', x, (Tat0 * smpl).eval(ns.u, **args), tri=smpl.tri, hull=smpl.hull, clim=[0, 1])
        # Update the arguments for the next iteration.
        args.update(u0=args['u1'], t0=args['t0'] + args['Δt'], Δt=2 * args['Δt'])

[SteinStoter]

That's a nice reference code that I'll definitely consider for future projects. For now, it's not feasible for me to switch to space-time though.

After having looked at it again, I believe this should do the trick:

time_stepper = solver.impliciteuler(dofs, res, inertia, timestep)
with treelog.iter.plain('timestep', time_stepper) as steps:
    for itime, lhs in enumerate(steps):
        time_stepper._wrapped.timestep *= 2

I find it a bit awkward that the actual impliciteuler instance is hiding in the _wrapped attribute of the "fake" impliciteuler instance, to be honest.

[joostvanzwieten]

That works for now. Please note that undocumented features may be changed in the future without deprecation period.