diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index d98ba1a51..1d9f8d66f 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -30,7 +30,7 @@ from psydac.api.fem import DiscreteBilinearForm from psydac.api.fem import DiscreteLinearForm from psydac.api.fem import DiscreteFunctional -from psydac.api.feec import DiscreteDeRham, DiscreteDeRhamMultipatch +from psydac.api.feec import DiscreteDeRham, MultipatchDiscreteDeRham from psydac.api.glt import DiscreteGltExpr from psydac.api.expr import DiscreteExpr from psydac.api.equation import DiscreteEquation @@ -199,7 +199,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): Create a discrete multipatch de Rham sequence from a symbolic one. This function creates the broken discrete spaces from the symbolic ones, and then - creates a DiscreteDeRhamMultipatch object from them. + creates a MultipatchDiscreteDeRham object from them. Parameters ---------- @@ -214,7 +214,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): Returns ------- - DiscreteDeRhamMultipatch + MultipatchDiscreteDeRham The discrete multipatch de Rham sequence containing the discrete spaces, differential operators and projectors. @@ -229,7 +229,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): spaces = [discretize_space(V, domain_h, basis=basis, **kwargs) \ for V, basis in zip(derham.spaces, bases)] - return DiscreteDeRhamMultipatch( + return MultipatchDiscreteDeRham( domain_h = domain_h, spaces = spaces ) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 2bd95ffbd..439c2f81f 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -26,9 +26,11 @@ from psydac.fem.basic import FemSpace, FemLinearOperator from psydac.fem.vector import VectorFemSpace +from psydac.fem.projectors import DirichletProjector, MultipatchDirichletProjector + from psydac.linalg.basic import IdentityOperator -__all__ = ('DiscreteDeRham', 'DiscreteDeRhamMultipatch',) +__all__ = ('DiscreteDeRham', 'MultipatchDiscreteDeRham',) #============================================================================== class DiscreteDeRham(BasicDiscrete): @@ -116,6 +118,7 @@ def __init__(self, domain_h, *spaces): self._hodge_operators = () self._conf_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property def dim(self): @@ -283,36 +286,77 @@ def projectors(self, *, kind='global', nquads=None): #-------------------------------------------------------------------------- def derivatives(self, kind='femlinop'): + assert kind in ('femlinop', 'linop') + if kind == 'femlinop': return self._derivatives elif kind == 'linop': return tuple(b_diff.linop for b_diff in self._derivatives) #-------------------------------------------------------------------------- - def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): + def dirichlet_projectors(self, kind='femlinop'): """ - return the conforming projectors of the broken multi-patch space + Returns operators that apply the correct homogeneous Dirichlet BCs. Parameters ---------- + kind : str + The kind of the projector, can be 'femlinop' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'linop' returns a psydac LinearOperator - p_moments : - The number of moments preserved by the projector. + Returns + ------- + d_projectors : tuple + Tuple of or + The Dirichlet projectors of each space and in desired form. - hom_bc: - Apply homogenous boundary conditions if True + Notes + ----- + See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form. + + """ + assert kind in ('femlinop', 'linop') + + if not self._dirichlet_proj: + d_projectors_linop = tuple(DirichletProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + d_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)) + self._dirichlet_proj = d_projectors_femlinop + + if kind == 'femlinop': + return self._dirichlet_proj + elif kind == 'linop': + return tuple(femlinop.linop for femlinop in self._dirichlet_proj) + + #-------------------------------------------------------------------------- + def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): + """ + return the conforming projectors of the broken multi-patch space + Parameters + ---------- kind : The kind of the projector, can be 'femlinop' or 'linop'. - 'femlinop' returns a psydac FemLinearOperator (default) - 'linop' returns a psydac LinearOperator + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. + + p_moments: + Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) + + hom_bc: + Apply homogenous boundary conditions if True + Returns ------- cP0, cP1, cP2 : Tuple of or The conforming projectors of each space and in desired form. """ + assert kind in ('femlinop', 'linop') if hom_bc is None: raise ValueError('please provide a value for "hom_bc" argument') @@ -442,6 +486,7 @@ def hodge_operator(self, space=None, dual=False, kind='femlinop', backend_langua H : or """ + assert kind in ('femlinop', 'linop') if not self._hodge_operators: self._init_hodge_operators(backend_language=backend_language) @@ -480,6 +525,7 @@ def hodge_operators(self, dual=False, kind='femlinop', backend_language='python' ------- The Hodge operators of all spaces and of the specified kind. """ + assert kind in ('femlinop', 'linop') if not self._hodge_operators: self._init_hodge_operators(backend_language=backend_language) @@ -488,7 +534,7 @@ def hodge_operators(self, dual=False, kind='femlinop', backend_language='python' #============================================================================== -class DiscreteDeRhamMultipatch(DiscreteDeRham): +class MultipatchDiscreteDeRham(DiscreteDeRham): """ Represents the discrete de Rham sequence for multipatch domains. It only works when the number of patches>1. @@ -538,6 +584,7 @@ def __init__(self, *, domain_h, spaces): self._hodge_operators = () self._conf_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property @@ -606,3 +653,38 @@ def projectors(self, *, kind='global', nquads=None): elif self.dim == 3: raise NotImplementedError("3D projectors are not available") + + #-------------------------------------------------------------------------- + def dirichlet_projectors(self, kind='femlinop'): + """ + Returns operators that apply the correct homogeneous Dirichlet BCs. + + Parameters + ---------- + kind : str + The kind of the projector, can be 'femlinop' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'linop' returns a psydac LinearOperator + + Returns + ------- + d_projectors : tuple + Tuple of or + The Dirichlet projectors of each space and in desired form. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form. + + """ + assert kind in ('femlinop', 'linop') + + if not self._dirichlet_proj: + d_projectors_linop = tuple(MultipatchDirichletProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + d_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)) + self._dirichlet_proj = d_projectors_femlinop + + if kind == 'femlinop': + return self._dirichlet_proj + elif kind == 'linop': + return tuple(femlinop.linop for femlinop in self._dirichlet_proj) diff --git a/psydac/feec/conforming_projectors.py b/psydac/feec/conforming_projectors.py index d045a702e..1ac1e52a8 100644 --- a/psydac/feec/conforming_projectors.py +++ b/psydac/feec/conforming_projectors.py @@ -1211,8 +1211,12 @@ def construct_h1_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=- def get_vertex_index(coords): - nbasis0 = Vh.spaces[coords[0]].nbasis - 1 - nbasis1 = Vh.spaces[coords[1]].nbasis - 1 + """ + Calculate the global index of the vertex basis function + from the geometric coordinates of a vertex in the domain. + """ + nbasis0 = Vh.spaces[0].nbasis - 1 + nbasis1 = Vh.spaces[1].nbasis - 1 # patch local index multi_index = [None] * ndim @@ -1223,6 +1227,11 @@ def get_vertex_index(coords): return l2g.get_index(0, 0, multi_index) def vertex_moment_indices(axis, coords, p_moments): + """ + Calculate the global indices of the basis functions + adjacent to the vertex basis function along axis + from the geometric coordinates of a vertex in the domain. + """ if coords[axis] == 0: return range(1, p_moments + 2) else: @@ -1233,6 +1242,9 @@ def vertex_moment_indices(axis, coords, p_moments): for co in [(0,0), (1,0), (0,1), (1,1)]: + if all(Vh.periodic): + break + # global index ig = get_vertex_index(co) @@ -1304,6 +1316,10 @@ def get_mu_minus(j, coarse_space, fine_space, R): # boundary condition for bn in domain.boundary: + + if Vh.periodic[bn.axis]: + continue + space_k = Vh axis = bn.axis @@ -1402,6 +1418,9 @@ def edge_moment_index(p, i, axis, ext): # boundary condition for bn in domain.boundary: + if Vh.periodic[bn.axis]: + continue + axis = bn.axis d = 1 - axis ext = bn.ext @@ -1431,9 +1450,13 @@ class ConformingProjectionV0(FemLinearOperator): ---------- V0h: The discrete space + + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. p_moments: Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) hom_bc : Apply homogenous boundary conditions if True @@ -1472,8 +1495,12 @@ class ConformingProjectionV1(FemLinearOperator): V1h: The discrete space + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. + p_moments: Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) hom_bc : Apply homogenous boundary conditions if True diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index e13185051..3aae76cb9 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,19 +1,24 @@ import numpy as np +from collections.abc import Iterable -from sympde.topology import element_of -from sympde.topology.space import ScalarFunction -from sympde.topology.mapping import Mapping +from scipy.sparse import diags + +from sympde.topology import element_of, Boundary from sympde.calculus import dot -from sympde.expr.expr import LinearForm, integral +from sympde.expr import LinearForm, integral, EssentialBC +from sympde.topology.datatype import SpaceType +from psydac.api.essential_bc import apply_essential_bc from psydac.api.settings import PSYDAC_BACKENDS - -from psydac.linalg.kron import KroneckerDenseMatrix from psydac.core.bsplines import hrefinement_matrix -from psydac.linalg.stencil import StencilVectorSpace from psydac.fem.basic import FemSpace +from psydac.linalg.basic import LinearOperator, Vector +from psydac.linalg.kron import KroneckerDenseMatrix +from psydac.linalg.stencil import StencilVectorSpace, StencilVector +from psydac.linalg.utilities import array_to_psydac -__all__ = ('knots_to_insert', 'knot_insertion_projection_operator') +__all__ = ('knots_to_insert', 'knot_insertion_projection_operator', 'get_dual_dofs', + 'DirichletProjector', 'MultipatchDirichletProjector') def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): """ Compute the point difference between the fine grid and coarse grid.""" @@ -158,3 +163,323 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste return tilde_f.toarray() else: return tilde_f + + +#=============================================================================== +class DirichletProjector(LinearOperator): + """ + A LinearOperator that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. + + bcs : Iterable | None + Iterable of sympde.topology.Boundary objects. + Allows the user to apply different kinds of BCs. + Must not be passed if a space_kind argument is passed. + + space_kind : str | SpaceType | None + Necessary only if fem_space.kind.name is undefined and no bcs are passed. + Must not be passed if a bcs argument is passed. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of such an operator. + + """ + def __init__(self, fem_space, *, bcs=None, space_kind=None): + + assert isinstance(fem_space, FemSpace) + assert bcs is None or isinstance(bcs, Iterable) + assert space_kind is None or isinstance(space_kind, (str, SpaceType)) + assert bcs is None or space_kind is None, \ + "Parameters bcs and space_kind are mutually exclusive" + + coeff_space = fem_space.coeff_space + self._domain = coeff_space + self._codomain = coeff_space + + if bcs is not None: + assert all(isinstance(bc, Boundary) for bc in bcs) + self._bcs = tuple(bcs) + else: + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + + #------------------------------------- + # Abstract interface + #------------------------------------- + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._domain + + @property + def dtype(self): + return None + + def tosparse(self): + np_ones = np.ones(shape=self.domain.dimension) + psy_ones = array_to_psydac(np_ones, self.domain) + diagonal = self @ psy_ones + sparse = diags(diagonal.toarray()) + return sparse + + def toarray(self): + np_ones = np.ones(shape=self.domain.dimension) + psy_ones = array_to_psydac(np_ones, self.domain) + diagonal = self @ psy_ones + array = np.diag(diagonal.toarray()) + return array + + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + if isinstance(self.domain, StencilVectorSpace): + apply_essential_bc(out, *self._bcs) + else: + for block, block_bcs in zip(out, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + + def transpose(self, conjugate=False): + return self + + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + @property + def bcs(self): + return self._bcs + + def _get_bcs(self, fem_space, *, space_kind=None): + """ + Returns a tuple of Boundaries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.kind.name determines the kind of BCs returned by this function. + + space_kind : str | SpaceType | None + Optional. Must match fem_space.kind.name if that value is different from "undefined". + Must be passed if fem_space.kind.name has value "undefined". + + Returns + ------- + bcs : tuple + tuple of sympde.topology.Boundary. + Coefficients corresponding to functions non-zero on these boundaries will be set to 0. + + """ + space = fem_space.symbolic_space + periodic = fem_space.periodic + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + kind = space_kind_str + dim = space.domain.dim + + if kind == 'l2': + return () + + u = element_of(space, name="u") + ebcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + if kind == "h1": + bcs = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim >= 2: + bcs += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + + elif kind == 'hcurl': + assert dim in (2, 3) + bcs_x = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_x += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs_y = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim == 3: + bcs_y += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + if dim == 3: + bcs_z = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_z += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + elif kind == 'hdiv': + assert dim in (2, 3) + bcs_x = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_y = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_z = [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return tuple(bcs) + +#=============================================================================== +class MultipatchDirichletProjector(DirichletProjector): + """ + A LinearOperator (for multipatch domains) that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. + + bcs : Iterable | None + Iterable of sympde.topology.Boundary objects. + Allows the user to apply different kinds of BCs. + Must not be passed if a space_kind argument is passed. + + space_kind : str | SpaceType | None + Necessary only if fem_space.kind.name is undefined and no bcs are passed. + Must not be passed if a bcs argument is passed. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of such an operator. + + """ + def __init__(self, fem_space, *, bcs=None, space_kind=None): + super().__init__(fem_space, bcs=bcs, space_kind=space_kind) + if fem_space.ldim != 2: + msg = f'The class {__class__.__name__} is implemented only in 2D.' + raise NotImplementedError(msg) + + #------------------------------------- + # Abstract interface + #------------------------------------- + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + + # apply bc on each patch + for p in out.blocks: + + if isinstance(p, StencilVector): + apply_essential_bc(p, *self._bcs) + else: + for block, block_bcs in zip(p, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + def _get_bcs(self, fem_space, *, space_kind=None): + """ + Returns a tuple of Boundaries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.kind.name determines the kind of BCs returned by this function. + + space_kind : str | SpaceType | None + Optional. Must match fem_space.kind.name if that value is different from "undefined". + Must be passed if fem_space.kind.name has value "undefined". + + Returns + ------- + bcs : tuple + tuple of sympde.topology.Boundary. + Coefficients corresponding to functions non-zero on these boundaries will be set to 0. + + """ + space = fem_space.symbolic_space + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + kind = space_kind_str + + if kind == 'l2': + return () + + u = element_of(space, name="u") + + if kind == "h1": + bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + elif kind == 'hcurl': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 0: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 1: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + elif kind == 'hdiv': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 1: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 0: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return tuple(bcs) diff --git a/psydac/fem/tests/test_dirichlet_projectors.py b/psydac/fem/tests/test_dirichlet_projectors.py new file mode 100644 index 000000000..044affb48 --- /dev/null +++ b/psydac/fem/tests/test_dirichlet_projectors.py @@ -0,0 +1,568 @@ +import numpy as np +import pytest +from mpi4py import MPI + +from sympy import sin, pi, sqrt, Tuple + +from sympde.calculus import inner, cross +from sympde.expr import integral, LinearForm, BilinearForm +from sympde.topology import elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace +from sympde.topology.datatype import H1Space, HcurlSpace + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.fem.projectors import DirichletProjector +from psydac.linalg.basic import LinearOperator, IdentityOperator +from psydac.linalg.block import BlockVectorSpace +from psydac.linalg.solvers import inverse + +class SquareTorus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + +class Annulus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)'} + + _ldim = 2 + _pdim = 2 + +class SinMapping1D(Mapping): + + _expressions = {'x': 'sin((pi/2)*x1)'} + + _ldim = 1 + _pdim = 1 + +def _test_LO_equality_using_rng(A, B): + """ + A simple tool to check with almost certainty that two linear operators are identical, + by applying them repeatedly to random vectors. + + """ + + assert isinstance(A, LinearOperator) + assert isinstance(B, LinearOperator) + assert A.domain is B.domain + assert A.codomain is B.codomain + + rng = np.random.default_rng(42) + + x = A.domain.zeros() + y1 = A.codomain.zeros() + y2 = y1.copy() + + n = 10 + + for _ in range(n): + + x *= 0. + + if isinstance(A.domain, BlockVectorSpace): + for block in x.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=x._data.shape, dtype="float64", out=x._data) + + A.dot(x, out=y1) + B.dot(x, out=y2) + + diff = y1 - y2 + err = diff.inner(diff) / diff.space.dimension**2 + tol = 1e-15 + + assert err < tol**2 + +#=============================================================================== +@pytest.mark.parametrize('dim', [1, 2]) + +def test_function_space_boundary_projector(dim): + + tol = 1e-15 + + ncells_3d = [8, 8, 8] + degree_3d = [2, 2, 2] + periodic_3d = [False, True, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + rng = np.random.default_rng(42) + + print() + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] + degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] + periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + # The function defined here satisfy the corresponding homogeneous Dirichlet BCs + if dim == 1: + x = domain.coordinates + V = ScalarFunctionSpace('V', domain, kind='H1') # testing various kind arguments + f = sin(2*pi*x) + if dim == 2: + x, y = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind=H1Space) # testing various kind arguments + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + else: + V = VectorFunctionSpace('V', domain, kind='hCuRl') # testing various kind arguments + f1 = x + f2 = y + f = Tuple(f1, f2) + if dim == 3: + x, y, z = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind='h1') # testing various kind arguments + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) + elif i == 1: + V = VectorFunctionSpace('V', domain, kind=HcurlSpace) # testing various kind arguments + f1 = z * (z - 1) * x + f2 = z * (z - 1) * y + f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f = Tuple(f1, f2, f3) + else: + V = VectorFunctionSpace('V', domain, kind='Hdiv') # testing various kind arguments + f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f3 = z * (z-1) * sin(x*y) + f = Tuple(f1, f2, f3) + + u, v = elements_of(V, names='u, v') + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + Vh = discretize(V, domain_h, degree=degree) + expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (Vh, Vh), backend=backend) + abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) + + I = IdentityOperator(Vh.coeff_space) + DP = DirichletProjector(Vh) + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc + lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f + l = LinearForm(v, integral(domain, lexpr)) + lh = discretize(l, domain_h, Vh, backend=backend) + rhs = lh.assemble() + A = M + 1e30*Mb + A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) + fc = A_inv @ rhs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err_sqr = diff.inner(diff) + print(f' || f - P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = Vh.coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DP @ rdm_coeffs + scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}') + assert scaled_boundary_int_proj_rdm_sqr < tol**2 + + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err_sqr = diff_arr.dot(diff_arr) + assert err_sqr < tol**2 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err_sqr = diff.inner(diff) + print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_sqr = M.dot_inner (fc, fc) + l2_norm2_sqr = M_0.dot_inner(fc, fc) + err_sqr = abs(l2_norm_sqr - l2_norm2_sqr) + print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to') + print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}') + # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15. + # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol. + assert err_sqr < tol + + print() + +#=============================================================================== +@pytest.mark.parametrize('dim', [1, 3]) +@pytest.mark.parallel + +def test_discrete_derham_boundary_projector(dim): + + tol = 1e-15 + + ncells = [8, 8, 8] + degree = [2, 2, 2] + periodic = [False, True, False] + + comm = MPI.COMM_WORLD + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + rng = np.random.default_rng(42) + + # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs + + f11 = lambda x : np.sin(2*np.pi*x) + + r2 = lambda x, y : np.sqrt(x**2 + y**2) + f21 = lambda x, y : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f22_1 = lambda x, y : x + f22_2 = lambda x, y : y + f22 = (f22_1, f22_2) + + f31 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) * z * (z - 1) + f32_1 = lambda x, y, z : z * (z - 1) * x + f32_2 = lambda x, y, z : z * (z - 1) * y + f32_3 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f32 = (f32_1, f32_2, f32_3) + f33_1 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_2 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_3 = lambda x, y, z : z * (z - 1) * np.sin(x*y) + f33 = (f33_1, f33_2, f33_3) + + funs = [[f11], [f21, f22], [f31, f32, f33]] + + print() + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] + degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] + periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] + + domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree_dim) + + d_projectors = derham_h.dirichlet_projectors(kind='linop') + + if dim == 2: + conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DP = d_projectors[i] + + if dim == 2: + CP = conf_projectors[i] + _test_LO_equality_using_rng(DP, CP) + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + f = funs[dim-1][i] + fc = derham_h.projectors()[i](f).coeffs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err_sqr = diff.inner(diff) + print(f' || f - P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DP @ rdm_coeffs + scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}') + assert scaled_boundary_int_proj_rdm_sqr < tol**2 + + # 2.2 + # Test tosparse(): (DP @ rdm_coeffs).toarray() should be equal to DP.tosparse().dot(rdm_coeffs.toarray()) + DP_spr = DP.tosparse() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_spr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err_sqr = diff_arr.dot(diff_arr) + assert err_sqr < tol**2 + + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err_sqr = diff.inner(diff) + print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_sqr = M.dot_inner (fc, fc) + l2_norm2_sqr = M_0.dot_inner(fc, fc) + err_sqr = abs(l2_norm_sqr - l2_norm2_sqr) + print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to') + print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}') + # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15. + # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol. + assert err_sqr < tol + + print() + +#=============================================================================== +def test_discrete_derham_boundary_projector_multipatch(): + + tol = 1e-15 + + ncells = [8, 8] + degree = [2, 2] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain + domain = build_multipatch_domain(domain_name='annulus_3') + + rng = np.random.default_rng(42) + + # The following are functions satisfying homogeneous Dirichlet BCs + r = lambda x, y : np.sqrt(x**2 + y**2) + f1 = lambda x, y : (r(x, y) - 0.5) * (r(x, y) - 1) + f2_1 = lambda x, y : x + f2_2 = lambda x, y : y + f2 = (f2_1, f2_2) + funs = [f1, f2] + print() + + boundary = domain.boundary + + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_h = {} + for D in domain.interior: + ncells_h[D.name] = ncells + + domain_h = discretize(domain, ncells=ncells_h, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree) + + projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) + + d_projectors = derham_h.dirichlet_projectors(kind='linop') + + nn = NormalVector('nn') + + for i in range(2): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + expr = u*v + if (i == 1): + boundary_expr = cross(nn, u) * cross(nn, v) + expr = inner(u,v) + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DP = d_projectors[i] + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + f = funs[i] + fc = projectors[i](f).coeffs + + # 1. + # The coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err_sqr = diff.inner(diff) + print(f' || f - P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + for patch in rdm_coeffs.blocks: + + if isinstance(patch.space, BlockVectorSpace): + for block in patch.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=patch._data.shape, dtype="float64", out=patch._data) + + rdm_coeffs2 = DP @ rdm_coeffs + scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}') + assert scaled_boundary_int_proj_rdm_sqr < tol**2 + + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err_sqr = diff_arr.dot(diff_arr) + assert err_sqr < tol**2 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err_sqr = diff.inner(diff) + print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}') + assert err_sqr < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_sqr = M.dot_inner (fc, fc) + l2_norm2_sqr = M_0.dot_inner(fc, fc) + err_sqr = abs(l2_norm_sqr - l2_norm2_sqr) + print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to') + print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}') + # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15. + # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol. + assert err_sqr < tol + + print() + + +# =============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== + +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 6d4b37e73..da7305d0d 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -1,11 +1,9 @@ +import numpy as np +import pytest -import numpy as np -import pytest -from psydac.linalg.solvers import inverse -from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector -from psydac.linalg.basic import LinearSolver -from psydac.ddm.cart import DomainDecomposition, CartDecomposition - +from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from psydac.linalg.solvers import inverse +from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector def define_data_hermitian(n, p, dtype=float): domain_decomposition = DomainDecomposition([n - p], [False]) @@ -54,9 +52,8 @@ def define_data(n, p, matrix_data, dtype=float): xe[s:e + 1] = np.random.random(e + 1 - s) return(V, A, xe) - #=============================================================================== -@pytest.mark.parametrize( 'n', [5, 10, 13] ) +@pytest.mark.parametrize('n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @pytest.mark.parametrize('dtype', [float, complex]) @pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) @@ -204,6 +201,7 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): assert errh_norm < tol assert solver == 'pcg' or errc_norm < tol + # =============================================================================== # SCRIPT FUNCTIONALITY #===============================================================================