pymor.timestepping

Extended functionality for pyMOR time steppers

Imports

Testing

from xarray import Coordinates
from pymor.basic import RectDomain, ConstantFunction, ExpressionFunction, StationaryProblem, \
    discretize_stationary_cg, thermal_block_problem, InstationaryProblem, discretize_instationary_cg, \
    NumpyMatrixOperator, ExpressionParameterFunctional, Parameters, InstationaryModel, Mu

A reference calculation using ImplicitEulerTimestepper:

import numpy as np
import matplotlib.pyplot as plt

from pymor.basic import (
    NumpyMatrixOperator,
    LincombOperator,
    ExpressionParameterFunctional,
    Parameters,
    InstationaryModel,
)
from pymor.vectorarrays.numpy import NumpyVectorSpace
from pymor.algorithms.timestepping import ImplicitEulerTimeStepper

# 2D state space [x, v]
space = NumpyVectorSpace(2)

A = NumpyMatrixOperator(np.array([
    [0.0, -1.0],
    [1.0,  0.0],
]))
B = NumpyMatrixOperator(np.array([
    [0.0, 0.0],
    [0.0, 1.0],
]))

operator = LincombOperator(
    [A, B],
    [
        ExpressionParameterFunctional("w0**2", Parameters({"w0": 1})),
        ExpressionParameterFunctional("c", Parameters({"c": 1})),
    ],
)

rhs = space.from_numpy(np.array([0.1, 0.0]))      # length-1 VectorArray
initial_data = space.from_numpy(np.array([1.0, 0.0]))
T = 100.0

model = InstationaryModel(
    T=T,
    initial_data=initial_data,
    operator=operator,
    rhs=rhs,
    time_stepper=ImplicitEulerTimeStepper(nt=10000),  # built-in pyMOR
    num_values=200,
)

mu = model.parameters.parse({"c": 0.1, "w0": "sin(t + .001)"})
reference_sol = model.solve(mu)
reference_sol_time = np.linspace(0, T, 200)
plt.plot(reference_sol_time, reference_sol.to_numpy()[0])
plt.plot(reference_sol_time, reference_sol.to_numpy()[1])

API

---

source

TimeStepper.solve

def solve(
    initial_time, # The time at which to begin time-stepping.
    end_time, # The time until which to perform time-stepping.
    initial_data, # The solution vector at initial_timeinitial_time.
    operator, # The |Operator| A.
    rhs:NoneType=None, # The right-hand side F (either |VectorArray| of length 1 or |Operator| with
    # source.dim == 1source.dim == 1). If NoneNone, zero right-hand side is assumed.
    mass:NoneType=None, # The |Operator| M. If NoneNone, the identity operator is assumed.
    mu:NoneType=None, # |Parameter values| for which operatoroperator and rhsrhs are evaluated. The current
    # time is added to mumu with key tt.
    num_values:NoneType=None, # The number of returned vectors of the solution trajectory. If NoneNone, each
    # intermediate vector that is calculated is returned.
):

Apply time-stepper to the equation.

The equation is of the form ::

M(mu) * d_t u + A(u, mu, t) = F(mu, t),
                 u(mu, t_0) = u_0(mu).

Parameters

initial_time The time at which to begin time-stepping. end_time The time until which to perform time-stepping. initial_data The solution vector at initial_timeinitial_time. operator The |Operator| A. rhs The right-hand side F (either |VectorArray| of length 1 or |Operator| with source.dim == 1source.dim == 1). If NoneNone, zero right-hand side is assumed. mass The |Operator| M. If NoneNone, the identity operator is assumed. mu |Parameter values| for which operatoroperator and rhsrhs are evaluated. The current time is added to mumu with key tt. num_values The number of returned vectors of the solution trajectory. If NoneNone, each intermediate vector that is calculated is returned.

Returns

|VectorArray| containing the solution trajectory.

---

source

init_progress_widget

def init_progress_widget(
    initial_time, end_time
):

---

source

