Coverage for src/pymor/operators/ei : 16%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
# -*- coding: utf-8 -*- # This file is part of the pyMOR project (http://www.pymor.org). # Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler # License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
'''Interpolate an |Operator| using Empirical Operator Interpolation.
Let `L` be an |Operator|, `0 <= c_1, ..., c_M <= L.dim_range` indices of interpolation DOFs and `b_1, ..., b_M in R^(L.dim_range)` collateral basis vectors. If moreover `ψ_j(U)` denotes the j-th component of `U`, the empirical interpolation `L_EI` of `L` w.r.t. the given data is given by ::
| M | L_EI(U, μ) = ∑ b_i⋅λ_i such that | i=1 | | ψ_(c_i)(L_EI(U, μ)) = ψ_(c_i)(L(U, μ)) for i=0,...,M
It is assumed that ::
ψ_(c_i)(b_j) = 0 for i < j
which means that the matrix of the interpolation problem is triangular.
Since the original operator only has to be evaluated at the given interpolation DOFs, |EmpiricalInterpolatedOperator| calls `operator.restricted(interpolation_dofs)` to obtain a restricted version of the operator which is stored and later used to quickly obtain the required evaluations. (The second return value of the `restricted` method has to be an array of source DOFs -- determined by the operator's stencil -- required to evaluate the restricted operator.)
The interpolation DOFs and the collateral basis can be generated using the algorithms provided in the :mod:`pymor.algorithms.ei` module.
Parameters ---------- operator The |Operator| to interpolate. The operator must implement a `restricted` method as described above. interpolation_dofs List or 1D |NumPy array| of the interpolation DOFs `c_1, ..., c_M`. collateral_basis |VectorArray| containing the collateral basis `b_1, ..., b_M`. name Name of the operator. '''
assert isinstance(operator, OperatorInterface) assert isinstance(collateral_basis, VectorArrayInterface) assert operator.dim_range == collateral_basis.dim assert operator.type_range == type(collateral_basis) assert hasattr(operator, 'restricted')
self.build_parameter_type(inherits=(operator,)) self.dim_source = operator.dim_source self.dim_range = operator.dim_range self.type_source = operator.type_source self.type_range = operator.type_range self.linear = operator.linear self.name = name or '{}_interpolated'.format(operator.name)
interpolation_dofs = np.array(interpolation_dofs, dtype=np.int32) self.interpolation_dofs = interpolation_dofs
if len(interpolation_dofs) > 0: restricted_operator, source_dofs = operator.restricted(interpolation_dofs) interpolation_matrix = collateral_basis.components(interpolation_dofs).T self.restricted_operator = restricted_operator self.source_dofs = source_dofs self.interpolation_matrix = interpolation_matrix self.collateral_basis = collateral_basis.copy()
mu = self.parse_parameter(mu) if len(self.interpolation_dofs) == 0: count = len(ind) if ind is not None else len(U) return self.type_range.zeros(dim=self.dim_range, count=count)
U_components = NumpyVectorArray(U.components(self.source_dofs, ind=ind), copy=False) AU = self.restricted_operator.apply(U_components, mu=mu) try: interpolation_coefficients = solve_triangular(self.interpolation_matrix, AU.data.T, lower=True, unit_diagonal=True).T except ValueError: # this exception occurs when AU contains NaNs ... interpolation_coefficients = np.empty((len(AU), len(self.collateral_basis))) + np.nan # interpolation_coefficients = np.linalg.solve(self.interpolation_matrix, AU._array.T).T return self.collateral_basis.lincomb(interpolation_coefficients)
assert source_basis is not None or self.dim_source == 0 assert range_basis is not None assert source_basis is None or self.dim_source == source_basis.dim assert self.dim_range == range_basis.dim assert source_basis is None or self.type_source == type(source_basis) assert self.type_range == type(range_basis)
if product is None: projected_collateral_basis = NumpyVectorArray(self.collateral_basis.dot(range_basis, pairwise=False)) else: projected_collateral_basis = NumpyVectorArray(product.apply2(self.collateral_basis, range_basis, pairwise=False))
return ProjectedEmpiciralInterpolatedOperator(self.restricted_operator, self.interpolation_matrix, NumpyVectorArray(source_basis.components(self.source_dofs), copy=False), projected_collateral_basis, name)
'''Project an |EmpiricalInterpolatedOperator|.
Not intended to be used directly. Instead use :meth:`~pymor.operators.interfaces.OperatorInterface.projected`. '''
projected_collateral_basis, name=None): self.dim_source = len(source_basis_dofs) self.dim_range = projected_collateral_basis.dim self.type_source = self.type_range = NumpyVectorArray self.linear = restricted_operator.linear self.build_parameter_type(inherits=(restricted_operator,)) self.restricted_operator = restricted_operator self.interpolation_matrix = interpolation_matrix self.source_basis_dofs = source_basis_dofs self.projected_collateral_basis = projected_collateral_basis self.name = name or '{}_projected'.format(restricted_operator.name)
mu = self.parse_parameter(mu) U_array = U._array if ind is None else U._array[ind] U_components = self.source_basis_dofs.lincomb(U_array) AU = self.restricted_operator.apply(U_components, mu=mu) try: interpolation_coefficients = solve_triangular(self.interpolation_matrix, AU.data.T, lower=True, unit_diagonal=True).T except ValueError: # this exception occurs when AU contains NaNs ... interpolation_coefficients = np.empty((len(AU), len(self.projected_collateral_basis))) + np.nan return self.projected_collateral_basis.lincomb(interpolation_coefficients)
assert dim_source is None or dim_source <= self.dim_source assert dim_range is None or dim_range <= self.dim_range assert dim_collateral is None or dim_collateral <= self.restricted_operator.dim_range name = name or '{}_projected_to_subbasis'.format(self.name)
interpolation_matrix = self.interpolation_matrix[:dim_collateral, :dim_collateral]
if dim_collateral is not None: restricted_operator, source_dofs = self.restricted_operator.restricted(np.arange(dim_collateral)) else: restricted_operator = self.restricted_operator
old_pcb = self.projected_collateral_basis projected_collateral_basis = NumpyVectorArray(old_pcb.data[:dim_collateral, :dim_range], copy=False)
old_sbd = self.source_basis_dofs source_basis_dofs = NumpyVectorArray(old_sbd.data[:dim_source], copy=False) if dim_collateral is None \ else NumpyVectorArray(old_sbd.data[:dim_source, source_dofs], copy=False)
return ProjectedEmpiciralInterpolatedOperator(restricted_operator, interpolation_matrix, source_basis_dofs, projected_collateral_basis, name=name) |