Unverified Commit 8e7aaa3e authored by Stephan Rave's avatar Stephan Rave Committed by GitHub
Browse files

Merge pull request #1410 from pymor/deaffinize

Add StationaryModel.deaffinize
parents 64bb8c5c 5e589040
Pipeline #102073 passed with stages
in 37 minutes and 35 seconds
......@@ -162,7 +162,14 @@ class ProjectRules(RuleTable):
# at least we can try to partially project the outer operators
projected_first = project(first, None, source_basis)
projected_last = project(last, range_basis, None)
return ConcatenationOperator((projected_last,) + op.operators[1:-1] + (projected_first,), name=op.name)
projected_op = ConcatenationOperator((projected_last,) + op.operators[1:-1] + (projected_first,), name=op.name)
# special handling for concatenations with ConstantOperators
# probably should be moved elsewhere
if not projected_op.parametric and any(isinstance(o, ConstantOperator) for o in projected_op.operators):
projected_op = ConstantOperator(projected_op.apply(projected_op.source.zeros()), projected_op.source)
return projected_op
@match_class(AdjointOperator)
def action_AdjointOperator(self, op):
......
......@@ -6,7 +6,7 @@ import numpy as np
from pymor.algorithms.timestepping import TimeStepper
from pymor.models.interface import Model
from pymor.operators.constructions import IdentityOperator, VectorOperator, ZeroOperator
from pymor.operators.constructions import IdentityOperator, VectorOperator, ZeroOperator, ConstantOperator
from pymor.vectorarrays.interface import VectorArray
from pymor.vectorarrays.numpy import NumpyVectorSpace
......@@ -146,6 +146,74 @@ class StationaryModel(Model):
else:
return gradients
def deaffinize(self, arg):
"""Build |Model| with linear solution space.
For many |Models| the solution manifold is contained in an
affine subspace of the :attr:`~pymor.models.interface.Model.solution_space`,
e.g. the affine space of functions with certain fixed boundary values.
Most MOR techniques, however, construct linear approximation spaces, which
are not fully contained in this affine subspace, even if these spaces are
created using snapshot vectors from the subspace. Depending on the
FOM, neglecting the affine structure of the solution space may lead to bad
approximations or even an ill-posed ROM. A standard approach to circumvent
this issue is to replace the FOM with an equivalent |Model| with linear
solution space. This method can be used to obtain such a |Model|.
Given a vector `u_0` from the affine solution space, the returned
:class:`StationaryModel` is equivalent to solving::
L(u(μ) + u_0, μ) = F(μ)
When :attr:`~StationaryModel.operator` is linear, the affine shift is
moved to the right-hand side to obtain::
L(u(μ), μ) = F(μ) - L(u_0, μ)
Solutions of the original |Model| can be obtained by adding `u_0` to the
solutions of the deaffinized |Model|.
The :attr:`~StationaryModel.output_functional` is adapted accordingly to
yield the same output for given |parameter values|.
Parameters
----------
arg
Either a |VectorArray| of length 1 containing the vector `u_0`.
Alternatively, |parameter values| can be provided, for which the
model is :meth:`solved <pymor.models.interface.Model.solve>` to
obtain `u_0`.
Returns
-------
The deaffinized |Model|.
"""
if not isinstance(arg, VectorArray):
mu = self.parameters.parse(arg)
arg = self.solve(mu)
assert arg in self.solution_space and len(arg) == 1
affine_shift = arg
if self.operator.linear:
new_operator = self.operator
new_rhs = self.rhs - self.operator @ VectorOperator(affine_shift)
else:
new_operator = \
self.operator @ (IdentityOperator(self.solution_space)
+ ConstantOperator(affine_shift, self.solution_space))
new_rhs = self.rhs
if self.output_functional is not None:
new_output_functional = \
self.output_functional @ (IdentityOperator(self.solution_space)
+ ConstantOperator(affine_shift, self.solution_space))
else:
new_output_functional = None
return self.with_(operator=new_operator, rhs=new_rhs, output_functional=new_output_functional,
error_estimator=None)
class InstationaryModel(Model):
"""Generic class for models of instationary problems.
......
......@@ -5,6 +5,9 @@
import numpy as np
from pymor.algorithms.basic import almost_equal
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 pymortests.base import runmodule
from pymortests.pickling import assert_picklable, assert_picklable_without_dumps_function
......@@ -27,5 +30,24 @@ def test_pickle_by_solving(model):
assert np.all(almost_equal(m.solve(mu), m2.solve(mu)))
def test_StationaryModel_deaffinize():
p = thermal_block_problem((2, 2)).with_(
dirichlet_data=ExpressionFunction('x[0]', 2),
outputs=[('l2', ConstantFunction(1., 2))]
)
m, _ = discretize_stationary_cg(p, diameter=1/10)
U_aff = m.solve([1, 1, 1, 1])
m_deaff = m.deaffinize(U_aff)
mu = m.parameters.parse([0.1, 10, 7, 1])
U = m.solve(mu)
U_deaff = m_deaff.solve(mu)
assert np.all(almost_equal(U, U_deaff + U_aff))
assert np.allclose(m.output(mu), m_deaff.output(mu))
if __name__ == "__main__":
runmodule(filename=__file__)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment