Merge pull request #1621 from pbuchfink/models-symplectic

Added basic models for symplectic framework

src/pymor/algorithms/timestepping.py

@@ -117,6 +117,35 @@ class ExplicitEulerTimeStepper(TimeStepper):
         return explicit_euler(operator, rhs, initial_data, initial_time, end_time, self.nt, mu, num_values)
 
 
+class ImplicitMidpointTimeStepper(TimeStepper):
+    """Implict midpoint rule time-stepper. Symplectic integrator + preserves quadratic invariants.
+
+    Solves equations of the form ::
+
+        M * d_t u + A(u, mu, t) = F(mu, t).
+
+    Parameters
+    ----------
+    nt
+        The number of time-steps the time-stepper will perform.
+    solver_options
+        The |solver_options| used to invert `M - dt/2*A`.
+        The special values `'mass'` and `'operator'` are
+        recognized, in which case the solver_options of
+        M (resp. A) are used.
+    """
+
+    def __init__(self, nt, solver_options='operator'):
+        self.__auto_init(locals())
+
+    def solve(self, initial_time, end_time, initial_data, operator, rhs=None, mass=None, mu=None,
+              num_values=None):
+        if not operator.linear:
+            raise NotImplementedError
+        return implicit_midpoint_rule(operator, rhs, mass, initial_data, initial_time, end_time,
+                                      self.nt, mu, num_values, solver_options=self.solver_options)
+
+
 def implicit_euler(A, F, M, U0, t0, t1, nt, mu=None, num_values=None, solver_options='operator'):
     assert isinstance(A, Operator)
     assert isinstance(F, (type(None), Operator, VectorArray))
@@ -230,6 +259,73 @@ def explicit_euler(A, F, U0, t0, t1, nt, mu=None, num_values=None):
     return R
 
 
+def implicit_midpoint_rule(A, F, M, U0, t0, t1, nt, mu=None, num_values=None, solver_options='operator'):
+    assert isinstance(A, Operator)
+    assert isinstance(F, (type(None), Operator, VectorArray))
+    assert isinstance(M, (type(None), Operator))
+    assert A.source == A.range
+    num_values = num_values or nt + 1
+    dt = (t1 - t0) / nt
+    DT = (t1 - t0) / (num_values - 1)
+
+    if F is None:
+        F_time_dep = False
+    elif isinstance(F, Operator):
+        assert F.source.dim == 1
+        assert F.range == A.range
+        F_time_dep = _depends_on_time(F, mu)
+        if not F_time_dep:
+            dt_F = F.as_vector(mu) * dt
+    else:
+        assert len(F) == 1
+        assert F in A.range
+        F_time_dep = False
+        dt_F = F * dt
+
+    if M is None:
+        from pymor.operators.constructions import IdentityOperator
+        M = IdentityOperator(A.source)
+
+    assert A.source == M.source == M.range
+    assert not M.parametric
+    assert U0 in A.source
+    assert len(U0) == 1
+
+    R = A.source.empty(reserve=nt+1)
+    R.append(U0)
+
+    if solver_options == 'operator':
+        options = A.solver_options
+    elif solver_options == 'mass':
+        options = M.solver_options
+    else:
+        options = solver_options
+
+    M_dt_A_impl = (M + A * (dt/2)).with_(solver_options=options)
+    if not _depends_on_time(M_dt_A_impl, mu):
+        M_dt_A_impl = M_dt_A_impl.assemble(mu)
+    M_dt_A_expl = (M - A * (dt/2)).with_(solver_options=options)
+    if not _depends_on_time(M_dt_A_expl, mu):
+        M_dt_A_expl = M_dt_A_expl.assemble(mu)
+
+    t = t0
+    U = U0.copy()
+
+    for n in range(nt):
+        mu = mu.with_(t=t + dt/2)
+        t += dt
+        rhs = M_dt_A_expl.apply(U, mu=mu)
+        if F_time_dep:
+            dt_F = F.as_vector(mu) * dt
+        if F:
+            rhs += dt_F
+        U = M_dt_A_impl.apply_inverse(rhs, mu=mu)
+        while t - t0 + (min(dt, DT) * 0.5) >= len(R) * DT:
+            R.append(U)
+
+    return R
+
+
 def _depends_on_time(obj, mu):
     if not mu:
         return False