tick

def tick():

---

source

update_progress_widget

def update_progress_widget(
    widget, tock, t, last_t
):

---

source

init_cvode

def init_cvode(
    cvode_rhs, num_values
):

---

source

solve_cvode

def solve_cvode(
    solver, t_list, initial_data
):

---

source

AdamsTimeStepper

def AdamsTimeStepper():

Interface for time-stepping algorithms.

Algorithms implementing this interface solve time-dependent initial value problems of the form ::

M(mu) * d_t u + A(u, mu, t) = F(mu, t),
                 u(mu, t_0) = u_0(mu).

Time-steppers used by |InstationaryModel| have to fulfill this interface.

Comparison with the reference calculation:

dof = Coordinates({'dof': ('dof', ['x', 'v'], {'long_name': "Degrees of freedom"})})

space = XarrayVectorSpace(dof)

A = XarrayMatrixOperator(np.array([[0., -1.], [1., 0]]), source=space, range=space)
B = XarrayMatrixOperator(np.array([[0., 0.], [0., 1.]]), source=space, range=space)
operator = LincombOperator(
    [A, B], 
    [
        ExpressionParameterFunctional("w0**2", Parameters({'w0': 1})), 
        ExpressionParameterFunctional("c", Parameters({'c': 1}))
    ]
)
rhs = space.from_numpy(np.array([0.1, 0.0]))
initial_data = space.from_numpy(np.array([1.0, 0.0]))
T = 100.0

model = InstationaryModel(T, initial_data, operator, rhs, time_stepper=AdamsTimeStepper(), num_values=200)

mu = dict(c=.1, w0="sin(t)")
mu = model.parameters.parse(mu)
adams_sol = model.solve(dict(c=.1, w0="sin(t+.001)"))
adams_sol.visualize()

Unable to display output for mime type(s): application/vnd.plotly.v1+json

plt.plot(adams_sol.array.Time, adams_sol.to_numpy()[0])
plt.plot(adams_sol.array.Time, adams_sol.to_numpy()[1])
plt.plot(reference_sol_time, reference_sol.to_numpy()[0])
plt.plot(reference_sol_time, reference_sol.to_numpy()[1])

The thermal block problem:

# domain = RectDomain([[0.,0.], [1.,1.]])
# diffusion = ConstantFunction(1, 2)
# rhs = ExpressionFunction('(sqrt( (x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 1.', 2)
# problem = StationaryProblem(
#    domain=domain,
#    diffusion=diffusion,
#    rhs=rhs,
# )
# m, data = discretize_stationary_cg(problem, diameter=1)

# p = thermal_block_problem([2,2])
# m, _ = discretize_stationary_cg(p, diameter=1)
# pp = InstationaryProblem(p, initial_data=ConstantFunction(0., 2), T=1.)
# mm, _ = discretize_instationary_cg(pp, nt=10, diameter=1)

# mm = mm.with_(time_stepper=AdamsTimeStepper(), num_values=200, operator=mm.operator, mass=None, output_functional=None)

# mm.solve({'diffusion': [.5, .6, .7, .8]});

# def numpy_to_xarray(op):
#     if isinstance(op, NumpyMatrixOperator): return XarrayMatrixOperator(op.matrix, source='elements', range='elements', name=op.name if op.name != 'NumpyMatrixOperator' else None)
#     if isinstance(op, LincombOperator): return LincombOperator([numpy_to_xarray(o) for o in op.operators], op.coefficients, name=op.name)
#     if isinstance(op, VectorOperator): return LincombOperator([numpy_to_xarray(o) for o in op.operators], op.coefficients, name=op.name)

# mm.operator

# range_space = numpy_to_xarray(mm.operator).range
# source_space = numpy_to_xarray(mm.operator).source

# mmm = mm.with_(
#     T=mm.T,
#     initial_data=source_space.from_numpy(mm.initial_data.array.to_numpy()),
#     operator=numpy_to_xarray(mm.operator),
#     rhs=range_space.from_numpy(mm.rhs.as_range_array().to_numpy()),
#     mass=None,
#     output_functional=None,
#     time_stepper=AdamsTimeStepper()
# )

