Unverified Commit b8937c65 authored by Felix Schindler's avatar Felix Schindler Committed by GitHub
Browse files

Merge pull request #1383 from pymor/dunegdt-bindings-only

bindings for dune-gdt (part 1: bindings)
parents 34b2c21b 3ecee16e
Pipeline #106338 passed with stages
in 36 minutes and 40 seconds
......@@ -8,6 +8,7 @@ import scipy.sparse as sps
import scipy.sparse.linalg as spsla
from pymor.algorithms.rules import RuleTable, match_class
from pymor.core.config import config
from pymor.operators.block import BlockOperatorBase
from pymor.operators.constructions import (AdjointOperator, ComponentProjectionOperator, ConcatenationOperator,
IdentityOperator, LincombOperator, LowRankOperator, LowRankUpdatedOperator,
......@@ -172,3 +173,23 @@ class ToMatrixRules(RuleTable):
return np.zeros((op.range.dim, op.source.dim))
else:
return getattr(sps, format + '_matrix')((op.range.dim, op.source.dim))
if config.HAVE_DUNEGDT:
from scipy.sparse import csr_matrix
from pymor.bindings.dunegdt import DuneXTMatrixOperator
@match_class(DuneXTMatrixOperator)
def action_DuneXTMatrixOperator(self, op):
data, row_ind, col_ind = op.matrix.to_csr()
mat = csr_matrix((data, (row_ind, col_ind)), shape=(op.matrix.rows, op.matrix.cols))
format = self.format
if format is None:
return mat
elif format == 'dense':
return mat.toarray()
else:
return mat.asformat(format)
ToMatrixRules.insert_rule(0, action_DuneXTMatrixOperator)
# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2016 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from pymor.core.config import config
if config.HAVE_DUNEGDT:
import numpy as np
from dune.xt.la import IstlVector
from pymor.operators.list import LinearComplexifiedListVectorArrayOperatorBase
from pymor.vectorarrays.interface import _create_random_values
from pymor.vectorarrays.list import (
ComplexifiedListVectorSpace, ComplexifiedVector, CopyOnWriteVector, NumpyVector)
class DuneXTVector(CopyOnWriteVector):
"""Wraps a vector from dune-xt to make it usable with ListVectorArray.
Parameters
----------
impl
The actual vector from dune.xt.la, usually IstlVector.
"""
def __init__(self, impl):
self.impl = impl
@classmethod
def from_instance(cls, instance):
return cls(instance.impl)
def _copy_data(self):
self.impl = self.impl.copy(True)
def _scal(self, alpha):
self.impl.scal(alpha)
def _axpy(self, alpha, x):
self.impl.axpy(alpha, x.impl)
def inner(self, other):
return self.impl.dot(other.impl)
def norm(self):
return self.impl.l2_norm()
def norm2(self):
return self.impl.l2_norm() ** 2
def sup_norm(self):
return self.impl.sup_norm()
def dofs(self, dof_indices):
impl = self.impl
return np.array([impl[i] for i in dof_indices])
def amax(self):
_amax = self.impl.amax()
return _amax[0], _amax[1]
def __add__(self, other):
return DuneXTVector(self.impl + other.impl)
def __iadd__(self, other):
self.impl += other.impl
return self
__radd__ = __add__
def __sub__(self, other):
return DuneXTVector(self.impl - other.impl)
def __isub__(self, other):
self.impl -= other.impl
return self
def __mul__(self, other):
return DuneXTVector(self.impl * other)
def __imul__(self, other):
self.impl *= other
return self
def __neg__(self):
return self * (-1)
def to_numpy(self, ensure_copy=False):
return np.array(self.impl, copy=ensure_copy)
class ComplexifiedDuneXTVector(ComplexifiedVector):
"""Required for DuneXTVectorSpace, Usually not to be used directly."""
def amax(self):
if self.imag_part is None:
return self.real_part.amax()
else:
real = np.array(self.real_part.impl, copy=False)
imag = np.array(self.imag_part.impl, copy=False)
return NumpyVector(real + imag * 1j).amax()
class DuneXTVectorSpace(ComplexifiedListVectorSpace):
"""A |VectorSpace| yielding DuneXTVector
Parameters
----------
dim
Dimension of the |VectorSpace|, i.e., length of the resulting vectors.
vector_type
Type of the actual vector from dune.xt.la, usually IstlVector.
id
Identifier of the |VectorSpace|.
"""
real_vector_type = DuneXTVector
vector_type = ComplexifiedDuneXTVector
def __init__(self, dim, dune_vector_type=IstlVector, id='STATE'):
self.__auto_init(locals())
def __eq__(self, other):
return type(other) is DuneXTVectorSpace \
and self.dune_vector_type == other.dune_vector_type \
and self.dim == other.dim \
and self.id == other.id
# since we implement __eq__, we also need to implement __hash__
def __hash__(self):
return id(self.dune_vector_type) + hash(self.dim)
def real_zero_vector(self):
return DuneXTVector(self.dune_vector_type(self.dim, 0.))
def real_full_vector(self, value):
return DuneXTVector(self.dune_vector_type(self.dim, value))
def real_random_vector(self, distribution, random_state, **kwargs):
values = _create_random_values(self.dim, distribution, random_state, **kwargs)
return self.real_vector_from_numpy(values)
def real_vector_from_numpy(self, data, ensure_copy=False):
v = self.real_zero_vector()
np_view = np.array(v.impl, copy=False)
np_view[:] = data
return v
def real_make_vector(self, obj):
return DuneXTVector(obj)
class DuneXTMatrixOperator(LinearComplexifiedListVectorArrayOperatorBase):
"""Wraps a dune-xt matrix as an |Operator|.
Parameters
----------
matrix
The actual matrix from dune.xt.la, usually IstlMatrix.
source_id
Identifier of the source |VectorSpace|.
range_id
Identifier of the source |VectorSpace|.
solver_options
If specified, either a string or a dict specifying the solver used in apply_inverse. See
https://zivgitlab.uni-muenster.de/ag-ohlberger/dune-community/dune-xt/-/tree/master/dune/xt/la/solver
for available options, depending on the type of `matrix`. E.g., for
dune.xt.la.IstlSparseMatrix, (as can be queried from dune.xt.la.IstlSparseMatrixSolver
via `types()` and `options(type)`):
- 'bicgstab.ssor'
- 'bicgstab.amg.ssor'
- 'bicgstab.amg.ilu0'
- 'bicgstab.ilut'
- 'bicgstab'
- 'cg'
name
Optional name of the resulting |Operator|.
"""
linear = True
def __init__(self, matrix, source_id='STATE', range_id='STATE', solver_options=None, name=None):
self.source = DuneXTVectorSpace(matrix.cols, matrix.vector_type(), source_id)
self.range = DuneXTVectorSpace(matrix.rows, matrix.vector_type(), range_id)
self.__auto_init(locals())
def _real_apply_one_vector(self, u, mu=None, prepare_data=None):
r = self.range.real_zero_vector()
self.matrix.mv(u.impl, r.impl)
return r
def _apply_adjoint_one_vector(self, v, mu=None, prepare_data=None):
r = self.source.real_zero_vector()
self.matrix.mtv(v.impl, r.impl)
return r
def _real_apply_inverse_one_vector(self, v, mu=None, initial_guess=None,
least_squares=False, prepare_data=None):
if least_squares:
raise NotImplementedError
r = (self.source.real_zero_vector() if initial_guess is None else
initial_guess.copy(deep=True))
options = self.solver_options.get('inverse') if self.solver_options else None
from dune.xt.la import make_solver
solver = make_solver(self.matrix)
if options:
solver.apply(v.impl, r.impl, options)
else:
solver.apply(v.impl, r.impl)
return r
def _assemble_lincomb(self, operators, coefficients, identity_shift=0., solver_options=None, name=None):
if not all(isinstance(op, DuneXTMatrixOperator) for op in operators):
return None
if identity_shift != 0:
return None
if np.iscomplexobj(coefficients):
return None
if coefficients[0] == 1:
matrix = operators[0].matrix.copy()
else:
matrix = operators[0].matrix * coefficients[0]
for op, c in zip(operators[1:], coefficients[1:]):
matrix.axpy(c, op.matrix) # TODO: Not guaranteed to work for all backends! For different
# sparsity patterns one would have to extract the patterns from the pruned
# matrices, merge them and create a new matrix.
return DuneXTMatrixOperator(matrix, self.source.id, self.range.id, solver_options=solver_options, name=name)
......@@ -43,15 +43,19 @@ def _get_dunegdt_version():
try:
version = dune.gdt.__version__
if parse(version) < parse('2021.1.2') or parse(version) >= parse('2021.2'):
warnings.warn(f'dune-gdt bindings have been tested for version 2021.1.2 (installed: {dune.gdt.__version__}).')
warnings.warn('dune-gdt bindings have been tested for version 2021.1.x (x >= 2) '
f'(installed: {dune.gdt.__version__}).')
except AttributeError:
warnings.warn(f'dune-gdt bindings have been tested for version 2021.1.2 (installed: unknown older than 2021.1.2).')
warnings.warn('dune-gdt bindings have been tested for version 2021.1.x (x >= 2) '
'(installed: unknown older than 2021.1.2).')
try:
xt_version = dune.xt.__version__
if parse(xt_version) < parse('2021.1.2') or parse(xt_version) >= parse('2021.2'):
warnings.warn(f'dune-gdt bindings have been tested for dune-xt 2021.1.2 (installed: {dune.xt.__version__}).')
warnings.warn('dune-gdt bindings have been tested for dune-xt 2021.1.x (x >= 2) '
f'(installed: {dune.xt.__version__}).')
except AttributeError:
warnings.warn(f'dune-gdt bindings have been tested for dune-xt version 2021.1.2 (installed: unknown older than 2021.1.2).')
warnings.warn('dune-gdt bindings have been tested for dune-xt version 2021.1.x (x >= 2) '
'(installed: unknown older than 2021.1.2).')
return version
......
......@@ -353,9 +353,13 @@ class ListVectorArray(VectorArray):
The associated |VectorSpace| is a subclass of
:class:`ListVectorSpace`.
For an example, see :class:`NumpyVector`, :class:`NumpyListVectorSpace`
or :class:`~pymor.bindings.fenics.FenicsVector`,
:class:`~pymor.bindings.fenics.FenicsVectorSpace`.
For an example, see :class:`NumpyVector` and :class:`NumpyListVectorSpace`,
:class:`~pymor.bindings.fenics.FenicsVector` and
:class:`~pymor.bindings.fenics.FenicsVectorSpace`,
:class:`~pymor.bindings.dunegdt.DuneXTVector` and
:class:`~pymor.bindings.dunegdt.DuneXTVectorSpace`,
:class:`~pymor.bindings.ngsolve.NGSolveVector` and
:class:`~pymor.bindings.ngsolve.NGSolveVectorSpace`.
"""
_NONE = ()
......
......@@ -9,6 +9,7 @@ from pymor.algorithms.basic import almost_equal
from pymor.algorithms.projection import project
from pymor.algorithms.to_matrix import to_matrix
from pymor.core.exceptions import InversionError, LinAlgError
from pymor.core.config import config
from pymor.operators.block import BlockDiagonalOperator
from pymor.operators.constructions import (SelectionOperator, InverseOperator, InverseAdjointOperator, IdentityOperator,
LincombOperator, VectorArrayOperator)
......@@ -478,3 +479,30 @@ def test_issue_1276():
v = B.source.ones()
B.apply_inverse(v)
if config.HAVE_DUNEGDT:
from dune.xt.la import IstlSparseMatrix, SparsityPatternDefault
from pymor.bindings.dunegdt import DuneXTMatrixOperator
def make_dunegdt_identity(N):
pattern = SparsityPatternDefault(N)
for n in range(N):
pattern.insert(n, n)
pattern.sort()
mat = IstlSparseMatrix(N, N, pattern)
for n in range(N):
mat.set_entry(n, n, 1.)
return DuneXTMatrixOperator(mat)
def test_dunegdt_identiy_apply():
op = make_dunegdt_identity(4)
U = op.source.ones(1)
V = op.apply(U)
assert (U - V).sup_norm() < 1e-14
def test_dunegdt_identiy_apply_inverse():
op = make_dunegdt_identity(4)
V = op.source.ones(1)
U = op.apply_inverse(V)
assert (U - V).sup_norm() < 1e-14
......@@ -7,6 +7,7 @@ from scipy.sparse import issparse
from types import FunctionType, MethodType
from pymor.core.base import BasicObject
from pymor.core.config import config
from pymor.core.pickle import dumps, loads, dumps_function, PicklingError
from pymor.discretizers.builtin.grids.subgrid import SubGrid
from pymor.operators.numpy import NumpyMatrixBasedOperator
......@@ -20,6 +21,21 @@ is_equal_ignored_attributes = \
is_equal_dispatch_table = {}
if config.HAVE_DUNEGDT:
from dune.xt.la import IstlVector
def _assert_IstlVector_equal(first, second):
if not isinstance(first, IstlVector):
return False
if not isinstance(second, IstlVector):
return False
if not len(first) == len(second):
return False
assert_is_equal(np.array(first, copy=False), np.array(second, copy=False))
return True
is_equal_dispatch_table[IstlVector] = _assert_IstlVector_equal
def func_with_closure_generator():
x = 42
......
......@@ -23,6 +23,9 @@ if config.HAVE_FENICS:
if config.HAVE_DEALII:
from pydealii.pymor.vectorarray import DealIIVectorSpace
if config.HAVE_DUNEGDT:
from pymor.bindings.dunegdt import DuneXTVectorSpace
if config.HAVE_NGSOLVE:
import ngsolve as ngs
import netgen.meshing as ngmsh
......@@ -141,6 +144,11 @@ if config.HAVE_DEALII:
return [(DealIIVectorSpace(d), ar) for d, ar in zip(dims, np_data_list)]
_other_vector_space_types.append('dealii')
if config.HAVE_DUNEGDT:
def _dunegdt_vector_spaces(draw, np_data_list, compatible, count, dims):
return [(DuneXTVectorSpace(d), ar) for d, ar in zip(dims, np_data_list)]
_other_vector_space_types.append('dunegdt')
_picklable_vector_space_types = ['numpy', 'numpy_list', 'block']
......
......@@ -7,6 +7,7 @@ import scipy.linalg as spla
import scipy.sparse as sps
from pymor.algorithms.to_matrix import to_matrix
from pymor.core.config import config
from pymor.operators.block import BlockOperator, BlockDiagonalOperator
from pymor.operators.constructions import (AdjointOperator, ComponentProjectionOperator, IdentityOperator,
LowRankOperator, LowRankUpdatedOperator, VectorArrayOperator, ZeroOperator)
......@@ -255,3 +256,25 @@ def test_to_matrix_ZeroOperator():
Zop = ZeroOperator(NumpyVectorSpace(n), NumpyVectorSpace(m))
assert_type_and_allclose(Z, Zop, 'sparse')
if config.HAVE_DUNEGDT:
from dune.xt.la import IstlSparseMatrix, SparsityPatternDefault
from pymor.bindings.dunegdt import DuneXTMatrixOperator
def test_to_matrix_DuneXTMatrixOperator():
np.random.seed(0)
A = np.random.randn(2, 2)
pattern = SparsityPatternDefault(2)
for ii in range(2):
for jj in range(2):
pattern.insert(ii, jj)
pattern.sort()
mat = IstlSparseMatrix(2, 2, pattern)
for ii in range(2):
for jj in range(2):
mat.set_entry(ii, jj, A[ii][jj])
Aop = DuneXTMatrixOperator(mat)
assert_type_and_allclose(A, Aop, 'sparse')
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