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

[models] remove the "dual" of a model but keep dual_operator

parent 02ac555f
......@@ -45,11 +45,8 @@ class StationaryModel(Model):
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(μ).
This equation is used for the adjoint approach to compute the gradient.
See Section 1.6.2 in [HPUU09]_ for more details.
dual_rhs
The right hand side vector G of the dual problem.
The dual |Operator| L* of L, which is e.g. used for computing the gradient
of the output_functional. See Section 1.6.2 in [HPUU09]_ for more details.
visualizer
A visualizer for the problem. This can be any object with
a `visualize(U, m, ...)` method. If `visualizer`
......@@ -61,7 +58,7 @@ class StationaryModel(Model):
"""
def __init__(self, operator, rhs, output_functional=None, products=None,
error_estimator=None, dual_operator=None, dual_rhs=None,
error_estimator=None, dual_operator=None,
visualizer=None, name=None):
if isinstance(rhs, VectorArray):
......@@ -91,18 +88,6 @@ class StationaryModel(Model):
def _compute_solution(self, mu=None, **kwargs):
return self.operator.apply_inverse(self.rhs.as_range_array(mu), mu=mu)
@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.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_single_direction(self, parameter, index, solution, mu):
lhs_d_mu = VectorOperator(self.operator.d_mu(parameter, index).apply(solution, mu=mu))
rhs_d_mu = self.rhs.d_mu(parameter, index)
......@@ -123,6 +108,7 @@ class StationaryModel(Model):
adjoint_approach
Use the adjoint approach for a more efficient way of computing the gradient.
See Section 1.6.2 in [HPUU09]_ for more details.
So far, this approach is only valid for linear output functionals and symmetric operators.
Returns
-------
......@@ -131,8 +117,12 @@ class StationaryModel(Model):
if adjoint_approach is False:
return super()._compute_output_d_mu(U, mu)
else:
assert self.output_functional is not None
assert self.output_functional.linear
assert self.dual_operator is not None
dual_problem = self.with_(operator=self.dual_operator, rhs=self.output_functional.H)
P = dual_problem.solve(mu)
gradient = []
P = self.dual.solve(mu)
for (parameter, size) in self.parameters.items():
for index in range(size):
output_partial_dmu = self.output_functional.d_mu(parameter, index).apply(U, mu=mu).to_numpy()[0,0]
......
......@@ -181,8 +181,6 @@ class StationaryRBReductor(ProjectionBasedReductor):
}
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):
......@@ -197,8 +195,6 @@ class StationaryRBReductor(ProjectionBasedReductor):
}
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,9 +51,9 @@ 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)
# define the dual operator
fom = fom.with_(dual_operator=fom.operator)
return fom, mu_bar
def record_results(function, parse, data, mu):
......
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