# mmm.solve({'diffusion': [.5, .6, .7, .8]})

---

source

BDFTimeStepper

def BDFTimeStepper(
    preconditioner:NoneType=None, spilu_drop_tol:float=0.0001, spilu_fill_factor:int=20,
    spilu_drop_rule:NoneType=None, spilu_permc_spec:str='COLAMD'
):

Interface for time-stepping algorithms.

Algorithms implementing this interface solve time-dependent initial value problems of the form ::

M(mu) * d_t u + A(u, mu, t) = F(mu, t),
                 u(mu, t_0) = u_0(mu).

Time-steppers used by |InstationaryModel| have to fulfill this interface.

Comparison with the reference calculation:

dof = Coordinates({'dof': ('dof', ['x', 'v'], {'long_name': "Degrees of freedom"})})
space = XarrayVectorSpace(dof)
A = XarrayMatrixOperator(np.array([[0., -1.], [1., 0]]), source=space, range=space)
B = XarrayMatrixOperator(np.array([[0., 0.], [0., 1.]]), source=space, range=space)
operator = LincombOperator(
    [A, B], 
    [
        ExpressionParameterFunctional("w0**2", Parameters({'w0': 1})), 
        ExpressionParameterFunctional("c", Parameters({'c': 1}))
    ]
)
rhs = space.from_numpy(np.array([[.1], [0.]]))
initial_data = space.from_numpy(np.array([1., 0.]))
T = 100
model = InstationaryModel(T, initial_data, operator, rhs, time_stepper=BDFTimeStepper(), num_values=100)
bdf_sol = model.solve(dict(c=.1, w0="sin(t)"))

bdf_sol.visualize()

Unable to display output for mime type(s): application/vnd.plotly.v1+json

plt.plot(bdf_sol.array.Time, bdf_sol.to_numpy()[0])
plt.plot(bdf_sol.array.Time, bdf_sol.to_numpy()[1])

plt.plot(reference_sol_time, reference_sol.to_numpy()[0])
plt.plot(reference_sol_time, reference_sol.to_numpy()[1])

The thermal block problem:

p = thermal_block_problem([2,2])
m, _ = discretize_stationary_cg(p)
pp = InstationaryProblem(p, initial_data=ConstantFunction(0., 2), T=1.)
mm, _ = discretize_instationary_cg(pp, nt=10)

solver_options={'inverse': {'type': 'scipy_lgmres_spilu', 'preconditioner_bandwidth': 2}}

mm = mm.with_(time_stepper=BDFTimeStepper(), num_values=200, operator=mm.operator)

mm.solve({'diffusion': [.5, .6, .7, .8]})

NumpyVectorArray(
    NumpyVectorSpace(20201),
    [[0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
      0.00000000e+00 0.00000000e+00]
     [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
      0.00000000e+00 0.00000000e+00]
     [0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
      0.00000000e+00 0.00000000e+00]
     ...
     [0.00000000e+00 2.75757447e-10 2.03167753e-06 ... 1.81201470e-05
      2.54533008e-05 3.17670706e-05]
     [0.00000000e+00 2.75757447e-10 2.03080205e-06 ... 1.79186198e-05
      2.48839961e-05 3.06432543e-05]
     [0.00000000e+00 2.75753797e-10 1.87420727e-06 ... 1.35152459e-05
      1.73482117e-05 2.00949822e-05]],
    _len=8)

---

source

StiffSwitchingTimeStepper

def StiffSwitchingTimeStepper():

Interface for time-stepping algorithms.

Algorithms implementing this interface solve time-dependent initial value problems of the form ::

M(mu) * d_t u + A(u, mu, t) = F(mu, t),
                 u(mu, t_0) = u_0(mu).

Time-steppers used by |InstationaryModel| have to fulfill this interface.