Source code for pymor.models.basic

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 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 TimeStepper
from pymor.models.interface import Model
from pymor.operators.block import BlockOperatorBase
from pymor.operators.constructions import IdentityOperator, VectorOperator, ZeroOperator
from pymor.vectorarrays.interface import VectorArray
from pymor.vectorarrays.numpy import NumpyVectorSpace


[docs]class StationaryModel(Model): """Generic class for models of stationary problems. This class describes discrete problems given by the equation:: L(u(μ), μ) = F(μ) with a vector-like right-hand side F and a (possibly non-linear) operator L. Note that even when solving a variational formulation where F is a functional and not a vector, F has to be specified as a vector-like |Operator| (mapping scalars to vectors). This ensures that in the complex case both L and F are anti-linear in the test variable. Parameters ---------- operator The |Operator| L. rhs The vector F. Either a |VectorArray| of length 1 or a vector-like |Operator|. 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. products A dict of inner product |Operators| defined on the discrete space the problem is posed on. For each product with key `'x'` a corresponding attribute `x_product`, as well as a norm method `x_norm` is added to the model. error_estimator An error estimator for the problem. This can be any object with 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)`. 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, operator, rhs, output_functional=None, products=None, error_estimator=None, visualizer=None, name=None): if isinstance(rhs, VectorArray): assert rhs in operator.range rhs = VectorOperator(rhs, name='rhs') assert rhs.range == operator.range and rhs.source.is_scalar and rhs.linear assert output_functional is None or output_functional.source == operator.source super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name) self.__auto_init(locals()) self.solution_space = operator.source self.linear = operator.linear and (output_functional is None or output_functional.linear) if output_functional is not None: self.dim_output = output_functional.range.dim
[docs] def __str__(self): return ( f'{self.name}\n' f' class: {self.__class__.__name__}\n' f' {"linear" if self.linear else "non-linear"}\n' f' solution_space: {self.solution_space}\n' f' dim_output: {self.dim_output}\n' )
[docs] def _compute_solution(self, mu=None, **kwargs): return self.operator.apply_inverse(self.rhs.as_range_array(mu), mu=mu)
[docs] def _compute_solution_d_mu_single_direction(self, parameter, index, solution, mu): lhs_d_mu = self.operator.d_mu(parameter, index).apply(solution, mu=mu) rhs_d_mu = self.rhs.d_mu(parameter, index).as_range_array(mu) rhs = rhs_d_mu - lhs_d_mu return self.operator.jacobian(solution, mu=mu).apply_inverse(rhs)
_compute_allowed_kwargs = frozenset({'use_adjoint'})
[docs] def _compute_output_d_mu(self, solution, mu, return_array=False, use_adjoint=None): """compute the gradient of the output functional w.r.t. the parameters Parameters ---------- solution Internal model state for the given |Parameter value| mu |Parameter value| for which to compute the gradient return_array if `True`, return the output gradient as a |NumPy array|. Otherwise, return a dict of gradients for each |Parameter|. use_adjoint if `None` use standard approach, if `True`, use the adjoint solution for a more efficient way of computing the gradient. See Section 1.6.2 in [HPUU09]_ for more details. So far, the adjoint approach is only valid for linear models. Returns ------- The gradient as a |NumPy array| or a dict of |NumPy arrays|. """ if use_adjoint is None: use_adjoint = True if (self.output_functional.linear and self.operator.linear) else False if not use_adjoint: return super()._compute_output_d_mu(solution, mu, return_array) else: assert self.output_functional is not None assert self.operator.linear assert self.output_functional.linear dual_solutions = self.operator.range.empty() for d in range(self.output_functional.range.dim): dual_problem = self.with_(operator=self.operator.H, rhs=self.output_functional.H.as_range_array(mu)[d]) dual_solutions.append(dual_problem.solve(mu)) gradients = [] if return_array else {} for (parameter, size) in self.parameters.items(): array = np.empty(shape=(size,self.output_functional.range.dim)) for index in range(size): output_partial_dmu = self.output_functional.d_mu(parameter, index).apply(solution, mu=mu).to_numpy()[0] lhs_d_mu = self.operator.d_mu(parameter, index).apply2(dual_solutions, solution, mu=mu)[:,0] rhs_d_mu = self.rhs.d_mu(parameter, index).apply_adjoint(dual_solutions, mu=mu).to_numpy()[:,0] array[index] = output_partial_dmu + rhs_d_mu - lhs_d_mu if return_array: gradients.extend(array) else: gradients[parameter] = array if return_array: return np.array(gradients) else: return gradients
[docs]class InstationaryModel(Model): """Generic class for models of instationary problems. This class describes instationary problems given by the equations:: M * ∂_t u(t, μ) + L(u(μ), t, μ) = F(t, μ) u(0, μ) = u_0(μ) for t in [0,T], where L is a (possibly non-linear) time-dependent |Operator|, F is a time-dependent vector-like |Operator|, and u_0 the initial data. The mass |Operator| M is assumed to be linear. Parameters ---------- T The final time T. initial_data The initial data `u_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 given |parameter values|. operator The |Operator| L. rhs The right-hand side F. mass The mass |Operator| `M`. If `None`, the identity is assumed. time_stepper The :class:`time-stepper <pymor.algorithms.timestepping.TimeStepper>` to be used by :meth:`~pymor.models.interface.Model.solve`. 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. products A dict of product |Operators| defined on the discrete space the problem is posed on. For each product with key `'x'` a corresponding attribute `x_product`, as well as a norm method `x_norm` is added to the model. error_estimator An error estimator for the problem. This can be any object with 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)`. 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, operator, rhs, mass=None, time_stepper=None, num_values=None, output_functional=None, products=None, error_estimator=None, visualizer=None, name=None): if isinstance(rhs, VectorArray): assert rhs in operator.range rhs = VectorOperator(rhs, name='rhs') if isinstance(initial_data, VectorArray): assert initial_data in operator.source initial_data = VectorOperator(initial_data, name='initial_data') mass = mass or IdentityOperator(operator.source) rhs = rhs or ZeroOperator(operator.source, NumpyVectorSpace(1)) assert isinstance(time_stepper, TimeStepper) assert initial_data.source.is_scalar assert operator.source == initial_data.range assert rhs.linear and rhs.range == operator.range and rhs.source.is_scalar assert mass.linear and mass.source == mass.range == operator.source assert output_functional is None or output_functional.source == operator.source super().__init__(products=products, error_estimator=error_estimator, visualizer=visualizer, name=name) self.parameters_internal = {'t': 1} self.__auto_init(locals()) self.solution_space = operator.source self.linear = operator.linear and (output_functional is None or output_functional.linear) if output_functional is not None: self.dim_output = output_functional.range.dim
[docs] def __str__(self): return ( f'{self.name}\n' f' class: {self.__class__.__name__}\n' f' {"linear" if self.linear else "non-linear"}\n' f' T: {self.T}\n' f' solution_space: {self.solution_space}\n' f' dim_output: {self.dim_output}\n' )
def with_time_stepper(self, **kwargs): return self.with_(time_stepper=self.time_stepper.with_(**kwargs))
[docs] def _compute_solution(self, mu=None, **kwargs): mu = mu.with_(t=0.) U0 = self.initial_data.as_range_array(mu) U = self.time_stepper.solve(operator=self.operator, rhs=None if isinstance(self.rhs, ZeroOperator) else self.rhs, initial_data=U0, mass=None if isinstance(self.mass, IdentityOperator) else self.mass, initial_time=0, end_time=self.T, mu=mu, num_values=self.num_values) return U
[docs] def to_lti(self): """Convert model to |LTIModel|. This method interprets the given model as an |LTIModel| in the following way:: - self.operator -> A self.rhs -> B self.output_functional -> C None -> D self.mass -> E """ if self.output_functional is None: raise ValueError('No output defined.') A = - self.operator B = self.rhs C = self.output_functional E = self.mass if not all(op.linear for op in [A, B, C, E]): raise ValueError('Operators not linear.') from pymor.models.iosys import LTIModel return LTIModel(A, B, C, E=E, visualizer=self.visualizer)