Commit dde4b3af authored by Tim Keil's avatar Tim Keil
Browse files

add dual_rhs and dual_operator to StationaryModel

parent bcbde26b
......@@ -44,6 +44,11 @@ class StationaryModel(Model):
an `estimate_error(U, mu, m)` method. If `error_estimator` is
not `None`, an `estimate_error(U, mu)` method is added to the
model which will call `error_estimator.estimate_error(U, mu, self)`.
dual_operator
The dual |Operator| L* of L such that p(μ) solves
L*(p(μ),μ) = G(μ)
dual_rhs
The right hand side vector G of the dual problem.
visualizer
A visualizer for the problem. This can be any object with
a `visualize(U, m, ...)` method. If `visualizer`
......@@ -55,7 +60,8 @@ class StationaryModel(Model):
"""
def __init__(self, operator, rhs, output_functional=None, products=None,
error_estimator=None, visualizer=None, name=None):
error_estimator=None, dual_operator=None, dual_rhs=None,
visualizer=None, name=None):
if isinstance(rhs, VectorArray):
assert rhs in operator.range
......@@ -86,13 +92,14 @@ class StationaryModel(Model):
@property
def dual(self):
"""Instantiate the dual model which is used to solve for a dual solution of the |StationaryModel|."""
try:
return self._dual
except AttributeError:
assert self.output_functional is not None
assert self.output_functional.linear
# TODO: assert that the operator is symmetric
self._dual = self.with_(rhs=self.output_functional.H)
assert self.dual_rhs is not None
assert self.dual_rhs.linear
assert self.dual_operator is not None
self._dual = self.with_(operator=self.dual_operator, rhs=self.dual_rhs)
return self._dual
def _compute_solution_d_mu(self, parameter, index, solution, mu):
......
......@@ -463,11 +463,6 @@ class Model(CacheableObject, ParametricObject):
)
return data['output_d_mu']
@property
def dual(self):
"""Instantiate the dual model which is used to solve for a dual solution of the |Model|."""
return NotImplemented
def estimate_error(self, mu=None, **kwargs):
"""Estimate the error for the computed internal state.
......
......@@ -179,6 +179,10 @@ class StationaryRBReductor(ProjectionBasedReductor):
'products': {k: project(v, RB, RB) for k, v in fom.products.items()},
'output_functional': project(fom.output_functional, None, RB) if fom.output_functional else None
}
if self.fom.dual_operator:
projected_operators['dual_operator'] = project(fom.dual_operator, RB, RB)
if self.fom.dual_rhs:
projected_operators['dual_rhs'] = project(fom.dual_rhs, RB, None)
return projected_operators
def project_operators_to_subbasis(self, dims):
......@@ -191,6 +195,10 @@ class StationaryRBReductor(ProjectionBasedReductor):
'output_functional': (project_to_subbasis(rom.output_functional, None, dim)
if rom.output_functional else None)
}
if self.fom.dual_operator:
projected_operators['dual_operator'] = project_to_subbasis(rom.dual_operator, dim, dim)
if self.fom.dual_rhs:
projected_operators['dual_rhs'] = project_to_subbasis(rom.dual_rhs, dim, None)
return projected_operators
def build_rom(self, projected_operators, error_estimator):
......
......@@ -51,6 +51,8 @@ def create_fom(args):
mu_bar = problem.parameters.parse([np.pi/2,np.pi/2])
fom, _ = discretize_stationary_cg(problem, diameter=1. / int(args['GRID_INTERVALS']),
mu_energy_product=mu_bar)
# define the dual quantities
fom = fom.with_(dual_operator=fom.operator, dual_rhs=fom.output_functional.H)
return fom, mu_bar
......
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