Time Stepping

5. Time Stepping#

st_ode_options(varargin)#

Options for the kinematic condition time stepping.

The admissible fields are

  • display: display additional information per timestep.

  • stepper: time stepping method: 'euler', 'imex', 'ssprk3'.

  • tmax: final time for the simulation.

  • maxit: maximum number of timesteps.

  • dt: time step. If negative, then it will be determined at runtime.

  • theta: Courant number used in a CFL condition if the time step is not fixed.

  • utol: tolerance on the max norm of the normal velocity field for steady state computations.

  • xtol: tolerance on the relative interface deformation for steady state computations.

  • odefreq: output frequency based on the time step.

  • odeplot(t, x, y, st): output function handle that takes the current time, the geometry information and the solution at the current time.

  • solver: integral representation used for solving for the velocity field and other variables. This is only used in st_ode_steady().

  • filter_strength: used to smooth the interface at each time step using st_ax_filter().

  • checkpoints: a cell array of files that can be used to retrieve forward solutions when solving the adjoint equations in st_ode_unsteady_adjoint().

Additional parameters can be passed in an unchecked 'user' field. The default values can be obtained by calling this function without any parameters.

Returns:

options struct will missing fields filled in with their default values.

st_ode_unsteady(x0, y0, bc, forcing, cost, checkpoint, varargin)#

Advance the unsteady kinematic condition in time.

This is similar to st_ode_steady(), but with a very different purpose. While st_ode_steady() is meant to drive the system to steady state as efficiently as possible, this function is meant for adjoint optimization.

As such, it uses a simpler time stepper and computes the full velocity field at each time step, instead of relying on a fixed point-type iteration. Furthermore, it allows checkpointing and other adjoint-focused features. However, it can still be used to drive the system to steady state in the same way, but we expect it to be slower (in computational time).

The system we are solving is the kinematic condition

\[\frac{\partial \mathbf{X}}{\partial t} = \mathbf{f}(\mathbf{u}, \mathbf{X}).\]

The right-hand side of the kinematic condition is described by forcing terms conforming to st_forcing_options(). If a cost struct is passed, with information about the cost functional, we computed it in the following way

\[\begin{split}\begin{aligned} J =\,\, & j_\Sigma(T) + \int_0^T j_\Sigma(t) \,\mathrm{d}[t] \\ \approx\,\, & j_\Sigma(T) + \sum_{n = 0}^N \Delta t^n j_\Sigma(t^n), \end{aligned}\end{split}\]

where \(j_\Sigma(t)\) is computed by cost.evaluate. This is a first-order approximation of the time integral, which matches the first-order forward Euler method we are using for the time integration.

Providing the cost functional is optional, but if it is provided, the following two fields are expected:

  • evaluate: a function handle taking (t, x, y, st), for evaluating the integrand \(j_\Sigma(t)\) of the cost functional at a given time \(t\), where st is a struct containing the solution to the Stokes problem.

  • variables: a cell array of variables necessary for the right-hand side and the cost functional, which is passed directly to st_repr_density_solution().

The additional ode are meant for the time stepper and the forcing term in the kinematic condition. They are defined in the common st_ode_options().

Parameters:

  • x0 – initial target geometry information.

  • y0 – source geometry information. This only requires information about the quadrature rules, as provided by st_point_quadrature() and st_layer_quadrature(). The actual geometry is updated from the target points during the evolution.

  • fn – boundary conditions for the Stokes problem.

  • cost – optional cost functional information.

  • checkpoint – created using st_cache_create().

  • ode – time stepper information.

Returns:

a tuple (x, y, state), where x and y contain the geometry information at tmax. The state contains information about the time stepping, like all the checkpoints that were saved, the times at which they were saved, the value of the cost functional, etc.

st_ode_steady(x0, y0, forcing, bc, checkpoint, varargin)#

Drive a drop to steady state.

The kinematic equation for the drop is given by:

\[\frac{\partial \mathbf{X}}{\partial t} = \mathbf{f}(\mathbf{u}, \mathbf{X}).\]

where \(f\) is given by forcing (see also st_forcing_options()). The equation above is solved using a second-order Runge-Kutta time integrator, where each step proceeds as follows:

  • First, the right-hand side of Equation 5.2.9 from [Pozrikidis1992] is evaluated for the current velocity. For the first time step, we initialize the velocity field to zero, so the double-layer is not evaluated.

  • Once an intermediary velocity has been obtained, we advance the solution using the kinematic condition.

  • At the end of the time step, the resulting geometry is smoothed, the points are equidistanced, the solution is plotted, etc., as requested by the user.

  • If not given explicitly, at every iteration, the timestep is recomputed using the formula (inspired by [Denner2015]):

    \[\Delta t = \theta \min_i \sqrt{Ca \Delta s_i^3}.\]

  • The time stepping is stopped if one of two conditions is satisfied:

    1. if \(\|\mathbf{u} \cdot \mathbf{n}\|_\infty < \epsilon_u\), which is the physical steady state condition, or

    2. if \(\|x^{n + 1} - x^{n}\|_\infty < \epsilon_x\), which checks if the time stepping has stagnated. As a rule, \(\epsilon_x \ll \epsilon_u\) to achieve desired results.

This pseudo-time integration is effectively a fixed point iteration for the pair \((\mathbf{x}, \mathbf{u})\). Note that, since this is the case, the velocity field \(\mathbf{u}\) only solves the Stokes equations once steady state was reached to a sufficiently small tolerance. All intermediate values cannot be considered valid solutions. The actual solution should be obtained only after steady state has been reached using, e.g., st_repr_velocity(). Exact solutions at each time step can be obtained by passing solver = 'density', but this is significantly slower.

Parameters:

  • x0 – initial target geometry information.

  • y0 – initial source geometry information.

  • bc – boundary conditions for the Stokes equations. See st_boundary_condition_options() for the format.

  • checkpoint – a cache created by st_cache_create(). We make use of its fileseries field to output the iterations.

  • varargin – options passed to st_ode_options(). Options can be provided by a struct followed by additional pairs of options.

Returns:

a tuple (x, y, state) of the geometry at the final time and a state struct with additional information about the simulation.

st_ode_unsteady_load(cache, x0, y0, subfolder, prefix)#

Load the final ODE state from a given cache.

Parameters:

  • cache – existing cache (see st_cache_create()) or a cache directory.

  • x0 – initial target point information used to reconstruct the geometry.

  • y0 – initial source point information.

  • subfolder – if given, a new cache is created at this subfolder with the same prefix as cache.

  • prefix – prefix for the cache created in subfolder, if given.

Returns:

same state as st_ode_unsteady().