src/pymor/models/symplectic.py

+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright pyMOR developers and contributors. All rights reserved.
+# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+import numpy as np
+from pymor.algorithms.timestepping import ImplicitMidpointTimeStepper
+from pymor.models.basic import InstationaryModel
+from pymor.operators.constructions import ConcatenationOperator, NumpyConversionOperator, VectorOperator
+from pymor.operators.interface import Operator
+from pymor.operators.numpy import NumpyVectorSpace
+from pymor.operators.symplectic import CanonicalSymplecticFormOperator
+from pymor.parameters.base import Mu
+from pymor.vectorarrays.block import BlockVectorSpace
+from pymor.vectorarrays.interface import VectorArray
+
+
+class QuadraticHamiltonianModel(InstationaryModel):
+    """Generic class for quadratic Hamiltonian systems.
+
+    This class describes Hamiltonian systems given by the equations::
+
+        ∂_t u(t, μ) = J * H_op(t, μ) * u(t, μ) + J * h(t, μ)
+            u(0, μ) = x_0(μ)
+
+    for t in [0,T], where H_op is a linear time-dependent |Operator|,
+    J is a canonical Poisson matrix, h is a (possibly) time-dependent
+    vector-like |Operator|, and x_0 the initial data.
+    The right-hand side of the Hamiltonian equation is J times the
+    gradient of the Hamiltonian
+
+        Ham(u, t, μ) = 1/2* u * H_op(t, μ) * u + u * h(t, μ)
+
+    Parameters
+    ----------
+    T
+        The final time T.
+    initial_data
+        The initial data `x_0`. Either a |VectorArray| of length 1 or
+        (for the |Parameter|-dependent case) a vector-like |Operator|
+        (i.e. a linear |Operator| with `source.dim == 1`) which
+        applied to `NumpyVectorArray(np.array([1]))` will yield the
+        initial data for a given |Parameter|.
+    H_op
+        The |Operator| H_op.
+    h
+        The state-independet part of the Hamiltonian h.
+    time_stepper
+        The :class:`time-stepper <pymor.algorithms.timestepping.TimeStepper>`
+        to be used by :meth:`~pymor.models.interface.Model.solve`.
+        Alternatively, the parameter nt can be specified to use the
+        :class:`implicit midpoint rule <pymor.algorithms.timestepping.ImplicitMidpointTimeStepper>`.
+    nt
+        If time_stepper is `None` and nt is specified, the
+        :class:`implicit midpoint rule <pymor.algorithms.timestepping.ImplicitMidpointTimeStepper>`
+        as time_stepper.
+    num_values
+        The number of returned vectors of the solution trajectory. If `None`, each
+        intermediate vector that is calculated is returned.
+    output_functional
+        |Operator| mapping a given solution to the model output. In many applications,
+        this will be a |Functional|, i.e. an |Operator| mapping to scalars.
+        This is not required, however.
+    visualizer
+        A visualizer for the problem. This can be any object with
+        a `visualize(U, m, ...)` method. If `visualizer`
+        is not `None`, a `visualize(U, *args, **kwargs)` method is added
+        to the model which forwards its arguments to the
+        visualizer's `visualize` method.
+    name
+        Name of the model.
+    """
+
+    def __init__(self, T, initial_data, H_op, h=None, time_stepper=None, nt=None, num_values=None,
+                 output_functional=None, visualizer=None, name=None):
+        assert isinstance(H_op, Operator) and H_op.linear and H_op.range == H_op.source \
+               and H_op.range.dim % 2 == 0
+
+        # interface to use ImplicitMidpointTimeStepper via parameter nt
+        if (time_stepper is None) == (nt is None):
+            raise ValueError('Either specify time_stepper or nt (not both)')
+        if time_stepper is None:
+            time_stepper = ImplicitMidpointTimeStepper(nt)
+
+        if isinstance(H_op.range, NumpyVectorSpace):
+            # make H_op compatible with blocked phase_space
+            assert H_op.range.dim % 2 == 0, 'H_op.range has to be even dimensional'
+            half_dim = H_op.range.dim // 2
+            phase_space = BlockVectorSpace([NumpyVectorSpace(half_dim)] * 2)
+            H_op = ConcatenationOperator([
+                NumpyConversionOperator(phase_space, direction='from_numpy'),
+                H_op,
+                NumpyConversionOperator(phase_space, direction='to_numpy'),
+            ])
+
+        if h is None:
+            h = H_op.range.zeros()
+
+        if isinstance(h, VectorArray):
+            h = VectorOperator(h, name='h')
+
+        if isinstance(h.range, NumpyVectorSpace):
+            # make h compatible with blocked phase_space
+            assert h.range.dim % 2 == 0, 'h.range has to be even dimensional'
+            half_dim = h.range.dim // 2
+            phase_space = H_op.range
+            h = ConcatenationOperator([
+                NumpyConversionOperator(phase_space, direction='from_numpy'),
+                h,
+            ])
+        assert h.range is H_op.range
+
+        if isinstance(initial_data.space, NumpyVectorSpace):
+            # make initial_data compatible with blocked phase_space
+            initial_data = H_op.source.from_numpy(initial_data.to_numpy())
+
+        # J based on blocked phase_space
+        self.J = CanonicalSymplecticFormOperator(H_op.source)
+
+        # minus is required since operator in an InstationaryModel is on the LHS
+        operator = -ConcatenationOperator([self.J, H_op])
+        rhs = ConcatenationOperator([self.J, h])
+        super().__init__(T, initial_data, operator, rhs,
+                         time_stepper=time_stepper,
+                         num_values=num_values,
+                         output_functional=output_functional,
+                         visualizer=visualizer,
+                         name=name)
+        self.__auto_init(locals())
+
+    def eval_hamiltonian(self, u, mu=None):
+        """Evaluate a quadratic Hamiltonian function
+
+        Ham(u, t, μ) = 1/2 * u * H_op(t, μ) * u + u * h(t, μ).
+        """
+        if not isinstance(mu, Mu):
+            mu = self.parameters.parse(mu)
+        assert self.parameters.assert_compatible(mu)
+        # compute linear part
+        ham_h = self.h.apply_adjoint(u, mu=mu)
+        # compute quadratic part
+        ham_H = ham_h.space.make_array(self.H_op.pairwise_apply2(u, u, mu=mu)[:, np.newaxis])
+
+        return 1/2 * ham_H + ham_h

