ode: implicit steppers

Defined in module: pressio4py.ode

Import as:       from pressio4py import ode

Overview

Provides functionalities to create steppers for implicit methods. Recall that implicit methods update the state of a system by solving a system of equations involving both the current and next state. An implicit stepper is an object that knows how to do one such implicit step.

Pressio implicit steppers are applicable to any system written in continuous-time form:

\[ \frac{d \boldsymbol{y}}{dt} = \boldsymbol{f}(\boldsymbol{y},t; ...) \]

and/or in a discrete-time form

\[ \boldsymbol{R}(\boldsymbol{y}, \boldsymbol{y_{n-1}}, ..., t_n, dt_n; ...) = \boldsymbol{0} \]

Here, $y$ is the state, $f$ the velocity, $t$ is time, and $R$ is the residual.

API, Parameters and Requirements

stepper = ode.create_implicit_stepper(scheme, state, system)
  • scheme:
    • value from the ode.stepscheme enum setting the desired stepping scheme:
enum valueMethodDiscrete Residual Formula
BDF1Backward Diff 1st order $R = y_{n+1}-y_{n}- hf(t_{n+1},y_{n+1})$
BDF2Backward Diff 2nd order $R = y_{n+1}-{\tfrac {4}{3}}y_{n}+{\tfrac {1}{3}}y_{n-1} - {\tfrac {2}{3}}hf(t_{n+1},y_{n+1})$
CrankNicolsonCrank-Nicolson $R = y_{n+1}- y_{n} - {\tfrac {1}{2}} h \left( f(t_{n+1},y_{n+1}) + f(t_{n},y_{n}) \right)$
  • state: -numpy.array storing your state
  • system:
    • object defining how to create an instance of the velocity $f$ and how to compute it.

    • Must conform to the following API:

      class MySys:
  def __init__(self):
    pass

  def createVelocity(self):
    return np.zeros(...)

  def velocity(self, stateIn, time, f):
    # compute f as needed
    # f[:] = ...

  def createJacobian(self):
    return np.zeros((...))

  def jacobian(self, stateIn, time, J):
      # compute J as needed
      # make sure to use J[:] to overwrite value

Note that currently, the implicit steppers are implemented only for dense Jacobians. This is on purpose, because pybind11 does not support passing by reference sparse types. Therefore, for the time being, we do not provide bindings for doing implicit stepping for systems with sparse Jacobians. Note, also, that this is not critical for the main purposes of this library because ROMs are inherently dense.

Stepper object

The returned stepper object exposes the following methods:

class Stepper:

  def order():
    return # order of the step scheme of this stepper instantiation

  def __call__(state, current_time, dt, step_number, solver)

  def createResidual()
    return # a residual instance

  def createJacobian()
    return # a Jacobian instance

  def residual(state, R)
  def jacobian(state, J)

What to do after a stepper is created?

Any stepper created using the functions above is guaranteed to satisfy the "steppable" concept discussed here. Therefore, once you create a stepper, you can use the advancers to step forward or you can use your own.
An example is below:

todo