diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 69dc89679..a0a6f6325 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -29,7 +29,7 @@ from psydac.api.fem import DiscreteLinearForm from psydac.api.fem import DiscreteFunctional from psydac.api.fem import DiscreteSumForm -from psydac.api.feec import DiscreteDerham +from psydac.api.feec import DiscreteDerham, DiscreteDerhamMultipatch from psydac.api.glt import DiscreteGltExpr from psydac.api.expr import DiscreteExpr from psydac.api.equation import DiscreteEquation @@ -47,6 +47,7 @@ __all__ = ( 'discretize', 'discretize_derham', + 'discretize_derham_multipatch', 'reduce_space_degrees', 'discretize_space', 'discretize_domain' @@ -175,11 +176,9 @@ def discretize_derham(derham, domain_h, *, get_H1vec_space=False, **kwargs): See Also -------- discretize_space - """ ldim = derham.shape - mapping = domain_h.domain.mapping # NOTE: assuming single-patch domain! bases = ['B'] + ldim * ['M'] spaces = [discretize_space(V, domain_h, basis=basis, **kwargs) for V, basis in zip(derham.spaces, bases)] @@ -192,7 +191,48 @@ def discretize_derham(derham, domain_h, *, get_H1vec_space=False, **kwargs): #We still need to specify the symbolic space because of "_recursive_element_of" not implemented in sympde spaces.append(Xh) - return DiscreteDerham(mapping, *spaces) + return DiscreteDerham(domain_h, *spaces) + +#============================================================================== +def discretize_derham_multipatch(derham, domain_h, *args, **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. + + Parameters + ---------- + derham : sympde.topology.space.Derham + The symbolic Derham sequence. + + domain_h : Geometry + Discrete domain where the spaces will be discretized. + + **kwargs : dict + Optional parameters for the space discretization. + + Returns + ------- + DiscreteDerhamMultipatch + The discrete multipatch De Rham sequence containing the discrete spaces, + differential operators and projectors. + + See Also + -------- + discretize_derham + discretize_space + """ + + ldim = derham.shape + bases = ['B'] + ldim * ['M'] + spaces = [discretize_space(V, domain_h, *args, basis=basis, **kwargs) \ + for V, basis in zip(derham.spaces, bases)] + + return DiscreteDerhamMultipatch( + domain_h = domain_h, + spaces = spaces + ) #============================================================================== def reduce_space_degrees(V, Vh, *, basis='B', sequence='DR'): @@ -329,7 +369,7 @@ def reduce_space_degrees(V, Vh, *, basis='B', sequence='DR'): #============================================================================== # TODO knots -def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, basis='B', sequence='DR'): +def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, basis='B', sequence='DR', pads=None): """ This function creates the discretized space starting from the symbolic space. @@ -457,6 +497,9 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, min_coords = interior.min_coords max_coords = interior.max_coords + if pads == None: + pads = [None] * len(degree_i) + assert len(ncells) == len(periodic) == len(degree_i) == len(multiplicity_i) == len(min_coords) == len(max_coords) if knots is None: @@ -465,10 +508,10 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] # Create 1D finite element spaces and precompute quadrature data - spaces[i] = [SplineSpace( p, multiplicity=m, grid=grid , periodic=P) for p,m,grid,P in zip(degree_i, multiplicity_i,grids, periodic)] + spaces[i] = [SplineSpace( p, multiplicity=m, grid=grid , periodic=P, pads=pd) for p,m,grid,P,pd in zip(degree_i, multiplicity_i,grids, periodic, pads)] else: # Create 1D finite element spaces and precompute quadrature data - spaces[i] = [SplineSpace( p, knots=T , periodic=P) for p,T, P in zip(degree_i, knots[interior.name], periodic)] + spaces[i] = [SplineSpace( p, knots=T , periodic=P, pads=pd) for p,T,P,pd in zip(degree_i, knots[interior.name], periodic, pads)] carts = create_cart(ddms, spaces) @@ -607,9 +650,12 @@ def discretize(a, *args, **kwargs): elif isinstance(a, BasicFunctionSpace): return discretize_space(a, *args, **kwargs) - elif isinstance(a, Derham): + elif isinstance(a, Derham) and not a.V0.is_broken: return discretize_derham(a, *args, **kwargs) + elif isinstance(a, Derham) and a.V0.is_broken: + return discretize_derham_multipatch(a, *args, **kwargs) + elif isinstance(a, Domain): return discretize_domain(a, *args, **kwargs) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index cfec097eb..11dcb66cb 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -1,18 +1,31 @@ -from sympde.topology.mapping import Mapping - -from psydac.api.basic import BasicDiscrete -from psydac.feec.derivatives import Derivative_1D, Gradient_2D, Gradient_3D -from psydac.feec.derivatives import ScalarCurl_2D, VectorCurl_2D, Curl_3D -from psydac.feec.derivatives import Divergence_2D, Divergence_3D -from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl, Projector_H1vec -from psydac.feec.global_projectors import Projector_Hdiv, Projector_L2 -from psydac.feec.pull_push import pull_1d_h1, pull_1d_l2 -from psydac.feec.pull_push import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2, pull_2d_h1vec -from psydac.feec.pull_push import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2, pull_3d_h1vec -from psydac.fem.basic import FemSpace -from psydac.fem.vector import VectorFemSpace - -__all__ = ('DiscreteDerham',) +from psydac.api.basic import BasicDiscrete + +from psydac.feec.derivatives import Derivative_1D, Gradient_2D, Gradient_3D +from psydac.feec.derivatives import ScalarCurl_2D, VectorCurl_2D, Curl_3D +from psydac.feec.derivatives import Divergence_2D, Divergence_3D +from psydac.feec.derivatives import BrokenGradient_2D +from psydac.feec.derivatives import BrokenScalarCurl_2D + +from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl, Projector_H1vec +from psydac.feec.global_projectors import Projector_Hdiv, Projector_L2 +from psydac.feec.global_projectors import Multipatch_Projector_H1 +from psydac.feec.global_projectors import Multipatch_Projector_Hcurl +from psydac.feec.global_projectors import Multipatch_Projector_L2 + +from psydac.feec.conforming_projectors import ConformingProjection_V0 +from psydac.feec.conforming_projectors import ConformingProjection_V1 + +from psydac.feec.hodge import HodgeOperator + +from psydac.feec.pull_push import pull_1d_h1, pull_1d_l2 +from psydac.feec.pull_push import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2, pull_2d_h1vec +from psydac.feec.pull_push import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2, pull_3d_h1vec + +from psydac.fem.basic import FemSpace, FemLinearOperator +from psydac.fem.vector import VectorFemSpace +from psydac.linalg.basic import IdentityOperator + +__all__ = ('DiscreteDerham', 'DiscreteDerhamMultipatch',) #============================================================================== class DiscreteDerham(BasicDiscrete): @@ -20,26 +33,19 @@ class DiscreteDerham(BasicDiscrete): Parameters ---------- - mapping : Mapping or None - Symbolic mapping from the logical space to the physical space, if any. + domain_h : Geometry + The discretized domain, which is a single-patch geometry. *spaces : list of FemSpace The discrete spaces of the de Rham sequence. Notes ----- - - The basic type Mapping is defined in module sympde.topology.mapping. - A discrete mapping (spline or NURBS) may be attached to it. - - This constructor should not be called directly, but rather from the `discretize_derham` function in `psydac.api.discretization`. - - - For the multipatch counterpart of this class please see - `MultipatchDiscreteDerham` in `psydac.feec.multipatch.api`. """ - def __init__(self, mapping, *spaces): + def __init__(self, domain_h, *spaces): - assert (mapping is None) or isinstance(mapping, Mapping) assert all(isinstance(space, FemSpace) for space in spaces) self.has_vec = isinstance(spaces[-1], VectorFemSpace) @@ -53,14 +59,19 @@ def __init__(self, mapping, *spaces): dim = len(spaces) - 1 self._spaces = spaces + self._domain_h = domain_h + self._sequence = tuple(space.symbolic_space.kind.name for space in spaces) self._dim = dim - self._mapping = mapping - self._callable_mapping = mapping.get_callable_mapping() if mapping else None + self._mapping = domain_h.domain.mapping + self._callable_mapping = self._mapping.get_callable_mapping() if self._mapping else None if dim == 1: D0 = Derivative_1D(spaces[0], spaces[1]) + spaces[0].diff = spaces[0].grad = D0 + self._derivatives = (D0,) + elif dim == 2: kind = spaces[1].symbolic_space.kind.name @@ -72,6 +83,8 @@ def __init__(self, mapping, *spaces): spaces[0].diff = spaces[0].grad = D0 spaces[1].diff = spaces[1].curl = D1 + self._derivatives = (D0, D1) + elif kind == 'hdiv': D0 = VectorCurl_2D(spaces[0], spaces[1]) @@ -80,6 +93,9 @@ def __init__(self, mapping, *spaces): spaces[0].diff = spaces[0].rot = D0 spaces[1].diff = spaces[1].div = D1 + self._derivatives = (D0, D1) + + elif dim == 3: D0 = Gradient_3D(spaces[0], spaces[1]) @@ -90,14 +106,28 @@ def __init__(self, mapping, *spaces): spaces[1].diff = spaces[1].curl = D1 spaces[2].diff = spaces[2].div = D2 + self._derivatives = (D0, D1, D2) + else: raise ValueError('Dimension {} is not available'.format(dim)) + self._Hodge_operators = () + self._conf_proj = () #-------------------------------------------------------------------------- @property def dim(self): """Dimension of the physical and logical domains, which are assumed to be the same.""" return self._dim + + @property + def domain_h(self): + """Discretized domain.""" + return self._domain_h + + @property + def spaces(self): + """Spaces of the proper de Rham sequence (excluding Hvec).""" + return self._spaces @property def V0(self): @@ -124,6 +154,10 @@ def V3(self): """Fourth space of the de Rham sequence : L2 space in 3d""" return self._spaces[3] + @property + def sequence(self): + return self._sequence + @property def H1vec(self): """Vector-valued H1 space built as the Cartesian product of N copies of V0, @@ -131,11 +165,6 @@ def H1vec(self): assert self.has_vec return self._H1vec - @property - def spaces(self): - """Spaces of the proper de Rham sequence (excluding Hvec).""" - return self._spaces - @property def mapping(self): """The mapping from the logical space to the physical space.""" @@ -146,21 +175,6 @@ def callable_mapping(self): """The mapping as a callable.""" return self._callable_mapping - @property - def derivatives_as_matrices(self): - """Differential operators of the De Rham sequence as LinearOperator objects.""" - return tuple(V.diff.matrix for V in self.spaces[:-1]) - - @property - def derivatives(self): - """Differential operators of the De Rham sequence as `DiffOperator` objects. - - Those are objects with `domain` and `codomain` properties that are `FemSpace`, - they act on `FemField` (they take a `FemField` of their `domain` as input and return - a `FemField` of their `codomain`. - """ - return tuple(V.diff for V in self.spaces[:-1]) - #-------------------------------------------------------------------------- def projectors(self, *, kind='global', nquads=None): """Projectors mapping callable functions of the physical coordinates to a @@ -264,3 +278,324 @@ def projectors(self, *, kind='global', nquads=None): else : return P0, P1, P2, P3 + #-------------------------------------------------------------------------- + def derivatives(self, kind='femlinop'): + if kind == 'femlinop': + return self._derivatives + elif kind == 'sparse': + return tuple(b_diff.tosparse for b_diff in self._derivatives) + elif kind == 'linop': + return tuple(b_diff.linop for b_diff in self._derivatives) + + #-------------------------------------------------------------------------- + def conforming_projectors(self, kind='femlinop', p_moments=-1, hom_bc=False): + """ + return the conforming projectors of the broken multi-patch space + + Parameters + ---------- + + p_moments : + The number of moments preserved by the projector. + + hom_bc: + Apply homogenous boundary conditions if True + + kind : + The kind of the projector, can be 'femlinop', 'sparse' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'sparse' returns a scipy sparse matrix + - 'linop' returns a psydac LinearOperator + + Returns + ------- + cP0, cP1, cP2 : Tuple of , or + The conforming projectors of each space and in desired form. + + """ + if hom_bc is None: + raise ValueError('please provide a value for "hom_bc" argument') + + if self.dim == 1: + raise NotImplementedError("1D projectors are not available") + + elif self.dim == 2: + if self.sequence[1] != 'hcurl': + raise NotImplementedError('2D sequence with H-div not available yet') + + else: + cP0 = ConformingProjection_V0(self.V0, p_moments=p_moments, hom_bc=hom_bc) + cP1 = ConformingProjection_V1(self.V1, p_moments=p_moments, hom_bc=hom_bc) + + I2 = IdentityOperator(self.V2.coeff_space) + cP2 = FemLinearOperator(fem_domain=self.V2, fem_codomain=self.V2, linop=I2, sparse_matrix=I2.tosparse()) + + self._conf_proj = (cP0, cP1, cP2) + + elif self.dim == 3: + raise NotImplementedError("3D projectors are not available") + + if kind == 'femlinop': + return cP0, cP1, cP2 + elif kind == 'sparse': + return cP0.tosparse, cP1.tosparse, cP2.tosparse + elif kind == 'linop': + return cP0.linop, cP1.linop, cP2.linop + + #-------------------------------------------------------------------------- + def _init_Hodge_operators(self, backend_language='python', load_dir=None): + """ + Initialize the Hodge operator for the multipatch de Rham sequence. + + Parameters + ---------- + + backend_language: + The backend used to accelerate the code + + load_dir: + Filename for storage in sparse matrix format + + """ + if not self._Hodge_operators: + + if self.dim == 1: + H0 = HodgeOperator(self.V0, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=0) + H1 = HodgeOperator(self.V1, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=1) + + self._Hodge_operators = (H0, H1) + + elif self.dim == 2: + + H0 = HodgeOperator(self.V0, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=0) + H1 = HodgeOperator(self.V1, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=1) + H2 = HodgeOperator(self.V2, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=2) + + self._Hodge_operators = (H0, H1, H2) + + elif self.dim == 3: + + H0 = HodgeOperator(self.V0, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=0) + H1 = HodgeOperator(self.V1, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=1) + H2 = HodgeOperator(self.V2, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=2) + H3 = HodgeOperator(self.V3, self.domain_h, backend_language=backend_language, load_dir=load_dir, load_space_index=3) + + self._Hodge_operators = (H0, H1, H2, H3) + + #-------------------------------------------------------------------------- + def _Hodge_kind(self, H, dual=False, kind='femlinop'): + """ + Helper function to return the Hodge operator in the specified form. + + Parameters + ---------- + H : + + dual : + If True, returns the dual Hodge operator + + kind : + The kind of the projector, can be 'femlinop', 'sparse' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'sparse' returns a scipy sparse matrix + - 'linop' returns a psydac LinearOperator + + Returns + ------- + Hodge operator in the specified form. + """ + + if not dual: + if kind == 'femlinop': + return H.Hodge + elif kind == 'sparse': + return H.tosparse + elif kind == 'linop': + return H.linop + else: + if kind == 'femlinop': + return H.dual_Hodge + elif kind == 'sparse': + return H.dual_tosparse + elif kind == 'linop': + return H.dual_linop + + #-------------------------------------------------------------------------- + def Hodge_operators(self, space=None, dual=False, kind='femlinop', backend_language='python', load_dir=None): + """ + Returns the Hodge operator for the given space and specified kind. + + Parameters + ---------- + space : str or None + The space for which to return the Hodge operator, can be 'V0', 'V1', 'V2' or None. + If None, returns a tuple with all three Hodge operators. + + dual : bool + If True, returns the dual Hodge operator. + + kind : + The kind of the projector, can be 'femlinop', 'sparse' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'sparse' returns a scipy sparse matrix + - 'linop' returns a psydac LinearOperator + + backend_language : str + The backend used to accelerate the code, default is 'python'. + + load_dir : str or None + Directory to load the Hodge operator from, if None the operator is computed on demand. + + Returns + ------- + Either one of the following Hodge operators of the specified kind and space or all of them if space is None. + + H0 : , or + + H1 : , or + + H2 : , or + """ + + if not self._Hodge_operators: + self._init_Hodge_operators(backend_language=backend_language, load_dir=load_dir) + + H0, H1, H2 = self._Hodge_operators + + if space == 'V0': + return self._Hodge_kind(self._Hodge_operators[0], dual=dual, kind=kind) + + elif space == 'V1': + return self._Hodge_kind(self._Hodge_operators[1], dual=dual, kind=kind) + + elif space == 'V2': + return self._Hodge_kind(self._Hodge_operators[2], dual=dual, kind=kind) + + elif space == 'V3': + return self._Hodge_kind(self._Hodge_operators[3], dual=dual, kind=kind) + + elif space is None: + return tuple(self._Hodge_kind(H, dual=dual, kind=kind) for H in self._Hodge_operators) + + +#============================================================================== +class DiscreteDerhamMultipatch(DiscreteDerham): + """ Represents the discrete De Rham sequence for multipatch domains. + It only works when the number of patches>1. + + Parameters + ---------- + domain_h: + The discrete domain + + spaces: + The discrete spaces that are contained in the De Rham sequence + """ + + def __init__(self, *, domain_h, spaces): + + dim = len(spaces) - 1 + self._spaces = tuple(spaces) + self._dim = dim + self._mapping = domain_h.domain.mapping + self._callable_mapping = [m.get_callable_mapping() for m in self._mapping.mappings.values()] if self._mapping else None + self._domain_h = domain_h + self._sequence = tuple(space.symbolic_space.kind.name for space in spaces) + + + if dim == 1: + raise NotImplementedError('1D FEEC multipatch non available yet') + + elif dim == 2: + + if self._sequence[1] == 'hcurl': + + self._derivatives = ( + BrokenGradient_2D(self.V0, self.V1), + BrokenScalarCurl_2D(self.V1, self.V2), # None, + ) + + elif self._sequence[1] == 'hdiv': + raise NotImplementedError('2D sequence with H-div not available yet') + + else: + raise ValueError('2D sequence not understood') + + elif dim == 3: + raise NotImplementedError('3D FEEC multipatch non available yet') + + else: + raise ValueError('Dimension {} is not available'.format(dim)) + + self._Hodge_operators = () + self._conf_proj = () + + #-------------------------------------------------------------------------- + @property + def H1vec(self): + raise NotImplementedError('Not implemented for Multipatch de Rham sequences.') + + #-------------------------------------------------------------------------- + def projectors(self, *, kind='global', nquads=None): + """ + This method returns the patch-wise commuting projectors on the broken multi-patch space + + Parameters + ---------- + kind: + The projectors kind, can be global or local + + nquads: + The number of quadrature points. + + Returns + ------- + P0: + Patch wise H1 projector + + P1: + Patch wise Hcurl projector + + P2: + Patch wise L2 projector + + Notes + ----- + - when applied to smooth functions they return conforming fields + - default 'global projectors' correspond to geometric interpolation/histopolation operators on Greville grids + - here 'global' is a patch-level notion, as the interpolation-type problems are solved on each patch independently + """ + if not (kind == 'global'): + raise NotImplementedError('only global projectors are available') + + if self.dim == 1: + raise NotImplementedError("1D projectors are not available") + + elif self.dim == 2: + P0 = Multipatch_Projector_H1(self.V0) + + if self.sequence[1] == 'hcurl': + P1 = Multipatch_Projector_Hcurl(self.V1, nquads=nquads) + else: + P1 = None # TODO: Multipatch_Projector_Hdiv(self.V1, nquads=nquads) + raise NotImplementedError('2D sequence with H-div not available yet') + + P2 = Multipatch_Projector_L2(self.V2, nquads=nquads) + + if self.mapping: + + P0_m = lambda f : P0([pull_2d_h1(f, m) for m in self.callable_mapping]) + + if self.sequence[1] == 'hcurl': + P1_m = lambda f : P1([pull_2d_hcurl(f, m) for m in self.callable_mapping]) + else: + raise NotImplementedError('2D sequence with H-div not available yet') + + P2_m = lambda f : P2([pull_2d_l2(f, m) for m in self.callable_mapping]) + + return P0_m, P1_m, P2_m + + return P0, P1, P2 + + elif self.dim == 3: + raise NotImplementedError("3D projectors are not available") diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index ad60a2c02..fd7e7f7dd 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -145,7 +145,7 @@ class CollelaMapping1D(Mapping): M1 = a1_h.assemble() # Differential operators - D0, = derham_h.derivatives_as_matrices + D0, = derham_h.derivatives(kind='linop') # Transpose of derivative matrix D0_T = D0.T diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 40b607615..c221af521 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -325,7 +325,7 @@ class CollelaMapping2D(Mapping): M2 = a2_h.assemble() # Differential operators (StencilMatrix or BlockLinearOperator objects) - D0, D1 = derham_h.derivatives_as_matrices + D0, D1 = derham_h.derivatives(kind='linop') # Discretize and assemble penalization matrix if not periodic: diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index d5220b5bc..90012ccab 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -123,7 +123,7 @@ def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, pe M2 = a2_h.assemble().tosparse().tocsr() # Diff operators - GRAD, CURL, DIV = derham_h.derivatives_as_matrices + GRAD, CURL, DIV = derham_h.derivatives(kind='linop') # Porjectors P0, P1, P2, P3 = derham_h.projectors(nquads=[5,5,5]) @@ -218,7 +218,7 @@ def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, M2 = a2_h.assemble() # Diff operators - GRAD, CURL, DIV = derham_h.derivatives_as_matrices + GRAD, CURL, DIV = derham_h.derivatives(kind='linop') # Porjectors P0, P1, P2, P3 = derham_h.projectors(nquads=[5,5,5]) diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index a84a41a2b..189b6c102 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -544,7 +544,7 @@ def test_assembly_no_synchr_args(backend): V1h = derham_h.V1 #differential operator - div, = derham_h.derivatives_as_matrices + div, = derham_h.derivatives(kind='linop') rho = element_of(V1h.symbolic_space, name='rho') g = element_of(V1h.symbolic_space, name='g') diff --git a/psydac/feec/multipatch/non_matching_operators.py b/psydac/feec/conforming_projectors.py similarity index 69% rename from psydac/feec/multipatch/non_matching_operators.py rename to psydac/feec/conforming_projectors.py index 4b172e7de..f9923f141 100644 --- a/psydac/feec/multipatch/non_matching_operators.py +++ b/psydac/feec/conforming_projectors.py @@ -1,20 +1,19 @@ -""" -This module provides utilities for constructing the conforming projections -for a H1-Hcurl-L2 broken FEEC de Rham sequence. -""" - +# coding: utf-8 +# Conga operators on piecewise (broken) de Rham sequences import os - import numpy as np -from scipy.sparse import eye as sparse_eye -from scipy.sparse import csr_matrix +from scipy.sparse import eye as sparse_eye +from scipy.sparse import csr_matrix +from scipy.special import comb from sympde.topology import Boundary, Interface -from psydac.fem.splines import SplineSpace -from psydac.utilities.quadratures import gauss_legendre -from psydac.core.bsplines import quadrature_grid, basis_ders_on_quad_grid, find_spans, elements_spans, cell_index, basis_ders_on_irregular_grid +from psydac.core.bsplines import quadrature_grid, basis_ders_on_quad_grid, find_spans, elements_spans, cell_index, basis_ders_on_irregular_grid +from psydac.fem.basic import FemLinearOperator +from psydac.fem.splines import SplineSpace +from psydac.utilities.quadratures import gauss_legendre +from psydac.linalg.utilities import SparseMatrixLinearOperator def get_patch_index_from_face(domain, face): @@ -118,7 +117,6 @@ def get_corners(domain, boundary_only): patches = domain.interior.args bd = domain.boundary - # corner_data[corner] = (patch_ind => local coordinates) corner_data = dict() if boundary_only: @@ -143,7 +141,6 @@ def get_corners(domain, boundary_only): else: for co in cos: corner_data[co] = dict() - for cb in co.corners: p_ind = patches.index(cb.domain) c_coord = cb.coordinates @@ -192,63 +189,48 @@ def construct_restriction_operator_1D( """ n_c = coarse_space_1d.nbasis n_f = fine_space_1d.nbasis - R = np.zeros((n_c, n_f)) if coarse_space_1d.basis == 'B': + #map V^+ to V^+_0 T = np.zeros((n_f, n_f)) - for i in range(1, n_f - 1): + for i in range(n_f): for j in range(n_f): - T[i, j] = int(i == j) - E[i, 0] * int(0 == j) - \ - E[i, -1] * int(n_f - 1 == j) + T[i, j] = int(i == j) - E[i, 0] * int(0 == j) - E[i, -1] * int(n_f - 1 == j) - cf_mass_mat = calculate_mixed_mass_matrix(coarse_space_1d, fine_space_1d)[ - 1:-1, 1:-1].transpose() - c_mass_mat = calculate_mass_matrix(coarse_space_1d)[1:-1, 1:-1] + cf_mass_mat = calculate_mixed_mass_matrix(coarse_space_1d, fine_space_1d).transpose() + c_mass_mat = calculate_mass_matrix(coarse_space_1d) if p_moments > 0: - - if not p_moments % 2 == 0: - p_moments += 1 - c_poly_mat = calculate_poly_basis_integral( - coarse_space_1d, p_moments=p_moments - 1)[:, 1:-1] - f_poly_mat = calculate_poly_basis_integral( - fine_space_1d, p_moments=p_moments - 1)[:, 1:-1] - - c_mass_mat[0:p_moments // 2, :] = c_poly_mat[0:p_moments // 2, :] - c_mass_mat[-p_moments // 2:, :] = c_poly_mat[-p_moments // 2:, :] - - cf_mass_mat[0:p_moments // 2, :] = f_poly_mat[0:p_moments // 2, :] - cf_mass_mat[-p_moments // 2:, :] = f_poly_mat[-p_moments // 2:, :] - - R0 = np.linalg.solve(c_mass_mat, cf_mass_mat) - R[1:-1, 1:-1] = R0 - R = R @ T - + # L^2 projection from V^+_0 to V^- + R[:, 1:-1] = np.linalg.solve(c_mass_mat, cf_mass_mat[:, 1:-1]) + gamma = get_1d_moment_correction(coarse_space_1d, p_moments=p_moments) + n = len(gamma) + + # maps V^- to V^+_0 in a moment preserving way + T2 = np.eye(n_c) + T2[0, 0] = T2[-1, -1] = 0 + T2[1:n+1, 0] += gamma + T2[-(n+1):-1, -1] += gamma[::-1] + + # maps V^+ to V^- in a moment preserving way + R = T2 @ R @ T + + else: + R[1:-1, 1:-1] = np.linalg.solve(c_mass_mat[1:-1, 1:-1], cf_mass_mat[1:-1, 1:-1]) + R = R @ T + + # add the degrees of freedom of T back R[0, 0] += 1 R[-1, -1] += 1 + else: - cf_mass_mat = calculate_mixed_mass_matrix( - coarse_space_1d, fine_space_1d).transpose() + cf_mass_mat = calculate_mixed_mass_matrix(coarse_space_1d, fine_space_1d).transpose() c_mass_mat = calculate_mass_matrix(coarse_space_1d) - if p_moments > 0: - - if not p_moments % 2 == 0: - p_moments += 1 - c_poly_mat = calculate_poly_basis_integral( - coarse_space_1d, p_moments=p_moments - 1) - f_poly_mat = calculate_poly_basis_integral( - fine_space_1d, p_moments=p_moments - 1) - - c_mass_mat[0:p_moments // 2, :] = c_poly_mat[0:p_moments // 2, :] - c_mass_mat[-p_moments // 2:, :] = c_poly_mat[-p_moments // 2:, :] - - cf_mass_mat[0:p_moments // 2, :] = f_poly_mat[0:p_moments // 2, :] - cf_mass_mat[-p_moments // 2:, :] = f_poly_mat[-p_moments // 2:, :] - + # The pure L^2 projection is already moment preserving R = np.linalg.solve(c_mass_mat, cf_mass_mat) return R @@ -287,8 +269,7 @@ def get_extension_restriction(coarse_space_1d, fine_space_1d, p_moments=-1): spl_type = coarse_space_1d.basis if not matching_interfaces: - grid = np.linspace( - fine_space_1d.breaks[0], fine_space_1d.breaks[-1], coarse_space_1d.ncells + 1) + grid = np.linspace(fine_space_1d.breaks[0], fine_space_1d.breaks[-1], coarse_space_1d.ncells + 1) coarse_space_1d_k_plus = SplineSpace( degree=fine_space_1d.degree, grid=grid, @@ -297,25 +278,18 @@ def get_extension_restriction(coarse_space_1d, fine_space_1d, p_moments=-1): E_1D = construct_extension_operator_1D( domain=coarse_space_1d_k_plus, codomain=fine_space_1d) + R_1D = construct_restriction_operator_1D( coarse_space_1d_k_plus, fine_space_1d, E_1D, p_moments) - + ER_1D = E_1D @ R_1D + assert np.allclose(R_1D @ E_1D, np.eye(coarse_space_1d.nbasis), 1e-12, 1e-12) + else: ER_1D = R_1D = E_1D = sparse_eye( fine_space_1d.nbasis, format="lil") - # TODO remove later - assert ( - np.allclose( - np.linalg.norm( - R_1D @ E_1D - - np.eye( - coarse_space_1d.nbasis)), - 0, - 1e-12, - 1e-12)) return E_1D, R_1D, ER_1D @@ -418,8 +392,7 @@ def calculate_mixed_mass_matrix(domain_space, codomain_space): fine_basis = basis_ders_on_quad_grid(fknots, fdeg, quad_x, 0, spl_type) coarse_basis = [ basis_ders_on_irregular_grid( - knots, deg, q, cell_index( - breaks, q), 0, spl_type) for q in quad_x] + knots, deg, q, cell_index(breaks, q), 0, spl_type) for q in quad_x] fine_spans = elements_spans(fknots, deg) coarse_spans = [find_spans(knots, deg, q[0])[0] for q in quad_x] @@ -471,7 +444,6 @@ def calculate_poly_basis_integral(space_1d, p_moments=-1): enddom = breaks[-1] begdom = breaks[0] denom = enddom - begdom - order = max(p_moments + 1, deg + 1) u, w = gauss_legendre(order) @@ -484,8 +456,7 @@ def calculate_poly_basis_integral(space_1d, p_moments=-1): Mass_mat = np.zeros((p_moments + 1, space_1d.nbasis)) for ie1 in range(Nel): # loop on cells - for pol in range( - p_moments + 1): # loops on basis function in each cell + for pol in range(p_moments + 1): # loops on basis function in each cell for il2 in range(deg + 1): # loops on basis function in each cell val = 0. @@ -494,7 +465,7 @@ def calculate_poly_basis_integral(space_1d, p_moments=-1): x = quad_x[ie1, q1] # val += quad_w[ie1, q1] * v0 * ((enddom-x)/denom)**pol val += quad_w[ie1, q1] * v0 * \ - ((enddom - x) / denom)**(p_moments - pol) * (x / denom)**pol + comb(p_moments, pol) * ((enddom - x) / denom)**(p_moments - pol) * ((x - begdom) / denom)**pol locind2 = il2 + spans[ie1] - deg Mass_mat[pol, locind2] += val @@ -531,16 +502,12 @@ def get_1d_moment_correction(space_1d, p_moments=-1): # local_shape[conf_axis]-2*(1+reg) such conforming functions p_max = space_1d.nbasis - 3 if p_max < p_moments: - print( - " ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** **") + print(" ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** **") print(" ** WARNING -- WARNING -- WARNING ") - print( - f" ** conf. projection imposing C0 smoothness on scalar space along this axis :") - print( - f" ** there are not enough dofs in a patch to preserve moments of degree {p_moments} !") + print(f" ** conf. projection imposing C0 smoothness on scalar space along this axis :") + print(f" ** there are not enough dofs in a patch to preserve moments of degree {p_moments} !") print(f" ** Only able to preserve up to degree --> {p_max} <-- ") - print( - " ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** **") + print(" ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** **") p_moments = p_max Mass_mat = calculate_poly_basis_integral(space_1d, p_moments) @@ -549,8 +516,10 @@ def get_1d_moment_correction(space_1d, p_moments=-1): return gamma -def construct_h1_conforming_projection( - Vh, reg_orders=0, p_moments=-1, hom_bc=False): +#============================================================================== +# Multipatch conforming projectors +#============================================================================== +def construct_h1_conforming_projection(Vh, reg_orders=0, p_moments=-1, hom_bc=False): """ Construct the conforming projection for a scalar space for a given regularity (0 continuous, -1 discontinuous). @@ -582,8 +551,7 @@ def construct_h1_conforming_projection( # moment corrections perpendicular to interfaces # assume same moments everywhere - gamma = get_1d_moment_correction( - Vh.patch_spaces[0].spaces[0], p_moments=p_moments) + gamma = get_1d_moment_correction(Vh.spaces[0].spaces[0], p_moments=p_moments) domain = Vh.symbolic_space.domain ndim = 2 @@ -592,22 +560,21 @@ def construct_h1_conforming_projection( l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) for k in range(n_patches): - Vk = Vh.patch_spaces[k] + Vk = Vh.spaces[k] # T is a TensorFemSpace and S is a 1D SplineSpace shapes = [S.nbasis for S in Vk.spaces] l2g.set_patch_shapes(k, shapes) # P vertex # vertex correction matrix - Proj_vertex = sparse_eye(dim_tot, format="lil") + Proj_vertex = sparse_eye(dim_tot, format="lil") corner_indices = set() corners = get_corners(domain, False) def get_vertex_index_from_patch(patch, coords): - # coords = co[patch] - nbasis0 = Vh.patch_spaces[patch].spaces[coords[0]].nbasis - 1 - nbasis1 = Vh.patch_spaces[patch].spaces[coords[1]].nbasis - 1 + nbasis0 = Vh.spaces[patch].spaces[coords[0]].nbasis - 1 + nbasis1 = Vh.spaces[patch].spaces[coords[1]].nbasis - 1 # patch local index multi_index = [None] * ndim @@ -621,12 +588,11 @@ def vertex_moment_indices(axis, coords, patch, p_moments): if coords[axis] == 0: return range(1, p_moments + 2) else: - return range(Vh.patch_spaces[patch].spaces[coords[axis]].nbasis - 1 - 1, - Vh.patch_spaces[patch].spaces[coords[axis]].nbasis - 1 - p_moments - 2, -1) + return range(Vh.spaces[patch].spaces[coords[axis]].nbasis - 1 - 1, + Vh.spaces[patch].spaces[coords[axis]].nbasis - 1 - p_moments - 2, -1) # loop over all vertices for (bd, co) in corners.items(): - # len(co)=#v is the number of adjacent patches at a vertex corr = len(co) @@ -639,9 +605,9 @@ def vertex_moment_indices(axis, coords, patch, p_moments): corner_indices.add(ig) for patch2 in co: - # local vertex coordinates in patch2 coords2 = co[patch2] + # global index jg = get_vertex_index_from_patch(patch2, coords2) @@ -701,7 +667,6 @@ def vertex_moment_indices(axis, coords, patch, p_moments): corners = get_corners(domain, True) if hom_bc: for (bd, co) in corners.items(): - for patch1 in co: # local vertex coordinates in patch2 @@ -714,6 +679,7 @@ def vertex_moment_indices(axis, coords, patch, p_moments): # local vertex coordinates in patch2 coords2 = co[patch2] + # global index jg = get_vertex_index_from_patch(patch2, coords2) @@ -830,8 +796,8 @@ def get_mu_minus(j, coarse_space, fine_space, R): k_minus = get_patch_index_from_face(domain, I.minus) k_plus = get_patch_index_from_face(domain, I.plus) - I_minus_ncells = Vh.patch_spaces[k_minus].ncells - I_plus_ncells = Vh.patch_spaces[k_plus].ncells + I_minus_ncells = Vh.spaces[k_minus].ncells + I_plus_ncells = Vh.spaces[k_plus].ncells # logical directions normal to interface if I_minus_ncells <= I_plus_ncells: @@ -848,8 +814,8 @@ def get_mu_minus(j, coarse_space, fine_space, R): d_fine = 1 - fine_axis d_coarse = 1 - coarse_axis - space_fine = Vh.patch_spaces[k_fine] - space_coarse = Vh.patch_spaces[k_coarse] + space_fine = Vh.spaces[k_fine] + space_coarse = Vh.spaces[k_coarse] coarse_space_1d = space_coarse.spaces[d_coarse] fine_space_1d = space_fine.spaces[d_fine] @@ -962,7 +928,7 @@ def get_mu_minus(j, coarse_space, fine_space, R): if hom_bc: for bn in domain.boundary: k = get_patch_index_from_face(domain, bn) - space_k = Vh.patch_spaces[k] + space_k = Vh.spaces[k] axis = bn.axis d = 1 - axis @@ -979,17 +945,19 @@ def get_mu_minus(j, coarse_space, fine_space, R): pg = edge_moment_index(p, i, axis, ext, space_k, k) Proj_edge[pg, ig] = gamma[p] else: - if corner_indices.issuperset({ig}): - mu_minus = get_mu_minus( - j, space_k_1d, space_k_1d, np.eye( - space_k_1d.nbasis)) + #if corner_indices.issuperset({ig}): + mu_minus = get_mu_minus( + i, space_k_1d, space_k_1d, np.eye( + space_k_1d.nbasis)) - for p in range(p_moments + 1): - for m in range(space_k_1d.nbasis): - pg = edge_moment_index( - p, m, axis, ext, space_k, k) - Proj_edge[pg, ig] = gamma[p] * mu_minus[m] - else: + for p in range(p_moments + 1): + for m in range(space_k_1d.nbasis): + pg = edge_moment_index( + p, m, axis, ext, space_k, k) + Proj_edge[pg, ig] = gamma[p] * mu_minus[m] + + if not corner_indices.issuperset({ig}): + corner_indices.add(ig) multi_index = [None] * ndim for p in range(p_moments + 1): @@ -997,16 +965,14 @@ def get_mu_minus(j, coarse_space, fine_space, R): 1 else space_k.spaces[axis].nbasis - 1 - p - 1 for pd in range(p_moments + 1): multi_index[1 - axis] = pd + \ - 1 if i == 0 else space_k.spaces[1 - - axis].nbasis - 1 - pd - 1 + 1 if i == 0 else space_k.spaces[1 - axis].nbasis - 1 - pd - 1 pg = l2g.get_index(k, 0, multi_index) Proj_edge[pg, ig] = gamma[p] * gamma[pd] return Proj_edge @ Proj_vertex -def construct_hcurl_conforming_projection( - Vh, reg_orders=0, p_moments=-1, hom_bc=False): +def construct_hcurl_conforming_projection(Vh, reg_orders=0, p_moments=-1, hom_bc=False): """ Construct the conforming projection for a vector Hcurl space for a given regularity (0 continuous, -1 discontinuous). @@ -1037,8 +1003,9 @@ def construct_hcurl_conforming_projection( return sparse_eye(dim_tot, format="lil") # moment corrections perpendicular to interfaces + # should be in the V^0 spaces gamma = [get_1d_moment_correction( - Vh.patch_spaces[0].spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] + Vh.spaces[0].spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] domain = Vh.symbolic_space.domain ndim = 2 @@ -1047,7 +1014,7 @@ def construct_hcurl_conforming_projection( l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) for k in range(n_patches): - Vk = Vh.patch_spaces[k] + Vk = Vh.spaces[k] # T is a TensorFemSpace and S is a 1D SplineSpace shapes = [[S.nbasis for S in T.spaces] for T in Vk.spaces] l2g.set_patch_shapes(k, *shapes) @@ -1073,7 +1040,7 @@ def edge_moment_index(p, i, axis, ext, space, k): multi_index[axis] = p + 1 if ext == - \ 1 else space.spaces[1 - axis].spaces[axis].nbasis - 1 - p - 1 return l2g.get_index(k, 1 - axis, multi_index) - + # loop over all interfaces for I in Interfaces: direction = I.ornt @@ -1086,8 +1053,8 @@ def edge_moment_index(p, i, axis, ext, space, k): minus_axis, plus_axis = I.minus.axis, I.plus.axis # logical directions along the interface d_minus, d_plus = 1 - minus_axis, 1 - plus_axis - I_minus_ncells = Vh.patch_spaces[k_minus].spaces[d_minus].ncells[d_minus] - I_plus_ncells = Vh.patch_spaces[k_plus].spaces[d_plus].ncells[d_plus] + I_minus_ncells = Vh.spaces[k_minus].spaces[d_minus].ncells[d_minus] + I_plus_ncells = Vh.spaces[k_plus].spaces[d_plus].ncells[d_plus] # logical directions normal to interface if I_minus_ncells <= I_plus_ncells: @@ -1104,8 +1071,8 @@ def edge_moment_index(p, i, axis, ext, space, k): d_fine = 1 - fine_axis d_coarse = 1 - coarse_axis - space_fine = Vh.patch_spaces[k_fine] - space_coarse = Vh.patch_spaces[k_coarse] + space_fine = Vh.spaces[k_fine] + space_coarse = Vh.spaces[k_coarse] coarse_space_1d = space_coarse.spaces[d_coarse].spaces[d_coarse] fine_space_1d = space_fine.spaces[d_fine].spaces[d_fine] @@ -1164,7 +1131,7 @@ def edge_moment_index(p, i, axis, ext, space, k): # boundary condition for bn in domain.boundary: k = get_patch_index_from_face(domain, bn) - space_k = Vh.patch_spaces[k] + space_k = Vh.spaces[k] axis = bn.axis if not hom_bc: @@ -1184,3 +1151,337 @@ def edge_moment_index(p, i, axis, ext, space, k): Proj_edge[pg, ig] = gamma[d][p] return Proj_edge + +#============================================================================== +# Singlepatch conforming projectors +#============================================================================== +def construct_hcurl_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=-1, hom_bc=False): + """ + Construct the conforming projection for a single patch vector Hcurl space for a given regularity (0 continuous, -1 discontinuous). + + Parameters + ---------- + Vh : TensorFemSpace + Finite Element Space coming from the discrete de Rham sequence. + + reg_orders : (int) + Regularity in each space direction -1 or 0. + + p_moments : (int) + Number of polynomial moments to be preserved. + + hom_bc : (bool) + Tangential homogeneous boundary conditions. + + Returns + ------- + cP : scipy.sparse.csr_array + Conforming projection as a sparse matrix. + """ + + dim_tot = Vh.nbasis + + # fully discontinuous space + if reg_orders < 0 or not hom_bc: + return sparse_eye(dim_tot, format="lil") + + # moment corrections perpendicular to interfaces + # should be in the V^0 spaces + + gamma = [get_1d_moment_correction(Vh.spaces[1 - d].spaces[d], p_moments=p_moments) for d in range(2)] + + domain = Vh.symbolic_space.domain + ndim = 2 + n_components = 2 + n_patches = len(domain) + + l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) + # T is a TensorFemSpace and S is a 1D SplineSpace + shapes = [[S.nbasis for S in T.spaces] for T in Vh.spaces] + l2g.set_patch_shapes(0, *shapes) + + # P edge + # edge correction matrix + Proj_edge = sparse_eye(dim_tot, format="lil") + + def get_edge_index(j, axis, ext): + multi_index = [None] * ndim + multi_index[axis] = 0 if ext == -1 else Vh.spaces[1 - axis].spaces[axis].nbasis - 1 + multi_index[1 - axis] = j + return l2g.get_index(0, 1 - axis, multi_index) + + def edge_moment_index(p, i, axis, ext): + multi_index = [None] * ndim + multi_index[1 - axis] = i + multi_index[axis] = p + 1 if ext == -1 else Vh.spaces[1 - axis].spaces[axis].nbasis - 1 - p - 1 + return l2g.get_index(0, 1 - axis, multi_index) + + + # boundary condition + for bn in domain.boundary: + + axis = bn.axis + d = 1 - axis + ext = bn.ext + space_1d = Vh.spaces[d].spaces[d] + + for i in range(0, space_1d.nbasis): + ig = get_edge_index(i, axis, ext) + Proj_edge[ig, ig] = 0 + + for p in range(p_moments + 1): + + pg = edge_moment_index(p, i, axis, ext) + Proj_edge[pg, ig] = gamma[d][p] + + return Proj_edge + + + + +def construct_h1_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=-1, hom_bc=False): + """ + Construct the conforming projection for a scalar space for a given regularity (0 continuous, -1 discontinuous). + + Parameters + ---------- + Vh : TensorFemSpace + Finite Element Space coming from the discrete de Rham sequence. + + reg_orders : (int) + Regularity in each space direction -1 or 0. + + p_moments : (int) + Number of moments to be preserved. + + hom_bc : (bool) + Homogeneous boundary conditions. + + Returns + ------- + cP : scipy.sparse.csr_array + Conforming projection as a sparse matrix. + """ + + dim_tot = Vh.nbasis + + # fully discontinuous space + if reg_orders < 0 or not hom_bc: + return sparse_eye(dim_tot, format="lil") + + # moment corrections perpendicular to interfaces + # assume same moments everywhere + gamma = get_1d_moment_correction(Vh.spaces[0], p_moments=p_moments) + + domain = Vh.symbolic_space.domain + ndim = 2 + n_components = 1 + n_patches = len(domain) + + l2g = Local2GlobalIndexMap(ndim, len(domain), n_components) + # T is a TensorFemSpace and S is a 1D SplineSpace + shapes = [S.nbasis for S in Vh.spaces] + l2g.set_patch_shapes(0, shapes) + + # P vertex + # vertex correction matrix + Proj_vertex = sparse_eye(dim_tot, format="lil") + + + def get_vertex_index(coords): + nbasis0 = Vh.spaces[coords[0]].nbasis - 1 + nbasis1 = Vh.spaces[coords[1]].nbasis - 1 + + # patch local index + multi_index = [None] * ndim + multi_index[0] = 0 if coords[0] == 0 else nbasis0 + multi_index[1] = 0 if coords[1] == 0 else nbasis1 + + # global index + return l2g.get_index(0, 0, multi_index) + + def vertex_moment_indices(axis, coords, p_moments): + if coords[axis] == 0: + return range(1, p_moments + 2) + else: + return range(Vh.spaces[coords[axis]].nbasis - 1 - 1, + Vh.spaces[coords[axis]].nbasis - 1 - p_moments - 2, -1) + + # boundary conditions + + for co in [(0,0), (1,0), (0,1), (1,1)]: + + # global index + ig = get_vertex_index(co) + + # conformity constraint + Proj_vertex[ig, ig] = 0 + + + if p_moments == -1: + continue + + # moment corrections from patch1 to patch1 + axis = 0 + d = 1 + multi_index_p = [None] * ndim + + d_moment_index = vertex_moment_indices(d, co, p_moments) + axis_moment_index = vertex_moment_indices(axis, co, p_moments) + + for pd in range(0, p_moments + 1): + multi_index_p[d] = d_moment_index[pd] + + for p in range(0, p_moments + 1): + multi_index_p[axis] = axis_moment_index[p] + + pg = l2g.get_index(0, 0, multi_index_p) + Proj_vertex[pg, ig] = gamma[p] * gamma[pd] + + # P edge + # edge correction matrix + Proj_edge = sparse_eye(dim_tot, format="lil") + + def get_edge_index(j, axis, ext): + multi_index = [None] * ndim + multi_index[axis] = 0 if ext == - 1 else Vh.spaces[axis].nbasis - 1 + multi_index[1 - axis] = j + return l2g.get_index(0, 0, multi_index) + + def edge_moment_index(p, i, axis, ext): + multi_index = [None] * ndim + multi_index[1 - axis] = i + multi_index[axis] = p + 1 if ext == -1 else Vh.spaces[axis].nbasis - 1 - p - 1 + return l2g.get_index(0, 0, multi_index) + + + def get_mu_minus(j, coarse_space, fine_space, R): + mu_plus = np.zeros(fine_space.nbasis) + mu_minus = np.zeros(coarse_space.nbasis) + + if j == 0: + mu_minus[0] = 1 + for p in range(p_moments + 1): + mu_plus[p + 1] = gamma[p] + else: + mu_minus[-1] = 1 + for p in range(p_moments + 1): + mu_plus[-1 - (p + 1)] = gamma[p] + + for m in range(coarse_space.nbasis): + for l in range(fine_space.nbasis): + mu_minus[m] += R[m, l] * mu_plus[l] + + if j == 0: + mu_minus[m] -= R[m, 0] + else: + mu_minus[m] -= R[m, -1] + + return mu_minus + + + # boundary condition + for bn in domain.boundary: + space_k = Vh + axis = bn.axis + + d = 1 - axis + ext = bn.ext + space_k_1d = space_k.spaces[d] + + for i in range(0, space_k_1d.nbasis): + ig = get_edge_index(i, axis, ext) + Proj_edge[ig, ig] = 0 + + if (i != 0 and i != space_k_1d.nbasis - 1): + for p in range(p_moments + 1): + + pg = edge_moment_index(p, i, axis, ext) + Proj_edge[pg, ig] = gamma[p] + else: + #if corner_indices.issuperset({ig}): + mu_minus = get_mu_minus( + i, space_k_1d, space_k_1d, np.eye( + space_k_1d.nbasis)) + + for p in range(p_moments + 1): + for m in range(space_k_1d.nbasis): + pg = edge_moment_index( + p, m, axis, ext) + Proj_edge[pg, ig] = gamma[p] * mu_minus[m] + + + return Proj_edge @ Proj_vertex + + +# =============================================================================== + +class ConformingProjection_V0(FemLinearOperator): + """ + Conforming projection from global broken V0 space to conforming global V0 space + Defined by averaging of interface dofs + + Parameters + ---------- + V0h: + The discrete space + + p_moments: + Number of polynomial moments to be preserved in the projection. + + hom_bc : + Apply homogenous boundary conditions if True + """ + # todo (MCP, 16.03.2021): + # - allow case without interfaces (single or multipatch) + + def __init__( + self, + V0h, + p_moments=-1, + hom_bc=False): + + FemLinearOperator.__init__(self, fem_domain=V0h, fem_codomain=V0h) + + if V0h.is_multipatch: + self._sparse_matrix = construct_h1_conforming_projection(V0h, reg_orders=0, p_moments=p_moments, hom_bc=hom_bc) + else: + self._sparse_matrix = construct_h1_singlepatch_conforming_projection(V0h, reg_orders=0, p_moments=p_moments, hom_bc=hom_bc) + + self._linop = SparseMatrixLinearOperator(self.linop_domain, self.linop_codomain, self._sparse_matrix) + + +class ConformingProjection_V1(FemLinearOperator): + """ + Conforming projection from global broken V1 space to conforming V1 global space + + proj.dot(v) returns the conforming projection of v, computed by solving linear system + + Parameters + ---------- + V1h: + The discrete space + + p_moments: + Number of polynomial moments to be preserved in the projection. + + hom_bc : + Apply homogenous boundary conditions if True + """ + # todo (MCP, 16.03.2021): + # - allow case without interfaces (single or multipatch) + + def __init__( + self, + V1h, + p_moments=-1, + hom_bc=False): + + FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V1h) + + if V1h.is_multipatch: + self._sparse_matrix = construct_hcurl_conforming_projection(V1h, reg_orders=0, p_moments=p_moments, hom_bc=hom_bc) + else: + self._sparse_matrix = construct_hcurl_singlepatch_conforming_projection(V1h, reg_orders=0, p_moments=p_moments, hom_bc=hom_bc) + + self._linop = SparseMatrixLinearOperator(self.linop_domain, self.linop_codomain, self._sparse_matrix) diff --git a/psydac/feec/derivatives.py b/psydac/feec/derivatives.py index f7f28ccbd..600db2d75 100644 --- a/psydac/feec/derivatives.py +++ b/psydac/feec/derivatives.py @@ -9,13 +9,12 @@ from psydac.fem.vector import VectorFemSpace from psydac.fem.tensor import TensorFemSpace from psydac.linalg.basic import IdentityOperator -from psydac.fem.basic import FemField, FemSpace +from psydac.fem.basic import FemField, FemSpace, FemLinearOperator from psydac.linalg.basic import LinearOperator from psydac.ddm.cart import DomainDecomposition, CartDecomposition __all__ = ( 'DirectionalDerivativeOperator', - 'DiffOperator', 'Derivative_1D', 'Gradient_2D', 'Gradient_3D', @@ -24,6 +23,10 @@ 'Curl_3D', 'Divergence_2D', 'Divergence_3D', + 'BrokenGradient_2D', + 'BrokenTransposedGradient_2D', + 'BrokenScalarCurl_2D', + 'BrokenTransposedScalarCurl_2D', 'block_tostencil' ) @@ -41,6 +44,8 @@ def block_tostencil(M): blocks[i1][i2] = mat.tostencil() return BlockLinearOperator(M.domain, M.codomain, blocks=blocks) +#==================================================================================================== +# Singlepatch derivative operators #==================================================================================================== class DirectionalDerivativeOperator(LinearOperator): """ @@ -361,40 +366,7 @@ def copy(self): self._diffdir, negative=self._negative, transposed=self._transposed) #==================================================================================================== -class DiffOperator: - def __init__(self, domain, codomain, matrix): - assert isinstance(domain, FemSpace) - assert isinstance(codomain, FemSpace) - assert isinstance(matrix, LinearOperator) - assert domain.coeff_space is matrix.domain - assert codomain.coeff_space is matrix.codomain - - self._domain = domain - self._codomain = codomain - self._matrix = matrix - - @property - def matrix(self): - return self._matrix - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._codomain - - def __call__(self, u): - assert isinstance(u, FemField) - assert u.space == self.domain - - coeffs = self.matrix.dot(u.coeffs) - - return FemField(self.codomain, coeffs=coeffs) - -#==================================================================================================== -class Derivative_1D(DiffOperator): +class Derivative_1D(FemLinearOperator): """ 1D derivative. @@ -415,10 +387,10 @@ def __init__(self, H1, L2): assert H1.degree[0] == L2.degree[0] + 1 # Store data in object - super().__init__(H1, L2, DirectionalDerivativeOperator(H1.coeff_space, L2.coeff_space, 0)) + super().__init__(fem_domain = H1, fem_codomain = L2, linop = DirectionalDerivativeOperator(H1.coeff_space, L2.coeff_space, 0)) #==================================================================================================== -class Gradient_2D(DiffOperator): +class Gradient_2D(FemLinearOperator): """ Gradient operator in 2D. @@ -434,7 +406,7 @@ class Gradient_2D(DiffOperator): def __init__(self, H1, Hcurl): assert isinstance( H1, TensorFemSpace); assert H1.ldim == 2 - assert isinstance(Hcurl, VectorFemSpace); assert Hcurl.ldim == 2 + assert isinstance(Hcurl, VectorFemSpace); assert Hcurl.ldim == 2 assert Hcurl.spaces[0].periodic == H1.periodic assert Hcurl.spaces[1].periodic == H1.periodic @@ -454,11 +426,10 @@ def __init__(self, H1, Hcurl): matrix = BlockLinearOperator(H1.coeff_space, Hcurl.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hcurl, matrix) - + super().__init__(fem_domain = H1, fem_codomain = Hcurl, linop = matrix) #==================================================================================================== -class Gradient_3D(DiffOperator): +class Gradient_3D(FemLinearOperator): """ Gradient operator in 3D. @@ -497,10 +468,10 @@ def __init__(self, H1, Hcurl): matrix = BlockLinearOperator(H1.coeff_space, Hcurl.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hcurl, matrix) + super().__init__(fem_domain = H1, fem_codomain = Hcurl, linop = matrix) #==================================================================================================== -class ScalarCurl_2D(DiffOperator): +class ScalarCurl_2D(FemLinearOperator): """ Scalar curl operator in 2D: computes a scalar field from a vector field. @@ -536,10 +507,10 @@ def __init__(self, Hcurl, L2): matrix = BlockLinearOperator(Hcurl.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hcurl, L2, matrix) + super().__init__(fem_domain = Hcurl, fem_codomain = L2, linop = matrix) #==================================================================================================== -class VectorCurl_2D(DiffOperator): +class VectorCurl_2D(FemLinearOperator): """ Vector curl operator in 2D: computes a vector field from a scalar field. This is sometimes called the 'rot' operator. @@ -576,10 +547,10 @@ def __init__(self, H1, Hdiv): matrix = BlockLinearOperator(H1.coeff_space, Hdiv.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hdiv, matrix) + super().__init__(fem_domain = H1, fem_codomain = Hdiv, linop = matrix) #==================================================================================================== -class Curl_3D(DiffOperator): +class Curl_3D(FemLinearOperator): """ Curl operator in 3D. @@ -626,10 +597,10 @@ def __init__(self, Hcurl, Hdiv): # ... # Store data in object - super().__init__(Hcurl, Hdiv, matrix) + super().__init__(fem_domain = Hcurl, fem_codomain = Hdiv, linop = matrix) #==================================================================================================== -class Divergence_2D(DiffOperator): +class Divergence_2D(FemLinearOperator): """ Divergence operator in 2D. @@ -665,10 +636,10 @@ def __init__(self, Hdiv, L2): matrix = BlockLinearOperator(Hdiv.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hdiv, L2, matrix) + super().__init__(fem_domain = Hdiv, fem_codomain = L2, linop = matrix) #==================================================================================================== -class Divergence_3D(DiffOperator): +class Divergence_3D(FemLinearOperator): """ Divergence operator in 3D. @@ -707,4 +678,70 @@ def __init__(self, Hdiv, L2): matrix = BlockLinearOperator(Hdiv.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hdiv, L2, matrix) + super().__init__(fem_domain = Hdiv, fem_codomain = L2, linop = matrix) + +#==================================================================================================== +# 2D Multipatch derivative operators +#==================================================================================================== +class BrokenGradient_2D(FemLinearOperator): + + def __init__(self, V0h, V1h): + + FemLinearOperator.__init__(self, fem_domain=V0h, fem_codomain=V1h) + + D0s = [Gradient_2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D0i.linop for i, D0i in enumerate(D0s)}) + + def transpose(self, conjugate=False): + # todo (MCP): define as the dual differential operator + return BrokenTransposedGradient_2D(self.fem_domain, self.fem_codomain) + +# ============================================================================== + +class BrokenTransposedGradient_2D(FemLinearOperator): + + def __init__(self, V0h, V1h): + + FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V0h) + + D0s = [Gradient_2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D0i.linop.T for i, D0i in enumerate(D0s)}) + + def transpose(self, conjugate=False): + # todo (MCP): discard + return BrokenGradient_2D(self.fem_codomain, self.fem_domain) + +# ============================================================================== +class BrokenScalarCurl_2D(FemLinearOperator): + def __init__(self, V1h, V2h): + + FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V2h) + + D1s = [ScalarCurl_2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D1i.linop for i, D1i in enumerate(D1s)}) + + def transpose(self, conjugate=False): + return BrokenTransposedScalarCurl_2D( + V1h=self.fem_domain, V2h=self.fem_codomain) + + +# ============================================================================== +class BrokenTransposedScalarCurl_2D(FemLinearOperator): + + def __init__(self, V1h, V2h): + + FemLinearOperator.__init__(self, fem_domain=V2h, fem_codomain=V1h) + + D1s = [ScalarCurl_2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D1i.linop.T for i, D1i in enumerate(D1s)}) + + def transpose(self, conjugate=False): + return BrokenScalarCurl_2D(V1h=self.fem_codomain, V2h=self.fem_domain) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index e5a61a0bb..e219c3fbf 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -4,7 +4,7 @@ from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix from psydac.linalg.stencil import StencilMatrix, StencilVectorSpace -from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.block import BlockLinearOperator, BlockVector from psydac.core.bsplines import quadrature_grid from psydac.utilities.quadratures import gauss_legendre from psydac.fem.basic import FemField @@ -18,6 +18,7 @@ from abc import ABCMeta, abstractmethod __all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2', + 'Multipatch_Projector_H1', 'Multipatch_Projector_Hcurl', 'Multipatch_Projector_L2', 'evaluate_dofs_1d_0form', 'evaluate_dofs_1d_1form', 'evaluate_dofs_2d_0form', 'evaluate_dofs_2d_1form_hcurl', 'evaluate_dofs_2d_1form_hdiv', 'evaluate_dofs_2d_2form', 'evaluate_dofs_3d_0form', 'evaluate_dofs_3d_1form', 'evaluate_dofs_3d_2form', 'evaluate_dofs_3d_3form') @@ -388,6 +389,8 @@ def __call__(self, fun): return FemField(self._space, coeffs=coeffs) +#============================================================================== +# SINGLEPATCH PROJECTORS #============================================================================== class Projector_H1(GlobalProjector): """ @@ -648,6 +651,7 @@ def __call__(self, fun): """ return super().__call__(fun) +#============================================================================== class Projector_H1vec(GlobalProjector): """ Projector from H1^3 = H1 x H1 x H1 to a conforming finite element space, i.e. @@ -712,6 +716,76 @@ def __call__(self, fun): """ return super().__call__(fun) +#============================================================================== +# MULTIPATCH PROJECTORS (2D) +#============================================================================== +class Multipatch_Projector_H1: + """ + to apply the H1 projection (2D) on every patch + """ + + def __init__(self, V0h): + + self._P0s = [Projector_H1(V) for V in V0h.spaces] + self._V0h = V0h # multipatch Fem Space + + def __call__(self, funs_log): + """ + project a list of functions given in the logical domain + """ + u0s = [P(fun) for P, fun, in zip(self._P0s, funs_log)] + + u0_coeffs = BlockVector(self._V0h.coeff_space, + blocks=[u0j.coeffs for u0j in u0s]) + + return FemField(self._V0h, coeffs=u0_coeffs) + +#============================================================================== +class Multipatch_Projector_Hcurl: + + """ + to apply the Hcurl projection (2D) on every patch + """ + + def __init__(self, V1h, nquads=None): + + self._P1s = [Projector_Hcurl(V, nquads=nquads) for V in V1h.spaces] + self._V1h = V1h # multipatch Fem Space + + def __call__(self, funs_log): + """ + project a list of functions given in the logical domain + """ + E1s = [P(fun) for P, fun, in zip(self._P1s, funs_log)] + + E1_coeffs = BlockVector(self._V1h.coeff_space, + blocks=[E1j.coeffs for E1j in E1s]) + + return FemField(self._V1h, coeffs=E1_coeffs) + +#============================================================================== +class Multipatch_Projector_L2: + + """ + to apply the L2 projection (2D) on every patch + """ + + def __init__(self, V2h, nquads=None): + + self._P2s = [Projector_L2(V, nquads=nquads) for V in V2h.spaces] + self._V2h = V2h # multipatch Fem Space + + def __call__(self, funs_log): + """ + project a list of functions given in the logical domain + """ + B2s = [P(fun) for P, fun, in zip(self._P2s, funs_log)] + + B2_coeffs = BlockVector(self._V2h.coeff_space, + blocks=[B2j.coeffs for B2j in B2s]) + + return FemField(self._V2h, coeffs=B2_coeffs) + #============================================================================== # 1D DEGREES OF FREEDOM #============================================================================== diff --git a/psydac/feec/hodge.py b/psydac/feec/hodge.py new file mode 100644 index 000000000..09b3f87ac --- /dev/null +++ b/psydac/feec/hodge.py @@ -0,0 +1,251 @@ +import os +import numpy as np + +from scipy.sparse import save_npz, load_npz +from scipy.sparse import block_diag +from scipy.sparse.linalg import inv + +from sympde.topology import element_of, elements_of +from sympde.topology.space import ScalarFunction +from sympde.calculus import dot +from sympde.expr.expr import BilinearForm +from sympde.expr.expr import integral + +from psydac.api.settings import PSYDAC_BACKENDS + +from psydac.linalg.utilities import SparseMatrixLinearOperator + +# =============================================================================== +class HodgeOperator: + """ + Change of basis operator: dual basis -> primal basis + + self._linop: matrix (LinearOperator) of the primal Hodge = this is the mass matrix ! + self.dual_linop: this is the INVERSE mass matrix (LinearOperator) + + Parameters + ---------- + Vh: + The discrete space + + domain_h: + The discrete domain of the projector + + metric : + the metric of the de Rham complex + + backend_language: + The backend used to accelerate the code + + load_dir: + storage files for the primal and dual Hodge sparse matrice + + load_space_index: + the space index in the derham sequence + + Notes + ----- + Either we use a storage, or these matrices are only computed on demand + # todo: we compute the sparse matrix when to_sparse_matrix is called -- but never the stencil matrix (should be fixed...) + We only support the identity metric, this implies that the dual Hodge is the inverse of the primal one. + # todo: allow for non-identity metrics + """ + + def __init__(self, Vh, domain_h, metric='identity', backend_language='python', load_dir=None, load_space_index=''): + + self._fem_domain = Vh + self._fem_codomain = Vh + + # FemLinearOperators + self._primal_Hodge = None + self._dual_Hodge = None + + # LinearOperators + self._linop = None + self._dual_linop = None + + # Sparse matrices + self._sparse_matrix = None + self._dual_sparse_matrix = None + + self._domain_h = domain_h + self._backend_language = backend_language + + assert metric == 'identity' + self._metric = metric + + if load_dir and isinstance(load_dir, str): + if not os.path.exists(load_dir): + os.makedirs(load_dir) + assert str(load_space_index) in ['0', '1', '2', '3'] + primal_Hodge_storage_fn = load_dir + '/H{}_m.npz'.format(load_space_index) + dual_Hodge_storage_fn = load_dir + '/dH{}_m.npz'.format(load_space_index) + + primal_Hodge_is_stored = os.path.exists(primal_Hodge_storage_fn) + dual_Hodge_is_stored = os.path.exists(dual_Hodge_storage_fn) + if dual_Hodge_is_stored: + assert primal_Hodge_is_stored + print(" ... loading dual Hodge sparse matrix from " + dual_Hodge_storage_fn) + self._dual_sparse_matrix = load_npz(dual_Hodge_storage_fn) + print("[HodgeOperator] loading primal Hodge sparse matrix from " + primal_Hodge_storage_fn) + self._sparse_matrix = load_npz(primal_Hodge_storage_fn) + else: + assert not primal_Hodge_is_stored + print("[HodgeOperator] assembling both sparse matrices for storage...") + self.assemble_matrix() + print("[HodgeOperator] storing primal Hodge sparse matrix in " + primal_Hodge_storage_fn) + save_npz(primal_Hodge_storage_fn, self._sparse_matrix) + self.assemble_dual_sparse_matrix() + print("[HodgeOperator] storing dual Hodge sparse matrix in " + dual_Hodge_storage_fn) + save_npz(dual_Hodge_storage_fn, self._dual_sparse_matrix) + else: + # matrices are not stored, we will probably compute them later + pass + + def assemble_matrix(self): + """ + the Hodge matrix is the patch-wise multi-patch mass matrix + it is not stored by default but assembled on demand + """ + from psydac.api.discretization import discretize + from psydac.fem.basic import FemLinearOperator + + if self._linop is None: + Vh = self._fem_domain + assert Vh == self._fem_codomain + + V = Vh.symbolic_space + domain = V.domain + # domain_h = V0h.domain: would be nice... + u, v = elements_of(V, names='u, v') + + if isinstance(u, ScalarFunction): + expr = u * v + else: + expr = dot(u, v) + + a = BilinearForm((u, v), integral(domain, expr)) + ah = discretize(a, self._domain_h, [Vh, Vh], backend=PSYDAC_BACKENDS[self._backend_language]) + + self._linop = ah.assemble() # Mass matrix in stencil format + self._sparse_matrix = self._linop.tosparse() + + self._primal_Hodge = FemLinearOperator(self._fem_domain, self._fem_codomain, linop=self._linop, sparse_matrix=self._sparse_matrix) + + # which of the two assemblys for the sparse dual matrix is better? For now use the exact one. + def assemble_dual_sparse_matrix(self): + """ + the dual Hodge sparse matrix is the patch-wise inverse of the multi-patch mass matrix + it is not stored by default but computed on demand, by local (patch-wise) exact inversion of the mass matrix + """ + from psydac.fem.basic import FemLinearOperator + + if self._dual_sparse_matrix is None: + if not self._linop: + self.assemble_matrix() + + M = self._linop # mass matrix of the (primal) basis + + if self._fem_domain.is_multipatch: + nrows = M.n_block_rows + ncols = M.n_block_cols + + inv_M_blocks = [] + for i in range(nrows): + Mii = M[i, i].tosparse() + inv_Mii = inv(Mii.tocsc()) + inv_Mii.eliminate_zeros() + inv_M_blocks.append(inv_Mii) + + inv_M = block_diag(inv_M_blocks) + + self._dual_sparse_matrix = inv_M + self._dual_linop = SparseMatrixLinearOperator(M.codomain, M.domain, inv_M) + self._dual_Hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop, sparse_matrix=self._dual_sparse_matrix) + + else: + M_m = M.tosparse() + inv_M = inv(M_m.tocsc()) + inv_M.eliminate_zeros() + + self._dual_sparse_matrix = inv_M + self._dual_linop = SparseMatrixLinearOperator(M.codomain, M.domain, inv_M) + self._dual_Hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop, sparse_matrix=self._dual_sparse_matrix) + + + def assemble_dual_matrix(self, solver ='gmres', **kwargs): + """ + the dual Hodge matrix is the patch-wise inverse of the multi-patch mass matrix + it is not stored by default but computed on demand, by approximate local (patch-wise) inversion of the mass matrix + """ + from psydac.linalg.solvers import inverse + from psydac.linalg.block import BlockLinearOperator + from psydac.fem.basic import FemLinearOperator + + if self._dual_linop is None: + if not self._linop: + self.assemble_matrix() + + M = self._linop # mass matrix of the (primal) basis + + if self._fem_domain.is_multipatch: + + nrows = M.n_block_rows + ncols = M.n_block_cols + + inv_M_blocks = [list(b) for b in M.blocks] + for i in range(nrows): + Mii = M[i, i] + inv_Mii = inverse(M[i,i], solver=solver, **kwargs) + inv_M_blocks[i][i] = inv_Mii + + self._dual_linop = BlockLinearOperator(M.codomain, M.domain, blocks=inv_M_blocks) + self._dual_sparse_matrix = None + self._dual_Hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop, sparse_matrix=self._dual_sparse_matrix) + + else: + inv_M = inverse(M, solver=solver, **kwargs) + self._dual_sparse_matrix = None + self._dual_Hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop, sparse_matrix=self._dual_sparse_matrix) + + @property + def linop(self): + if self._linop is None: + self.assemble_matrix() + + return self._linop + + @property + def tosparse(self): + if self._sparse_matrix is None: + self.assemble_matrix() + + return self._sparse_matrix + + @property + def dual_linop(self): + if self._dual_linop is None: + self.assemble_dual_sparse_matrix() + + return self._dual_linop + + @property + def dual_tosparse(self): + if self._dual_sparse_matrix is None: + self.assemble_dual_sparse_matrix() + + return self._dual_sparse_matrix + + @property + def Hodge(self): + if self._linop is None: + self.assemble_matrix() + + return self._primal_Hodge + + @property + def dual_Hodge(self): + if self._dual_linop is None: + self.assemble_dual_sparse_matrix() + + return self._dual_Hodge \ No newline at end of file diff --git a/psydac/feec/multipatch/api.py b/psydac/feec/multipatch/api.py deleted file mode 100644 index 566815421..000000000 --- a/psydac/feec/multipatch/api.py +++ /dev/null @@ -1,301 +0,0 @@ -# coding: utf-8 -import os - -from sympde.topology import Derham -from sympde.topology import element_of, elements_of -from sympde.topology.space import ScalarFunction -from sympde.calculus import grad, dot, inner, rot, div -from sympde.calculus import laplace, bracket, convect -from sympde.calculus import jump, avg, Dn, minus, plus -from sympde.expr.expr import LinearForm, BilinearForm, integral - -from psydac.api.settings import PSYDAC_BACKENDS - -from psydac.api.discretization import discretize as discretize_single_patch -from psydac.api.discretization import discretize_space -from psydac.api.discretization import DiscreteDerham -from psydac.feec.multipatch.operators import BrokenGradient_2D -from psydac.feec.multipatch.operators import BrokenScalarCurl_2D -from psydac.feec.multipatch.operators import Multipatch_Projector_H1 -from psydac.feec.multipatch.operators import Multipatch_Projector_Hcurl -from psydac.feec.multipatch.operators import Multipatch_Projector_L2 -from psydac.feec.multipatch.operators import ConformingProjection_V0 -from psydac.feec.multipatch.operators import ConformingProjection_V1 -from psydac.feec.multipatch.fem_linear_operators import IdLinearOperator - - -__all__ = ('DiscreteDerhamMultipatch', 'discretize', 'discretize_derham_multipatch') - -#============================================================================== -class DiscreteDerhamMultipatch(DiscreteDerham): - """ Represents the discrete De Rham sequence for multipatch domains. - It only works when the number of patches>1 - - Parameters - ---------- - mapping: - The mapping of the multipatch domain, the multipatch mapping contains the mapping of each patch - - domain_h: - The discrete domain - - spaces: - The discrete spaces that are contained in the De Rham sequence - - sequence: - The space kind of each space in the De Rham sequence - """ - - def __init__(self, *, mapping, domain_h, spaces, sequence=None): - - - dim = len(spaces) - 1 - self._dim = dim - self._mapping = mapping - self._spaces = tuple(spaces) - self._domain_h = domain_h - - if sequence: - if len(sequence) != dim + 1: - raise ValueError('Expected len(sequence) = {}, got {} instead'. - format(dim + 1, len(sequence))) - - if dim == 1: - self._sequence = ('h1', 'l2') - raise NotImplementedError('1D FEEC multipatch non available yet') - - elif dim == 2: - if sequence is None: - raise ValueError('Sequence must be specified in 2D case') - - elif tuple(sequence) == ('h1', 'hcurl', 'l2'): - self._sequence = tuple(sequence) - self._broken_diff_ops = ( - BrokenGradient_2D(self.V0, self.V1), - BrokenScalarCurl_2D(self.V1, self.V2), # None, - ) - - elif tuple(sequence) == ('h1', 'hdiv', 'l2'): - self._sequence = tuple(sequence) - raise NotImplementedError('2D sequence with H-div not available yet') - - else: - raise ValueError('2D sequence not understood') - - elif dim == 3: - self._sequence = ('h1', 'hcurl', 'hdiv', 'l2') - raise NotImplementedError('3D FEEC multipatch non available yet') - - else: - raise ValueError('Dimension {} is not available'.format(dim)) - - #-------------------------------------------------------------------------- - @property - def sequence(self): - return self._sequence - - # ... - @property - def broken_derivatives_as_operators(self): - return self._broken_diff_ops - - # ... - @property - def broken_derivatives_as_matrices(self): - return tuple(b_diff.matrix for b_diff in self._broken_diff_ops) - - #-------------------------------------------------------------------------- - def projectors(self, *, kind='global', nquads=None): - """ - This method returns the patch-wise commuting projectors on the broken multi-patch space - - Parameters - ---------- - kind: - The projectors kind, can be global or local - - nquads: - The number of quadrature points. - - Returns - ------- - P0: - Patch wise H1 projector - - P1: - Patch wise Hcurl projector - - P2: - Patch wise L2 projector - - Notes - ----- - - when applied to smooth functions they return conforming fields - - default 'global projectors' correspond to geometric interpolation/histopolation operators on Greville grids - - here 'global' is a patch-level notion, as the interpolation-type problems are solved on each patch independently - """ - if not (kind == 'global'): - raise NotImplementedError('only global projectors are available') - - if self.dim == 1: - raise NotImplementedError("1D projectors are not available") - - elif self.dim == 2: - P0 = Multipatch_Projector_H1(self.V0) - - if self.sequence[1] == 'hcurl': - P1 = Multipatch_Projector_Hcurl(self.V1, nquads=nquads) - else: - P1 = None # TODO: Multipatch_Projector_Hdiv(self.V1, nquads=nquads) - raise NotImplementedError('2D sequence with H-div not available yet') - - P2 = Multipatch_Projector_L2(self.V2, nquads=nquads) - return P0, P1, P2 - - elif self.dim == 3: - raise NotImplementedError("3D projectors are not available") - - #-------------------------------------------------------------------------- - def conforming_projection(self, space, hom_bc=False, backend_language="python", load_dir=None): - """ - return the conforming projectors of the broken multi-patch space - - Parameters - ---------- - space : - The space of the projector - - hom_bc: - Apply homogenous boundary conditions if True - - backend_language: - The backend used to accelerate the code - - load_dir: - Filename for storage in sparse matrix format - - Returns - ------- - Cp: - The conforming projector - - """ - if hom_bc is None: - raise ValueError('please provide a value for "hom_bc" argument') - - if isinstance(load_dir, str): - if not os.path.exists(load_dir): - os.makedirs(load_dir) - if space == 'V0': - P_name = 'cP0' - elif space == 'V1': - P_name = 'cP1' - elif space == 'V2': - P_name = 'cP2' - else: - raise ValueError(space) - - if hom_bc: - storage_fn = load_dir + '/{}_hom_m.npz'.format(P_name) - else: - storage_fn = load_dir + '/{}_m.npz'.format(P_name) - else: - storage_fn = None - - cP = None - if self.dim == 1: - raise NotImplementedError("1D projectors are not available") - - elif self.dim == 2: - if space == 'V0': - cP = ConformingProjection_V0(self.V0, self._domain_h, hom_bc=hom_bc, backend_language=backend_language, storage_fn=storage_fn) - elif space == 'V1': - if self.sequence[1] == 'hcurl': - cP = ConformingProjection_V1(self.V1, self._domain_h, hom_bc=hom_bc, backend_language=backend_language, storage_fn=storage_fn) - else: - raise NotImplementedError('2D sequence with H-div not available yet') - - elif space == 'V2': - cP = IdLinearOperator(self.V2) # no storage needed! - else: - raise ValueError('Invalid value for "space" argument: {}'.format(space)) - - elif self.dim == 3: - raise NotImplementedError("3D projectors are not available") - - return cP - - def get_dual_dofs(self, space, f, backend_language="python", return_format='stencil_array'): - """ - return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(V^k)) of a given function f, as a stencil array or numpy array - - Parameters - ---------- - space : - The space of the dual dofs - - f : - The function used for evaluation - - backend_language: - The backend used to accelerate the code - - return_format: - The format of the dofs, can be 'stencil_array' or 'numpy_array' - - Returns - ------- - tilde_f: - The dual dofs - """ - if space == 'V0': - Vh = self.V0 - elif space == 'V1': - Vh = self.V1 - elif space == 'V2': - Vh = self.V2 - else: - raise NotImplementedError("The space of kind {} is not available".format(space)) - - V = Vh.symbolic_space - v = element_of(V, name='v') - - if isinstance(v, ScalarFunction): - expr = f*v - else: - expr = dot(f,v) - - l = LinearForm(v, integral( V.domain, expr)) - lh = discretize(l, self._domain_h, Vh, backend=PSYDAC_BACKENDS[backend_language]) - tilde_f = lh.assemble() - - if return_format == 'numpy_array': - return tilde_f.toarray() - else: - return tilde_f - -#============================================================================== -def discretize_derham_multipatch(derham, domain_h, *args, **kwargs): - - ldim = derham.shape - mapping = derham.spaces[0].domain.mapping - - bases = ['B'] + ldim * ['M'] - spaces = [discretize_space(V, domain_h, *args, basis=basis, **kwargs) \ - for V, basis in zip(derham.spaces, bases)] - - return DiscreteDerhamMultipatch( - mapping = mapping, - domain_h = domain_h, - spaces = spaces, - sequence = [V.kind.name for V in derham.spaces] - ) - -#============================================================================== -def discretize(expr, *args, **kwargs): - - if isinstance(expr, Derham) and expr.V0.is_broken: - return discretize_derham_multipatch(expr, *args, **kwargs) - - else: - return discretize_single_patch(expr, *args, **kwargs) diff --git a/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py b/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py index 7b1e33792..656d91143 100644 --- a/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py +++ b/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py @@ -12,38 +12,21 @@ V0h --grad-> V1h -—curl-> V2h """ - -from mpi4py import MPI - import os import numpy as np -from collections import OrderedDict - -from sympy import lambdify -from scipy.sparse.linalg import spsolve -from sympde.calculus import dot -from sympde.expr.expr import LinearForm -from sympde.expr.expr import integral, Norm from sympde.topology import Derham -from sympde.topology import element_of -from psydac.api.settings import PSYDAC_BACKENDS -from psydac.feec.multipatch.api import discretize -from psydac.feec.pull_push import pull_2d_h1 -from psydac.feec.multipatch.utils_conga_2d import P0_phys +from psydac.api.discretization import discretize +from psydac.linalg.basic import IdentityOperator +from psydac.linalg.solvers import inverse -from psydac.feec.multipatch.fem_linear_operators import IdLinearOperator -from psydac.feec.multipatch.operators import HodgeOperator from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_h1 -from psydac.feec.multipatch.utilities import time_count -from psydac.feec.multipatch.non_matching_operators import construct_h1_conforming_projection, construct_hcurl_conforming_projection -from psydac.api.postprocessing import OutputManager, PostProcessManager +from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_h1 -from psydac.linalg.utilities import array_to_psydac -from psydac.fem.basic import FemField +from psydac.fem.projectors import get_dual_dofs +from psydac.fem.basic import FemField from psydac.api.postprocessing import OutputManager, PostProcessManager @@ -98,9 +81,6 @@ def solve_h1_source_pbm( print('building the multipatch domain...') domain = build_multipatch_domain(domain_name=domain_name) - mappings = OrderedDict([(P.logical_domain, P.mapping) - for P in domain.interior]) - mappings_list = list(mappings.values()) if isinstance(nc, int): ncells = [nc, nc] @@ -114,102 +94,89 @@ def solve_h1_source_pbm( derham = Derham(domain, ["H1", "Hcurl", "L2"]) derham_h = discretize(derham, domain_h, degree=degree) - # multi-patch (broken) spaces - V0h = derham_h.V0 - V1h = derham_h.V1 - V2h = derham_h.V2 + V0h, V1h, V2h = derham_h.spaces print('dim(V0h) = {}'.format(V0h.nbasis)) print('dim(V1h) = {}'.format(V1h.nbasis)) print('dim(V2h) = {}'.format(V2h.nbasis)) print('broken differential operators...') # broken (patch-wise) differential operators - bD0, bD1 = derham_h.broken_derivatives_as_operators - bD0_m = bD0.to_sparse_matrix() - - print('building the discrete operators:') - print('commuting projection operators...') - nquads = [4 * (d + 1) for d in degree] - P0, P1, P2 = derham_h.projectors(nquads=nquads) - - I0 = IdLinearOperator(V0h) - I0_m = I0.to_sparse_matrix() + bD0, bD1 = derham_h.derivatives(kind='linop') print('Hodge operators...') - # multi-patch (broken) linear operators / matrices - H0 = HodgeOperator(V0h, domain_h, backend_language=backend_language) - H1 = HodgeOperator(V1h, domain_h, backend_language=backend_language) - - H0_m = H0.to_sparse_matrix() # = mass matrix of V0 - dH0_m = H0.get_dual_Hodge_sparse_matrix() # = inverse mass matrix of V0 - H1_m = H1.to_sparse_matrix() # = mass matrix of V1 + # multi-patch (broken) linear operators + H0 = derham_h.Hodge_operators(space='V0', kind='linop', backend_language=backend_language) + H1 = derham_h.Hodge_operators(space='V1', kind='linop', backend_language=backend_language) + dH0 = derham_h.Hodge_operators(space='V0', kind='linop', dual=True, backend_language=backend_language) print('conforming projection operators...') # conforming Projections (should take into account the boundary conditions # of the continuous deRham sequence) - cP0_m = construct_h1_conforming_projection(V0h, hom_bc=True) - - def lift_u_bc(u_bc): - if u_bc is not None: - print('lifting the boundary condition in V0h... [warning: Not Tested Yet!]') - d_ubc_c = derham_h.get_dual_dofs(space='V0', f=u_bc, backend_language=backend_language, return_format='numpy_array') - ubc_c = dH0_m.dot(d_ubc_c) + cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc = True) - ubc_c = ubc_c - cP0_m.dot(ubc_c) - else: - ubc_c = None - return ubc_c + print('building the discrete operators:') - # Conga (projection-based) stiffness matrices: - # div grad: - pre_DG_m = - bD0_m.transpose() @ H1_m @ bD0_m + I0 = IdentityOperator(V0h.coeff_space) + + # div grad + DG = - bD0.T @ H1 @ bD0 # jump penalization: - jump_penal_m = I0_m - cP0_m - JP0_m = jump_penal_m.transpose() @ H0_m @ jump_penal_m + JP0 = (I0 - cP0).T @ H0 @ (I0 - cP0) # useful for the boundary condition (if present) - pre_A_m = cP0_m.transpose() @ (eta * H0_m - mu * pre_DG_m) - A_m = pre_A_m @ cP0_m + gamma_h * JP0_m + pre_A = cP0.T @ (eta * H0 - mu * DG) + + A = pre_A @ cP0 + gamma_h * JP0 + + + f_scal, u_bc, u_ex = get_source_and_solution_h1(source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,) - print('getting the source and ref solution...') - f_scal, u_bc, u_ex = get_source_and_solution_h1( - source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, - ) + df = get_dual_dofs(Vh=V0h, f=f_scal, domain_h=domain_h, backend_language=backend_language) + f = dH0.dot(df) + df = cP0.T @ df - # compute approximate source f_h - b_c = derham_h.get_dual_dofs(space='V0', f=f_scal, backend_language=backend_language, return_format='numpy_array') - # source in primal sequence for plotting - f_c = dH0_m.dot(b_c) - b_c = cP0_m.transpose() @ b_c + def lift_u_bc(u_bc): + if u_bc is not None: + du_bc = get_dual_dofs(Vh=V0h, f=u_bc, domain_h = domain_h, backend_language=backend_language) + u_bc = dH0.dot(du_bc) + u_bc = u_bc - cP0.dot(u_bc) - ubc_c = lift_u_bc(u_bc) + else: + u_bc = None + + return u_bc - if ubc_c is not None: + ubc = lift_u_bc(u_bc) + + if ubc is not None: # modified source for the homogeneous pbm print('modifying the source with lifted bc solution...') - b_c = b_c - pre_A_m.dot(ubc_c) + df = df - pre_A.dot(ubc) # direct solve with scipy spsolve - print('solving source problem with scipy.spsolve...') - uh_c = spsolve(A_m, b_c) + print('solving source problem with conjugate gradient...') + solver = inverse(A, solver='cg', tol=1e-8) + u = solver.solve(df) # project the homogeneous solution on the conforming problem space print('projecting the homogeneous solution on the conforming problem space...') - uh_c = cP0_m.dot(uh_c) + u = cP0.dot(u) - if ubc_c is not None: + if ubc is not None: # adding the lifted boundary condition print('adding the lifted boundary condition...') - uh_c += ubc_c + u += ubc - print('getting and plotting the FEM solution from numpy coefs array...') if u_ex: - u_ex_c = derham_h.get_dual_dofs(space='V0', f=u_ex, backend_language=backend_language, return_format='numpy_array') - u_ex_c = dH0_m.dot(u_ex_c) + u_ex = get_dual_dofs(Vh=V0h, f=u_ex, domain_h=domain_h, backend_language=backend_language) + u_ex = dH0.dot(u_ex) + if plot_dir is not None: + print('plotting the FEM solution...') + if not os.path.exists(plot_dir): os.makedirs(plot_dir) @@ -217,17 +184,14 @@ def lift_u_bc(u_bc): OM.add_spaces(V0h=V0h) OM.set_static() - stencil_coeffs = array_to_psydac(uh_c, V0h.coeff_space) - vh = FemField(V0h, coeffs=stencil_coeffs) - OM.export_fields(vh=vh) + uh = FemField(V0h, coeffs=u) + OM.export_fields(uh=uh) - stencil_coeffs = array_to_psydac(f_c, V0h.coeff_space) - fh = FemField(V0h, coeffs=stencil_coeffs) + fh = FemField(V0h, coeffs=f) OM.export_fields(fh=fh) if u_ex: - stencil_coeffs = array_to_psydac(u_ex_c, V0h.coeff_space) - uh_ex = FemField(V0h, coeffs=stencil_coeffs) + uh_ex = FemField(V0h, coeffs=u_ex) OM.export_fields(uh_ex=uh_ex) OM.export_space_info() @@ -243,7 +207,7 @@ def lift_u_bc(u_bc): grid=None, npts_per_cell=[6] * 2, snapshots='all', - fields='vh') + fields='uh') PM.export_to_vtk( plot_dir + "/f_h", @@ -263,9 +227,9 @@ def lift_u_bc(u_bc): PM.close() if u_ex: - err = uh_c - u_ex_c - rel_err = np.sqrt(np.dot(err, H0_m.dot(err)))/np.sqrt(np.dot(u_ex_c,H0_m.dot(u_ex_c))) - + err = u - u_ex + rel_err = np.sqrt( err.inner(H0.dot(err))) / np.sqrt(u_ex.inner(H0.dot(u_ex))) + return rel_err @@ -273,19 +237,21 @@ def lift_u_bc(u_bc): omega = np.sqrt(170) # source eta = -omega**2 + mu=0 + gamma_h = 10 source_type = 'manu_poisson_elliptic' domain_name = 'pretzel_f' - nc = 10 + nc = 4 deg = 2 run_dir = '{}_{}_nc={}_deg={}/'.format(domain_name, source_type, nc, deg) solve_h1_source_pbm( nc=nc, deg=deg, eta=eta, - mu=1, # 1, + mu=mu, # 1, domain_name=domain_name, source_type=source_type, backend_language='pyccel-gcc', diff --git a/psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py b/psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py index db54565c7..2009941eb 100644 --- a/psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py +++ b/psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py @@ -2,34 +2,25 @@ Solve the eigenvalue problem for the curl-curl operator in 2D with a FEEC discretization """ import os -from mpi4py import MPI - import numpy as np -import matplotlib.pyplot as plt -from collections import OrderedDict + from sympde.topology import Derham -from psydac.feec.multipatch.api import discretize +from psydac.api.discretization import discretize from psydac.api.settings import PSYDAC_BACKENDS -from psydac.feec.multipatch.fem_linear_operators import IdLinearOperator -from psydac.feec.multipatch.operators import HodgeOperator from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -from psydac.feec.multipatch.utilities import time_count, get_run_dir, get_plot_dir, get_mat_dir, get_sol_dir, diag_fn -from psydac.feec.multipatch.utils_conga_2d import write_diags_to_file +from psydac.feec.multipatch.utilities import time_count -from sympde.topology import Square -from sympde.topology import IdentityMapping, PolarMapping from scipy.sparse.linalg import spilu, lgmres from scipy.sparse.linalg import LinearOperator, eigsh, minres -from scipy.sparse import csr_matrix from scipy.linalg import norm +from psydac.linalg.basic import IdentityOperator from psydac.linalg.utilities import array_to_psydac from psydac.fem.basic import FemField from psydac.feec.multipatch.multipatch_domain_utilities import build_cartesian_multipatch_domain -from psydac.feec.multipatch.non_matching_operators import construct_h1_conforming_projection, construct_hcurl_conforming_projection from psydac.api.postprocessing import OutputManager, PostProcessManager @@ -110,10 +101,6 @@ def hcurl_solve_eigen_pbm(ncells=np.array([[8, 4], [4, 4]]), degree=(3, 3), doma elif ncells.ndim == 2: ncells = {patch.name: [ncells[int(patch.name[2])][int(patch.name[4])], ncells[int(patch.name[2])][int(patch.name[4])]] for patch in domain.interior} - - mappings = OrderedDict([(P.logical_domain, P.mapping) - for P in domain.interior]) - mappings_list = list(mappings.values()) t_stamp = time_count(t_stamp) print(' .. discrete domain...') @@ -128,9 +115,7 @@ def hcurl_solve_eigen_pbm(ncells=np.array([[8, 4], [4, 4]]), degree=(3, 3), doma print(' .. discrete derham sequence...') derham_h = discretize(derham, domain_h, degree=degree) - V0h = derham_h.V0 - V1h = derham_h.V1 - V2h = derham_h.V2 + V0h, V1h, V2h = derham_h.spaces print('dim(V0h) = {}'.format(V0h.nbasis)) print('dim(V1h) = {}'.format(V1h.nbasis)) print('dim(V2h) = {}'.format(V2h.nbasis)) @@ -142,64 +127,30 @@ def hcurl_solve_eigen_pbm(ncells=np.array([[8, 4], [4, 4]]), degree=(3, 3), doma print('building the discrete operators:') print('commuting projection operators...') - I1 = IdLinearOperator(V1h) - I1_m = I1.to_sparse_matrix() + I1 = IdentityOperator(V1h.coeff_space) t_stamp = time_count(t_stamp) print('Hodge operators...') # multi-patch (broken) linear operators / matrices - H0 = HodgeOperator( - V0h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=0) - H1 = HodgeOperator( - V1h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=1) - H2 = HodgeOperator( - V2h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=2) - - H0_m = H0.to_sparse_matrix() # = mass matrix of V0 - dH0_m = H0.get_dual_Hodge_sparse_matrix() # = inverse mass matrix of V0 - H1_m = H1.to_sparse_matrix() # = mass matrix of V1 - dH1_m = H1.get_dual_Hodge_sparse_matrix() # = inverse mass matrix of V1 - H2_m = H2.to_sparse_matrix() # = mass matrix of V2 - dH2_m = H2.get_dual_Hodge_sparse_matrix() # = inverse mass matrix of V2 + H0, H1, H2 = derham_h.Hodge_operators(kind='linop', backend_language=backend_language, load_dir=m_load_dir) + dH0, dH1, dH2 = derham_h.Hodge_operators(kind='linop', dual=True, backend_language=backend_language, load_dir=m_load_dir) t_stamp = time_count(t_stamp) print('conforming projection operators...') # conforming Projections (should take into account the boundary conditions # of the continuous deRham sequence) - cP0_m = construct_h1_conforming_projection(V0h, hom_bc=True) - cP1_m = construct_hcurl_conforming_projection(V1h, hom_bc=True) + cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc = True) - t_stamp = time_count(t_stamp) - print('broken differential operators...') - bD0, bD1 = derham_h.broken_derivatives_as_operators - bD0_m = bD0.to_sparse_matrix() - bD1_m = bD1.to_sparse_matrix() t_stamp = time_count(t_stamp) - print('converting some matrices to csr format...') + print('broken differential operators...') + bD0, bD1 = derham_h.derivatives(kind='linop') - H1_m = H1_m.tocsr() - dH1_m = dH1_m.tocsr() - H2_m = H2_m.tocsr() - bD1_m = bD1_m.tocsr() if not os.path.exists(plot_dir): os.makedirs(plot_dir) print('computing the full operator matrix...') - A_m = np.zeros_like(H1_m) # Conga (projection-based) stiffness matrices if mu != 0: @@ -208,33 +159,29 @@ def hcurl_solve_eigen_pbm(ncells=np.array([[8, 4], [4, 4]]), degree=(3, 3), doma print('mu = {}'.format(mu)) print('curl-curl stiffness matrix...') - pre_CC_m = bD1_m.transpose() @ H2_m @ bD1_m - CC_m = cP1_m.transpose() @ pre_CC_m @ cP1_m # Conga stiffness matrix - A_m += mu * CC_m + CC = cP1.T @ bD1.T @ H2 @ bD1 @ cP1 # Conga stiffness matrix + A = mu * CC if nu != 0: - pre_GD_m = - H1_m @ bD0_m @ cP0_m @ dH0_m @ cP0_m.transpose() @ bD0_m.transpose() @ H1_m - GD_m = cP1_m.transpose() @ pre_GD_m @ cP1_m # Conga stiffness matrix - A_m -= nu * GD_m + GD = - cP1.T @ H1 @ bD0 @ cP0 @ dH0 @ cP0.T @ bD0.T @ H1 @ cP1 + A -= nu * GD # jump stabilization in V1h: if gamma_h != 0 or generalized_pbm: t_stamp = time_count(t_stamp) print('jump stabilization matrix...') - jump_stab_m = I1_m - cP1_m - JS_m = jump_stab_m.transpose() @ H1_m @ jump_stab_m - A_m += gamma_h * JS_m + JS = (I1 - cP1).T @ H1 @ (I1 - cP1) + A += gamma_h * JS if generalized_pbm: print('adding jump stabilization to RHS of generalized eigenproblem...') - B_m = cP1_m.transpose() @ H1_m @ cP1_m + JS_m + B = cP1.T @ H1 @ cP1 + JS else: - B_m = H1_m + B = H1 t_stamp = time_count(t_stamp) print('solving matrix eigenproblem...') - all_eigenvalues, all_eigenvectors_transp = get_eigenvalues( - nb_eigs_solve, sigma, A_m, B_m) + all_eigenvalues, all_eigenvectors_transp = get_eigenvalues(nb_eigs_solve, sigma, A.tosparse(), B.tosparse()) # Eigenvalue processing t_stamp = time_count(t_stamp) print('sorting out eigenvalues...') @@ -267,6 +214,9 @@ def hcurl_solve_eigen_pbm(ncells=np.array([[8, 4], [4, 4]]), degree=(3, 3), doma OM.export_space_info() nb_eigs = len(eigenvalues) + H1_m = H1.tosparse() + cP1_m = cP1.tosparse() + for i in range(min(nb_eigs_plot, nb_eigs)): print('looking at emode i = {}... '.format(i)) diff --git a/psydac/feec/multipatch/examples/hcurl_eigen_pbms_dg_2d.py b/psydac/feec/multipatch/examples/hcurl_eigen_pbms_dg_2d.py index 58c6a2bb6..710c2c1b4 100644 --- a/psydac/feec/multipatch/examples/hcurl_eigen_pbms_dg_2d.py +++ b/psydac/feec/multipatch/examples/hcurl_eigen_pbms_dg_2d.py @@ -3,14 +3,13 @@ A. Buffa and I. Perugia, “Discontinuous Galerkin Approximation of the Maxwell Eigenproblem” SIAM Journal on Numerical Analysis 44 (2006) """ - import os from mpi4py import MPI from collections import OrderedDict import numpy as np import matplotlib.pyplot -from scipy.sparse.linalg import spsolve, inv + from scipy.sparse.linalg import LinearOperator, eigsh, minres from sympde.calculus import grad, dot, curl, cross @@ -26,14 +25,12 @@ from sympde.expr.equation import find, EssentialBC from psydac.linalg.utilities import array_to_psydac -from psydac.api.tests.build_domain import build_pretzel from psydac.fem.basic import FemField -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.feec.pull_push import pull_2d_hcurl from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -from psydac.feec.multipatch.utilities import time_count, get_run_dir, get_plot_dir, get_mat_dir, get_sol_dir, diag_fn -from psydac.feec.multipatch.api import discretize +from psydac.feec.multipatch.utilities import time_count +from psydac.api.discretization import discretize from psydac.feec.multipatch.multipatch_domain_utilities import build_cartesian_multipatch_domain from psydac.api.postprocessing import OutputManager, PostProcessManager diff --git a/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py b/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py index b7442e054..c7f8481cf 100644 --- a/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py +++ b/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py @@ -14,40 +14,29 @@ """ import os -from mpi4py import MPI import numpy as np -from collections import OrderedDict -from sympy import lambdify, Matrix - -from scipy.sparse.linalg import spsolve - -from sympde.calculus import dot -from sympde.topology import element_of -from sympde.expr.expr import LinearForm -from sympde.expr.expr import integral, Norm from sympde.topology import Derham -from psydac.api.settings import PSYDAC_BACKENDS -from psydac.feec.pull_push import pull_2d_hcurl -from psydac.feec.multipatch.api import discretize -from psydac.feec.multipatch.fem_linear_operators import IdLinearOperator -from psydac.feec.multipatch.operators import HodgeOperator +from psydac.api.discretization import discretize from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_hcurl -from psydac.feec.multipatch.utils_conga_2d import DiagGrid, P0_phys, P1_phys, P2_phys, get_Vh_diags_for +from psydac.feec.multipatch.utils_conga_2d import P1_phys from psydac.feec.multipatch.utilities import time_count -from psydac.linalg.utilities import array_to_psydac +# from psydac.linalg.utilities import array_to_psydac from psydac.fem.basic import FemField -from psydac.feec.multipatch.non_matching_operators import construct_h1_conforming_projection, construct_hcurl_conforming_projection from psydac.api.postprocessing import OutputManager, PostProcessManager +from psydac.linalg.basic import IdentityOperator +from psydac.fem.projectors import get_dual_dofs +from psydac.linalg.solvers import inverse + def solve_hcurl_source_pbm( - nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='P_geom', source_type='manu_J', + nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='tilde_Pi', source_type='manu_J', eta=-10., mu=1., nu=1., gamma_h=10., - project_sol=False, plot_dir=None, + project_sol=True, plot_dir=None, m_load_dir=None, ): """ @@ -107,9 +96,9 @@ def solve_hcurl_source_pbm( t_stamp = time_count() print(' .. multi-patch domain...') domain = build_multipatch_domain(domain_name=domain_name) - mappings = OrderedDict([(P.logical_domain, P.mapping) - for P in domain.interior]) - mappings_list = list(mappings.values()) + # mappings = OrderedDict([(P.logical_domain, P.mapping) + # for P in domain.interior]) + # mappings_list = list(mappings.values()) if isinstance(nc, int): ncells = [nc, nc] @@ -117,8 +106,6 @@ def solve_hcurl_source_pbm( ncells = {patch.name: [nc[i], nc[i]] for (i, patch) in enumerate(domain.interior)} - # for diagnosttics - diag_grid = DiagGrid(mappings=mappings, N_diag=100) t_stamp = time_count(t_stamp) print(' .. derham sequence...') @@ -134,14 +121,14 @@ def solve_hcurl_source_pbm( t_stamp = time_count(t_stamp) print(' .. commuting projection operators...') - nquads = [4 * (d + 1) for d in degree] + nquads = [10 * (d + 1) for d in degree] P0, P1, P2 = derham_h.projectors(nquads=nquads) t_stamp = time_count(t_stamp) print(' .. multi-patch spaces...') - V0h = derham_h.V0 - V1h = derham_h.V1 - V2h = derham_h.V2 + V0h, V1h, V2h = derham_h.spaces + mappings = derham_h.callable_mapping + print('dim(V0h) = {}'.format(V0h.nbasis)) print('dim(V1h) = {}'.format(V1h.nbasis)) print('dim(V2h) = {}'.format(V2h.nbasis)) @@ -151,101 +138,44 @@ def solve_hcurl_source_pbm( t_stamp = time_count(t_stamp) print(' .. Id operator and matrix...') - I1 = IdLinearOperator(V1h) - I1_m = I1.to_sparse_matrix() + I1 = IdentityOperator(V1h.coeff_space) t_stamp = time_count(t_stamp) print(' .. Hodge operators...') # multi-patch (broken) linear operators / matrices # other option: define as Hodge Operators: - H0 = HodgeOperator( - V0h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=0) - H1 = HodgeOperator( - V1h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=1) - H2 = HodgeOperator( - V2h, - domain_h, - backend_language=backend_language, - load_dir=m_load_dir, - load_space_index=2) - - t_stamp = time_count(t_stamp) - print(' .. Hodge matrix H0_m = M0_m ...') - H0_m = H0.to_sparse_matrix() - t_stamp = time_count(t_stamp) - print(' .. dual Hodge matrix dH0_m = inv_M0_m ...') - dH0_m = H0.get_dual_Hodge_sparse_matrix() - - t_stamp = time_count(t_stamp) - print(' .. Hodge matrix H1_m = M1_m ...') - H1_m = H1.to_sparse_matrix() - t_stamp = time_count(t_stamp) - print(' .. dual Hodge matrix dH1_m = inv_M1_m ...') - dH1_m = H1.get_dual_Hodge_sparse_matrix() + H0, H1, H2 = derham_h.Hodge_operators(kind='linop', backend_language=backend_language) + dH0, dH1, dH2 = derham_h.Hodge_operators(kind='linop', dual=True, backend_language=backend_language) - t_stamp = time_count(t_stamp) - print(' .. Hodge matrix H2_m = M2_m ...') - H2_m = H2.to_sparse_matrix() - dH2_m = H2.get_dual_Hodge_sparse_matrix() t_stamp = time_count(t_stamp) print(' .. conforming Projection operators...') # conforming Projections (should take into account the boundary conditions # of the continuous deRham sequence) - cP0_m = construct_h1_conforming_projection(V0h, hom_bc=True) - cP1_m = construct_hcurl_conforming_projection(V1h, hom_bc=True) + cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc = True) + t_stamp = time_count(t_stamp) print(' .. broken differential operators...') # broken (patch-wise) differential operators - bD0, bD1 = derham_h.broken_derivatives_as_operators - bD0_m = bD0.to_sparse_matrix() - bD1_m = bD1.to_sparse_matrix() - - if plot_dir is not None and not os.path.exists(plot_dir): - os.makedirs(plot_dir) - - def lift_u_bc(u_bc): - if u_bc is not None: - print('lifting the boundary condition in V1h...') - # note: for simplicity we apply the full P1 on u_bc, but we only - # need to set the boundary dofs - uh_bc = P1_phys(u_bc, P1, domain, mappings_list) - ubc_c = uh_bc.coeffs.toarray() - # removing internal dofs (otherwise ubc_c may already be a very - # good approximation of uh_c ...) - ubc_c = ubc_c - cP1_m.dot(ubc_c) - else: - ubc_c = None - return ubc_c + bD0, bD1 = derham_h.derivatives(kind='linop') # Conga (projection-based) stiffness matrices # curl curl: t_stamp = time_count(t_stamp) print(' .. curl-curl stiffness matrix...') - print(bD1_m.shape, H2_m.shape) - pre_CC_m = bD1_m.transpose() @ H2_m @ bD1_m - # CC_m = cP1_m.transpose() @ pre_CC_m @ cP1_m # Conga stiffness matrix + pre_CC = bD1.T @ H2 @ bD1 # grad div: t_stamp = time_count(t_stamp) print(' .. grad-div stiffness matrix...') - pre_GD_m = - H1_m @ bD0_m @ cP0_m @ dH0_m @ cP0_m.transpose() @ bD0_m.transpose() @ H1_m - # GD_m = cP1_m.transpose() @ pre_GD_m @ cP1_m # Conga stiffness matrix + pre_GD = - H1 @ bD0 @ cP0 @ dH0 @ cP0.T @ bD0.T @ H1 # jump stabilization: t_stamp = time_count(t_stamp) print(' .. jump stabilization matrix...') - jump_penal_m = I1_m - cP1_m - JP_m = jump_penal_m.transpose() @ H1_m @ jump_penal_m + JS = (I1 - cP1).T @ H1 @ (I1 - cP1) + t_stamp = time_count(t_stamp) print(' .. full operator matrix...') @@ -254,69 +184,81 @@ def lift_u_bc(u_bc): print('nu = {}'.format(nu)) print('STABILIZATION: gamma_h = {}'.format(gamma_h)) # useful for the boundary condition (if present) - pre_A_m = cP1_m.transpose() @ (eta * H1_m + mu * pre_CC_m - nu * pre_GD_m) - A_m = pre_A_m @ cP1_m + gamma_h * JP_m + pre_A = eta * cP1.T @ H1 + if mu != 0: + pre_A += mu * cP1.T @ pre_CC + if nu != 0: + pre_A -= nu * cP1.T @ pre_GD + + A = pre_A @ cP1 + gamma_h * JS t_stamp = time_count(t_stamp) print() print(' -- getting source --') - f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl( - source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,) + f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl(source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name,) # compute approximate source f_h t_stamp = time_count(t_stamp) # f_h = L2 projection of f_vect, with filtering if tilde_Pi - print(' .. projecting the source with ' + - source_proj +' projection...') - - tilde_f_c = derham_h.get_dual_dofs( - space='V1', - f=f_vect, - backend_language=backend_language, - return_format='numpy_array') + print(' .. projecting the source with ' + source_proj +' projection...') + + tilde_f = get_dual_dofs(Vh=V1h, f=f_vect, domain_h=domain_h, backend_language=backend_language) + if source_proj == 'tilde_Pi': - print(' .. filtering the discrete source with P0.T ...') - tilde_f_c = cP1_m.transpose() @ tilde_f_c + print(' .. filtering the discrete source with P1.T ...') + tilde_f = cP1.T @ tilde_f + + def lift_u_bc(u_bc): + if u_bc is not None: + ubc = P1_phys(u_bc, P1, domain).coeffs + ubc = ubc - cP1.dot(ubc) + + else: + ubc = None + + return ubc - ubc_c = lift_u_bc(u_bc) - if ubc_c is not None: + ubc = lift_u_bc(u_bc) + + if ubc is not None: # modified source for the homogeneous pbm t_stamp = time_count(t_stamp) print(' .. modifying the source with lifted bc solution...') - tilde_f_c = tilde_f_c - pre_A_m.dot(ubc_c) + tilde_f = tilde_f - pre_A.dot(ubc) # direct solve with scipy spsolve t_stamp = time_count(t_stamp) - print() - print(' -- solving source problem with scipy.spsolve...') - uh_c = spsolve(A_m, tilde_f_c) + print('solving source problem with conjugate gradient...') + solver = inverse(A, solver='cg', tol=1e-8) + u = solver.solve(tilde_f) # project the homogeneous solution on the conforming problem space + t_stamp = time_count(t_stamp) if project_sol: - t_stamp = time_count(t_stamp) print(' .. projecting the homogeneous solution on the conforming problem space...') - uh_c = cP1_m.dot(uh_c) - else: - print(' .. NOT projecting the homogeneous solution on the conforming problem space') + u = cP1.dot(u) - if ubc_c is not None: + if ubc is not None: # adding the lifted boundary condition t_stamp = time_count(t_stamp) print(' .. adding the lifted boundary condition...') - uh_c += ubc_c + u += ubc - uh = FemField(V1h, coeffs=array_to_psydac(uh_c, V1h.coeff_space)) + uh = FemField(V1h, coeffs=u) #need cp1 here? - f_c = dH1_m.dot(tilde_f_c) - jh = FemField(V1h, coeffs=array_to_psydac(f_c, V1h.coeff_space)) + f = dH1.dot(tilde_f) + jh = FemField(V1h, coeffs=f) t_stamp = time_count(t_stamp) print(' -- plots and diagnostics --') if plot_dir: + if not os.path.exists(plot_dir): + os.makedirs(plot_dir) + OM = OutputManager(plot_dir + '/spaces.yml', plot_dir + '/fields.h5') OM.add_spaces(V1h=V1h) OM.set_static() @@ -349,11 +291,12 @@ def lift_u_bc(u_bc): time_count(t_stamp) if u_ex: - u_ex_c = P1_phys(u_ex, P1, domain, mappings_list).coeffs.toarray() - err = u_ex_c - uh_c - l2_error = np.sqrt(np.dot(err, H1_m.dot(err)))/np.sqrt(np.dot(u_ex_c,H1_m.dot(u_ex_c))) + u_ex_p = P1_phys(u_ex, P1, domain).coeffs + + err = u_ex_p - u + print(err.inner(H1.dot(err))) + l2_error = np.sqrt( err.inner(H1.dot(err))) / np.sqrt(u_ex_p.inner(H1.dot(u_ex_p))) print(l2_error) - #return l2_error diags['err'] = l2_error return diags diff --git a/psydac/feec/multipatch/examples/ppc_test_cases.py b/psydac/feec/multipatch/examples/ppc_test_cases.py index 94b772f4d..2fd889c49 100644 --- a/psydac/feec/multipatch/examples/ppc_test_cases.py +++ b/psydac/feec/multipatch/examples/ppc_test_cases.py @@ -1,22 +1,11 @@ # coding: utf-8 - -from sympy.functions.special.error_functions import erf -from mpi4py import MPI - import os import numpy as np from sympy import pi, cos, sin, Tuple, exp, atan, atan2 +from sympy.functions.special.error_functions import erf -from sympde.topology import Derham - -from psydac.fem.basic import FemField -from psydac.feec.multipatch.api import discretize -from psydac.feec.multipatch.operators import HodgeOperator -from psydac.fem.plotting_utilities import get_plotting_grid, my_small_plot, my_small_streamplot -from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -comm = MPI.COMM_WORLD # todo [MCP, 12/02/2022]: add an 'equation' argument to be able to return diff --git a/psydac/feec/multipatch/fem_linear_operators.py b/psydac/feec/multipatch/fem_linear_operators.py deleted file mode 100644 index aafb77931..000000000 --- a/psydac/feec/multipatch/fem_linear_operators.py +++ /dev/null @@ -1,218 +0,0 @@ -# coding: utf-8 - -from mpi4py import MPI - -from scipy.sparse import eye as sparse_id - -from psydac.linalg.basic import LinearOperator -from psydac.fem.basic import FemField - -#=============================================================================== -class FemLinearOperator( LinearOperator ): - """ - Linear operators with an additional Fem layer - """ - - def __init__( self, fem_domain=None, fem_codomain=None, matrix=None, sparse_matrix=None): - """ - we may store the matrix of the linear operator with different formats - :param matrix: stencil format - :param sparse_matrix: scipy sparse format - """ - assert fem_domain - self._fem_domain = fem_domain - if fem_codomain: - self._fem_codomain = fem_codomain - else: - self._fem_codomain = fem_domain - self._domain = self._fem_domain.coeff_space - self._codomain = self._fem_codomain.coeff_space - - self._matrix = matrix - self._sparse_matrix = sparse_matrix - - @property - def domain( self ): - return self._domain - - @property - def codomain( self ): - return self._codomain - - @property - def fem_domain( self ): - return self._fem_domain - - @property - def fem_codomain( self ): - return self._fem_codomain - - @property - def matrix( self ): - return self._matrix - - @property - def T(self): - return self.transpose() - - @property - def dtype( self ): - return self.domain.dtype - - def toarray(self): - return self._matrix.toarray() - #raise NotImplementedError('toarray() is not defined for FEMLinearOperators.') - - def tosparse(self): - return self._matrix.tosparse() - #raise NotImplementedError('tosparse() is not defined for FEMLinearOperators.') - - # ... - def transpose(self, conjugate=False): - raise NotImplementedError('Class does not provide a transpose() method') - - # ... - def to_sparse_matrix( self , **kwargs): - if self._sparse_matrix is not None: - return self._sparse_matrix - elif self._matrix is not None: - return self._matrix.tosparse() - else: - raise NotImplementedError('Class does not provide a get_sparse_matrix() method without a matrix') - - # ... - def __call__( self, f ): - if self._matrix is not None: - coeffs = self._matrix.dot(f.coeffs) - return FemField(self.fem_codomain, coeffs=coeffs) - else: - raise NotImplementedError('Class does not provide a __call__ method without a matrix') - - # ... - def dot( self, f_coeffs, out=None ): - # coeffs layer - if self._matrix is not None: - f = FemField(self.fem_domain, coeffs=f_coeffs) - return self(f).coeffs - else: - raise NotImplementedError('Class does not provide a dot method without a matrix') - - # ... - def __mul__(self, c): - return MultLinearOperator(c, self) - - # ... - def __add__(self, C): - assert isinstance(C, FemLinearOperator) - return SumLinearOperator(C, self) - - # ... - def __sub__(self, C): - assert isinstance(C, FemLinearOperator) - return SumLinearOperator(C, -self) - - # ... - def __neg__(self): - return MultLinearOperator(-1, self) - - -#============================================================================== -class ComposedLinearOperator( FemLinearOperator ): - """ - operator L = L_1 .. L_n - with L_i = self._operators[i-1] - (so, the last one is applied first, like in a product) - """ - def __init__( self, operators ): - n = len(operators) - assert all([isinstance(operators[i], FemLinearOperator) for i in range(n)]) - assert all([operators[i].fem_domain == operators[i+1].fem_codomain for i in range(n-1)]) - FemLinearOperator.__init__( - self, fem_domain=operators[-1].fem_domain, fem_codomain=operators[0].fem_codomain - ) - self._operators = operators - self._n = n - - # matrix not defined by matrix product because it could break the Stencil Matrix structure - - def to_sparse_matrix( self, **kwargs): - mat = self._operators[-1].to_sparse_matrix() - for i in range(2, self._n+1): - mat = self._operators[-i].to_sparse_matrix() * mat - return mat - - def __call__( self, f ): - v = self._operators[-1](f) - for i in range(2, self._n+1): - v = self._operators[-i](v) - return v - - def dot( self, f_coeffs, out=None ): - v_coeffs = self._operators[-1].dot(f_coeffs) - for i in range(2, self._n+1): - v_coeffs = self._operators[-i].dot(v_coeffs) - return v_coeffs - - -#============================================================================== -class IdLinearOperator( FemLinearOperator ): - - def __init__( self, V ): - FemLinearOperator.__init__(self, fem_domain=V) - - def to_sparse_matrix( self , **kwargs): - return sparse_id( self.fem_domain.nbasis ) - - def __call__( self, f ): - return f - - def dot( self, f_coeffs, out=None ): - return f_coeffs - -#============================================================================== -class SumLinearOperator( FemLinearOperator ): - - def __init__( self, B, A ): - assert isinstance(A, FemLinearOperator) - assert isinstance(B, FemLinearOperator) - assert B.fem_domain == A.fem_domain - assert B.fem_codomain == A.fem_codomain - FemLinearOperator.__init__( - self, fem_domain=A.fem_domain, fem_codomain=A.fem_codomain - ) - self._A = A - self._B = B - - def to_sparse_matrix( self, **kwargs): - return self._A.to_sparse_matrix() + self._B.to_sparse_matrix() - - def __call__( self, f ): - # fem layer - return self._B(f) + self._A(f) - - def dot( self, f_coeffs, out=None ): - # coeffs layer - return self._B.dot(f_coeffs) + self._A.dot(f_coeffs) - -#============================================================================== -class MultLinearOperator( FemLinearOperator ): - - def __init__( self, c, A ): - assert isinstance(A, FemLinearOperator) - FemLinearOperator.__init__( - self, fem_domain=A.fem_domain, fem_codomain=A.fem_codomain - ) - self._A = A - self._c = c - - def to_sparse_matrix( self, **kwargs): - return self._c * self._A.to_sparse_matrix() - - def __call__( self, f ): - # fem layer - return self._c * self._A(f) - - def dot( self, f_coeffs, out=None ): - # coeffs layer - return self._c * self._A.dot(f_coeffs) - diff --git a/psydac/feec/multipatch/operators.py b/psydac/feec/multipatch/operators.py deleted file mode 100644 index 0ee00210c..000000000 --- a/psydac/feec/multipatch/operators.py +++ /dev/null @@ -1,1248 +0,0 @@ -# coding: utf-8 - -# Conga operators on piecewise (broken) de Rham sequences - -from sympy import Tuple -from mpi4py import MPI -import os -import numpy as np - -from scipy.sparse import save_npz, load_npz -from scipy.sparse import kron, block_diag -from scipy.sparse.linalg import inv - -from sympde.topology import Boundary, Interface, Union -from sympde.topology import element_of, elements_of -from sympde.topology.space import ScalarFunction -from sympde.calculus import grad, dot, inner, rot, div -from sympde.calculus import laplace, bracket, convect -from sympde.calculus import jump, avg, Dn, minus, plus -from sympde.expr.expr import LinearForm, BilinearForm -from sympde.expr.expr import integral - -from psydac.core.bsplines import collocation_matrix, histopolation_matrix - -from psydac.api.discretization import discretize -from psydac.api.essential_bc import apply_essential_bc_stencil -from psydac.api.settings import PSYDAC_BACKENDS -from psydac.linalg.block import BlockVectorSpace, BlockVector, BlockLinearOperator -from psydac.linalg.stencil import StencilVector, StencilMatrix, StencilInterfaceMatrix -from psydac.linalg.solvers import inverse -from psydac.fem.basic import FemField - - -from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl, Projector_L2 -from psydac.feec.derivatives import Gradient_2D, ScalarCurl_2D -from psydac.feec.multipatch.fem_linear_operators import FemLinearOperator - - -def get_patch_index_from_face(domain, face): - """ Return the patch index of subdomain/boundary - - Parameters - ---------- - domain : - The Symbolic domain - - face : - A patch or a boundary of a patch - - Returns - ------- - i : - The index of a subdomain/boundary in the multipatch domain - """ - - if domain.mapping: - domain = domain.logical_domain - if face.mapping: - face = face.logical_domain - - domains = domain.interior.args - if isinstance(face, Interface): - raise NotImplementedError( - "This face is an interface, it has several indices -- I am a machine, I cannot choose. Help.") - elif isinstance(face, Boundary): - i = domains.index(face.domain) - else: - i = domains.index(face) - return i - - -def get_interface_from_corners(corner1, corner2, domain): - """ Return the interface between two corners from two different patches that correspond to a single (physical) vertex. - - Parameters - ---------- - corner1 : - The first corner of the 2D interface - - corner2 : - The second corner of the 2D interface - - domain : - The Symbolic domain - - Returns - ------- - interface: - The interface between two vertices - - """ - - interface = [] - interfaces = domain.interfaces - - if not isinstance(interfaces, Union): - interfaces = (interfaces,) - - for i in interfaces: - if i.plus.domain in [corner1.domain, corner2.domain]: - if i.minus.domain in [corner1.domain, corner2.domain]: - interface.append(i) - - bd1 = corner1.boundaries - bd2 = corner2.boundaries - - new_interface = [] - - for i in interface: - if i.minus in bd1 + bd2: - if i.plus in bd2 + bd1: - new_interface.append(i) - - if len(new_interface) == 1: - return new_interface[0] - if len(new_interface) > 1: - raise ValueError( - 'found more than one interface for the corners {} and {}'.format( - corner1, corner2)) - return None - - -def get_row_col_index(corner1, corner2, interface, axis, V1, V2): - """ Return the row and column index of a corner in the StencilInterfaceMatrix - for dofs of H1 type spaces - - Parameters - ---------- - corner1 : - The first corner of the 2D interface - - corner2 : - The second corner of the 2D interface - - interface : - The interface between the two corners - - axis : - Axis of the interface - - V1 : - Test Space - - V2 : - Trial Space - - Returns - ------- - index: - The StencilInterfaceMatrix index of the corner, it has the form (i1, i2, k1, k2) in 2D, - where (i1, i2) identifies the row and (k1, k2) the diagonal. - """ - start = V1.coeff_space.starts - end = V1.coeff_space.ends - degree = V2.degree - start_end = (start, end) - - row = [None] * len(start) - col = [0] * len(start) - - assert corner1.boundaries[0].axis == corner2.boundaries[0].axis - - for bd in corner1.boundaries: - row[bd.axis] = start_end[(bd.ext + 1) // 2][bd.axis] - - if interface is None and corner1.domain != corner2.domain: - bd = [i for i in corner1.boundaries if i.axis == axis][0] - if bd.ext == 1: - row[bd.axis] = degree[bd.axis] - - if interface is None: - return row + col - - axis = interface.axis - - if interface.minus.domain == corner1.domain: - if interface.minus.ext == -1: - row[axis] = 0 - else: - row[axis] = degree[axis] - else: - if interface.plus.ext == -1: - row[axis] = 0 - else: - row[axis] = degree[axis] - - if interface.minus.ext == interface.plus.ext: - pass - elif interface.minus.domain == corner1.domain: - if interface.minus.ext == -1: - col[axis] = degree[axis] - else: - col[axis] = -degree[axis] - else: - if interface.plus.ext == -1: - col[axis] = degree[axis] - else: - col[axis] = -degree[axis] - - return row + col - - -# =============================================================================== -def allocate_interface_matrix(corners, test_space, trial_space): - """ Allocate the interface matrix for a vertex shared by two patches - - Parameters - ---------- - corners: - The patch corners corresponding to the common shared vertex - - test_space: - The test space - - trial_space: - The trial space - - Returns - ------- - mat: - The interface matrix shared by two patches - """ - bi, bj = list(zip(*corners)) - permutation = np.arange(bi[0].domain.dim) - - flips = [] - k = 0 - while k < len(bi): - c1 = np.array(bi[k].coordinates) - c2 = np.array(bj[k].coordinates)[permutation] - flips.append( - np.array([-1 if d1 != d2 else 1 for d1, d2 in zip(c1, c2)])) - - if np.sum(abs(flips[0] - flips[-1])) != 0: - prod = [f1 * f2 for f1, f2 in zip(flips[0], flips[-1])] - while -1 in prod: - i1 = prod.index(-1) - if -1 in prod[i1 + 1:]: - i2 = i1 + 1 + prod[i1 + 1:].index(-1) - prod = prod[i2 + 1:] - permutation[i1], permutation[i2] = permutation[i2], permutation[i1] - k = -1 - flips = [] - else: - break - - k += 1 - - assert all(abs(flips[0] - i).sum() == 0 for i in flips) - cs = list(zip(*[i.coordinates for i in bi])) - axis = [all(i[0] == j for j in i) for i in cs].index(True) - ext = 1 if cs[axis][0] == 1 else -1 - s = test_space.get_assembly_grids( - )[axis].spans[-1 if ext == 1 else 0] - test_space.degree[axis] - - mat = StencilInterfaceMatrix( - trial_space.coeff_space, - test_space.coeff_space, - s, - s, - axis, - flip=flips[0], - permutation=list(permutation)) - return mat - -# =============================================================================== -# The following operators are not compatible with the changes in the Stencil format -# and their datatype does not allow for non-matching interfaces, but they might be -# useful for future implementations -# =============================================================================== - - -class ConformingProjection_V0(FemLinearOperator): - """ - Conforming projection from global broken V0 space to conforming global V0 space - Defined by averaging of interface dofs - - Parameters - ---------- - V0h: - The discrete space - - domain_h: - The discrete domain of the projector - - hom_bc : - Apply homogenous boundary conditions if True - - backend_language: - The backend used to accelerate the code - - storage_fn: - filename to store/load the operator sparse matrix - """ - # todo (MCP, 16.03.2021): - # - avoid discretizing a bilinear form - # - allow case without interfaces (single or multipatch) - - def __init__( - self, - V0h, - domain_h, - hom_bc=False, - backend_language='python', - storage_fn=None): - - FemLinearOperator.__init__(self, fem_domain=V0h) - - V0 = V0h.symbolic_space - domain = V0.domain - self.symbolic_domain = domain - - if storage_fn and os.path.exists(storage_fn): - print( - "[ConformingProjection_V0] loading operator sparse matrix from " + - storage_fn) - self._sparse_matrix = load_npz(storage_fn) - - else: - # assemble the operator matrix - u, v = elements_of(V0, names='u, v') - expr = u * v # dot(u,v) - - Interfaces = domain.interfaces # note: interfaces does not include the boundary - # this penalization is for an H1-conforming space - expr_I = (plus(u) - minus(u)) * (plus(v) - minus(v)) - - a = BilinearForm((u, v), integral(domain, expr) + - integral(Interfaces, expr_I)) - # print('[[ forcing python backend for ConformingProjection_V0]] ') - # backend_language = 'python' - ah = discretize( - a, domain_h, [ - V0h, V0h], backend=PSYDAC_BACKENDS[backend_language]) - - # self._A = ah.assemble() - self._A = ah.forms[0]._matrix - - spaces = self._A.domain.spaces - - if isinstance(Interfaces, Interface): - Interfaces = (Interfaces, ) - - for b1 in self._A.blocks: - for A in b1: - if A is None: - continue - A[:, :, :, :] = 0 - - indices = [slice(None, None)] * domain.dim + [0] * domain.dim - - for i in range(len(self._A.blocks)): - self._A[i, i][tuple(indices)] = 1 - - for I in Interfaces: - - axis = I.axis - i_minus = get_patch_index_from_face(domain, I.minus) - i_plus = get_patch_index_from_face(domain, I.plus) - - sp_minus = spaces[i_minus] - sp_plus = spaces[i_plus] - - s_minus = sp_minus.starts[axis] - e_minus = sp_minus.ends[axis] - - s_plus = sp_plus.starts[axis] - e_plus = sp_plus.ends[axis] - - d_minus = V0h.spaces[i_minus].degree[axis] - d_plus = V0h.spaces[i_plus].degree[axis] - - indices = [slice(None, None)] * domain.dim + [0] * domain.dim - - minus_ext = I.minus.ext - plus_ext = I.plus.ext - - if minus_ext == 1: - indices[axis] = e_minus - else: - indices[axis] = s_minus - self._A[i_minus, i_minus][tuple(indices)] = 1 / 2 - - if plus_ext == 1: - indices[axis] = e_plus - else: - indices[axis] = s_plus - - self._A[i_plus, i_plus][tuple(indices)] = 1 / 2 - - if plus_ext == minus_ext: - if minus_ext == 1: - indices[axis] = d_minus - else: - indices[axis] = s_minus - - self._A[i_minus, i_plus][tuple(indices)] = 1 / 2 - - if plus_ext == 1: - indices[axis] = d_plus - else: - indices[axis] = s_plus - - self._A[i_plus, i_minus][tuple(indices)] = 1 / 2 - - else: - if minus_ext == 1: - indices[axis] = d_minus - else: - indices[axis] = s_minus - - if plus_ext == 1: - indices[domain.dim + axis] = d_plus - else: - indices[domain.dim + axis] = -d_plus - - self._A[i_minus, i_plus][tuple(indices)] = 1 / 2 - - if plus_ext == 1: - indices[axis] = d_plus - else: - indices[axis] = s_plus - - if minus_ext == 1: - indices[domain.dim + axis] = d_minus - else: - indices[domain.dim + axis] = -d_minus - - self._A[i_plus, i_minus][tuple(indices)] = 1 / 2 - - domain = domain.logical_domain - corner_blocks = {} - for c in domain.corners: - for b1 in c.corners: - i = get_patch_index_from_face(domain, b1.domain) - for b2 in c.corners: - j = get_patch_index_from_face(domain, b2.domain) - if (i, j) in corner_blocks: - corner_blocks[i, j] += [(b1, b2)] - else: - corner_blocks[i, j] = [(b1, b2)] - - for c in domain.corners: - if len(c) == 2: - continue - for b1 in c.corners: - i = get_patch_index_from_face(domain, b1.domain) - for b2 in c.corners: - j = get_patch_index_from_face(domain, b2.domain) - interface = get_interface_from_corners(b1, b2, domain) - axis = None - if self._A[i, j] is None: - self._A[i, j] = allocate_interface_matrix( - corner_blocks[i, j], V0h.spaces[i], V0h.spaces[j]) - - if i != j and self._A[i, j]: - axis = self._A[i, j]._dim - index = get_row_col_index( - b1, b2, interface, axis, V0h.spaces[i], V0h.spaces[j]) - self._A[i, j][tuple(index)] = 1 / len(c) - - if hom_bc: - for bn in domain.boundary: - self.set_homogenous_bc(bn) - - self._matrix = self._A - self._sparse_matrix = self._matrix.tosparse() # self._sparse_matrix - - if storage_fn: - print( - "[ConformingProjection_V0] storing operator sparse matrix in " + - storage_fn) - save_npz(storage_fn, self._sparse_matrix) - - def set_homogenous_bc(self, boundary, rhs=None): - domain = self.symbolic_domain - Vh = self.fem_domain - if domain.mapping: - domain = domain.logical_domain - if boundary.mapping: - boundary = boundary.logical_domain - - corners = domain.corners - i = get_patch_index_from_face(domain, boundary) - if rhs: - apply_essential_bc_stencil( - rhs[i], axis=boundary.axis, ext=boundary.ext, order=0) - for j in range(len(domain)): - if self._A[i, j] is None: - continue - apply_essential_bc_stencil( - self._A[i, j], axis=boundary.axis, ext=boundary.ext, order=0) - - for c in corners: - faces = [f for b in c.corners for f in b.boundaries] - if len(c) == 2: - continue - if boundary in faces: - for b1 in c.corners: - i = get_patch_index_from_face(domain, b1.domain) - for b2 in c.corners: - j = get_patch_index_from_face(domain, b2.domain) - interface = get_interface_from_corners(b1, b2, domain) - axis = None - if i != j: - axis = self._A[i, j].dim - index = get_row_col_index( - b1, b2, interface, axis, Vh.spaces[i], Vh.spaces[j]) - self._A[i, j][tuple(index)] = 0. - - if i == j and rhs: - rhs[i][tuple(index[:2])] = 0. - -# =============================================================================== - - -class ConformingProjection_V1(FemLinearOperator): - """ - Conforming projection from global broken V1 space to conforming V1 global space - - proj.dot(v) returns the conforming projection of v, computed by solving linear system - - Parameters - ---------- - V1h: - The discrete space - - domain_h: - The discrete domain of the projector - - hom_bc : - Apply homogenous boundary conditions if True - - backend_language: - The backend used to accelerate the code - - storage_fn: - filename to store/load the operator sparse matrix - """ - # todo (MCP, 16.03.2021): - # - avoid discretizing a bilinear form - # - allow case without interfaces (single or multipatch) - - def __init__( - self, - V1h, - domain_h, - hom_bc=False, - backend_language='python', - storage_fn=None): - - FemLinearOperator.__init__(self, fem_domain=V1h) - - V1 = V1h.symbolic_space - domain = V1.domain - self.symbolic_domain = domain - - if storage_fn and os.path.exists(storage_fn): - print( - "[ConformingProjection_V1] loading operator sparse matrix from " + - storage_fn) - self._sparse_matrix = load_npz(storage_fn) - - else: - # assemble the operator matrix - u, v = elements_of(V1, names='u, v') - expr = dot(u, v) - # - Interfaces = domain.interfaces # note: interfaces does not include the boundary - # this penalization is for an H1-conforming space - expr_I = dot(plus(u) - minus(u), plus(v) - minus(v)) - - a = BilinearForm((u, v), integral(domain, expr) + - integral(Interfaces, expr_I)) - # print('[[ forcing python backend for ConformingProjection_V1]] ') - # backend_language = 'python' - ah = discretize( - a, domain_h, [ - V1h, V1h], backend=PSYDAC_BACKENDS[backend_language]) - # - # # self._A = ah.assemble() - self._A = ah.forms[0]._matrix - # C1 = V1h.coeff_space - # self._A = BlockLinearOperator(C1, C1) - - for b1 in self._A.blocks: - for b2 in b1: - if b2 is None: - continue - for b3 in b2.blocks: - for A in b3: - if A is None: - continue - A[:, :, :, :] = 0 - - spaces = self._A.domain.spaces - - if isinstance(Interfaces, Interface): - Interfaces = (Interfaces, ) - - indices = [slice(None, None)] * domain.dim + [0] * domain.dim - - for i in range(len(self._A.blocks)): - self._A[i, i][0, 0][tuple(indices)] = 1 - self._A[i, i][1, 1][tuple(indices)] = 1 - - # empty list if no interfaces ? - if Interfaces is not None: - - for I in Interfaces: - - i_minus = get_patch_index_from_face(domain, I.minus) - i_plus = get_patch_index_from_face(domain, I.plus) - - indices = [slice(None, None)] * \ - domain.dim + [0] * domain.dim - - sp1 = spaces[i_minus] - sp2 = spaces[i_plus] - - s11 = sp1.spaces[0].starts[I.axis] - e11 = sp1.spaces[0].ends[I.axis] - s12 = sp1.spaces[1].starts[I.axis] - e12 = sp1.spaces[1].ends[I.axis] - - s21 = sp2.spaces[0].starts[I.axis] - e21 = sp2.spaces[0].ends[I.axis] - s22 = sp2.spaces[1].starts[I.axis] - e22 = sp2.spaces[1].ends[I.axis] - - d11 = V1h.spaces[i_minus].spaces[0].degree[I.axis] - d12 = V1h.spaces[i_minus].spaces[1].degree[I.axis] - - d21 = V1h.spaces[i_plus].spaces[0].degree[I.axis] - d22 = V1h.spaces[i_plus].spaces[1].degree[I.axis] - - s_minus = [s11, s12] - e_minus = [e11, e12] - - s_plus = [s21, s22] - e_plus = [e21, e22] - - d_minus = [d11, d12] - d_plus = [d21, d22] - - minus_ext = I.minus.ext - plus_ext = I.plus.ext - - axis = I.axis - for k in range(domain.dim): - if k == I.axis: - continue - - if minus_ext == 1: - indices[axis] = e_minus[k] - else: - indices[axis] = s_minus[k] - self._A[i_minus, i_minus][k, k][tuple(indices)] = 1 / 2 - - if plus_ext == 1: - indices[axis] = e_plus[k] - else: - indices[axis] = s_plus[k] - - self._A[i_plus, i_plus][k, k][tuple(indices)] = 1 / 2 - - if plus_ext == minus_ext: - if minus_ext == 1: - indices[axis] = d_minus[k] - else: - indices[axis] = s_minus[k] - - self._A[i_minus, i_plus][k, k][tuple( - indices)] = 1 / 2 * I.direction - - if plus_ext == 1: - indices[axis] = d_plus[k] - else: - indices[axis] = s_plus[k] - - self._A[i_plus, i_minus][k, k][tuple( - indices)] = 1 / 2 * I.direction - - else: - if minus_ext == 1: - indices[axis] = d_minus[k] - else: - indices[axis] = s_minus[k] - - if plus_ext == 1: - indices[domain.dim + axis] = d_plus[k] - else: - indices[domain.dim + axis] = -d_plus[k] - - self._A[i_minus, i_plus][k, k][tuple( - indices)] = 1 / 2 * I.direction - - if plus_ext == 1: - indices[axis] = d_plus[k] - else: - indices[axis] = s_plus[k] - - if minus_ext == 1: - indices[domain.dim + axis] = d_minus[k] - else: - indices[domain.dim + axis] = -d_minus[k] - - self._A[i_plus, i_minus][k, k][tuple( - indices)] = 1 / 2 * I.direction - - if hom_bc: - for bn in domain.boundary: - self.set_homogenous_bc(bn) - - self._matrix = self._A - self._sparse_matrix = self._matrix.tosparse() - - if storage_fn: - print( - "[ConformingProjection_V1] storing operator sparse matrix in " + - storage_fn) - save_npz(storage_fn, self._sparse_matrix) - - def set_homogenous_bc(self, boundary): - domain = self.symbolic_domain - Vh = self.fem_domain - - i = get_patch_index_from_face(domain, boundary) - axis = boundary.axis - ext = boundary.ext - for j in range(len(domain)): - if self._A[i, j] is None: - continue - apply_essential_bc_stencil( - self._A[i, j][1 - axis, 1 - axis], axis=axis, ext=ext, order=0) - - -# =============================================================================== -def get_K0_and_K0_inv(V0h, uniform_patches=False): - """ - Compute the change of basis matrices K0 and K0^{-1} in V0h. - - With - K0_ij = sigma^0_i(B_j) = B_jx(n_ix) * B_jy(n_iy) - where sigma_i is the geometric (interpolation) dof - and B_j is the tensor-product B-spline - """ - if uniform_patches: - print(' [[WARNING -- hack in get_K0_and_K0_inv: using copies of 1st-patch matrices in every patch ]] ') - - V0 = V0h.symbolic_space # VOh is FemSpace - domain = V0.domain - K0_blocks = [] - K0_inv_blocks = [] - for k, D in enumerate(domain.interior): - if uniform_patches and k > 0: - K0_k = K0_blocks[0].copy() - K0_inv_k = K0_inv_blocks[0].copy() - - else: - V0_k = V0h.spaces[k] # fem space on patch k: (TensorFemSpace) - K0_k_factors = [None, None] - for d in [0, 1]: - # 1d fem space alond dim d (SplineSpace) - V0_kd = V0_k.spaces[d] - K0_k_factors[d] = collocation_matrix( - knots=V0_kd.knots, - degree=V0_kd.degree, - periodic=V0_kd.periodic, - normalization=V0_kd.basis, - xgrid=V0_kd.greville - ) - K0_k = kron(*K0_k_factors) - K0_k.eliminate_zeros() - K0_inv_k = inv(K0_k.tocsc()) - K0_inv_k.eliminate_zeros() - - K0_blocks.append(K0_k) - K0_inv_blocks.append(K0_inv_k) - K0 = block_diag(K0_blocks) - K0_inv = block_diag(K0_inv_blocks) - return K0, K0_inv - - -# =============================================================================== -def get_K1_and_K1_inv(V1h, uniform_patches=False): - """ - Compute the change of basis matrices K1 and K1^{-1} in Hcurl space V1h. - - With - K1_ij = sigma^1_i(B_j) = int_{e_ix}(M_jx) * B_jy(n_iy) - if i = horizontal edge [e_ix, n_iy] and j = (M_jx o B_jy) x-oriented MoB spline - or - = B_jx(n_ix) * int_{e_iy}(M_jy) - if i = vertical edge [n_ix, e_iy] and j = (B_jx o M_jy) y-oriented BoM spline - (above, 'o' denotes tensor-product for functions) - """ - if uniform_patches: - print(' [[WARNING -- hack in get_K1_and_K1_inv: using copies of 1st-patch matrices in every patch ]] ') - - V1 = V1h.symbolic_space # V1h is FemSpace - domain = V1.domain - K1_blocks = [] - K1_inv_blocks = [] - for k, D in enumerate(domain.interior): - if uniform_patches and k > 0: - K1_k = K1_blocks[0].copy() - K1_inv_k = K1_inv_blocks[0].copy() - - else: - # fem space on patch k: - V1_k = V1h.spaces[k] - K1_k_blocks = [] - for c in [0, 1]: # dim of component - # fem space for comp. dc (TensorFemSpace) - V1_kc = V1_k.spaces[c] - K1_kc_factors = [None, None] - for d in [0, 1]: # dim of variable - # 1d fem space for comp c alond dim d (SplineSpace) - V1_kcd = V1_kc.spaces[d] - if c == d: - K1_kc_factors[d] = histopolation_matrix( - knots=V1_kcd.knots, - degree=V1_kcd.degree, - periodic=V1_kcd.periodic, - normalization=V1_kcd.basis, - xgrid=V1_kcd.ext_greville - ) - else: - K1_kc_factors[d] = collocation_matrix( - knots=V1_kcd.knots, - degree=V1_kcd.degree, - periodic=V1_kcd.periodic, - normalization=V1_kcd.basis, - xgrid=V1_kcd.greville - ) - K1_kc = kron(*K1_kc_factors) - K1_kc.eliminate_zeros() - K1_k_blocks.append(K1_kc) - K1_k = block_diag(K1_k_blocks) - K1_k.eliminate_zeros() - K1_inv_k = inv(K1_k.tocsc()) - K1_inv_k.eliminate_zeros() - - K1_blocks.append(K1_k) - K1_inv_blocks.append(K1_inv_k) - - K1 = block_diag(K1_blocks) - K1_inv = block_diag(K1_inv_blocks) - return K1, K1_inv - - -# #=============================================================================== -# def get_M_and_M_inv(Vh, subdomains_h, is_scalar, backend_language='python'): -# """ -# compute the mass matrix M and M^{-1} in multipatch space Vh -# DOES NOT WORK -- SHOULD WE HAVE THE POSSIBILITY OF DOING THAT ? -# """ -# from pprint import pprint -# -# V = Vh.symbolic_space # VOh is FemSpace -# domain = V.domain -# M_blocks = [] -# M_inv_blocks = [] -# -# # print('type(domain_h) = ', type(domain_h)) -# # -# # print('type(domain_h._patches) = ', type(domain_h._patches)) -# # print('len(domain_h._patches) = ', len(domain_h._patches)) -# # -# # mappings = domain_h.mappings -# # print('type(mappings) = ', type(mappings)) -# # print('len(mappings) = ', len(mappings)) -# # -# # mappings_list = list(mappings.values()) -# # print('len(mappings_list) = ', len(mappings_list)) -# # -# # print('type(mappings_list[0]) = ', type(mappings_list[0])) -# -# for k, Dh_k in enumerate(subdomains_h): -# -# print('k = ', k) -# print('type(Dh_k) = ', type(Dh_k)) -# # print('Dh = ', Dh) -# D_k = domain.interior[k] -# -# # exit() -# -# # for k, D in enumerate(domain.interior): -# -# V_k = V.spaces[k] -# Vh_k = Vh.spaces[k] -# -# # print(type(domain_h)) -# # -# # pprint(dir(domain_h)) -# # -# # -# # print(len(domain_h._patches)) -# # exit() -# # Dh_k = domain_h.spaces[k] # fem space on patch k: (TensorFemSpace) -# u, v = elements_of(V_k, names='u, v') -# if is_scalar: -# expr = u*v -# else: -# expr = dot(u,v) -# a_k = BilinearForm((u,v), integral(D_k, expr)) -# a_kh = discretize(a_k, Dh_k, [Vh_k, Vh_k], backend=PSYDAC_BACKENDS[backend_language]) # 'pyccel-gcc']) -# -# M_k = a_kh.assemble().toarray() -# M_k.eliminate_zeros() -# M_inv_k = inv(M_k.tocsc()) -# M_inv_k.eliminate_zeros() -# -# M_blocks.append(M_k) -# M_inv_blocks.append(M_inv_k) -# M = block_diag(M_blocks) -# M_inv = block_diag(M_inv_blocks) -# return M, M_inv - -# =============================================================================== -class HodgeOperator(FemLinearOperator): - """ - Change of basis operator: dual basis -> primal basis - - self._matrix: matrix of the primal Hodge = this is the mass matrix ! - self.dual_Hodge_matrix: this is the INVERSE mass matrix - - Parameters - ---------- - Vh: - The discrete space - - domain_h: - The discrete domain of the projector - - metric : - the metric of the de Rham complex - - backend_language: - The backend used to accelerate the code - - load_dir: - storage files for the primal and dual Hodge sparse matrice - - load_space_index: - the space index in the derham sequence - - Notes - ----- - Either we use a storage, or these matrices are only computed on demand - # todo: we compute the sparse matrix when to_sparse_matrix is called -- but never the stencil matrix (should be fixed...) - We only support the identity metric, this implies that the dual Hodge is the inverse of the primal one. - # todo: allow for non-identity metrics - """ - - def __init__( - self, - Vh, - domain_h, - metric='identity', - backend_language='python', - load_dir=None, - load_space_index=''): - - FemLinearOperator.__init__(self, fem_domain=Vh) - self._domain_h = domain_h - self._backend_language = backend_language - self._dual_Hodge_sparse_matrix = None - - assert metric == 'identity' - self._metric = metric - - if load_dir and isinstance(load_dir, str): - if not os.path.exists(load_dir): - os.makedirs(load_dir) - assert str(load_space_index) in ['0', '1', '2', '3'] - primal_Hodge_storage_fn = load_dir + \ - '/H{}_m.npz'.format(load_space_index) - dual_Hodge_storage_fn = load_dir + \ - '/dH{}_m.npz'.format(load_space_index) - - primal_Hodge_is_stored = os.path.exists(primal_Hodge_storage_fn) - dual_Hodge_is_stored = os.path.exists(dual_Hodge_storage_fn) - if dual_Hodge_is_stored: - assert primal_Hodge_is_stored - print( - " ... loading dual Hodge sparse matrix from " + - dual_Hodge_storage_fn) - self._dual_Hodge_sparse_matrix = load_npz( - dual_Hodge_storage_fn) - print( - "[HodgeOperator] loading primal Hodge sparse matrix from " + - primal_Hodge_storage_fn) - self._sparse_matrix = load_npz(primal_Hodge_storage_fn) - else: - assert not primal_Hodge_is_stored - print( - "[HodgeOperator] assembling both sparse matrices for storage...") - self.assemble_primal_Hodge_matrix() - print( - "[HodgeOperator] storing primal Hodge sparse matrix in " + - primal_Hodge_storage_fn) - save_npz(primal_Hodge_storage_fn, self._sparse_matrix) - self.assemble_dual_Hodge_matrix() - print( - "[HodgeOperator] storing dual Hodge sparse matrix in " + - dual_Hodge_storage_fn) - save_npz(dual_Hodge_storage_fn, self._dual_Hodge_sparse_matrix) - else: - # matrices are not stored, we will probably compute them later - pass - - def to_sparse_matrix(self): - """ - the Hodge matrix is the patch-wise multi-patch mass matrix - it is not stored by default but assembled on demand - """ - - if (self._sparse_matrix is not None) or (self._matrix is not None): - return FemLinearOperator.to_sparse_matrix(self) - - self.assemble_primal_Hodge_matrix() - - return self._sparse_matrix - - def assemble_primal_Hodge_matrix(self): - """ - the Hodge matrix is the patch-wise multi-patch mass matrix - it is not stored by default but assembled on demand - """ - - if self._matrix is None: - Vh = self.fem_domain - assert Vh == self.fem_codomain - - V = Vh.symbolic_space - domain = V.domain - # domain_h = V0h.domain: would be nice... - u, v = elements_of(V, names='u, v') - - if isinstance(u, ScalarFunction): - expr = u * v - else: - expr = dot(u, v) - - a = BilinearForm((u, v), integral(domain, expr)) - ah = discretize(a, self._domain_h, [ - Vh, Vh], backend=PSYDAC_BACKENDS[self._backend_language]) - - self._matrix = ah.assemble() # Mass matrix in stencil format - self._sparse_matrix = self._matrix.tosparse() - - def get_dual_Hodge_sparse_matrix(self): - if self._dual_Hodge_sparse_matrix is None: - self.assemble_dual_Hodge_matrix() - - return self._dual_Hodge_sparse_matrix - - def assemble_dual_Hodge_matrix(self): - """ - the dual Hodge matrix is the patch-wise inverse of the multi-patch mass matrix - it is not stored by default but computed on demand, by local (patch-wise) inversion of the mass matrix - """ - - if self._dual_Hodge_sparse_matrix is None: - if not self._matrix: - self.assemble_primal_Hodge_matrix() - - M = self._matrix # mass matrix of the (primal) basis - nrows = M.n_block_rows - ncols = M.n_block_cols - - inv_M_blocks = [] - for i in range(nrows): - Mii = M[i, i].tosparse() - inv_Mii = inv(Mii.tocsc()) - inv_Mii.eliminate_zeros() - inv_M_blocks.append(inv_Mii) - - inv_M = block_diag(inv_M_blocks) - self._dual_Hodge_sparse_matrix = inv_M - -# ============================================================================== - - -class BrokenGradient_2D(FemLinearOperator): - - def __init__(self, V0h, V1h): - - FemLinearOperator.__init__(self, fem_domain=V0h, fem_codomain=V1h) - - D0s = [Gradient_2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] - - self._matrix = BlockLinearOperator(self.domain, self.codomain, blocks={ - (i, i): D0i._matrix for i, D0i in enumerate(D0s)}) - - def transpose(self, conjugate=False): - # todo (MCP): define as the dual differential operator - return BrokenTransposedGradient_2D(self.fem_domain, self.fem_codomain) - -# ============================================================================== - - -class BrokenTransposedGradient_2D(FemLinearOperator): - - def __init__(self, V0h, V1h): - - FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V0h) - - D0s = [Gradient_2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] - - self._matrix = BlockLinearOperator(self.domain, self.codomain, blocks={ - (i, i): D0i._matrix.T for i, D0i in enumerate(D0s)}) - - def transpose(self, conjugate=False): - # todo (MCP): discard - return BrokenGradient_2D(self.fem_codomain, self.fem_domain) - - -# ============================================================================== -class BrokenScalarCurl_2D(FemLinearOperator): - def __init__(self, V1h, V2h): - - FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V2h) - - D1s = [ScalarCurl_2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] - - self._matrix = BlockLinearOperator(self.domain, self.codomain, blocks={ - (i, i): D1i._matrix for i, D1i in enumerate(D1s)}) - - def transpose(self, conjugate=False): - return BrokenTransposedScalarCurl_2D( - V1h=self.fem_domain, V2h=self.fem_codomain) - - -# ============================================================================== -class BrokenTransposedScalarCurl_2D(FemLinearOperator): - - def __init__(self, V1h, V2h): - - FemLinearOperator.__init__(self, fem_domain=V2h, fem_codomain=V1h) - - D1s = [ScalarCurl_2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] - - self._matrix = BlockLinearOperator(self.domain, self.codomain, blocks={ - (i, i): D1i._matrix.T for i, D1i in enumerate(D1s)}) - - def transpose(self, conjugate=False): - return BrokenScalarCurl_2D(V1h=self.fem_codomain, V2h=self.fem_domain) - - -# ============================================================================== - -# def multipatch_Moments_Hcurl(f, V1h, domain_h): - -def ortho_proj_Hcurl(EE, V1h, domain_h, M1, backend_language='python'): - """ - return orthogonal projection of E on V1h, given M1 the mass matrix - """ - assert isinstance(EE, Tuple) - V1 = V1h.symbolic_space - v = element_of(V1, name='v') - l = LinearForm(v, integral(V1.domain, dot(v, EE))) - lh = discretize( - l, - domain_h, - V1h, - backend=PSYDAC_BACKENDS[backend_language]) - b = lh.assemble() - M1_inv = inverse(M1.mat(), 'pcg', pc='jacobi', tol=1e-10) - sol_coeffs = M1_inv @ b - - return FemField(V1h, coeffs=sol_coeffs) - -# ============================================================================== - - -class Multipatch_Projector_H1: - """ - to apply the H1 projection (2D) on every patch - """ - - def __init__(self, V0h): - - self._P0s = [Projector_H1(V) for V in V0h.spaces] - self._V0h = V0h # multipatch Fem Space - - def __call__(self, funs_log): - """ - project a list of functions given in the logical domain - """ - u0s = [P(fun) for P, fun, in zip(self._P0s, funs_log)] - - u0_coeffs = BlockVector(self._V0h.coeff_space, - blocks=[u0j.coeffs for u0j in u0s]) - - return FemField(self._V0h, coeffs=u0_coeffs) - -# ============================================================================== - - -class Multipatch_Projector_Hcurl: - - """ - to apply the Hcurl projection (2D) on every patch - """ - - def __init__(self, V1h, nquads=None): - - self._P1s = [Projector_Hcurl(V, nquads=nquads) for V in V1h.spaces] - self._V1h = V1h # multipatch Fem Space - - def __call__(self, funs_log): - """ - project a list of functions given in the logical domain - """ - E1s = [P(fun) for P, fun, in zip(self._P1s, funs_log)] - - E1_coeffs = BlockVector(self._V1h.coeff_space, - blocks=[E1j.coeffs for E1j in E1s]) - - return FemField(self._V1h, coeffs=E1_coeffs) - -# ============================================================================== - - -class Multipatch_Projector_L2: - - """ - to apply the L2 projection (2D) on every patch - """ - - def __init__(self, V2h, nquads=None): - - self._P2s = [Projector_L2(V, nquads=nquads) for V in V2h.spaces] - self._V2h = V2h # multipatch Fem Space - - def __call__(self, funs_log): - """ - project a list of functions given in the logical domain - """ - B2s = [P(fun) for P, fun, in zip(self._P2s, funs_log)] - - B2_coeffs = BlockVector(self._V2h.coeff_space, - blocks=[B2j.coeffs for B2j in B2s]) - - return FemField(self._V2h, coeffs=B2_coeffs) diff --git a/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py b/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py index bb1b09004..a64c3a190 100644 --- a/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py +++ b/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py @@ -1,12 +1,12 @@ # coding: utf-8 import numpy as np +import pytest from psydac.feec.multipatch.examples.hcurl_source_pbms_conga_2d import solve_hcurl_source_pbm from psydac.feec.multipatch.examples.hcurl_eigen_pbms_conga_2d import hcurl_solve_eigen_pbm from psydac.feec.multipatch.examples.hcurl_eigen_pbms_dg_2d import hcurl_solve_eigen_pbm_dg -from psydac.feec.multipatch.examples.timedomain_maxwell import solve_td_maxwell_pbm - +# from psydac.feec.multipatch.examples.timedomain_maxwell import solve_td_maxwell_pbm def test_time_harmonic_maxwell_pretzel_f(): nc = 4 @@ -31,7 +31,6 @@ def test_time_harmonic_maxwell_pretzel_f(): assert abs(diags["err"] - 0.007201508128407582) < 1e-10 - def test_time_harmonic_maxwell_pretzel_f_nc(): deg = 2 nc = np.array([8, 8, 8, 8, 8, 4, 4, 4, 4, @@ -54,8 +53,7 @@ def test_time_harmonic_maxwell_pretzel_f_nc(): source_proj=source_proj, backend_language='pyccel-gcc') - assert abs(diags["err"] - 0.004849165663310541) < 1e-10 - + assert abs(diags["err"] - 0.004849225522124346) < 1e-7 def test_maxwell_eigen_curved_L_shape(): domain_name = 'curved_L_shape' @@ -99,7 +97,6 @@ def test_maxwell_eigen_curved_L_shape(): assert abs(error - 0.01291539899483907) < 1e-10 - def test_maxwell_eigen_curved_L_shape_nc(): domain_name = 'curved_L_shape' domain = [[1, 3], [0, np.pi / 4]] @@ -144,7 +141,6 @@ def test_maxwell_eigen_curved_L_shape_nc(): assert abs(error - 0.010504876643873904) < 1e-10 - def test_maxwell_eigen_curved_L_shape_dg(): domain_name = 'curved_L_shape' domain = [[1, 3], [0, np.pi / 4]] @@ -187,7 +183,7 @@ def test_maxwell_eigen_curved_L_shape_dg(): assert abs(error - 0.035139029534570064) < 1e-10 - +@pytest.mark.skip(reason="need to adapt notation") def test_maxwell_timedomain(): solve_td_maxwell_pbm(nc = 4, deg = 2, final_time = 2, domain_name = 'square_2') diff --git a/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py b/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py index 7a0d94cbb..fe6aed82f 100644 --- a/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py +++ b/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py @@ -2,7 +2,6 @@ from psydac.feec.multipatch.examples.h1_source_pbms_conga_2d import solve_h1_source_pbm - def test_poisson_pretzel_f(): source_type = 'manu_poisson_2' @@ -21,7 +20,6 @@ def test_poisson_pretzel_f(): assert abs(l2_error - 1.0585687717792318e-05) < 1e-10 - def test_poisson_pretzel_f_nc(): source_type = 'manu_poisson_2' diff --git a/psydac/feec/multipatch/utils_conga_2d.py b/psydac/feec/multipatch/utils_conga_2d.py index 351511d5e..b457f0232 100644 --- a/psydac/feec/multipatch/utils_conga_2d.py +++ b/psydac/feec/multipatch/utils_conga_2d.py @@ -5,70 +5,41 @@ from sympy import lambdify from sympde.topology import Derham +from sympde.topology.callable_mapping import BasicCallableMapping from psydac.api.settings import PSYDAC_BACKENDS from psydac.feec.pull_push import pull_2d_h1, pull_2d_hcurl, pull_2d_l2 -from psydac.feec.multipatch.api import discretize -from psydac.feec.multipatch.utilities import time_count # , export_sol, import_sol +from psydac.api.discretization import discretize +from psydac.feec.multipatch.utilities import time_count from psydac.linalg.utilities import array_to_psydac from psydac.fem.basic import FemField from psydac.fem.plotting_utilities import get_plotting_grid, get_grid_quad_weights, get_grid_vals +from scipy.sparse import kron, block_diag +from psydac.core.bsplines import collocation_matrix, histopolation_matrix +from psydac.linalg.solvers import inverse -# commuting projections on the physical domain (should probably be in the -# interface) -def P0_phys(f_phys, P0, domain, mappings_list): - f = lambdify(domain.coordinates, f_phys) - f_log = [pull_2d_h1(f, m.get_callable_mapping()) for m in mappings_list] - return P0(f_log) - - -def P1_phys(f_phys, P1, domain, mappings_list): - f_x = lambdify(domain.coordinates, f_phys[0]) - f_y = lambdify(domain.coordinates, f_phys[1]) - f_log = [pull_2d_hcurl([f_x, f_y], m.get_callable_mapping()) - for m in mappings_list] - return P1(f_log) - - -def P2_phys(f_phys, P2, domain, mappings_list): - f = lambdify(domain.coordinates, f_phys) - f_log = [pull_2d_l2(f, m.get_callable_mapping()) for m in mappings_list] - return P2(f_log) # commuting projections on the physical domain (should probably be in the # interface) - - -def P_phys_h1(f_phys, P0, domain, mappings_list): +def P0_phys(f_phys, P0, domain): f = lambdify(domain.coordinates, f_phys) - if len(mappings_list) == 1: - m = mappings_list[0] - f_log = pull_2d_h1(f, m) - else: - f_log = [pull_2d_h1(f, m) for m in mappings_list] - return P0(f_log) + return P0(f) -def P_phys_hcurl(f_phys, P1, domain, mappings_list): - f_x = lambdify(domain.coordinates, f_phys[0]) - f_y = lambdify(domain.coordinates, f_phys[1]) - f_log = [pull_2d_hcurl([f_x, f_y], m) for m in mappings_list] - return P1(f_log) - -def P_phys_hdiv(f_phys, P1, domain, mappings_list): +def P1_phys(f_phys, P1, domain): f_x = lambdify(domain.coordinates, f_phys[0]) f_y = lambdify(domain.coordinates, f_phys[1]) - f_log = [pull_2d_hdiv([f_x, f_y], m) for m in mappings_list] - return P1(f_log) + + return P1([f_x, f_y]) -def P_phys_l2(f_phys, P2, domain, mappings_list): +def P2_phys(f_phys, P2, domain): f = lambdify(domain.coordinates, f_phys) - f_log = [pull_2d_l2(f, m) for m in mappings_list] - return P2(f_log) + + return P2(f) def get_kind(space='V*'): @@ -83,6 +54,142 @@ def get_kind(space='V*'): raise ValueError(space) return kind +# =============================================================================== +def get_K0_and_K0_inv(V0h, uniform_patches=False): + """ + Compute the change of basis matrices K0 and K0^{-1} in V0h. + + With + K0_ij = sigma^0_i(B_j) = B_jx(n_ix) * B_jy(n_iy) + where sigma_i is the geometric (interpolation) dof + and B_j is the tensor-product B-spline + """ + if uniform_patches: + print(' [[WARNING -- hack in get_K0_and_K0_inv: using copies of 1st-patch matrices in every patch ]] ') + + V0 = V0h.symbolic_space # VOh is FemSpace + domain = V0.domain + K0_blocks = [] + K0_inv_blocks = [] + for k, D in enumerate(domain.interior): + if uniform_patches and k > 0: + K0_k = K0_blocks[0].copy() + K0_inv_k = K0_inv_blocks[0].copy() + + else: + V0_k = V0h.spaces[k] # fem space on patch k: (TensorFemSpace) + K0_k_factors = [None, None] + for d in [0, 1]: + # 1d fem space alond dim d (SplineSpace) + V0_kd = V0_k.spaces[d] + K0_k_factors[d] = collocation_matrix( + knots=V0_kd.knots, + degree=V0_kd.degree, + periodic=V0_kd.periodic, + normalization=V0_kd.basis, + xgrid=V0_kd.greville + ) + K0_k = kron(*K0_k_factors) + K0_k.eliminate_zeros() + K0_inv_k = inv(K0_k.tocsc()) + K0_inv_k.eliminate_zeros() + + K0_blocks.append(K0_k) + K0_inv_blocks.append(K0_inv_k) + K0 = block_diag(K0_blocks) + K0_inv = block_diag(K0_inv_blocks) + return K0, K0_inv + + +# =============================================================================== +def get_K1_and_K1_inv(V1h, uniform_patches=False): + """ + Compute the change of basis matrices K1 and K1^{-1} in Hcurl space V1h. + + With + K1_ij = sigma^1_i(B_j) = int_{e_ix}(M_jx) * B_jy(n_iy) + if i = horizontal edge [e_ix, n_iy] and j = (M_jx o B_jy) x-oriented MoB spline + or + = B_jx(n_ix) * int_{e_iy}(M_jy) + if i = vertical edge [n_ix, e_iy] and j = (B_jx o M_jy) y-oriented BoM spline + (above, 'o' denotes tensor-product for functions) + """ + if uniform_patches: + print(' [[WARNING -- hack in get_K1_and_K1_inv: using copies of 1st-patch matrices in every patch ]] ') + + V1 = V1h.symbolic_space # V1h is FemSpace + domain = V1.domain + K1_blocks = [] + K1_inv_blocks = [] + for k, D in enumerate(domain.interior): + if uniform_patches and k > 0: + K1_k = K1_blocks[0].copy() + K1_inv_k = K1_inv_blocks[0].copy() + + else: + # fem space on patch k: + V1_k = V1h.spaces[k] + K1_k_blocks = [] + for c in [0, 1]: # dim of component + # fem space for comp. dc (TensorFemSpace) + V1_kc = V1_k.spaces[c] + K1_kc_factors = [None, None] + for d in [0, 1]: # dim of variable + # 1d fem space for comp c alond dim d (SplineSpace) + V1_kcd = V1_kc.spaces[d] + if c == d: + K1_kc_factors[d] = histopolation_matrix( + knots=V1_kcd.knots, + degree=V1_kcd.degree, + periodic=V1_kcd.periodic, + normalization=V1_kcd.basis, + xgrid=V1_kcd.ext_greville + ) + else: + K1_kc_factors[d] = collocation_matrix( + knots=V1_kcd.knots, + degree=V1_kcd.degree, + periodic=V1_kcd.periodic, + normalization=V1_kcd.basis, + xgrid=V1_kcd.greville + ) + K1_kc = kron(*K1_kc_factors) + K1_kc.eliminate_zeros() + K1_k_blocks.append(K1_kc) + K1_k = block_diag(K1_k_blocks) + K1_k.eliminate_zeros() + K1_inv_k = inv(K1_k.tocsc()) + K1_inv_k.eliminate_zeros() + + K1_blocks.append(K1_k) + K1_inv_blocks.append(K1_inv_k) + + K1 = block_diag(K1_blocks) + K1_inv = block_diag(K1_inv_blocks) + return K1, K1_inv + +# =============================================================================== + + +def ortho_proj_Hcurl(EE, V1h, domain_h, M1, backend_language='python'): + """ + return orthogonal projection of E on V1h, given M1 the mass matrix + """ + assert isinstance(EE, Tuple) + V1 = V1h.symbolic_space + v = element_of(V1, name='v') + l = LinearForm(v, integral(V1.domain, dot(v, EE))) + lh = discretize( + l, + domain_h, + V1h, + backend=PSYDAC_BACKENDS[backend_language]) + b = lh.assemble() + M1_inv = inverse(M1.mat(), 'pcg', pc='jacobi', tol=1e-10) + sol_coeffs = M1_inv @ b + + return FemField(V1h, coeffs=sol_coeffs) + # =============================================================================== class DiagGrid(): diff --git a/psydac/feec/tests/test_axis_projection.py b/psydac/feec/tests/test_axis_projection.py index 61e1dfb4f..7a81f0953 100644 --- a/psydac/feec/tests/test_axis_projection.py +++ b/psydac/feec/tests/test_axis_projection.py @@ -1,6 +1,6 @@ -from sympde.topology import Square, Derham, element_of -from sympde.expr.expr import BilinearForm, integral -from psydac.feec.multipatch.api import discretize +from sympde.topology import Square, Derham, element_of +from sympde.expr.expr import BilinearForm, integral +from psydac.api.discretization import discretize from psydac.api.settings import PSYDAC_BACKENDS def test_axis_projection(): diff --git a/psydac/feec/tests/test_commuting_projections_dual.py b/psydac/feec/tests/test_commuting_projections_dual.py index 944915bc4..5bc1a68a7 100644 --- a/psydac/feec/tests/test_commuting_projections_dual.py +++ b/psydac/feec/tests/test_commuting_projections_dual.py @@ -47,7 +47,7 @@ def test_transpose_div_3d(Nel, Nq, p, bc, m): u2 = discretize(f2, domain_h, derham_h.V2, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() u3 = discretize(f3, domain_h, derham_h.V3, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() - divT_u3 = - div.matrix.T.dot(u3) + divT_u3 = - div.linop.T.dot(u3) error = abs((u2-divT_u3).toarray()).max() assert error < 2e-10 @@ -105,7 +105,7 @@ def test_transpose_curl_3d(Nel, Nq, p, bc, m): u1 = discretize(f1, domain_h, derham_h.V1, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() u2 = discretize(f2, domain_h, derham_h.V2, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() - curlT_u2 = curl.matrix.T.dot(u2) + curlT_u2 = curl.linop.T.dot(u2) error = abs((u1-curlT_u2).toarray()).max() assert error < 2e-9 @@ -152,7 +152,7 @@ def test_transpose_grad_3d(Nel, Nq, p, bc, m): u0 = discretize(f0, domain_h, derham_h.V0, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() u1 = discretize(f1, domain_h, derham_h.V1, nquads=Nq, backend=PSYDAC_BACKENDS['pyccel-gcc']).assemble() - gradT_u1 = -grad.matrix.T.dot(u1) + gradT_u1 = -grad.linop.T.dot(u1) error = abs((u0-gradT_u1).toarray()).max() assert error < 5e-10 diff --git a/psydac/feec/multipatch/tests/test_feec_conf_projectors_cart_2d.py b/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py similarity index 54% rename from psydac/feec/multipatch/tests/test_feec_conf_projectors_cart_2d.py rename to psydac/feec/tests/test_feec_conf_projectors_cart_2d.py index 3d0dd85ff..fca2845e6 100644 --- a/psydac/feec/multipatch/tests/test_feec_conf_projectors_cart_2d.py +++ b/psydac/feec/tests/test_feec_conf_projectors_cart_2d.py @@ -2,16 +2,15 @@ from collections import OrderedDict import numpy as np -from sympy import Tuple +from sympy import Tuple, lambdify from scipy.sparse.linalg import norm as sp_norm from sympde.topology.domain import Domain from sympde.topology import Derham, Square, IdentityMapping -from psydac.feec.multipatch.api import discretize -from psydac.feec.multipatch.operators import HodgeOperator -from psydac.feec.multipatch.non_matching_operators import construct_h1_conforming_projection, construct_hcurl_conforming_projection -from psydac.feec.multipatch.utils_conga_2d import P_phys_l2, P_phys_hdiv, P_phys_hcurl, P_phys_h1 +from psydac.api.discretization import discretize + +from psydac.fem.projectors import get_dual_dofs def get_polynomial_function(degree, hom_bc_axes, domain): @@ -37,7 +36,10 @@ def get_polynomial_function(degree, hom_bc_axes, domain): # else: g0_y = (y - 0.75)**degree[1] - return g0_x * g0_y + expr = g0_x * g0_y + callable_function = lambdify(domain.coordinates, expr) + + return expr, callable_function # ============================================================================== @@ -46,10 +48,10 @@ def get_polynomial_function(degree, hom_bc_axes, domain): @pytest.mark.parametrize('nc', [5]) @pytest.mark.parametrize('reg', [0]) @pytest.mark.parametrize('hom_bc', [False, True]) -@pytest.mark.parametrize('domain_name', ["4patch_nc", "2patch_nc"]) +@pytest.mark.parametrize('domain_name', ["1patch", "4patch_nc", "2patch_nc"]) @pytest.mark.parametrize("nonconforming, full_mom_pres", [(True, True), (False, True)]) -# NOTE (MCP march 2025): momentum conservation fails for nc = 4 and degree = 3, why? + def test_conf_projectors_2d( V1_type, degree, @@ -60,8 +62,12 @@ def test_conf_projectors_2d( domain_name, nonconforming ): + if domain_name == '1patch': + log_domain = Square('Omega', bounds1=(0, 1), bounds2=(0, 1)) + mapping = IdentityMapping('M1', dim=2) + domain = mapping(log_domain) - if domain_name == '2patch_nc': + elif domain_name == '2patch_nc': A = Square('A', bounds1=(0, 0.5), bounds2=(0, 1)) B = Square('B', bounds1=(0.5, 1.), bounds2=(0, 1)) @@ -95,8 +101,10 @@ def test_conf_projectors_2d( ((0, 1, 1), (2, 1, -1), 1), ((1, 1, 1), (3, 1, -1), 1)], name='domain') + if domain_name == '1patch': + ncells_h = {domain.name: [nc, nc]} - if nonconforming: + elif nonconforming: if len(domain) == 2: ncells_h = { 'M1(A)': [nc, nc], @@ -115,24 +123,14 @@ def test_conf_projectors_2d( for k, D in enumerate(domain.interior): ncells_h[D.name] = [nc, nc] - derham = Derham(domain, ["H1", "Hcurl", "L2"]) domain_h = discretize(domain, ncells=ncells_h) # Vh space - derham_h = discretize(derham, domain_h, degree=degree) - V0h = derham_h.V0 - V1h = derham_h.V1 - V2h = derham_h.V2 - - mappings = OrderedDict([(P.logical_domain, P.mapping) - for P in domain.interior]) - mappings_list = [m.get_callable_mapping() for m in mappings.values()] - p_derham = Derham(domain, ["H1", V1_type, "L2"]) + derham = Derham(domain, ["H1", V1_type, "L2"]) nquads = [(d + 1) for d in degree] - p_derham_h = discretize(p_derham, domain_h, degree=degree) - p_V0h = p_derham_h.V0 - p_V1h = p_derham_h.V1 - p_V2h = p_derham_h.V2 + derham_h = discretize(derham, domain_h, degree=degree) + V0h, V1h, V2h = derham_h.spaces + # full moment preservation only possible if enough interior functions in a # patch (<=> enough cells) @@ -145,38 +143,22 @@ def test_conf_projectors_2d( # moment preservation... # geometric projections (operators) - p_geomP0, p_geomP1, p_geomP2 = p_derham_h.projectors(nquads=nquads) + geomP0, geomP1, geomP2 = derham_h.projectors(nquads=nquads) # conforming projections (scipy matrices) - cP0 = construct_h1_conforming_projection(V0h, reg, mom_pres, hom_bc) - cP1 = construct_hcurl_conforming_projection(V1h, reg, mom_pres, hom_bc) - cP2 = construct_h1_conforming_projection(V2h, reg - 1, mom_pres, hom_bc) + cP0, cP1, cP2 = derham_h.conforming_projectors(kind='sparse', p_moments=mom_pres, hom_bc=hom_bc) - HOp0 = HodgeOperator(p_V0h, domain_h) - M0 = HOp0.to_sparse_matrix() # mass matrix - M0_inv = HOp0.get_dual_Hodge_sparse_matrix() # inverse mass matrix + M0, M1, M2 = derham_h.Hodge_operators(kind='sparse') + M0_inv, M1_inv, M2_inv = derham_h.Hodge_operators(kind='sparse', dual=True) - HOp1 = HodgeOperator(p_V1h, domain_h) - M1 = HOp1.to_sparse_matrix() # mass matrix - M1_inv = HOp1.get_dual_Hodge_sparse_matrix() # inverse mass matrix + bD0, bD1 = derham_h.derivatives(kind='sparse') - HOp2 = HodgeOperator(p_V2h, domain_h) - M2 = HOp2.to_sparse_matrix() # mass matrix - M2_inv = HOp2.get_dual_Hodge_sparse_matrix() # inverse mass matrix - - bD0, bD1 = p_derham_h.broken_derivatives_as_operators - - bD0 = bD0.to_sparse_matrix() # broken grad - bD1 = bD1.to_sparse_matrix() # broken curl or div D0 = bD0 @ cP0 # Conga grad D1 = bD1 @ cP1 # Conga curl or div - assert np.allclose(sp_norm(cP0 - cP0 @ cP0), 0, 1e-12, - 1e-12) # cP0 is a projection - assert np.allclose(sp_norm(cP1 - cP1 @ cP1), 0, 1e-12, - 1e-12) # cP1 is a projection - assert np.allclose(sp_norm(cP2 - cP2 @ cP2), 0, 1e-12, - 1e-12) # cP2 is a projection + assert np.allclose(sp_norm(cP0 - cP0 @ cP0), 0, 1e-12, 1e-12) # cP0 is a projection + assert np.allclose(sp_norm(cP1 - cP1 @ cP1), 0, 1e-12, 1e-12) # cP1 is a projection + assert np.allclose(sp_norm(cP2 - cP2 @ cP2), 0, 1e-12, 1e-12) # cP2 is a projection # D0 maps in the conforming V1 space (where cP1 coincides with Id) assert np.allclose(sp_norm(D0 - cP1 @ D0), 0, 1e-12, 1e-12) @@ -186,14 +168,11 @@ def test_conf_projectors_2d( # comparing projections of polynomials which should be exact # tests on cP0: - g0 = get_polynomial_function( - degree=degree, hom_bc_axes=[ - hom_bc, hom_bc], domain=domain) - g0h = P_phys_h1(g0, p_geomP0, domain, mappings_list) + g0, g0_fun = get_polynomial_function(degree=degree, hom_bc_axes=[hom_bc, hom_bc], domain=domain) + g0h = geomP0(g0_fun) g0_c = g0h.coeffs.toarray() - tilde_g0_c = p_derham_h.get_dual_dofs( - space='V0', f=g0, return_format='numpy_array') + tilde_g0_c = get_dual_dofs(Vh=V0h, f=g0, domain_h = domain_h, return_format='numpy_array') g0_L2_c = M0_inv @ tilde_g0_c # (P0_geom - P0_L2) polynomial = 0 @@ -206,46 +185,27 @@ def test_conf_projectors_2d( # the following projection should be exact for polynomials of proper degree (no bc) # conf_P0* : L2 -> V0 defined by := # for all phi in V0 - g0 = get_polynomial_function(degree=degree, hom_bc_axes=[ - False, False], domain=domain) - g0h = P_phys_h1(g0, p_geomP0, domain, mappings_list) + g0, g0_fun = get_polynomial_function(degree=degree, hom_bc_axes=[False, False], domain=domain) + g0h = geomP0(g0_fun) g0_c = g0h.coeffs.toarray() - tilde_g0_c = p_derham_h.get_dual_dofs( - space='V0', f=g0, return_format='numpy_array') + tilde_g0_c = get_dual_dofs(Vh=V0h, f=g0, domain_h = domain_h, return_format='numpy_array') + g0_star_c = M0_inv @ cP0.transpose() @ tilde_g0_c # (P10_geom - P0_star) polynomial = 0 assert np.allclose(g0_c, g0_star_c, 1e-12, 1e-12) # tests on cP1: + G1_x, G1_x_fun = get_polynomial_function(degree=[degree[0] - 1,degree[1]], hom_bc_axes=[False, hom_bc], domain=domain) + G1_y, G1_y_fun = get_polynomial_function(degree=[degree[0], degree[1] - 1], hom_bc_axes=[hom_bc, False], domain=domain) + + G1 = Tuple(G1_x, G1_y) + G1_fun = [G1_x_fun, G1_y_fun] - G1 = Tuple( - get_polynomial_function( - degree=[ - degree[0] - 1, - degree[1]], - hom_bc_axes=[ - False, - hom_bc], - domain=domain), - get_polynomial_function( - degree=[ - degree[0], - degree[1] - 1], - hom_bc_axes=[ - hom_bc, - False], - domain=domain) - ) - - if V1_type == "Hcurl": - G1h = P_phys_hcurl(G1, p_geomP1, domain, mappings_list) - elif V1_type == "Hdiv": - G1h = P_phys_hdiv(G1, p_geomP1, domain, mappings_list) - + G1h = geomP1(G1_fun) G1_c = G1h.coeffs.toarray() - tilde_G1_c = p_derham_h.get_dual_dofs( - space='V1', f=G1, return_format='numpy_array') + tilde_G1_c = get_dual_dofs(Vh=V1h, f=G1, domain_h=domain_h, return_format='numpy_array') + G1_L2_c = M1_inv @ tilde_G1_c assert np.allclose(G1_c, G1_L2_c, 1e-12, 1e-12) @@ -254,48 +214,29 @@ def test_conf_projectors_2d( if full_mom_pres: # as above - G1 = Tuple( - get_polynomial_function( - degree=[ - degree[0] - 1, - degree[1]], - hom_bc_axes=[ - False, - False], - domain=domain), - get_polynomial_function( - degree=[ - degree[0], - degree[1] - 1], - hom_bc_axes=[ - False, - False], - domain=domain) - ) - - G1h = P_phys_hcurl(G1, p_geomP1, domain, mappings_list) - G1_c = G1h.coeffs.toarray() - - tilde_G1_c = p_derham_h.get_dual_dofs( - space='V1', f=G1, return_format='numpy_array') - G1_star_c = M1_inv @ cP1.transpose() @ tilde_G1_c - # (P1_geom - P1_star) polynomial = 0 - assert np.allclose(G1_c, G1_star_c, 1e-12, 1e-12) + G1_x, G1_x_fun = get_polynomial_function(degree=[degree[0] - 1,degree[1]], hom_bc_axes=[False, False], domain=domain) + G1_y, G1_y_fun = get_polynomial_function(degree=[degree[0], degree[1] - 1], hom_bc_axes=[False, False], domain=domain) + + G1 = Tuple(G1_x, G1_y) + G1_fun = [G1_x_fun, G1_y_fun] + + G1h = geomP1(G1_fun) + G1_c = G1h.coeffs.toarray() + + G1h = geomP1(G1_fun) + G1_c = G1h.coeffs.toarray() + + tilde_G1_c = get_dual_dofs(Vh=V1h, f=G1, domain_h=domain_h, return_format='numpy_array') + G1_star_c = M1_inv @ cP1.transpose() @ tilde_G1_c + # (P1_geom - P1_star) polynomial = 0 + assert np.allclose(G1_c, G1_star_c, 1e-12, 1e-12) # tests on cP2 (non trivial for reg = 1): - g2 = get_polynomial_function( - degree=[ - degree[0] - 1, - degree[1] - 1], - hom_bc_axes=[ - False, - False], - domain=domain) - g2h = P_phys_l2(g2, p_geomP2, domain, mappings_list) + g2, g2_fun = get_polynomial_function(degree=[degree[0] - 1, degree[1] - 1], hom_bc_axes=[False, False], domain=domain) + g2h = geomP2(g2_fun) g2_c = g2h.coeffs.toarray() - tilde_g2_c = p_derham_h.get_dual_dofs( - space='V2', f=g2, return_format='numpy_array') + tilde_g2_c = get_dual_dofs(Vh=V2h, f=g2, domain_h=domain_h, return_format='numpy_array') g2_L2_c = M2_inv @ tilde_g2_c # (P2_geom - P2_L2) polynomial = 0 @@ -305,7 +246,6 @@ def test_conf_projectors_2d( if full_mom_pres: # as above, here with same degree and bc as - # tilde_g2_c = p_derham_h.get_dual_dofs(space='V2', f=g2, return_format='numpy_array', nquads=nquads) g2_star_c = M2_inv @ cP2.transpose() @ tilde_g2_c # (P2_geom - P2_star) polynomial = 0 assert np.allclose(g2_c, g2_star_c, 1e-12, 1e-12) diff --git a/psydac/fem/basic.py b/psydac/fem/basic.py index 09a1466fb..17192c086 100644 --- a/psydac/fem/basic.py +++ b/psydac/fem/basic.py @@ -7,9 +7,9 @@ """ from abc import ABCMeta, abstractmethod -from psydac.linalg.basic import Vector +from psydac.linalg.basic import Vector, LinearOperator -__all__ = ('FemSpace', 'FemField') +__all__ = ('FemSpace', 'FemField', 'FemLinearOperator') #=============================================================================== # ABSTRACT BASE CLASS: FINITE ELEMENT SPACE @@ -380,3 +380,97 @@ def __isub__(self, other): assert self._space is other._space self._coeffs -= other._coeffs return self + +#=============================================================================== +# CONCRETE CLASS: Linear Operator acting on a FEM field +#=============================================================================== +class FemLinearOperator: + """ + Linear operators with an additional Fem layer. + There is also a shorthand access to sparse matrices as they are sometimes used in the FEEC interfaces. + + Parameters + ---------- + fem_domain : psydac.fem.basic.FemSpace + The discrete space + + fem_codomain : psydac.fem.basic.FemSpace + Number of polynomial moments to be preserved in the projection. + + linop : (optional) + Linear Operator. + + sparse_matrix : (optional) + Sparse matrix representation of the linear operator. + """ + + def __init__(self, fem_domain, fem_codomain, linop=None, sparse_matrix=None): + assert isinstance(fem_domain, FemSpace) + assert isinstance(fem_codomain, FemSpace) + + self._fem_domain = fem_domain + self._fem_codomain = fem_codomain + + self._linop_domain = fem_domain.coeff_space + self._linop_codomain = fem_codomain.coeff_space + + self._linop = linop + self._sparse_matrix = sparse_matrix + + @property + def fem_domain(self): + return self._fem_domain + + @property + def fem_codomain(self): + return self._fem_codomain + + @property + def linop_domain(self): + return self._linop_domain + + @property + def linop_codomain(self): + return self._linop_codomain + + @property + def linop(self): + return self._linop + + @property + def toarray(self): + if self._sparse_matrix is not None: + return self._sparse_matrix.toarray() + else: + return self._linop.toarray() + + @property + def tosparse(self): + if self._sparse_matrix is None: + self._sparse_matrix = self._linop.tosparse() + + return self._sparse_matrix + + + def __call__(self, u): + assert isinstance(u, FemField) + assert u.space == self.fem_domain + + if self._linop is not None: + coeffs = self._linop.dot(u.coeffs) + elif self._sparse_matrix is not None: + coeffs = self._sparse_matrix.dot(u.coeffs) + else: + raise NotImplementedError('Class does not provide a __call__ method without a matrix or linear operator') + + return FemField(self.fem_codomain, coeffs=coeffs) + + def dot( self, f_coeffs, out=None ): + assert isinstance(f_coeffs, Vector) + assert f_coeffs.space is self._linop_domain + + if self._linop is not None or self._sparse_matrix is not None: + f = FemField(self.fem_domain, coeffs=f_coeffs) + return self(f).coeffs + else: + raise NotImplementedError('Class does not provide a dot method without a matrix or linear operator') diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index d40db6a8c..e13185051 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,8 +1,17 @@ import numpy as np -from psydac.linalg.kron import KroneckerDenseMatrix -from psydac.core.bsplines import hrefinement_matrix -from psydac.linalg.stencil import StencilVectorSpace +from sympde.topology import element_of +from sympde.topology.space import ScalarFunction +from sympde.topology.mapping import Mapping +from sympde.calculus import dot +from sympde.expr.expr import LinearForm, integral + +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 __all__ = ('knots_to_insert', 'knot_insertion_projection_operator') @@ -100,3 +109,52 @@ def knot_insertion_projection_operator(domain, codomain): ops.append(np.eye(d.nbasis)) return KroneckerDenseMatrix(domain.coeff_space, codomain.coeff_space, *ops) + + +def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='stencil_array'): + """ + return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(Vh)) of a given function f, as a stencil array or numpy array + + Parameters + ---------- + Vh : FemSpace + The discrete space for the dual dofs + + f : + The function used for evaluation + + domain_h : + The discrete domain corresponding to Vh + + backend_language: + The backend used to accelerate the code + + return_format: + The format of the dofs, can be 'stencil_array' or 'numpy_array' + + Returns + ------- + tilde_f: + The dual dofs + """ + + from psydac.api.discretization import discretize + + assert isinstance(Vh, FemSpace) + + V = Vh.symbolic_space + v = element_of(V, name='v') + + if Vh.is_vector_valued: + expr = dot(f,v) + else: + expr = f*v + + l = LinearForm(v, integral( V.domain, expr)) + lh = discretize(l, domain_h, Vh, backend=PSYDAC_BACKENDS[backend_language]) + tilde_f = lh.assemble() + + if return_format == 'numpy_array': + return tilde_f.toarray() + else: + return tilde_f diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index f9f63406d..062723e35 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -700,7 +700,7 @@ def toarray(self): def tosparse(self): from scipy.sparse import csr_matrix - return self._scalar*csr_matrix(self._operator.toarray()) + return self._scalar*csr_matrix(self._operator.tosparse()) def transpose(self, conjugate=False): return ScaledLinearOperator(domain=self.codomain, codomain=self.domain, c=self._scalar if not conjugate else np.conjugate(self._scalar), A=self._operator.transpose(conjugate=conjugate)) diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py index 57a9a6b86..d0d406133 100644 --- a/psydac/linalg/utilities.py +++ b/psydac/linalg/utilities.py @@ -3,7 +3,7 @@ import numpy as np from math import sqrt -from psydac.linalg.basic import Vector +from psydac.linalg.basic import Vector, MatrixFreeLinearOperator from psydac.linalg.stencil import StencilVectorSpace, StencilVector from psydac.linalg.block import BlockVector, BlockVectorSpace from psydac.linalg.topetsc import petsc_local_to_psydac, get_npts_per_block @@ -11,7 +11,8 @@ __all__ = ( 'array_to_psydac', 'petsc_to_psydac', - '_sym_ortho' + '_sym_ortho', + 'SparseMatrixLinearOperator' ) #============================================================================== @@ -198,3 +199,28 @@ def _sym_ortho(a, b): s = c * tau r = a / c return c, s, r + +#============================================================================== +class SparseMatrixLinearOperator(MatrixFreeLinearOperator): + """ + Wrap a sparse matrix into a MatrixFreeLinearOperator + """ + + def __init__(self, domain, codomain, sparse_matrix): + self._matrix = sparse_matrix + + def dot_sparse(v): + res = self._matrix @ v.toarray() + Mv = array_to_psydac(res, self.codomain) + return Mv + + super().__init__(domain, codomain, dot_sparse) + + def tosparse(self): + return self._matrix + + def toarray(self): + return self._matrix.toarray() + + def transpose(self, conjugate=False): + return SparseMatrixLinearOperator(self.codomain, self.domain, self._matrix.T) \ No newline at end of file