src/pymortests/model.py

@@ -3,12 +3,18 @@
 # License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
 
 import numpy as np
+import pytest
 
 from pymor.algorithms.basic import almost_equal
+from pymor.algorithms.timestepping import ImplicitMidpointTimeStepper
 from pymor.analyticalproblems.functions import ExpressionFunction, ConstantFunction
 from pymor.analyticalproblems.thermalblock import thermal_block_problem
 from pymor.discretizers.builtin import discretize_stationary_cg
 from pymor.core.pickle import dumps, loads
+from pymor.models.symplectic import QuadraticHamiltonianModel
+from pymor.operators.block import BlockDiagonalOperator
+from pymor.operators.constructions import IdentityOperator
+from pymor.vectorarrays.numpy import NumpyVectorSpace
 from pymortests.base import runmodule
 from pymortests.pickling import assert_picklable, assert_picklable_without_dumps_function
 
@@ -49,5 +55,31 @@ def test_StationaryModel_deaffinize():
     assert np.allclose(m.output(mu), m_deaff.output(mu))
 
 
+@pytest.mark.parametrize('block_phase_space', (False, True))
+def test_quadratic_hamiltonian_model(block_phase_space):
+    """Check QuadraticHamiltonianModel with implicit midpoint rule."""
+    if block_phase_space:
+        H_op = BlockDiagonalOperator([IdentityOperator(NumpyVectorSpace(1))] * 2)
+    else:
+        phase_space = NumpyVectorSpace(2)
+        H_op = IdentityOperator(phase_space)
+
+    model = QuadraticHamiltonianModel(
+        3,
+        H_op.source.ones(),
+        H_op,
+        h=H_op.range.from_numpy(np.array([1, 0])),
+        time_stepper=ImplicitMidpointTimeStepper(3),
+        name='test_mass_spring_model'
+    )
+
+    res = model.solve()
+
+    ham = model.eval_hamiltonian(res).to_numpy()
+
+    # check preservation of the Hamiltonian
+    assert np.allclose(ham, ham[0])
+
+
 if __name__ == "__main__":
     runmodule(filename=__file__)