From fe9b0d074ccfcdd23f190be6ec8f812b2122725f Mon Sep 17 00:00:00 2001 From: jowezarek Date: Thu, 30 Oct 2025 15:04:28 +0100 Subject: [PATCH 01/25] add dirichlet projectors and tests --- psydac/api/feec.py | 47 ++- psydac/linalg/tests/test_solvers.py | 542 +++++++++++++++++++++++++++- 2 files changed, 580 insertions(+), 9 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 2bd95ffbd..510a4f510 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -1,3 +1,10 @@ +import numpy as np + +from scipy.sparse import dia_matrix + +from sympde.expr import integral, BilinearForm +from sympde.topology import elements_of, Line, Derham + from psydac.api.basic import BasicDiscrete from psydac.feec.derivatives import Derivative1D, Gradient2D, Gradient3D @@ -26,7 +33,12 @@ from psydac.fem.basic import FemSpace, FemLinearOperator from psydac.fem.vector import VectorFemSpace -from psydac.linalg.basic import IdentityOperator + +from psydac.linalg.basic import LinearOperator, IdentityOperator +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.direct_solvers import BandedSolver +from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix +from psydac.linalg.stencil import StencilVectorSpace __all__ = ('DiscreteDeRham', 'DiscreteDeRhamMultipatch',) @@ -288,6 +300,39 @@ def derivatives(self, kind='femlinop'): elif kind == 'linop': return tuple(b_diff.linop for b_diff in self._derivatives) + #-------------------------------------------------------------------------- + def dirichlet_projectors(self, kind='femlinop'): + """ + Returns operators that apply the correct Dirichlet boundary conditions. + + Parameters + ---------- + kind : str + The kind of the projector, can be 'femlinop' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'linop' returns a psydac LinearOperator + + Returns + ------- + d_projectors : list + List of or + The Dirichlet boundary projectors of each space and in desired form. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form. + + """ + assert kind in ('femlinop', 'linop') + + from psydac.linalg.tests.test_solvers import DirichletBoundaryProjector + d_projectors = [DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + + if kind == 'femlinop': + d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)] + + return d_projectors + #-------------------------------------------------------------------------- def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): """ diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 6d4b37e73..b649b5cb9 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -1,11 +1,28 @@ - -import numpy as np -import pytest -from psydac.linalg.solvers import inverse -from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector -from psydac.linalg.basic import LinearSolver -from psydac.ddm.cart import DomainDecomposition, CartDecomposition - +import time +import numpy as np +import pytest + +from sympy import sin, pi, sqrt, Tuple + +from scipy.sparse import dia_matrix + +from sympde.calculus import inner, cross +from sympde.expr import integral, LinearForm, BilinearForm, EssentialBC +from sympde.topology import element_of, elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace +from sympde.topology.datatype import SpaceType, H1Space, HcurlSpace + +from psydac.api.discretization import discretize +from psydac.api.essential_bc import apply_essential_bc +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from psydac.fem.basic import FemSpace +from psydac.linalg.basic import LinearOperator, Vector, IdentityOperator +from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator +from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix +from psydac.linalg.solvers import inverse +from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector +from psydac.linalg.tests.test_kron_direct_solver import matrix_to_bandsolver +from psydac.linalg.direct_solvers import BandedSolver def define_data_hermitian(n, p, dtype=float): domain_decomposition = DomainDecomposition([n - p], [False]) @@ -54,6 +71,193 @@ def define_data(n, p, matrix_data, dtype=float): xe[s:e + 1] = np.random.random(e + 1 - s) return(V, A, xe) +class SquareTorus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + +class Annulus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)'} + + _ldim = 2 + _pdim = 2 + +class SinMapping1D(Mapping): + + _expressions = {'x': 'sin((pi/2)*x1)'} + + _ldim = 1 + _pdim = 1 + +def _test_LO_equality_using_rng(A, B): + """ + A simple tool to check with almost certainty that two linear operators are identical, + by applying them repeatedly to random vectors. + + """ + + assert isinstance(A, LinearOperator) + assert isinstance(B, LinearOperator) + assert A.domain is B.domain + assert A.codomain is B.codomain + + rng = np.random.default_rng(42) + + x = A.domain.zeros() + y1 = A.codomain.zeros() + y2 = y1.copy() + + n = 10 + + for _ in range(n): + + x *= 0. + + if isinstance(A.domain, BlockVectorSpace): + for block in x.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=x._data.shape, dtype="float64", out=x._data) + + A.dot(x, out=y1) + B.dot(x, out=y2) + + diff = y1 - y2 + err = A.codomain.inner(diff, diff) + + assert err < 1e-15 + +class DirichletBoundaryProjector(LinearOperator): + + def __init__(self, fem_space, bcs=None, space_kind=None): + + assert isinstance(fem_space, FemSpace) + + coeff_space = fem_space.coeff_space + self._domain = coeff_space + self._codomain = coeff_space + + if bcs is not None: + self._bcs = bcs + else: + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._domain + + @property + def dtype(self): + return None + + @property + def bcs(self): + return self._bcs + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def _get_bcs(self, fem_space, space_kind=None): + """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + space = fem_space.symbolic_space + periodic = fem_space.periodic + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + + kind = space_kind_str + dim = space.domain.dim + + if kind == 'l2': + return None + + u = element_of(space, name="u") + ebcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + if kind == "h1": + bcs = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim >= 2: + bcs += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + + elif kind == 'hcurl': + assert dim in (2, 3) + bcs_x = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_x += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs_y = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim == 3: + bcs_y += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + if dim == 3: + bcs_z = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_z += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + elif kind == 'hdiv': + assert dim in (2, 3) + bcs_x = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_y = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_z = [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return bcs + + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + if isinstance(self.domain, StencilVectorSpace): + apply_essential_bc(out, *self._bcs) + else: + for block, block_bcs in zip(out, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out #=============================================================================== @pytest.mark.parametrize( 'n', [5, 10, 13] ) @@ -204,6 +408,328 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): assert errh_norm < tol assert solver == 'pcg' or errc_norm < tol +#=============================================================================== +def test_function_space_dirichlet_projector(): + + ncells_3d = [8, 8, 8] + degree_3d = [2, 2, 2] + periodic_3d = [False, True, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + dims = [1, 2, 3] + rng = np.random.default_rng(42) + + print() + for dim in dims: + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] + degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] + periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + # The function defined here satisfy the corresponding homogeneous Dirichlet BCs + if dim == 1: + x = domain.coordinates + V = ScalarFunctionSpace('V', domain, kind='H1') + f = sin(2*pi*x) + if dim == 2: + x, y = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind=H1Space) + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + else: + V = VectorFunctionSpace('V', domain, kind='hCuRl') + f1 = x + f2 = y + f = Tuple(f1, f2) + if dim == 3: + x, y, z = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind='h1') + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) + elif i == 1: + V = VectorFunctionSpace('V', domain, kind=HcurlSpace) + f1 = z * (z - 1) * x + f2 = z * (z - 1) * y + f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f = Tuple(f1, f2, f3) + else: + V = VectorFunctionSpace('V', domain, kind='Hdiv') + f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f3 = z * (z-1) * sin(x*y) + f = Tuple(f1, f2, f3) + + u, v = elements_of(V, names='u, v') + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + Vh = discretize(V, domain_h, degree=degree) + expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (Vh, Vh), backend=backend) + abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) + + I = IdentityOperator(Vh.coeff_space) + DBP = DirichletBoundaryProjector(Vh) + + M = ah.assemble() + M_0 = DBP @ M @ DBP + (I - DBP) + Mb = abh.assemble() + + # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc + lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f + l = LinearForm(v, integral(domain, lexpr)) + lh = discretize(l, domain_h, Vh, backend=backend) + rhs = lh.assemble() + A = M + 1e30*Mb + A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) + fc = A_inv @ rhs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DBP @ fc + diff = fc - fc2 + err = np.linalg.norm(diff.toarray()) + print(f' | f - P @ f | = {err}') + assert err < 1e-15 + + # 2. + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = Vh.coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DBP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < 1e-15 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DBP @ fc2 + diff = fc2 - fc3 + err = np.linalg.norm(diff.toarray()) + print(f' | P @ f - P @ P @ f | = {err}') + assert err == 0. + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_squared = M.dot_inner (fc, fc) + l2_norm_squared2 = M_0.dot_inner(fc, fc) + diff = l2_norm_squared - l2_norm_squared2 + print(f' || f ||^2 = {l2_norm_squared} should be equal to') + print(f' || P @ f ||^2 = {l2_norm_squared2}') + assert diff < 1e-15 + + print() + + print() + +#=============================================================================== +def test_discrete_derham_dirichlet_projector(): + + ncells = [8, 8, 8] + degree = [2, 2, 2] + periodic = [False, True, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + dims = [1, 2, 3] + rng = np.random.default_rng(42) + + # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs + + f11 = lambda x : np.sin(2*np.pi*x) + + r2 = lambda x, y : np.sqrt(x**2 + y**2) + f21 = lambda x, y : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f22_1 = lambda x, y : x + f22_2 = lambda x, y : y + f22 = (f22_1, f22_2) + + f31 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) * z * (z - 1) + f32_1 = lambda x, y, z : z * (z - 1) * x + f32_2 = lambda x, y, z : z * (z - 1) * y + f32_3 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f32 = (f32_1, f32_2, f32_3) + f33_1 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_2 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_3 = lambda x, y, z : z * (z - 1) * np.sin(x*y) + f33 = (f33_1, f33_2, f33_3) + + funs = [[f11], [f21, f22], [f31, f32, f33]] + + print() + for dim in dims: + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] + degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] + periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] + + domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree_dim) + + db_projectors = derham_h.dirichlet_projectors(kind='linop') + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DBP = db_projectors[i] + + M = ah.assemble() + M_0 = DBP @ M @ DBP + (I - DBP) + Mb = abh.assemble() + + f = funs[dim-1][i] + fc = derham_h.projectors()[i](f).coeffs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DBP @ fc + diff = fc - fc2 + err = np.linalg.norm(diff.toarray()) + print(f' | f - P @ f | = {err}') + assert err < 1e-15 + + # 2. + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DBP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < 1e-15 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DBP @ fc2 + diff = fc2 - fc3 + err = np.linalg.norm(diff.toarray()) + print(f' | P @ f - P @ P @ f | = {err}') + assert err == 0. + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_squared = M.dot_inner (fc, fc) + l2_norm_squared2 = M_0.dot_inner(fc, fc) + diff = l2_norm_squared - l2_norm_squared2 + print(f' || f ||^2 = {l2_norm_squared} should be equal to') + print(f' || P @ f ||^2 = {l2_norm_squared2}') + assert diff < 1e-15 + + print() + + print() + # =============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== From 9ec94bd17f0b9b3285610beff783090c177e4e60 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Thu, 30 Oct 2025 17:08:17 +0100 Subject: [PATCH 02/25] add diagonal method to KroneckerStencilMatrix, add tests --- psydac/linalg/kron.py | 53 ++++++++- psydac/linalg/stencil.py | 30 ++++-- .../linalg/tests/test_kron_stencil_matrix.py | 102 ++++++++++++++++++ 3 files changed, 176 insertions(+), 9 deletions(-) diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index e4abfdd33..77ddf91ba 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -6,7 +6,7 @@ from scipy.sparse import coo_matrix from psydac.linalg.basic import LinearOperator, LinearSolver -from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix +from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix, StencilDiagonalMatrix __all__ = ('KroneckerStencilMatrix', 'KroneckerLinearSolver', @@ -218,6 +218,57 @@ def toarray(self): def transpose(self, conjugate=False): mats_tr = [Mi.transpose(conjugate=conjugate) for Mi in self.mats] return KroneckerStencilMatrix(self.codomain, self.domain, *mats_tr) + + def diagonal(self, *, inverse=False, sqrt=False, out=None): + """ + Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. + + Parameters + ---------- + inverse : bool + If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. + + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root. + + out : StencilDiagonalMatrix + If provided, write the diagonal entries into this matrix. (Default: None). + + Returns + ------- + StencilDiagonalMatrix + The matrix which contains the main diagonal of self (or its inverse (square root)). + + """ + # Check `inverse` and `sqrt` argument + assert isinstance(inverse, bool) + assert isinstance(sqrt, bool) + + # Determine domain and codomain of the StencilDiagonalMatrix + V, W = self.domain, self.codomain + if inverse: + V, W = W, V + + # Check `out` argument + if out is not None: + assert isinstance(out, StencilDiagonalMatrix) + assert out.domain is V + assert out.codomain is W + + # Obtain nested numpy array of diagonal entries (or their inverse (square root)) + # by using the `diagonal` method of StencilMatrices + diag = 1. + for mat in self.mats[::-1]: + diag = np.array([d*diag for d in mat.diagonal(inverse=inverse, sqrt=sqrt)._data], dtype='float64') + + if out is not None: + np.copyto(diag, out._data) + else: + out = StencilDiagonalMatrix(V, W, diag) + + return out #============================================================================== class KroneckerDenseMatrix(LinearOperator): diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 73599655a..3f5ec8d87 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1419,11 +1419,12 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): Returns ------- StencilDiagonalMatrix - The matrix which contains the main diagonal of self (or its inverse). + The matrix which contains the main diagonal of self (or its inverse (square root)). """ - # Check `inverse` argument + # Check `inverse` and `sqrt` argument assert isinstance(inverse, bool) + assert isinstance(sqrt, bool) # Determine domain and codomain of the StencilDiagonalMatrix V, W = self.domain, self.codomain @@ -1436,7 +1437,6 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): assert out.domain is V assert out.codomain is W - # Extract diagonal data from self and identify output array diagonal_indices = self._get_diagonal_indices() diag = self._data[diagonal_indices] @@ -2049,16 +2049,21 @@ def copy(self, *, out=None): return out - def diagonal(self, *, inverse = False, out = None): + def diagonal(self, *, inverse = False, sqrt = False, out = None): """ Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. - In the default case (inverse=False, out=None) self is returned. + In the default case (inverse=False, sqrt=False, out=None) self is returned. Parameters ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. + + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root. out : StencilDiagonalMatrix If provided, write the diagonal entries into this matrix. (Default: None). @@ -2066,11 +2071,12 @@ def diagonal(self, *, inverse = False, out = None): Returns ------- StencilDiagonalMatrix - Either self, or another StencilDiagonalMatrix with the diagonal inverse. + Either self, or another StencilDiagonalMatrix with the diagonal (or its inverse (square root)). """ - # Check `inverse` argument + # Check `inverse` and `sqrt` argument assert isinstance(inverse, bool) + assert isinstance(sqrt, bool) # Determine domain and codomain of the `out` matrix V, W = self.domain, self.codomain @@ -2086,12 +2092,20 @@ def diagonal(self, *, inverse = False, out = None): assert out.codomain is W data = out._data + diag = self._data + # Calculate entries, or set `out=self` in default case if inverse: data = np.divide(1, diag, out=data) elif out: np.copyto(data, diag) - else: + + if sqrt: + if (not inverse) and (out is None): + data = diag.copy() + np.sqrt(data, out=data) + + if (not inverse) and (not sqrt) and (out is None): out = self # If needed create a new StencilDiagonalMatrix object diff --git a/psydac/linalg/tests/test_kron_stencil_matrix.py b/psydac/linalg/tests/test_kron_stencil_matrix.py index 75f712e4e..d4b76d6da 100644 --- a/psydac/linalg/tests/test_kron_stencil_matrix.py +++ b/psydac/linalg/tests/test_kron_stencil_matrix.py @@ -9,6 +9,12 @@ from psydac.linalg.stencil import StencilVector from psydac.linalg.stencil import StencilMatrix from psydac.linalg.kron import KroneckerStencilMatrix + +from sympde.topology import Square, Line, Derham, elements_of +from sympde.expr import BilinearForm, integral +from sympde.calculus import inner +from psydac.linalg.block import BlockLinearOperator +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL #=============================================================================== def compute_global_starts_ends(domain_decomposition, npts): ndims = len(npts) @@ -112,3 +118,99 @@ def test_KroneckerStencilMatrix(dtype, npts, pads, periodic): # Test dot product assert np.array_equal(M_sp.dot(w.toarray()), M.dot(w).toarray()) + +#============================================================================== +def test_KroneckerStencilMatrix_diagonal(): + """We create three mass matrices (Stencil/Block and Kronecker) belonging to a 2D de Rham sequence, and compare their diagonals.""" + + from psydac.api.discretization import discretize + + ncells = [6, 7] + degree = [3, 2] + mult = [1, 2] + periodic = [False, True] + + backend = PSYDAC_BACKEND_GPYCCEL + + # 1. Obtain StencilMatrix / BlockLinearOperator (of StencilMatrices) mass matrices + + domain = Square('S', bounds1=(0,1), bounds2=(0,2)) + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity=mult) + + V0, V1, V2 = derham.spaces + V0h, V1h, V2h = derham_h.spaces + V0cs, V1cs, V2cs = [Vh.coeff_space for Vh in derham_h.spaces] + + u0, v0 = elements_of(V0, names='u0, v0') + u1, v1 = elements_of(V1, names='u1, v1') + u2, v2 = elements_of(V2, names='u2, v2') + + a0 = BilinearForm((u0, v0), integral(domain, u0*v0)) + a1 = BilinearForm((u1, v1), integral(domain, inner(u1, v1))) + a2 = BilinearForm((u2, v2), integral(domain, u2*v2)) + + a0h = discretize(a0, domain_h, (V0h, V0h), backend=backend) + a1h = discretize(a1, domain_h, (V1h, V1h), backend=backend) + a2h = discretize(a2, domain_h, (V2h, V2h), backend=backend) + + M0 = a0h.assemble() + M1 = a1h.assemble() + M2 = a2h.assemble() + + # 2. Obtain KroneckerStencilMatrix / BlockLinearOperator (of KroneckerStencilMatrices) mass matrices + + domain_1d_x = Line('L', bounds=(0,1)) + domain_1d_y = Line('L', bounds=(0,2)) + domains_1d = [domain_1d_x, domain_1d_y] + + M0s_1d = [] + M1s_1d = [] + + for n, d, m, p, domain_1d in zip(ncells, degree, mult, periodic, domains_1d): + derham_1d = Derham(domain_1d) + + domain_1d_h = discretize(domain_1d, ncells=[n, ], periodic=[p, ]) + derham_1d_h = discretize(derham_1d, domain_1d_h, degree=[d, ], multiplicity=[m, ]) + + V0_1d, V1_1d = derham_1d.spaces + V0h_1d, V1h_1d = derham_1d_h.spaces + + u0_1d, v0_1d = elements_of(V0_1d, names='u0, v0') + u1_1d, v1_1d = elements_of(V1_1d, names='u1, v1') + + a0_1d = BilinearForm((u0_1d, v0_1d), integral(domain_1d, u0_1d*v0_1d)) + a1_1d = BilinearForm((u1_1d, v1_1d), integral(domain_1d, u1_1d*v1_1d)) + + a0h_1d = discretize(a0_1d, domain_1d_h, (V0h_1d, V0h_1d)) + a1h_1d = discretize(a1_1d, domain_1d_h, (V1h_1d, V1h_1d)) + + M0s_1d.append(a0h_1d.assemble()) + M1s_1d.append(a1h_1d.assemble()) + + M0_kron = KroneckerStencilMatrix(V0cs, V0cs, *M0s_1d) + M1_kron = BlockLinearOperator(V1cs, V1cs, [[KroneckerStencilMatrix(V1cs[0], V1cs[0], M1s_1d[0], M0s_1d[1]), None], + [None, KroneckerStencilMatrix(V1cs[1], V1cs[1], M0s_1d[0], M1s_1d[1])]]) + M2_kron = KroneckerStencilMatrix(V2cs, V2cs, *M1s_1d) + + # 3. Test! + + for M, M_kron in zip([M0, M1, M2], [M0_kron, M1_kron, M2_kron]): + M_diag_arr = [] + M_kron_diag_arr = [] + + options = [True, False] + + for inverse in options: + for sqrt in options: + M_diag = M.diagonal(inverse=inverse, sqrt=sqrt) + M_kron_diag = M_kron.diagonal(inverse=inverse, sqrt=sqrt) + M_diag_arr.append(M_diag) + M_kron_diag_arr.append(M_kron_diag) + + for M_diag, M_kron_diag in zip(M_diag_arr, M_kron_diag_arr): + diff = M_diag.toarray() - M_kron_diag.toarray() + err = np.linalg.norm(diff) + assert err < 1e-10 # arbitrary bound (largest occuring error = 6.7e-13) From a898caed9ca34db79684e0c50b5de8fab432a8f9 Mon Sep 17 00:00:00 2001 From: Frederik Schnack Date: Fri, 31 Oct 2025 14:12:11 +0100 Subject: [PATCH 03/25] fix periodic BC in conforming projectors and add test --- psydac/feec/conforming_projectors.py | 14 ++++++++++++-- psydac/linalg/tests/test_solvers.py | 9 ++++++++- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/psydac/feec/conforming_projectors.py b/psydac/feec/conforming_projectors.py index d045a702e..638c33e15 100644 --- a/psydac/feec/conforming_projectors.py +++ b/psydac/feec/conforming_projectors.py @@ -1211,8 +1211,8 @@ def construct_h1_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=- def get_vertex_index(coords): - nbasis0 = Vh.spaces[coords[0]].nbasis - 1 - nbasis1 = Vh.spaces[coords[1]].nbasis - 1 + nbasis0 = Vh.spaces[0].nbasis - 1 + nbasis1 = Vh.spaces[1].nbasis - 1 # patch local index multi_index = [None] * ndim @@ -1233,6 +1233,9 @@ def vertex_moment_indices(axis, coords, p_moments): for co in [(0,0), (1,0), (0,1), (1,1)]: + if all(Vh.periodic): + break + # global index ig = get_vertex_index(co) @@ -1304,6 +1307,10 @@ def get_mu_minus(j, coarse_space, fine_space, R): # boundary condition for bn in domain.boundary: + + if Vh.periodic[bn.axis]: + continue + space_k = Vh axis = bn.axis @@ -1402,6 +1409,9 @@ def edge_moment_index(p, i, axis, ext): # boundary condition for bn in domain.boundary: + if Vh.periodic[bn.axis]: + continue + axis = bn.axis d = 1 - axis ext = bn.ext diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index b649b5cb9..6c7b7ac41 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -649,6 +649,9 @@ def test_discrete_derham_dirichlet_projector(): db_projectors = derham_h.dirichlet_projectors(kind='linop') + if dim == 2: + conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + nn = NormalVector('nn') for i in range(dim): @@ -676,6 +679,10 @@ def test_discrete_derham_dirichlet_projector(): I = IdentityOperator(derham_h.spaces[i].coeff_space) DBP = db_projectors[i] + if dim == 2: + CP = conf_projectors[i] + _test_LO_equality_using_rng(DBP, CP) + M = ah.assemble() M_0 = DBP @ M @ DBP + (I - DBP) Mb = abh.assemble() @@ -736,4 +743,4 @@ def test_discrete_derham_dirichlet_projector(): if __name__ == "__main__": import sys - pytest.main( sys.argv ) + pytest.main( sys.argv ) \ No newline at end of file From b8d618cb9752665bb085dba4cea52f2ef5f1f7a3 Mon Sep 17 00:00:00 2001 From: Frederik Schnack Date: Fri, 31 Oct 2025 15:36:45 +0100 Subject: [PATCH 04/25] add preliminary DirichletMultipatchBoundaryProjector --- psydac/api/feec.py | 33 ++++ psydac/linalg/tests/test_solvers.py | 248 +++++++++++++++++++++++++++- 2 files changed, 280 insertions(+), 1 deletion(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 510a4f510..17c76e37b 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -651,3 +651,36 @@ def projectors(self, *, kind='global', nquads=None): elif self.dim == 3: raise NotImplementedError("3D projectors are not available") + + #-------------------------------------------------------------------------- + def dirichlet_projectors(self, kind='femlinop'): + """ + Returns operators that apply the correct Dirichlet boundary conditions. + + Parameters + ---------- + kind : str + The kind of the projector, can be 'femlinop' or 'linop'. + - 'femlinop' returns a psydac FemLinearOperator (default) + - 'linop' returns a psydac LinearOperator + + Returns + ------- + d_projectors : list + List of or + The Dirichlet boundary projectors of each space and in desired form. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form. + + """ + assert kind in ('femlinop', 'linop') + + from psydac.linalg.tests.test_solvers import DirichletMultipatchBoundaryProjector + d_projectors = [DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + + if kind == 'femlinop': + d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)] + + return d_projectors diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 6c7b7ac41..6eb3d62cf 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -737,10 +737,256 @@ def test_discrete_derham_dirichlet_projector(): print() +#=============================================================================== +class DirichletMultipatchBoundaryProjector(LinearOperator): + + def __init__(self, fem_space, bcs=None, space_kind=None): + + assert isinstance(fem_space, FemSpace) + assert fem_space.is_multipatch + + coeff_space = fem_space.coeff_space + self._domain = coeff_space + self._codomain = coeff_space + + if bcs is not None: + self._bcs = bcs + else: + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._domain + + @property + def dtype(self): + return None + + @property + def bcs(self): + return self._bcs + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def _get_bcs(self, fem_space, space_kind=None): + """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + space = fem_space.symbolic_space + periodic = fem_space.periodic + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + kind = space_kind_str + dim = space.domain.dim + assert dim==2 + + if kind == 'l2': + return None + + u = element_of(space, name="u") + + if kind == "h1": + bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + + elif kind == 'hcurl': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 0: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 1: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + + elif kind == 'hdiv': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 1: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 0: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return bcs + + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + + # apply bc on each patch + for p in out.blocks: + + if isinstance(p, StencilVector): + apply_essential_bc(p, *self._bcs) + else: + for block, block_bcs in zip(p, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + +#=============================================================================== +def test_discrete_derham_dirichlet_projector_multipatch(): + + ncells = [8, 8] + degree = [2, 2] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain + domain = build_multipatch_domain(domain_name='annulus_3') + + rng = np.random.default_rng(42) + + # The following are functions satisfying homogeneous Dirichlet BCs + r = lambda x, y : np.sqrt(x**2 + y**2) + f1 = lambda x, y : (r(x, y) - 0.5) * (r(x, y) - 1) + f2_1 = lambda x, y : x + f2_2 = lambda x, y : y + f2 = (f2_1, f2_2) + funs = [f1, f2] + print() + + boundary = domain.boundary + + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_h = {} + for k, D in enumerate(domain.interior): + ncells_h[D.name] = ncells + + domain_h = discretize(domain, ncells=ncells_h, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree) + + projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) + + db_projectors = derham_h.dirichlet_projectors(kind='linop') + + conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + + nn = NormalVector('nn') + + for i in range(2): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + expr = u*v + if (i == 1): + boundary_expr = cross(nn, u) * cross(nn, v) + expr = inner(u,v) + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DBP = db_projectors[i] + + M = ah.assemble() + M_0 = DBP @ M @ DBP + (I - DBP) + Mb = abh.assemble() + + f = funs[i] + fc = projectors[i](f).coeffs + + # 1. + # The coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DBP @ fc + diff = fc - fc2 + err = np.linalg.norm(diff.toarray()) + print(f' | f - P @ f | = {err}') + assert err < 1e-15 + + # 2. + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + for patch in rdm_coeffs.blocks: + + if isinstance(patch.space, BlockVectorSpace): + for block in patch.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=patch._data.shape, dtype="float64", out=patch._data) + + rdm_coeffs2 = DBP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < 1e-15 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DBP @ fc2 + diff = fc2 - fc3 + err = np.linalg.norm(diff.toarray()) + print(f' | P @ f - P @ P @ f | = {err}') + assert err == 0. + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm_squared = M.dot_inner (fc, fc) + l2_norm_squared2 = M_0.dot_inner(fc, fc) + diff = l2_norm_squared - l2_norm_squared2 + print(f' || f ||^2 = {l2_norm_squared} should be equal to') + print(f' || P @ f ||^2 = {l2_norm_squared2}') + assert diff < 1e-15 + + print() + # =============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) \ No newline at end of file + pytest.main( sys.argv ) From f35d98e569397deb4c1be62e05799ab3e553f24d Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 31 Oct 2025 16:13:59 +0100 Subject: [PATCH 05/25] move Dirichlet projectors to fem.projectors --- psydac/api/feec.py | 4 +- psydac/fem/projectors.py | 272 ++++++++++++++++++++++++++- psydac/linalg/tests/test_solvers.py | 273 +--------------------------- 3 files changed, 271 insertions(+), 278 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 17c76e37b..8b3e56283 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -325,7 +325,7 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - from psydac.linalg.tests.test_solvers import DirichletBoundaryProjector + from psydac.fem.projectors import DirichletBoundaryProjector d_projectors = [DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]] if kind == 'femlinop': @@ -677,7 +677,7 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - from psydac.linalg.tests.test_solvers import DirichletMultipatchBoundaryProjector + from psydac.fem.projectors import DirichletMultipatchBoundaryProjector d_projectors = [DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] if kind == 'femlinop': diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index e13185051..42758981b 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,19 +1,20 @@ import numpy as np 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 sympde.expr import LinearForm, integral, EssentialBC +from sympde.topology.datatype import SpaceType +from psydac.api.essential_bc import apply_essential_bc from psydac.api.settings import PSYDAC_BACKENDS - -from psydac.linalg.kron import KroneckerDenseMatrix from psydac.core.bsplines import hrefinement_matrix -from psydac.linalg.stencil import StencilVectorSpace from psydac.fem.basic import FemSpace +from psydac.linalg.basic import LinearOperator, Vector +from psydac.linalg.kron import KroneckerDenseMatrix +from psydac.linalg.stencil import StencilVectorSpace, StencilVector -__all__ = ('knots_to_insert', 'knot_insertion_projection_operator') +__all__ = ('knots_to_insert', 'knot_insertion_projection_operator', 'get_dual_dofs', + 'DirichletBoundaryProjector', 'DirichletMultipatchBoundaryProjector') def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): """ Compute the point difference between the fine grid and coarse grid.""" @@ -27,7 +28,6 @@ def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): assert abs(intersection-coarse_grid).max()_{L2} i = 1, .. dim(Vh)) of a given function f, as a stencil array or numpy array @@ -158,3 +157,258 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste return tilde_f.toarray() else: return tilde_f + +#=============================================================================== +class DirichletBoundaryProjector(LinearOperator): + + def __init__(self, fem_space, bcs=None, space_kind=None): + + assert isinstance(fem_space, FemSpace) + + coeff_space = fem_space.coeff_space + self._domain = coeff_space + self._codomain = coeff_space + + if bcs is not None: + self._bcs = bcs + else: + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._domain + + @property + def dtype(self): + return None + + @property + def bcs(self): + return self._bcs + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def _get_bcs(self, fem_space, space_kind=None): + """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + space = fem_space.symbolic_space + periodic = fem_space.periodic + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + + kind = space_kind_str + dim = space.domain.dim + + if kind == 'l2': + return None + + u = element_of(space, name="u") + ebcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + if kind == "h1": + bcs = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim >= 2: + bcs += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + + elif kind == 'hcurl': + assert dim in (2, 3) + bcs_x = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_x += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs_y = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + if dim == 3: + bcs_y += [ebcs[4], ebcs[5]] if periodic[2] == False else [] + if dim == 3: + bcs_z = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_z += [ebcs[2], ebcs[3]] if periodic[1] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + elif kind == 'hdiv': + assert dim in (2, 3) + bcs_x = [ebcs[0], ebcs[1]] if periodic[0] == False else [] + bcs_y = [ebcs[2], ebcs[3]] if periodic[1] == False else [] + if dim == 3: + bcs_z = [ebcs[4], ebcs[5]] if periodic[2] == False else [] + bcs = [bcs_x, bcs_y] + if dim == 3: + bcs.append(bcs_z) + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return bcs + + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + if isinstance(self.domain, StencilVectorSpace): + apply_essential_bc(out, *self._bcs) + else: + for block, block_bcs in zip(out, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + +#=============================================================================== +class DirichletMultipatchBoundaryProjector(LinearOperator): + + def __init__(self, fem_space, bcs=None, space_kind=None): + + assert isinstance(fem_space, FemSpace) + assert fem_space.is_multipatch + + coeff_space = fem_space.coeff_space + self._domain = coeff_space + self._codomain = coeff_space + + if bcs is not None: + self._bcs = bcs + else: + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._domain + + @property + def dtype(self): + return None + + @property + def bcs(self): + return self._bcs + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + raise NotImplementedError + + def transpose(self, conjugate=False): + return self + + def _get_bcs(self, fem_space, space_kind=None): + """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + space = fem_space.symbolic_space + periodic = fem_space.periodic + + space_kind_str = space.kind.name + if space_kind is not None: + # Check whether kind is a valid input + if isinstance(space_kind, str): + kind_str = space_kind.lower() + assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) + elif isinstance(space_kind, SpaceType): + kind_str = space_kind.name + else: + raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') + + # If fem_space has a kind, it must be compatible with kind + if space_kind_str != 'undefined': + assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' + else: + # If space_kind_str = 'undefined': Update the variable using kind + space_kind_str = kind_str + + kind = space_kind_str + dim = space.domain.dim + assert dim==2 + + if kind == 'l2': + return None + + u = element_of(space, name="u") + + if kind == "h1": + bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + + elif kind == 'hcurl': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 0: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 1: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + + elif kind == 'hdiv': + bcs_x = [] + bcs_y = [] + + for bn in space.domain.boundary: + if bn.axis == 1: + bcs_y.append(EssentialBC(u, 0, bn, position=0)) + elif bn.axis == 0: + bcs_x.append(EssentialBC(u, 0, bn, position=0)) + + bcs = [bcs_x, bcs_y] + + else: + raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') + + return bcs + + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + + # apply bc on each patch + for p in out.blocks: + + if isinstance(p, StencilVector): + apply_essential_bc(p, *self._bcs) + else: + for block, block_bcs in zip(p, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 6eb3d62cf..481423ee9 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -1,28 +1,21 @@ -import time import numpy as np import pytest from sympy import sin, pi, sqrt, Tuple -from scipy.sparse import dia_matrix - from sympde.calculus import inner, cross -from sympde.expr import integral, LinearForm, BilinearForm, EssentialBC -from sympde.topology import element_of, elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace -from sympde.topology.datatype import SpaceType, H1Space, HcurlSpace +from sympde.expr import integral, LinearForm, BilinearForm +from sympde.topology import elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace +from sympde.topology.datatype import H1Space, HcurlSpace from psydac.api.discretization import discretize -from psydac.api.essential_bc import apply_essential_bc from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.fem.basic import FemSpace -from psydac.linalg.basic import LinearOperator, Vector, IdentityOperator -from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator -from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix +from psydac.fem.projectors import DirichletBoundaryProjector +from psydac.linalg.basic import LinearOperator, IdentityOperator +from psydac.linalg.block import BlockVectorSpace from psydac.linalg.solvers import inverse from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector -from psydac.linalg.tests.test_kron_direct_solver import matrix_to_bandsolver -from psydac.linalg.direct_solvers import BandedSolver def define_data_hermitian(n, p, dtype=float): domain_decomposition = DomainDecomposition([n - p], [False]) @@ -133,132 +126,6 @@ def _test_LO_equality_using_rng(A, B): assert err < 1e-15 -class DirichletBoundaryProjector(LinearOperator): - - def __init__(self, fem_space, bcs=None, space_kind=None): - - assert isinstance(fem_space, FemSpace) - - coeff_space = fem_space.coeff_space - self._domain = coeff_space - self._codomain = coeff_space - - if bcs is not None: - self._bcs = bcs - else: - self._bcs = self._get_bcs(fem_space, space_kind=space_kind) - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._domain - - @property - def dtype(self): - return None - - @property - def bcs(self): - return self._bcs - - def tosparse(self): - raise NotImplementedError - - def toarray(self): - raise NotImplementedError - - def transpose(self, conjugate=False): - return self - - def _get_bcs(self, fem_space, space_kind=None): - """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" - space = fem_space.symbolic_space - periodic = fem_space.periodic - - space_kind_str = space.kind.name - if space_kind is not None: - # Check whether kind is a valid input - if isinstance(space_kind, str): - kind_str = space_kind.lower() - assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) - elif isinstance(space_kind, SpaceType): - kind_str = space_kind.name - else: - raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') - - # If fem_space has a kind, it must be compatible with kind - if space_kind_str != 'undefined': - assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' - else: - # If space_kind_str = 'undefined': Update the variable using kind - space_kind_str = kind_str - - - kind = space_kind_str - dim = space.domain.dim - - if kind == 'l2': - return None - - u = element_of(space, name="u") - ebcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] - - if kind == "h1": - bcs = [ebcs[0], ebcs[1]] if periodic[0] == False else [] - if dim >= 2: - bcs += [ebcs[2], ebcs[3]] if periodic[1] == False else [] - if dim == 3: - bcs += [ebcs[4], ebcs[5]] if periodic[2] == False else [] - - elif kind == 'hcurl': - assert dim in (2, 3) - bcs_x = [ebcs[2], ebcs[3]] if periodic[1] == False else [] - if dim == 3: - bcs_x += [ebcs[4], ebcs[5]] if periodic[2] == False else [] - bcs_y = [ebcs[0], ebcs[1]] if periodic[0] == False else [] - if dim == 3: - bcs_y += [ebcs[4], ebcs[5]] if periodic[2] == False else [] - if dim == 3: - bcs_z = [ebcs[0], ebcs[1]] if periodic[0] == False else [] - bcs_z += [ebcs[2], ebcs[3]] if periodic[1] == False else [] - bcs = [bcs_x, bcs_y] - if dim == 3: - bcs.append(bcs_z) - - elif kind == 'hdiv': - assert dim in (2, 3) - bcs_x = [ebcs[0], ebcs[1]] if periodic[0] == False else [] - bcs_y = [ebcs[2], ebcs[3]] if periodic[1] == False else [] - if dim == 3: - bcs_z = [ebcs[4], ebcs[5]] if periodic[2] == False else [] - bcs = [bcs_x, bcs_y] - if dim == 3: - bcs.append(bcs_z) - - else: - raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') - - return bcs - - def dot(self, v, out=None): - if out is not None: - assert isinstance(out, Vector) - assert out.space is self.codomain - else: - out = self.codomain.zeros() - - v.copy(out=out) - if isinstance(self.domain, StencilVectorSpace): - apply_essential_bc(out, *self._bcs) - else: - for block, block_bcs in zip(out, self._bcs): - apply_essential_bc(block, *block_bcs) - - return out - #=============================================================================== @pytest.mark.parametrize( 'n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @@ -737,134 +604,6 @@ def test_discrete_derham_dirichlet_projector(): print() -#=============================================================================== -class DirichletMultipatchBoundaryProjector(LinearOperator): - - def __init__(self, fem_space, bcs=None, space_kind=None): - - assert isinstance(fem_space, FemSpace) - assert fem_space.is_multipatch - - coeff_space = fem_space.coeff_space - self._domain = coeff_space - self._codomain = coeff_space - - if bcs is not None: - self._bcs = bcs - else: - self._bcs = self._get_bcs(fem_space, space_kind=space_kind) - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._domain - - @property - def dtype(self): - return None - - @property - def bcs(self): - return self._bcs - - def tosparse(self): - raise NotImplementedError - - def toarray(self): - raise NotImplementedError - - def transpose(self, conjugate=False): - return self - - def _get_bcs(self, fem_space, space_kind=None): - """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" - space = fem_space.symbolic_space - periodic = fem_space.periodic - - space_kind_str = space.kind.name - if space_kind is not None: - # Check whether kind is a valid input - if isinstance(space_kind, str): - kind_str = space_kind.lower() - assert(kind_str in ['h1', 'hcurl', 'hdiv', 'l2', 'undefined']) - elif isinstance(space_kind, SpaceType): - kind_str = space_kind.name - else: - raise TypeError(f'Expecting space_kind {space_kind} to be a str or of SpaceType') - - # If fem_space has a kind, it must be compatible with kind - if space_kind_str != 'undefined': - assert space_kind_str == kind_str, f'fem_space and space_kind are not compatible.' - else: - # If space_kind_str = 'undefined': Update the variable using kind - space_kind_str = kind_str - - kind = space_kind_str - dim = space.domain.dim - assert dim==2 - - if kind == 'l2': - return None - - u = element_of(space, name="u") - - if kind == "h1": - bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] - - - elif kind == 'hcurl': - bcs_x = [] - bcs_y = [] - - for bn in space.domain.boundary: - if bn.axis == 0: - bcs_y.append(EssentialBC(u, 0, bn, position=0)) - elif bn.axis == 1: - bcs_x.append(EssentialBC(u, 0, bn, position=0)) - - bcs = [bcs_x, bcs_y] - - - elif kind == 'hdiv': - bcs_x = [] - bcs_y = [] - - for bn in space.domain.boundary: - if bn.axis == 1: - bcs_y.append(EssentialBC(u, 0, bn, position=0)) - elif bn.axis == 0: - bcs_x.append(EssentialBC(u, 0, bn, position=0)) - - bcs = [bcs_x, bcs_y] - - else: - raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') - - return bcs - - def dot(self, v, out=None): - if out is not None: - assert isinstance(out, Vector) - assert out.space is self.codomain - else: - out = self.codomain.zeros() - - v.copy(out=out) - - # apply bc on each patch - for p in out.blocks: - - if isinstance(p, StencilVector): - apply_essential_bc(p, *self._bcs) - else: - for block, block_bcs in zip(p, self._bcs): - apply_essential_bc(block, *block_bcs) - - return out - #=============================================================================== def test_discrete_derham_dirichlet_projector_multipatch(): From 9f6231e8cc13a4c381f52ad86981d91a9c789947 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 15:23:51 +0100 Subject: [PATCH 06/25] move import, add dirichlet_proj as attribute, return tuple, add identity proj --- psydac/api/feec.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 8b3e56283..fbeb96d61 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -33,6 +33,7 @@ from psydac.fem.basic import FemSpace, FemLinearOperator from psydac.fem.vector import VectorFemSpace +from psydac.fem.projectors import DirichletBoundaryProjector, DirichletMultipatchBoundaryProjector from psydac.linalg.basic import LinearOperator, IdentityOperator from psydac.linalg.block import BlockLinearOperator @@ -128,6 +129,7 @@ def __init__(self, domain_h, *spaces): self._hodge_operators = () self._conf_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property def dim(self): @@ -325,13 +327,15 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - from psydac.fem.projectors import DirichletBoundaryProjector - d_projectors = [DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + if not self._dirichlet_proj: + d_projectors_linop = tuple([DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + [IdentityOperator(self.spaces[-1].coeff_space)]) + d_projectors_femlinop = tuple([FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)]) + self._dirichlet_proj = d_projectors_femlinop if kind == 'femlinop': - d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)] - - return d_projectors + return self._dirichlet_proj + elif kind == 'linop': + return tuple([femlinop.linop for femlinop in self._dirichlet_proj]) #-------------------------------------------------------------------------- def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): @@ -583,6 +587,7 @@ def __init__(self, *, domain_h, spaces): self._hodge_operators = () self._conf_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property @@ -677,10 +682,12 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - from psydac.fem.projectors import DirichletMultipatchBoundaryProjector - d_projectors = [DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + if not self._dirichlet_proj: + d_projectors_linop = tuple([DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + IdentityOperator(self.spaces[-1].coeff_space)) + d_projectors_femlinop = tuple([FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)]) + self._dirichlet_proj = d_projectors_femlinop if kind == 'femlinop': - d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)] - - return d_projectors + return self._dirichlet_proj + elif kind == 'linop': + return tuple([femlinop.linop for femlinop in self._dirichlet_proj]) From 8bbb6b6affec3c5478d9dec9afe9c438ad64e6d9 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 15:25:33 +0100 Subject: [PATCH 07/25] remove periodic from _get_bcs in multipatch case --- psydac/fem/projectors.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 42758981b..0bb6cf4d6 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -330,7 +330,6 @@ def transpose(self, conjugate=False): def _get_bcs(self, fem_space, space_kind=None): """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" space = fem_space.symbolic_space - periodic = fem_space.periodic space_kind_str = space.kind.name if space_kind is not None: From 08f552536a3c9ad2202f4abf4bb199a5687878fc Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 15:37:13 +0100 Subject: [PATCH 08/25] change tolerance in _test_LO_equality... --- psydac/linalg/tests/test_solvers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 481423ee9..b638b7b13 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -122,7 +122,7 @@ def _test_LO_equality_using_rng(A, B): B.dot(x, out=y2) diff = y1 - y2 - err = A.codomain.inner(diff, diff) + err = np.sqrt(diff.inner(diff) / diff.space.dimension) assert err < 1e-15 From 07b14cb41169aaeb35da16c0092966aef1aba9e6 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 16:14:35 +0100 Subject: [PATCH 09/25] make dim a parameter, add sqrts, compute errs better --- psydac/linalg/tests/test_solvers.py | 498 ++++++++++++++-------------- 1 file changed, 247 insertions(+), 251 deletions(-) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index b638b7b13..237de35a7 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -127,7 +127,7 @@ def _test_LO_equality_using_rng(A, B): assert err < 1e-15 #=============================================================================== -@pytest.mark.parametrize( 'n', [5, 10, 13] ) +@pytest.mark.parametrize('n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @pytest.mark.parametrize('dtype', [float, complex]) @pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) @@ -276,7 +276,9 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): assert solver == 'pcg' or errc_norm < tol #=============================================================================== -def test_function_space_dirichlet_projector(): +@pytest.mark.parametrize('dim', [1, 2, 3]) + +def test_function_space_dirichlet_projector(dim): ncells_3d = [8, 8, 8] degree_3d = [2, 2, 2] @@ -295,155 +297,153 @@ def test_function_space_dirichlet_projector(): mapping_3d = SquareTorus ('ST') mappings = [mapping_1d, mapping_2d, mapping_3d] - dims = [1, 2, 3] rng = np.random.default_rng(42) print() - for dim in dims: - print(f' ----- Test projectors in dimension {dim} -----') - print() + print(f' ----- Test projectors in dimension {dim} -----') + print() - domain = mappings[dim-1](logical_domains[dim-1]) - from sympde.utilities.utils import plot_domain - #plot_domain(domain, draw=True, isolines=True) + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) - # Obtain "true" boundary, i.e., remove periodic y-direction boundary - if dim == 1: - boundary = domain.boundary - elif dim == 2: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) - else: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), - domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) - - ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] - degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] - periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] - - domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) - - nn = NormalVector('nn') - - for i in range(dim): - print(f' - Test DBP{i}') - - # The function defined here satisfy the corresponding homogeneous Dirichlet BCs - if dim == 1: - x = domain.coordinates - V = ScalarFunctionSpace('V', domain, kind='H1') - f = sin(2*pi*x) - if dim == 2: - x, y = domain.coordinates - if i == 0: - V = ScalarFunctionSpace('V', domain, kind=H1Space) - f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - else: - V = VectorFunctionSpace('V', domain, kind='hCuRl') - f1 = x - f2 = y - f = Tuple(f1, f2) - if dim == 3: - x, y, z = domain.coordinates - if i == 0: - V = ScalarFunctionSpace('V', domain, kind='h1') - f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) - elif i == 1: - V = VectorFunctionSpace('V', domain, kind=HcurlSpace) - f1 = z * (z - 1) * x - f2 = z * (z - 1) * y - f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f = Tuple(f1, f2, f3) - else: - V = VectorFunctionSpace('V', domain, kind='Hdiv') - f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f3 = z * (z-1) * sin(x*y) - f = Tuple(f1, f2, f3) + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] + degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] + periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + + nn = NormalVector('nn') - u, v = elements_of(V, names='u, v') + for i in range(dim): + print(f' - Test DBP{i}') + + # The function defined here satisfy the corresponding homogeneous Dirichlet BCs + if dim == 1: + x = domain.coordinates + V = ScalarFunctionSpace('V', domain, kind='H1') + f = sin(2*pi*x) + if dim == 2: + x, y = domain.coordinates if i == 0: - boundary_expr = u*v - if (i == 1) and (dim == 2): - boundary_expr = cross(nn, u) * cross(nn, v) - if (i == 1) and (dim == 3): - boundary_expr = inner(cross(nn, u), cross(nn, v)) - if i == 2: - boundary_expr = inner(nn, u) * inner(nn, v) - - Vh = discretize(V, domain_h, degree=degree) - expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v - - a = BilinearForm((u, v), integral(domain, expr)) - ab = BilinearForm((u, v), integral(boundary, boundary_expr)) - - ah = discretize(a, domain_h, (Vh, Vh), backend=backend) - abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) - - I = IdentityOperator(Vh.coeff_space) - DBP = DirichletBoundaryProjector(Vh) - - M = ah.assemble() - M_0 = DBP @ M @ DBP + (I - DBP) - Mb = abh.assemble() - - # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc - lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f - l = LinearForm(v, integral(domain, lexpr)) - lh = discretize(l, domain_h, Vh, backend=backend) - rhs = lh.assemble() - A = M + 1e30*Mb - A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) - fc = A_inv @ rhs - - # 1. - # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet - # boundary conditions should not change under application of the corresponding projector - fc2 = DBP @ fc - diff = fc - fc2 - err = np.linalg.norm(diff.toarray()) - print(f' | f - P @ f | = {err}') - assert err < 1e-15 - - # 2. - # After applying a projector to a random vector, we want to verify that the - # corresponding boundary integral vanishes - rdm_coeffs = Vh.coeff_space.zeros() - print(' Random boundary integrals:') - for _ in range(3): - if isinstance(rdm_coeffs.space, BlockVectorSpace): - for block in rdm_coeffs.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) - print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < 1e-15 - - # 3. - # We want to verify that applying a projector twice does not change the vector twice - fc3 = DBP @ fc2 - diff = fc2 - fc3 - err = np.linalg.norm(diff.toarray()) - print(f' | P @ f - P @ P @ f | = {err}') - assert err == 0. - - # 4. - # Finally, the modified mass matrix should still compute inner products correctly - l2_norm_squared = M.dot_inner (fc, fc) - l2_norm_squared2 = M_0.dot_inner(fc, fc) - diff = l2_norm_squared - l2_norm_squared2 - print(f' || f ||^2 = {l2_norm_squared} should be equal to') - print(f' || P @ f ||^2 = {l2_norm_squared2}') - assert diff < 1e-15 - - print() + V = ScalarFunctionSpace('V', domain, kind=H1Space) + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + else: + V = VectorFunctionSpace('V', domain, kind='hCuRl') + f1 = x + f2 = y + f = Tuple(f1, f2) + if dim == 3: + x, y, z = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind='h1') + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) + elif i == 1: + V = VectorFunctionSpace('V', domain, kind=HcurlSpace) + f1 = z * (z - 1) * x + f2 = z * (z - 1) * y + f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f = Tuple(f1, f2, f3) + else: + V = VectorFunctionSpace('V', domain, kind='Hdiv') + f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f3 = z * (z-1) * sin(x*y) + f = Tuple(f1, f2, f3) + + u, v = elements_of(V, names='u, v') + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + Vh = discretize(V, domain_h, degree=degree) + expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (Vh, Vh), backend=backend) + abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) + + I = IdentityOperator(Vh.coeff_space) + DBP = DirichletBoundaryProjector(Vh) + + M = ah.assemble() + M_0 = DBP @ M @ DBP + (I - DBP) + Mb = abh.assemble() + + # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc + lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f + l = LinearForm(v, integral(domain, lexpr)) + lh = discretize(l, domain_h, Vh, backend=backend) + rhs = lh.assemble() + A = M + 1e30*Mb + A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) + fc = A_inv @ rhs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DBP @ fc + diff = fc - fc2 + err = np.sqrt(diff.inner(diff)) + print(f' | f - P @ f | = {err}') + assert err < 1e-15 + + # 2. + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = Vh.coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DBP @ rdm_coeffs + boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) + boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < 1e-15 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DBP @ fc2 + diff = fc2 - fc3 + err = np.sqrt(diff.inner(diff)) + print(f' | P @ f - P @ P @ f | = {err}') + assert err == 0. + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm = np.sqrt(M.dot_inner (fc, fc)) + l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) + diff = l2_norm - l2_norm2 + print(f' || f || = {l2_norm} should be equal to') + print(f' || P @ f || = {l2_norm2}') + assert diff < 1e-15 print() #=============================================================================== -def test_discrete_derham_dirichlet_projector(): +@pytest.mark.parametrize('dim', [1, 2, 3]) + +def test_discrete_derham_dirichlet_projector(dim): ncells = [8, 8, 8] degree = [2, 2, 2] @@ -462,7 +462,6 @@ def test_discrete_derham_dirichlet_projector(): mapping_3d = SquareTorus ('ST') mappings = [mapping_1d, mapping_2d, mapping_3d] - dims = [1, 2, 3] rng = np.random.default_rng(42) # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs @@ -488,119 +487,116 @@ def test_discrete_derham_dirichlet_projector(): funs = [[f11], [f21, f22], [f31, f32, f33]] print() - for dim in dims: - print(f' ----- Test projectors in dimension {dim} -----') - print() + print(f' ----- Test projectors in dimension {dim} -----') + print() - domain = mappings[dim-1](logical_domains[dim-1]) - from sympde.utilities.utils import plot_domain - #plot_domain(domain, draw=True, isolines=True) + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) - # Obtain "true" boundary, i.e., remove periodic y-direction boundary - if dim == 1: - boundary = domain.boundary - elif dim == 2: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) - else: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), - domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] + degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] + periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] + + domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree_dim) + + db_projectors = derham_h.dirichlet_projectors(kind='linop') + + if dim == 2: + conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') - derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) + u, v = elements_of(derham.spaces[i], names='u, v') - ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] - degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] - periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) - domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) - derham_h = discretize(derham, domain_h, degree=degree_dim) + expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v - db_projectors = derham_h.dirichlet_projectors(kind='linop') + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DBP = db_projectors[i] if dim == 2: - conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + CP = conf_projectors[i] + _test_LO_equality_using_rng(DBP, CP) - nn = NormalVector('nn') + M = ah.assemble() + M_0 = DBP @ M @ DBP + (I - DBP) + Mb = abh.assemble() - for i in range(dim): - print(f' - Test DBP{i}') + f = funs[dim-1][i] + fc = derham_h.projectors()[i](f).coeffs - u, v = elements_of(derham.spaces[i], names='u, v') + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DBP @ fc + diff = fc - fc2 + err = np.sqrt(diff.inner(diff)) + print(f' | f - P @ f | = {err}') + assert err < 1e-15 - if i == 0: - boundary_expr = u*v - if (i == 1) and (dim == 2): - boundary_expr = cross(nn, u) * cross(nn, v) - if (i == 1) and (dim == 3): - boundary_expr = inner(cross(nn, u), cross(nn, v)) - if i == 2: - boundary_expr = inner(nn, u) * inner(nn, v) - - expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v - - a = BilinearForm((u, v), integral(domain, expr)) - ab = BilinearForm((u, v), integral(boundary, boundary_expr)) - - ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) - abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) - - I = IdentityOperator(derham_h.spaces[i].coeff_space) - DBP = db_projectors[i] - - if dim == 2: - CP = conf_projectors[i] - _test_LO_equality_using_rng(DBP, CP) - - M = ah.assemble() - M_0 = DBP @ M @ DBP + (I - DBP) - Mb = abh.assemble() - - f = funs[dim-1][i] - fc = derham_h.projectors()[i](f).coeffs - - # 1. - # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet - # boundary conditions should not change under application of the corresponding projector - fc2 = DBP @ fc - diff = fc - fc2 - err = np.linalg.norm(diff.toarray()) - print(f' | f - P @ f | = {err}') - assert err < 1e-15 - - # 2. - # After applying a projector to a random vector, we want to verify that the - # corresponding boundary integral vanishes - rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() - print(' Random boundary integrals:') - for _ in range(3): - if isinstance(rdm_coeffs.space, BlockVectorSpace): - for block in rdm_coeffs.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) - print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < 1e-15 - - # 3. - # We want to verify that applying a projector twice does not change the vector twice - fc3 = DBP @ fc2 - diff = fc2 - fc3 - err = np.linalg.norm(diff.toarray()) - print(f' | P @ f - P @ P @ f | = {err}') - assert err == 0. - - # 4. - # Finally, the modified mass matrix should still compute inner products correctly - l2_norm_squared = M.dot_inner (fc, fc) - l2_norm_squared2 = M_0.dot_inner(fc, fc) - diff = l2_norm_squared - l2_norm_squared2 - print(f' || f ||^2 = {l2_norm_squared} should be equal to') - print(f' || P @ f ||^2 = {l2_norm_squared2}') - assert diff < 1e-15 - - print() + # 2. + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DBP @ rdm_coeffs + boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) + boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < 1e-15 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DBP @ fc2 + diff = fc2 - fc3 + err = np.sqrt(diff.inner(diff)) + print(f' | P @ f - P @ P @ f | = {err}') + assert err == 0. + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm = np.sqrt(M.dot_inner (fc, fc)) + l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) + diff = l2_norm - l2_norm2 + print(f' || f || = {l2_norm} should be equal to') + print(f' || P @ f || = {l2_norm2}') + assert diff < 1e-15 print() @@ -679,7 +675,7 @@ def test_discrete_derham_dirichlet_projector_multipatch(): # boundary conditions should not change under application of the corresponding projector fc2 = DBP @ fc diff = fc - fc2 - err = np.linalg.norm(diff.toarray()) + err = np.sqrt(diff.inner(diff)) print(f' | f - P @ f | = {err}') assert err < 1e-15 @@ -698,8 +694,8 @@ def test_discrete_derham_dirichlet_projector_multipatch(): rng.random(size=patch._data.shape, dtype="float64", out=patch._data) rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) + boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) + boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') assert boundary_int_proj_rdm < 1e-15 @@ -707,17 +703,17 @@ def test_discrete_derham_dirichlet_projector_multipatch(): # We want to verify that applying a projector twice does not change the vector twice fc3 = DBP @ fc2 diff = fc2 - fc3 - err = np.linalg.norm(diff.toarray()) + err = np.sqrt(diff.inner(diff)) print(f' | P @ f - P @ P @ f | = {err}') assert err == 0. # 4. # Finally, the modified mass matrix should still compute inner products correctly - l2_norm_squared = M.dot_inner (fc, fc) - l2_norm_squared2 = M_0.dot_inner(fc, fc) - diff = l2_norm_squared - l2_norm_squared2 - print(f' || f ||^2 = {l2_norm_squared} should be equal to') - print(f' || P @ f ||^2 = {l2_norm_squared2}') + l2_norm = np.sqrt(M.dot_inner (fc, fc)) + l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) + diff = l2_norm - l2_norm2 + print(f' || f || = {l2_norm} should be equal to') + print(f' || P @ f || = {l2_norm2}') assert diff < 1e-15 print() From 8e1466ebfbb4c2f66897fe16af51001c2face565 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 16:14:52 +0100 Subject: [PATCH 10/25] add empty lines --- psydac/fem/projectors.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 0bb6cf4d6..9702c3253 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -28,6 +28,7 @@ def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): assert abs(intersection-coarse_grid).max()_{L2} i = 1, .. dim(Vh)) of a given function f, as a stencil array or numpy array @@ -158,6 +160,7 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste else: return tilde_f + #=============================================================================== class DirichletBoundaryProjector(LinearOperator): From 19b54caf0bbfc90307c64e3ec4fbddab67e958e2 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 16:21:21 +0100 Subject: [PATCH 11/25] include type checks --- psydac/fem/projectors.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 9702c3253..044ec2a4f 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,6 +1,7 @@ import numpy as np +from collections.abc import Iterable -from sympde.topology import element_of +from sympde.topology import element_of, Boundary from sympde.calculus import dot from sympde.expr import LinearForm, integral, EssentialBC from sympde.topology.datatype import SpaceType @@ -164,18 +165,21 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste #=============================================================================== class DirichletBoundaryProjector(LinearOperator): - def __init__(self, fem_space, bcs=None, space_kind=None): + def __init__(self, fem_space, *, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) + assert bcs is None or isinstance(bcs, Iterable) + assert space_kind is None or isinstance(space_kind, (str, SpaceType)) coeff_space = fem_space.coeff_space self._domain = coeff_space self._codomain = coeff_space if bcs is not None: - self._bcs = bcs + assert all(isinstance(bc, Boundary) for bc in bcs) + self._bcs = tuple(bcs) else: - self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + self._bcs = tuple(self._get_bcs(fem_space, space_kind=space_kind)) @property def domain(self): @@ -225,7 +229,6 @@ def _get_bcs(self, fem_space, space_kind=None): # If space_kind_str = 'undefined': Update the variable using kind space_kind_str = kind_str - kind = space_kind_str dim = space.domain.dim @@ -295,15 +298,18 @@ def __init__(self, fem_space, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) assert fem_space.is_multipatch + assert bcs is None or isinstance(bcs, Iterable) + assert space_kind is None or isinstance(space_kind, (str, SpaceType)) coeff_space = fem_space.coeff_space self._domain = coeff_space self._codomain = coeff_space if bcs is not None: - self._bcs = bcs + assert all(isinstance(bc, Boundary) for bc in bcs) + self._bcs = tuple(bcs) else: - self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + self._bcs = tuple(self._get_bcs(fem_space, space_kind=space_kind)) @property def domain(self): @@ -364,7 +370,6 @@ def _get_bcs(self, fem_space, space_kind=None): if kind == "h1": bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] - elif kind == 'hcurl': bcs_x = [] bcs_y = [] @@ -377,7 +382,6 @@ def _get_bcs(self, fem_space, space_kind=None): bcs = [bcs_x, bcs_y] - elif kind == 'hdiv': bcs_x = [] bcs_y = [] From ed03e2311237ad504bf538d222edea5408e83e5b Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 16:33:03 +0100 Subject: [PATCH 12/25] bug fix --- psydac/api/feec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index fbeb96d61..dfa286977 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -683,7 +683,7 @@ def dirichlet_projectors(self, kind='femlinop'): assert kind in ('femlinop', 'linop') if not self._dirichlet_proj: - d_projectors_linop = tuple([DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + IdentityOperator(self.spaces[-1].coeff_space)) + d_projectors_linop = tuple([DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + [IdentityOperator(self.spaces[-1].coeff_space)]) d_projectors_femlinop = tuple([FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)]) self._dirichlet_proj = d_projectors_femlinop From 051ebf6b230e1b5d869687e7347b1b83de36d49a Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 3 Nov 2025 16:54:55 +0100 Subject: [PATCH 13/25] add docstrings --- psydac/fem/projectors.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 044ec2a4f..fe08538cf 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -164,7 +164,25 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste #=============================================================================== class DirichletBoundaryProjector(LinearOperator): + """ + A LinearOperator that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. + bcs : Iterable + Iterable of sympde.topology.Boundary objects. Necessary only if fem_space.kind.name is undefined. + + space_kind : str | SpaceType + Necessary only if fem_space.kind.name is undefined. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of such an operator. + + """ def __init__(self, fem_space, *, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) @@ -293,7 +311,25 @@ def dot(self, v, out=None): #=============================================================================== class DirichletMultipatchBoundaryProjector(LinearOperator): + """ + A LinearOperator (for multipatch domains) that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. + bcs : Iterable + Iterable of sympde.topology.Boundary objects. Necessary only if fem_space.kind.name is undefined. + + space_kind : str | SpaceType + Necessary only if fem_space.kind.name is undefined. + + Notes + ----- + See examples/vector_potential_3d.py for a use case of such an operator. + + """ def __init__(self, fem_space, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) From e1de56aef84c72e7c918d9acd73c70352b7311c2 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Tue, 4 Nov 2025 13:41:03 +0100 Subject: [PATCH 14/25] addressing comments --- psydac/api/feec.py | 50 +++---- psydac/fem/projectors.py | 216 ++++++++++++++++------------ psydac/linalg/tests/test_solvers.py | 134 ++++++++--------- 3 files changed, 222 insertions(+), 178 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index dfa286977..c5b470e96 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -33,7 +33,7 @@ from psydac.fem.basic import FemSpace, FemLinearOperator from psydac.fem.vector import VectorFemSpace -from psydac.fem.projectors import DirichletBoundaryProjector, DirichletMultipatchBoundaryProjector +from psydac.fem.projectors import BoundaryProjector, MultipatchBoundaryProjector from psydac.linalg.basic import LinearOperator, IdentityOperator from psydac.linalg.block import BlockLinearOperator @@ -129,7 +129,7 @@ def __init__(self, domain_h, *spaces): self._hodge_operators = () self._conf_proj = () - self._dirichlet_proj = () + self._boundary_proj = () #-------------------------------------------------------------------------- @property def dim(self): @@ -303,9 +303,9 @@ def derivatives(self, kind='femlinop'): return tuple(b_diff.linop for b_diff in self._derivatives) #-------------------------------------------------------------------------- - def dirichlet_projectors(self, kind='femlinop'): + def boundary_projectors(self, kind='femlinop'): """ - Returns operators that apply the correct Dirichlet boundary conditions. + Returns operators that apply the correct homogeneous Dirichlet BCs. Parameters ---------- @@ -316,9 +316,9 @@ def dirichlet_projectors(self, kind='femlinop'): Returns ------- - d_projectors : list - List of or - The Dirichlet boundary projectors of each space and in desired form. + b_projectors : tuple + Tuple of or + The boundary projectors of each space and in desired form. Notes ----- @@ -327,15 +327,15 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - if not self._dirichlet_proj: - d_projectors_linop = tuple([DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + [IdentityOperator(self.spaces[-1].coeff_space)]) - d_projectors_femlinop = tuple([FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)]) - self._dirichlet_proj = d_projectors_femlinop + if not self._boundary_proj: + b_projectors_linop = tuple(BoundaryProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + b_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=b_projector) for Vh, b_projector in zip(self.spaces, b_projectors_linop)) + self._boundary_proj = b_projectors_femlinop if kind == 'femlinop': - return self._dirichlet_proj + return self._boundary_proj elif kind == 'linop': - return tuple([femlinop.linop for femlinop in self._dirichlet_proj]) + return tuple(femlinop.linop for femlinop in self._boundary_proj) #-------------------------------------------------------------------------- def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): @@ -587,7 +587,7 @@ def __init__(self, *, domain_h, spaces): self._hodge_operators = () self._conf_proj = () - self._dirichlet_proj = () + self._boundary_proj = () #-------------------------------------------------------------------------- @property @@ -658,9 +658,9 @@ def projectors(self, *, kind='global', nquads=None): raise NotImplementedError("3D projectors are not available") #-------------------------------------------------------------------------- - def dirichlet_projectors(self, kind='femlinop'): + def boundary_projectors(self, kind='femlinop'): """ - Returns operators that apply the correct Dirichlet boundary conditions. + Returns operators that apply the correct homogeneous Dirichlet BCs. Parameters ---------- @@ -671,9 +671,9 @@ def dirichlet_projectors(self, kind='femlinop'): Returns ------- - d_projectors : list - List of or - The Dirichlet boundary projectors of each space and in desired form. + b_projectors : tuple + Tuple of or + The boundary projectors of each space and in desired form. Notes ----- @@ -682,12 +682,12 @@ def dirichlet_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - if not self._dirichlet_proj: - d_projectors_linop = tuple([DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]] + [IdentityOperator(self.spaces[-1].coeff_space)]) - d_projectors_femlinop = tuple([FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)]) - self._dirichlet_proj = d_projectors_femlinop + if not self._boundary_proj: + b_projectors_linop = tuple(MultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + b_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=b_projector) for Vh, b_projector in zip(self.spaces, b_projectors_linop)) + self._boundary_proj = b_projectors_femlinop if kind == 'femlinop': - return self._dirichlet_proj + return self._boundary_proj elif kind == 'linop': - return tuple([femlinop.linop for femlinop in self._dirichlet_proj]) + return tuple(femlinop.linop for femlinop in self._boundary_proj) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index fe08538cf..823ff0b29 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -163,7 +163,7 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste #=============================================================================== -class DirichletBoundaryProjector(LinearOperator): +class BoundaryProjector(LinearOperator): """ A LinearOperator that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. @@ -172,11 +172,14 @@ class DirichletBoundaryProjector(LinearOperator): fem_space : psydac.fem.basic.FemSpace fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. - bcs : Iterable - Iterable of sympde.topology.Boundary objects. Necessary only if fem_space.kind.name is undefined. + bcs : Iterable | None + Iterable of sympde.topology.Boundary objects. + Allows the user to apply different kinds of BCs. + Must not be passed if a space_kind argument is passed. - space_kind : str | SpaceType - Necessary only if fem_space.kind.name is undefined. + space_kind : str | SpaceType | None + Necessary only if fem_space.kind.name is undefined and no bcs are passed. + Must not be passed if a space_kind argument is passed. Notes ----- @@ -188,6 +191,10 @@ def __init__(self, fem_space, *, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) assert bcs is None or isinstance(bcs, Iterable) assert space_kind is None or isinstance(space_kind, (str, SpaceType)) + if bcs is not None: + assert space_kind is None + if space_kind is not None: + assert bcs is None coeff_space = fem_space.coeff_space self._domain = coeff_space @@ -195,10 +202,13 @@ def __init__(self, fem_space, *, bcs=None, space_kind=None): if bcs is not None: assert all(isinstance(bc, Boundary) for bc in bcs) - self._bcs = tuple(bcs) + self._bcs = bcs else: - self._bcs = tuple(self._get_bcs(fem_space, space_kind=space_kind)) + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + #------------------------------------- + # Abstract interface + #------------------------------------- @property def domain(self): return self._domain @@ -210,22 +220,59 @@ def codomain(self): @property def dtype(self): return None - - @property - def bcs(self): - return self._bcs - + def tosparse(self): raise NotImplementedError def toarray(self): raise NotImplementedError + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + if isinstance(self.domain, StencilVectorSpace): + apply_essential_bc(out, *self._bcs) + else: + for block, block_bcs in zip(out, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + def transpose(self, conjugate=False): return self - def _get_bcs(self, fem_space, space_kind=None): - """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + @property + def bcs(self): + return self._bcs + + def _get_bcs(self, fem_space, *, space_kind=None): + """ + Returns a tuple of Boundries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.kind.name determines the kind of BCs returned by this function. + + space_kind : str | SpaceType | None + Optional. Must match fem_space.kind.name if that value is different from "undefined". + Must be passed if fem_space.kind.name has value "undefined". + + Returns + ------- + bcs : tuple + tuple of sympde.topology.Boundary. + Coefficients corresponding to functions non-zero on these boundaries will be set to 0. + + """ space = fem_space.symbolic_space periodic = fem_space.periodic @@ -291,26 +338,10 @@ def _get_bcs(self, fem_space, space_kind=None): else: raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') - return bcs - - def dot(self, v, out=None): - if out is not None: - assert isinstance(out, Vector) - assert out.space is self.codomain - else: - out = self.codomain.zeros() - - v.copy(out=out) - if isinstance(self.domain, StencilVectorSpace): - apply_essential_bc(out, *self._bcs) - else: - for block, block_bcs in zip(out, self._bcs): - apply_essential_bc(block, *block_bcs) - - return out + return tuple(bcs) #=============================================================================== -class DirichletMultipatchBoundaryProjector(LinearOperator): +class MultipatchBoundaryProjector(BoundaryProjector): """ A LinearOperator (for multipatch domains) that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. @@ -319,23 +350,33 @@ class DirichletMultipatchBoundaryProjector(LinearOperator): fem_space : psydac.fem.basic.FemSpace fem_space.coeff_space is domain and codomain of this LO. fem_space.kind.name determines the BCs that get applied. - bcs : Iterable - Iterable of sympde.topology.Boundary objects. Necessary only if fem_space.kind.name is undefined. + bcs : Iterable | None + Iterable of sympde.topology.Boundary objects. + Allows the user to apply different kinds of BCs. + Must not be passed if a space_kind argument is passed. - space_kind : str | SpaceType - Necessary only if fem_space.kind.name is undefined. + space_kind : str | SpaceType | None + Necessary only if fem_space.kind.name is undefined and no bcs are passed. + Must not be passed if a space_kind argument is passed. Notes ----- See examples/vector_potential_3d.py for a use case of such an operator. """ - def __init__(self, fem_space, bcs=None, space_kind=None): + def __init__(self, fem_space, *, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) assert fem_space.is_multipatch assert bcs is None or isinstance(bcs, Iterable) assert space_kind is None or isinstance(space_kind, (str, SpaceType)) + if bcs is not None: + assert space_kind is None + if space_kind is not None: + assert bcs is None + + dim = fem_space.symbolic_space.domain.dim + assert dim == 2, 'The class MultipatchBoundaryProjector is implemented only in 2D.' coeff_space = fem_space.coeff_space self._domain = coeff_space @@ -343,37 +384,56 @@ def __init__(self, fem_space, bcs=None, space_kind=None): if bcs is not None: assert all(isinstance(bc, Boundary) for bc in bcs) - self._bcs = tuple(bcs) + self._bcs = bcs else: - self._bcs = tuple(self._get_bcs(fem_space, space_kind=space_kind)) + self._bcs = self._get_bcs(fem_space, space_kind=space_kind) - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._domain - - @property - def dtype(self): - return None - - @property - def bcs(self): - return self._bcs - - def tosparse(self): - raise NotImplementedError - - def toarray(self): - raise NotImplementedError - - def transpose(self, conjugate=False): - return self - - def _get_bcs(self, fem_space, space_kind=None): - """Returns the correct Dirichlet boundary conditions for the passed fem_space.""" + #------------------------------------- + # Abstract interface + #------------------------------------- + def dot(self, v, out=None): + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + + # apply bc on each patch + for p in out.blocks: + + if isinstance(p, StencilVector): + apply_essential_bc(p, *self._bcs) + else: + for block, block_bcs in zip(p, self._bcs): + apply_essential_bc(block, *block_bcs) + + return out + + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + def _get_bcs(self, fem_space, *, space_kind=None): + """ + Returns a tuple of Boundries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + + Parameters + ---------- + fem_space : psydac.fem.basic.FemSpace + fem_space.kind.name determines the kind of BCs returned by this function. + + space_kind : str | SpaceType | None + Optional. Must match fem_space.kind.name if that value is different from "undefined". + Must be passed if fem_space.kind.name has value "undefined". + + Returns + ------- + bcs : tuple + tuple of sympde.topology.Boundary. + Coefficients corresponding to functions non-zero on these boundaries will be set to 0. + + """ space = fem_space.symbolic_space space_kind_str = space.kind.name @@ -395,8 +455,6 @@ def _get_bcs(self, fem_space, space_kind=None): space_kind_str = kind_str kind = space_kind_str - dim = space.domain.dim - assert dim==2 if kind == 'l2': return None @@ -433,24 +491,4 @@ def _get_bcs(self, fem_space, space_kind=None): else: raise ValueError(f'{kind} must be either "h1", "hcurl" or "hdiv"') - return bcs - - def dot(self, v, out=None): - if out is not None: - assert isinstance(out, Vector) - assert out.space is self.codomain - else: - out = self.codomain.zeros() - - v.copy(out=out) - - # apply bc on each patch - for p in out.blocks: - - if isinstance(p, StencilVector): - apply_essential_bc(p, *self._bcs) - else: - for block, block_bcs in zip(p, self._bcs): - apply_essential_bc(block, *block_bcs) - - return out + return tuple(bcs) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 237de35a7..946931ea8 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -11,7 +11,7 @@ from psydac.api.discretization import discretize from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.fem.projectors import DirichletBoundaryProjector +from psydac.fem.projectors import BoundaryProjector from psydac.linalg.basic import LinearOperator, IdentityOperator from psydac.linalg.block import BlockVectorSpace from psydac.linalg.solvers import inverse @@ -122,9 +122,10 @@ def _test_LO_equality_using_rng(A, B): B.dot(x, out=y2) diff = y1 - y2 - err = np.sqrt(diff.inner(diff) / diff.space.dimension) + err = diff.inner(diff) / diff.space.dimension**2 + tol = 1e-15 - assert err < 1e-15 + assert err < tol**2 #=============================================================================== @pytest.mark.parametrize('n', [5, 10, 13] ) @@ -278,7 +279,9 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): #=============================================================================== @pytest.mark.parametrize('dim', [1, 2, 3]) -def test_function_space_dirichlet_projector(dim): +def test_function_space_boundary_projector(dim): + + tol = 1e-15 ncells_3d = [8, 8, 8] degree_3d = [2, 2, 2] @@ -330,31 +333,31 @@ def test_function_space_dirichlet_projector(dim): # The function defined here satisfy the corresponding homogeneous Dirichlet BCs if dim == 1: x = domain.coordinates - V = ScalarFunctionSpace('V', domain, kind='H1') + V = ScalarFunctionSpace('V', domain, kind='H1') # testing various kind arguments f = sin(2*pi*x) if dim == 2: x, y = domain.coordinates if i == 0: - V = ScalarFunctionSpace('V', domain, kind=H1Space) + V = ScalarFunctionSpace('V', domain, kind=H1Space) # testing various kind arguments f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) else: - V = VectorFunctionSpace('V', domain, kind='hCuRl') + V = VectorFunctionSpace('V', domain, kind='hCuRl') # testing various kind arguments f1 = x f2 = y f = Tuple(f1, f2) if dim == 3: x, y, z = domain.coordinates if i == 0: - V = ScalarFunctionSpace('V', domain, kind='h1') + V = ScalarFunctionSpace('V', domain, kind='h1') # testing various kind arguments f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) elif i == 1: - V = VectorFunctionSpace('V', domain, kind=HcurlSpace) + V = VectorFunctionSpace('V', domain, kind=HcurlSpace) # testing various kind arguments f1 = z * (z - 1) * x f2 = z * (z - 1) * y f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) f = Tuple(f1, f2, f3) else: - V = VectorFunctionSpace('V', domain, kind='Hdiv') + V = VectorFunctionSpace('V', domain, kind='Hdiv') # testing various kind arguments f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) f3 = z * (z-1) * sin(x*y) @@ -380,7 +383,7 @@ def test_function_space_dirichlet_projector(dim): abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) I = IdentityOperator(Vh.coeff_space) - DBP = DirichletBoundaryProjector(Vh) + DBP = BoundaryProjector(Vh) M = ah.assemble() M_0 = DBP @ M @ DBP + (I - DBP) @@ -400,9 +403,9 @@ def test_function_space_dirichlet_projector(dim): # boundary conditions should not change under application of the corresponding projector fc2 = DBP @ fc diff = fc - fc2 - err = np.sqrt(diff.inner(diff)) + err = diff.inner(diff) print(f' | f - P @ f | = {err}') - assert err < 1e-15 + assert err < tol**2 # 2. # After applying a projector to a random vector, we want to verify that the @@ -416,34 +419,35 @@ def test_function_space_dirichlet_projector(dim): else: rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) - boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < 1e-15 + assert boundary_int_proj_rdm < tol**2 # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DBP @ fc2 diff = fc2 - fc3 - err = np.sqrt(diff.inner(diff)) - print(f' | P @ f - P @ P @ f | = {err}') - assert err == 0. + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 # 4. # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = np.sqrt(M.dot_inner (fc, fc)) - l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) - diff = l2_norm - l2_norm2 - print(f' || f || = {l2_norm} should be equal to') - print(f' || P @ f || = {l2_norm2}') - assert diff < 1e-15 + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol print() #=============================================================================== @pytest.mark.parametrize('dim', [1, 2, 3]) -def test_discrete_derham_dirichlet_projector(dim): +def test_discrete_derham_boundary_projector(dim): + + tol = 1e-15 ncells = [8, 8, 8] degree = [2, 2, 2] @@ -512,7 +516,7 @@ def test_discrete_derham_dirichlet_projector(dim): domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) derham_h = discretize(derham, domain_h, degree=degree_dim) - db_projectors = derham_h.dirichlet_projectors(kind='linop') + b_projectors = derham_h.boundary_projectors(kind='linop') if dim == 2: conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) @@ -542,7 +546,7 @@ def test_discrete_derham_dirichlet_projector(dim): abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) I = IdentityOperator(derham_h.spaces[i].coeff_space) - DBP = db_projectors[i] + DBP = b_projectors[i] if dim == 2: CP = conf_projectors[i] @@ -560,9 +564,9 @@ def test_discrete_derham_dirichlet_projector(dim): # boundary conditions should not change under application of the corresponding projector fc2 = DBP @ fc diff = fc - fc2 - err = np.sqrt(diff.inner(diff)) - print(f' | f - P @ f | = {err}') - assert err < 1e-15 + err = diff.inner(diff) + print(f' | f - P @ f |^2 = {err}') + assert err < tol**2 # 2. # After applying a projector to a random vector, we want to verify that the @@ -576,32 +580,35 @@ def test_discrete_derham_dirichlet_projector(dim): else: rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) - boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < 1e-15 + assert boundary_int_proj_rdm < tol**2 # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DBP @ fc2 diff = fc2 - fc3 - err = np.sqrt(diff.inner(diff)) - print(f' | P @ f - P @ P @ f | = {err}') - assert err == 0. + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 # 4. # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = np.sqrt(M.dot_inner (fc, fc)) - l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) - diff = l2_norm - l2_norm2 - print(f' || f || = {l2_norm} should be equal to') - print(f' || P @ f || = {l2_norm2}') - assert diff < 1e-15 + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + print(f' || f ||^2 = {l2_norm} should be equal to') + print(f' || P @ f ||^2 = {l2_norm2}') + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol print() #=============================================================================== -def test_discrete_derham_dirichlet_projector_multipatch(): +def test_discrete_derham_boundary_projector_multipatch(): + + tol = 1e-15 ncells = [8, 8] degree = [2, 2] @@ -628,7 +635,7 @@ def test_discrete_derham_dirichlet_projector_multipatch(): derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) ncells_h = {} - for k, D in enumerate(domain.interior): + for D in domain.interior: ncells_h[D.name] = ncells domain_h = discretize(domain, ncells=ncells_h, comm=comm) @@ -636,9 +643,7 @@ def test_discrete_derham_dirichlet_projector_multipatch(): projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) - db_projectors = derham_h.dirichlet_projectors(kind='linop') - - conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + b_projectors = derham_h.boundary_projectors(kind='linop') nn = NormalVector('nn') @@ -661,7 +666,7 @@ def test_discrete_derham_dirichlet_projector_multipatch(): abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) I = IdentityOperator(derham_h.spaces[i].coeff_space) - DBP = db_projectors[i] + DBP = b_projectors[i] M = ah.assemble() M_0 = DBP @ M @ DBP + (I - DBP) @@ -675,9 +680,9 @@ def test_discrete_derham_dirichlet_projector_multipatch(): # boundary conditions should not change under application of the corresponding projector fc2 = DBP @ fc diff = fc - fc2 - err = np.sqrt(diff.inner(diff)) - print(f' | f - P @ f | = {err}') - assert err < 1e-15 + err = diff.inner(diff) + print(f' | f - P @ f |^2 = {err}') + assert err < tol**2 # 2. # After applying a projector to a random vector, we want to verify that the @@ -694,27 +699,28 @@ def test_discrete_derham_dirichlet_projector_multipatch(): rng.random(size=patch._data.shape, dtype="float64", out=patch._data) rdm_coeffs2 = DBP @ rdm_coeffs - boundary_int_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension) - boundary_int_proj_rdm = np.sqrt(Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension) + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < 1e-15 + assert boundary_int_proj_rdm < tol**2 # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DBP @ fc2 diff = fc2 - fc3 - err = np.sqrt(diff.inner(diff)) - print(f' | P @ f - P @ P @ f | = {err}') - assert err == 0. + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 # 4. # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = np.sqrt(M.dot_inner (fc, fc)) - l2_norm2 = np.sqrt(M_0.dot_inner(fc, fc)) - diff = l2_norm - l2_norm2 - print(f' || f || = {l2_norm} should be equal to') - print(f' || P @ f || = {l2_norm2}') - assert diff < 1e-15 + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + print(f' || f ||^2 = {l2_norm} should be equal to') + print(f' || P @ f ||^2 = {l2_norm2}') + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol print() From bf1d441444af494beaf4c4f70954d680ccf8f65c Mon Sep 17 00:00:00 2001 From: jowezarek Date: Tue, 4 Nov 2025 13:46:26 +0100 Subject: [PATCH 15/25] Codacity fix --- psydac/fem/projectors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 823ff0b29..48acce72d 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -15,7 +15,7 @@ from psydac.linalg.stencil import StencilVectorSpace, StencilVector __all__ = ('knots_to_insert', 'knot_insertion_projection_operator', 'get_dual_dofs', - 'DirichletBoundaryProjector', 'DirichletMultipatchBoundaryProjector') + 'BoundaryProjector', 'MultipatchBoundaryProjector') def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): """ Compute the point difference between the fine grid and coarse grid.""" From daaabd679ebf003929189c3a77e3aad0ba5eb7d1 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Tue, 4 Nov 2025 16:10:06 +0100 Subject: [PATCH 16/25] addressing comments pt. 2 --- psydac/api/feec.py | 59 ++++++++++++----------------- psydac/fem/projectors.py | 10 ++--- psydac/linalg/tests/test_solvers.py | 38 +++++++++---------- 3 files changed, 48 insertions(+), 59 deletions(-) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index c5b470e96..8f68e2eee 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -1,10 +1,3 @@ -import numpy as np - -from scipy.sparse import dia_matrix - -from sympde.expr import integral, BilinearForm -from sympde.topology import elements_of, Line, Derham - from psydac.api.basic import BasicDiscrete from psydac.feec.derivatives import Derivative1D, Gradient2D, Gradient3D @@ -33,13 +26,9 @@ from psydac.fem.basic import FemSpace, FemLinearOperator from psydac.fem.vector import VectorFemSpace -from psydac.fem.projectors import BoundaryProjector, MultipatchBoundaryProjector +from psydac.fem.projectors import DirichletProjector, MultipatchDirichletProjector -from psydac.linalg.basic import LinearOperator, IdentityOperator -from psydac.linalg.block import BlockLinearOperator -from psydac.linalg.direct_solvers import BandedSolver -from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix -from psydac.linalg.stencil import StencilVectorSpace +from psydac.linalg.basic import IdentityOperator __all__ = ('DiscreteDeRham', 'DiscreteDeRhamMultipatch',) @@ -129,7 +118,7 @@ def __init__(self, domain_h, *spaces): self._hodge_operators = () self._conf_proj = () - self._boundary_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property def dim(self): @@ -303,7 +292,7 @@ def derivatives(self, kind='femlinop'): return tuple(b_diff.linop for b_diff in self._derivatives) #-------------------------------------------------------------------------- - def boundary_projectors(self, kind='femlinop'): + def dirichlet_projectors(self, kind='femlinop'): """ Returns operators that apply the correct homogeneous Dirichlet BCs. @@ -316,9 +305,9 @@ def boundary_projectors(self, kind='femlinop'): Returns ------- - b_projectors : tuple - Tuple of or - The boundary projectors of each space and in desired form. + d_projectors : tuple + Tuple of or + The Dirichlet projectors of each space and in desired form. Notes ----- @@ -327,15 +316,15 @@ def boundary_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - if not self._boundary_proj: - b_projectors_linop = tuple(BoundaryProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) - b_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=b_projector) for Vh, b_projector in zip(self.spaces, b_projectors_linop)) - self._boundary_proj = b_projectors_femlinop + if not self._dirichlet_proj: + d_projectors_linop = tuple(DirichletProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + d_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)) + self._dirichlet_proj = d_projectors_femlinop if kind == 'femlinop': - return self._boundary_proj + return self._dirichlet_proj elif kind == 'linop': - return tuple(femlinop.linop for femlinop in self._boundary_proj) + return tuple(femlinop.linop for femlinop in self._dirichlet_proj) #-------------------------------------------------------------------------- def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): @@ -587,7 +576,7 @@ def __init__(self, *, domain_h, spaces): self._hodge_operators = () self._conf_proj = () - self._boundary_proj = () + self._dirichlet_proj = () #-------------------------------------------------------------------------- @property @@ -658,7 +647,7 @@ def projectors(self, *, kind='global', nquads=None): raise NotImplementedError("3D projectors are not available") #-------------------------------------------------------------------------- - def boundary_projectors(self, kind='femlinop'): + def dirichlet_projectors(self, kind='femlinop'): """ Returns operators that apply the correct homogeneous Dirichlet BCs. @@ -671,9 +660,9 @@ def boundary_projectors(self, kind='femlinop'): Returns ------- - b_projectors : tuple - Tuple of or - The boundary projectors of each space and in desired form. + d_projectors : tuple + Tuple of or + The Dirichlet projectors of each space and in desired form. Notes ----- @@ -682,12 +671,12 @@ def boundary_projectors(self, kind='femlinop'): """ assert kind in ('femlinop', 'linop') - if not self._boundary_proj: - b_projectors_linop = tuple(MultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) - b_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=b_projector) for Vh, b_projector in zip(self.spaces, b_projectors_linop)) - self._boundary_proj = b_projectors_femlinop + if not self._dirichlet_proj: + d_projectors_linop = tuple(MultipatchDirichletProjector(Vh) for Vh in self.spaces[:-1]) + (IdentityOperator(self.spaces[-1].coeff_space),) + d_projectors_femlinop = tuple(FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces, d_projectors_linop)) + self._dirichlet_proj = d_projectors_femlinop if kind == 'femlinop': - return self._boundary_proj + return self._dirichlet_proj elif kind == 'linop': - return tuple(femlinop.linop for femlinop in self._boundary_proj) + return tuple(femlinop.linop for femlinop in self._dirichlet_proj) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 48acce72d..c473f9fb4 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -15,7 +15,7 @@ from psydac.linalg.stencil import StencilVectorSpace, StencilVector __all__ = ('knots_to_insert', 'knot_insertion_projection_operator', 'get_dual_dofs', - 'BoundaryProjector', 'MultipatchBoundaryProjector') + 'DirichletProjector', 'MultipatchDirichletProjector') def knots_to_insert(coarse_grid, fine_grid, tol=1e-14): """ Compute the point difference between the fine grid and coarse grid.""" @@ -163,7 +163,7 @@ def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='ste #=============================================================================== -class BoundaryProjector(LinearOperator): +class DirichletProjector(LinearOperator): """ A LinearOperator that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. @@ -298,7 +298,7 @@ def _get_bcs(self, fem_space, *, space_kind=None): dim = space.domain.dim if kind == 'l2': - return None + return () u = element_of(space, name="u") ebcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] @@ -341,7 +341,7 @@ def _get_bcs(self, fem_space, *, space_kind=None): return tuple(bcs) #=============================================================================== -class MultipatchBoundaryProjector(BoundaryProjector): +class MultipatchDirichletProjector(DirichletProjector): """ A LinearOperator (for multipatch domains) that applies homogeneous (unless manually given different bcs) Dirichlet boundary conditions. @@ -457,7 +457,7 @@ def _get_bcs(self, fem_space, *, space_kind=None): kind = space_kind_str if kind == 'l2': - return None + return () u = element_of(space, name="u") diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 946931ea8..b6eb8f539 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -11,7 +11,7 @@ from psydac.api.discretization import discretize from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.fem.projectors import BoundaryProjector +from psydac.fem.projectors import DirichletProjector from psydac.linalg.basic import LinearOperator, IdentityOperator from psydac.linalg.block import BlockVectorSpace from psydac.linalg.solvers import inverse @@ -383,10 +383,10 @@ def test_function_space_boundary_projector(dim): abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) I = IdentityOperator(Vh.coeff_space) - DBP = BoundaryProjector(Vh) + DP = DirichletProjector(Vh) M = ah.assemble() - M_0 = DBP @ M @ DBP + (I - DBP) + M_0 = DP @ M @ DP + (I - DP) Mb = abh.assemble() # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc @@ -401,7 +401,7 @@ def test_function_space_boundary_projector(dim): # 1. # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet # boundary conditions should not change under application of the corresponding projector - fc2 = DBP @ fc + fc2 = DP @ fc diff = fc - fc2 err = diff.inner(diff) print(f' | f - P @ f | = {err}') @@ -418,7 +418,7 @@ def test_function_space_boundary_projector(dim): rng.random(size=block._data.shape, dtype="float64", out=block._data) else: rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DBP @ rdm_coeffs + rdm_coeffs2 = DP @ rdm_coeffs boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') @@ -426,7 +426,7 @@ def test_function_space_boundary_projector(dim): # 3. # We want to verify that applying a projector twice does not change the vector twice - fc3 = DBP @ fc2 + fc3 = DP @ fc2 diff = fc2 - fc3 err = diff.inner(diff) print(f' | P @ f - P @ P @ f |^2 = {err}') @@ -516,7 +516,7 @@ def test_discrete_derham_boundary_projector(dim): domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) derham_h = discretize(derham, domain_h, degree=degree_dim) - b_projectors = derham_h.boundary_projectors(kind='linop') + d_projectors = derham_h.dirichlet_projectors(kind='linop') if dim == 2: conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) @@ -546,14 +546,14 @@ def test_discrete_derham_boundary_projector(dim): abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) I = IdentityOperator(derham_h.spaces[i].coeff_space) - DBP = b_projectors[i] + DP = d_projectors[i] if dim == 2: CP = conf_projectors[i] - _test_LO_equality_using_rng(DBP, CP) + _test_LO_equality_using_rng(DP, CP) M = ah.assemble() - M_0 = DBP @ M @ DBP + (I - DBP) + M_0 = DP @ M @ DP + (I - DP) Mb = abh.assemble() f = funs[dim-1][i] @@ -562,7 +562,7 @@ def test_discrete_derham_boundary_projector(dim): # 1. # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet # boundary conditions should not change under application of the corresponding projector - fc2 = DBP @ fc + fc2 = DP @ fc diff = fc - fc2 err = diff.inner(diff) print(f' | f - P @ f |^2 = {err}') @@ -579,7 +579,7 @@ def test_discrete_derham_boundary_projector(dim): rng.random(size=block._data.shape, dtype="float64", out=block._data) else: rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DBP @ rdm_coeffs + rdm_coeffs2 = DP @ rdm_coeffs boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') @@ -587,7 +587,7 @@ def test_discrete_derham_boundary_projector(dim): # 3. # We want to verify that applying a projector twice does not change the vector twice - fc3 = DBP @ fc2 + fc3 = DP @ fc2 diff = fc2 - fc3 err = diff.inner(diff) print(f' | P @ f - P @ P @ f |^2 = {err}') @@ -643,7 +643,7 @@ def test_discrete_derham_boundary_projector_multipatch(): projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) - b_projectors = derham_h.boundary_projectors(kind='linop') + d_projectors = derham_h.dirichlet_projectors(kind='linop') nn = NormalVector('nn') @@ -666,10 +666,10 @@ def test_discrete_derham_boundary_projector_multipatch(): abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) I = IdentityOperator(derham_h.spaces[i].coeff_space) - DBP = b_projectors[i] + DP = d_projectors[i] M = ah.assemble() - M_0 = DBP @ M @ DBP + (I - DBP) + M_0 = DP @ M @ DP + (I - DP) Mb = abh.assemble() f = funs[i] @@ -678,7 +678,7 @@ def test_discrete_derham_boundary_projector_multipatch(): # 1. # The coefficients of functions satisfying homogeneous Dirichlet # boundary conditions should not change under application of the corresponding projector - fc2 = DBP @ fc + fc2 = DP @ fc diff = fc - fc2 err = diff.inner(diff) print(f' | f - P @ f |^2 = {err}') @@ -698,7 +698,7 @@ def test_discrete_derham_boundary_projector_multipatch(): else: rng.random(size=patch._data.shape, dtype="float64", out=patch._data) - rdm_coeffs2 = DBP @ rdm_coeffs + rdm_coeffs2 = DP @ rdm_coeffs boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') @@ -706,7 +706,7 @@ def test_discrete_derham_boundary_projector_multipatch(): # 3. # We want to verify that applying a projector twice does not change the vector twice - fc3 = DBP @ fc2 + fc3 = DP @ fc2 diff = fc2 - fc3 err = diff.inner(diff) print(f' | P @ f - P @ P @ f |^2 = {err}') From 2cf7a3b83c9215f655f34026636b49456a915534 Mon Sep 17 00:00:00 2001 From: Frederik Schnack Date: Wed, 5 Nov 2025 12:55:11 +0100 Subject: [PATCH 17/25] add docstrings to conforming projections, change DiscreteDeRhamMultipatch to MultipatchDiscreteDeRham --- psydac/api/discretization.py | 8 ++++---- psydac/api/feec.py | 26 +++++++++++++++++--------- psydac/feec/conforming_projectors.py | 17 +++++++++++++++++ 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index d98ba1a51..1d9f8d66f 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -30,7 +30,7 @@ from psydac.api.fem import DiscreteBilinearForm from psydac.api.fem import DiscreteLinearForm from psydac.api.fem import DiscreteFunctional -from psydac.api.feec import DiscreteDeRham, DiscreteDeRhamMultipatch +from psydac.api.feec import DiscreteDeRham, MultipatchDiscreteDeRham from psydac.api.glt import DiscreteGltExpr from psydac.api.expr import DiscreteExpr from psydac.api.equation import DiscreteEquation @@ -199,7 +199,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): Create a discrete multipatch de Rham sequence from a symbolic one. This function creates the broken discrete spaces from the symbolic ones, and then - creates a DiscreteDeRhamMultipatch object from them. + creates a MultipatchDiscreteDeRham object from them. Parameters ---------- @@ -214,7 +214,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): Returns ------- - DiscreteDeRhamMultipatch + MultipatchDiscreteDeRham The discrete multipatch de Rham sequence containing the discrete spaces, differential operators and projectors. @@ -229,7 +229,7 @@ def discretize_derham_multipatch(derham, domain_h, **kwargs): spaces = [discretize_space(V, domain_h, basis=basis, **kwargs) \ for V, basis in zip(derham.spaces, bases)] - return DiscreteDeRhamMultipatch( + return MultipatchDiscreteDeRham( domain_h = domain_h, spaces = spaces ) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 8f68e2eee..439c2f81f 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -30,7 +30,7 @@ from psydac.linalg.basic import IdentityOperator -__all__ = ('DiscreteDeRham', 'DiscreteDeRhamMultipatch',) +__all__ = ('DiscreteDeRham', 'MultipatchDiscreteDeRham',) #============================================================================== class DiscreteDeRham(BasicDiscrete): @@ -286,6 +286,8 @@ def projectors(self, *, kind='global', nquads=None): #-------------------------------------------------------------------------- def derivatives(self, kind='femlinop'): + assert kind in ('femlinop', 'linop') + if kind == 'femlinop': return self._derivatives elif kind == 'linop': @@ -333,24 +335,28 @@ def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, h 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' or 'linop'. - 'femlinop' returns a psydac FemLinearOperator (default) - 'linop' returns a psydac LinearOperator + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. + + p_moments: + Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) + + hom_bc: + Apply homogenous boundary conditions if True + Returns ------- cP0, cP1, cP2 : Tuple of or The conforming projectors of each space and in desired form. """ + assert kind in ('femlinop', 'linop') if hom_bc is None: raise ValueError('please provide a value for "hom_bc" argument') @@ -480,6 +486,7 @@ def hodge_operator(self, space=None, dual=False, kind='femlinop', backend_langua H : or """ + assert kind in ('femlinop', 'linop') if not self._hodge_operators: self._init_hodge_operators(backend_language=backend_language) @@ -518,6 +525,7 @@ def hodge_operators(self, dual=False, kind='femlinop', backend_language='python' ------- The Hodge operators of all spaces and of the specified kind. """ + assert kind in ('femlinop', 'linop') if not self._hodge_operators: self._init_hodge_operators(backend_language=backend_language) @@ -526,7 +534,7 @@ def hodge_operators(self, dual=False, kind='femlinop', backend_language='python' #============================================================================== -class DiscreteDeRhamMultipatch(DiscreteDeRham): +class MultipatchDiscreteDeRham(DiscreteDeRham): """ Represents the discrete de Rham sequence for multipatch domains. It only works when the number of patches>1. diff --git a/psydac/feec/conforming_projectors.py b/psydac/feec/conforming_projectors.py index 638c33e15..1ac1e52a8 100644 --- a/psydac/feec/conforming_projectors.py +++ b/psydac/feec/conforming_projectors.py @@ -1211,6 +1211,10 @@ def construct_h1_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=- def get_vertex_index(coords): + """ + Calculate the global index of the vertex basis function + from the geometric coordinates of a vertex in the domain. + """ nbasis0 = Vh.spaces[0].nbasis - 1 nbasis1 = Vh.spaces[1].nbasis - 1 @@ -1223,6 +1227,11 @@ def get_vertex_index(coords): return l2g.get_index(0, 0, multi_index) def vertex_moment_indices(axis, coords, p_moments): + """ + Calculate the global indices of the basis functions + adjacent to the vertex basis function along axis + from the geometric coordinates of a vertex in the domain. + """ if coords[axis] == 0: return range(1, p_moments + 2) else: @@ -1441,9 +1450,13 @@ class ConformingProjectionV0(FemLinearOperator): ---------- V0h: The discrete space + + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. p_moments: Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) hom_bc : Apply homogenous boundary conditions if True @@ -1482,8 +1495,12 @@ class ConformingProjectionV1(FemLinearOperator): V1h: The discrete space + mom_pres: + If True, preserve polynomial moments of maximal order in the projection. + p_moments: Number of polynomial moments to be preserved in the projection. + (Gets overwritten if the parameter mom_pres equals True) hom_bc : Apply homogenous boundary conditions if True From c803120c4d99fe44ab7bf1513af492daa52fc584 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 7 Nov 2025 11:33:14 +0100 Subject: [PATCH 18/25] add toarray, tosparse, and tests --- psydac/fem/projectors.py | 15 +++++++++++++-- psydac/linalg/tests/test_solvers.py | 30 ++++++++++++++++++++++++++--- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index c473f9fb4..80dbd6b28 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,6 +1,8 @@ import numpy as np from collections.abc import Iterable +from scipy.sparse import diags + from sympde.topology import element_of, Boundary from sympde.calculus import dot from sympde.expr import LinearForm, integral, EssentialBC @@ -13,6 +15,7 @@ from psydac.linalg.basic import LinearOperator, Vector from psydac.linalg.kron import KroneckerDenseMatrix from psydac.linalg.stencil import StencilVectorSpace, StencilVector +from psydac.linalg.utilities import array_to_psydac __all__ = ('knots_to_insert', 'knot_insertion_projection_operator', 'get_dual_dofs', 'DirichletProjector', 'MultipatchDirichletProjector') @@ -222,10 +225,18 @@ def dtype(self): return None def tosparse(self): - raise NotImplementedError + np_ones = np.ones(shape=self.domain.dimension) + psy_ones = array_to_psydac(np_ones, self.domain) + diagonal = self @ psy_ones + sparse = diags(diagonal.toarray()) + return sparse def toarray(self): - raise NotImplementedError + np_ones = np.ones(shape=self.domain.dimension) + psy_ones = array_to_psydac(np_ones, self.domain) + diagonal = self @ psy_ones + array = np.diag(diagonal.toarray()) + return array def dot(self, v, out=None): if out is not None: diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index b6eb8f539..9316d37a6 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -407,7 +407,7 @@ def test_function_space_boundary_projector(dim): print(f' | f - P @ f | = {err}') assert err < tol**2 - # 2. + # 2.1 # After applying a projector to a random vector, we want to verify that the # corresponding boundary integral vanishes rdm_coeffs = Vh.coeff_space.zeros() @@ -424,6 +424,14 @@ def test_function_space_boundary_projector(dim): print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') assert boundary_int_proj_rdm < tol**2 + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DP @ fc2 @@ -568,7 +576,7 @@ def test_discrete_derham_boundary_projector(dim): print(f' | f - P @ f |^2 = {err}') assert err < tol**2 - # 2. + # 2.1 # After applying a projector to a random vector, we want to verify that the # corresponding boundary integral vanishes rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() @@ -585,6 +593,14 @@ def test_discrete_derham_boundary_projector(dim): print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') assert boundary_int_proj_rdm < tol**2 + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DP @ fc2 @@ -684,7 +700,7 @@ def test_discrete_derham_boundary_projector_multipatch(): print(f' | f - P @ f |^2 = {err}') assert err < tol**2 - # 2. + # 2.1 # After applying a projector to a random vector, we want to verify that the # corresponding boundary integral vanishes rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() @@ -704,6 +720,14 @@ def test_discrete_derham_boundary_projector_multipatch(): print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') assert boundary_int_proj_rdm < tol**2 + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + # 3. # We want to verify that applying a projector twice does not change the vector twice fc3 = DP @ fc2 From 01beb17d13503c9d2e712f729f02ac5aaf97b1aa Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 7 Nov 2025 11:43:58 +0100 Subject: [PATCH 19/25] move dirichlet projector tests to new file --- psydac/fem/tests/test_dirichlet_projectors.py | 560 ++++++++++++++++++ psydac/linalg/tests/test_solvers.py | 546 ----------------- 2 files changed, 560 insertions(+), 546 deletions(-) create mode 100644 psydac/fem/tests/test_dirichlet_projectors.py diff --git a/psydac/fem/tests/test_dirichlet_projectors.py b/psydac/fem/tests/test_dirichlet_projectors.py new file mode 100644 index 000000000..7da37f787 --- /dev/null +++ b/psydac/fem/tests/test_dirichlet_projectors.py @@ -0,0 +1,560 @@ +import numpy as np +import pytest + +from sympy import sin, pi, sqrt, Tuple + +from sympde.calculus import inner, cross +from sympde.expr import integral, LinearForm, BilinearForm +from sympde.topology import elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace +from sympde.topology.datatype import H1Space, HcurlSpace + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.fem.projectors import DirichletProjector +from psydac.linalg.basic import LinearOperator, IdentityOperator +from psydac.linalg.block import BlockVectorSpace +from psydac.linalg.solvers import inverse + +class SquareTorus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + +class Annulus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)'} + + _ldim = 2 + _pdim = 2 + +class SinMapping1D(Mapping): + + _expressions = {'x': 'sin((pi/2)*x1)'} + + _ldim = 1 + _pdim = 1 + +def _test_LO_equality_using_rng(A, B): + """ + A simple tool to check with almost certainty that two linear operators are identical, + by applying them repeatedly to random vectors. + + """ + + assert isinstance(A, LinearOperator) + assert isinstance(B, LinearOperator) + assert A.domain is B.domain + assert A.codomain is B.codomain + + rng = np.random.default_rng(42) + + x = A.domain.zeros() + y1 = A.codomain.zeros() + y2 = y1.copy() + + n = 10 + + for _ in range(n): + + x *= 0. + + if isinstance(A.domain, BlockVectorSpace): + for block in x.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=x._data.shape, dtype="float64", out=x._data) + + A.dot(x, out=y1) + B.dot(x, out=y2) + + diff = y1 - y2 + err = diff.inner(diff) / diff.space.dimension**2 + tol = 1e-15 + + assert err < tol**2 + +#=============================================================================== +@pytest.mark.parametrize('dim', [1, 2, 3]) + +def test_function_space_boundary_projector(dim): + + tol = 1e-15 + + ncells_3d = [8, 8, 8] + degree_3d = [2, 2, 2] + periodic_3d = [False, True, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + rng = np.random.default_rng(42) + + print() + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] + degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] + periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + # The function defined here satisfy the corresponding homogeneous Dirichlet BCs + if dim == 1: + x = domain.coordinates + V = ScalarFunctionSpace('V', domain, kind='H1') # testing various kind arguments + f = sin(2*pi*x) + if dim == 2: + x, y = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind=H1Space) # testing various kind arguments + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + else: + V = VectorFunctionSpace('V', domain, kind='hCuRl') # testing various kind arguments + f1 = x + f2 = y + f = Tuple(f1, f2) + if dim == 3: + x, y, z = domain.coordinates + if i == 0: + V = ScalarFunctionSpace('V', domain, kind='h1') # testing various kind arguments + f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) + elif i == 1: + V = VectorFunctionSpace('V', domain, kind=HcurlSpace) # testing various kind arguments + f1 = z * (z - 1) * x + f2 = z * (z - 1) * y + f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f = Tuple(f1, f2, f3) + else: + V = VectorFunctionSpace('V', domain, kind='Hdiv') # testing various kind arguments + f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) + f3 = z * (z-1) * sin(x*y) + f = Tuple(f1, f2, f3) + + u, v = elements_of(V, names='u, v') + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + Vh = discretize(V, domain_h, degree=degree) + expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (Vh, Vh), backend=backend) + abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) + + I = IdentityOperator(Vh.coeff_space) + DP = DirichletProjector(Vh) + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc + lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f + l = LinearForm(v, integral(domain, lexpr)) + lh = discretize(l, domain_h, Vh, backend=backend) + rhs = lh.assemble() + A = M + 1e30*Mb + A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) + fc = A_inv @ rhs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err = diff.inner(diff) + print(f' | f - P @ f | = {err}') + assert err < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = Vh.coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < tol**2 + + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol + + print() + +#=============================================================================== +@pytest.mark.parametrize('dim', [1, 2, 3]) + +def test_discrete_derham_boundary_projector(dim): + + tol = 1e-15 + + ncells = [8, 8, 8] + degree = [2, 2, 2] + periodic = [False, True, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain_1d = Line ('L', bounds= (0, 1)) + logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] + + mapping_1d = SinMapping1D('LM') + mapping_2d = Annulus ('A' ) + mapping_3d = SquareTorus ('ST') + mappings = [mapping_1d, mapping_2d, mapping_3d] + + rng = np.random.default_rng(42) + + # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs + + f11 = lambda x : np.sin(2*np.pi*x) + + r2 = lambda x, y : np.sqrt(x**2 + y**2) + f21 = lambda x, y : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f22_1 = lambda x, y : x + f22_2 = lambda x, y : y + f22 = (f22_1, f22_2) + + f31 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) * z * (z - 1) + f32_1 = lambda x, y, z : z * (z - 1) * x + f32_2 = lambda x, y, z : z * (z - 1) * y + f32_3 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f32 = (f32_1, f32_2, f32_3) + f33_1 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_2 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) + f33_3 = lambda x, y, z : z * (z - 1) * np.sin(x*y) + f33 = (f33_1, f33_2, f33_3) + + funs = [[f11], [f21, f22], [f31, f32, f33]] + + print() + print(f' ----- Test projectors in dimension {dim} -----') + print() + + domain = mappings[dim-1](logical_domains[dim-1]) + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + # Obtain "true" boundary, i.e., remove periodic y-direction boundary + if dim == 1: + boundary = domain.boundary + elif dim == 2: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) + else: + boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), + domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) + + derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] + degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] + periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] + + domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree_dim) + + d_projectors = derham_h.dirichlet_projectors(kind='linop') + + if dim == 2: + conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) + + nn = NormalVector('nn') + + for i in range(dim): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + if (i == 1) and (dim == 2): + boundary_expr = cross(nn, u) * cross(nn, v) + if (i == 1) and (dim == 3): + boundary_expr = inner(cross(nn, u), cross(nn, v)) + if i == 2: + boundary_expr = inner(nn, u) * inner(nn, v) + + expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DP = d_projectors[i] + + if dim == 2: + CP = conf_projectors[i] + _test_LO_equality_using_rng(DP, CP) + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + f = funs[dim-1][i] + fc = derham_h.projectors()[i](f).coeffs + + # 1. + # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err = diff.inner(diff) + print(f' | f - P @ f |^2 = {err}') + assert err < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + if isinstance(rdm_coeffs.space, BlockVectorSpace): + for block in rdm_coeffs.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) + rdm_coeffs2 = DP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < tol**2 + + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + print(f' || f ||^2 = {l2_norm} should be equal to') + print(f' || P @ f ||^2 = {l2_norm2}') + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol + + print() + +#=============================================================================== +def test_discrete_derham_boundary_projector_multipatch(): + + tol = 1e-15 + + ncells = [8, 8] + degree = [2, 2] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain + domain = build_multipatch_domain(domain_name='annulus_3') + + rng = np.random.default_rng(42) + + # The following are functions satisfying homogeneous Dirichlet BCs + r = lambda x, y : np.sqrt(x**2 + y**2) + f1 = lambda x, y : (r(x, y) - 0.5) * (r(x, y) - 1) + f2_1 = lambda x, y : x + f2_2 = lambda x, y : y + f2 = (f2_1, f2_2) + funs = [f1, f2] + print() + + boundary = domain.boundary + + derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) + + ncells_h = {} + for D in domain.interior: + ncells_h[D.name] = ncells + + domain_h = discretize(domain, ncells=ncells_h, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree) + + projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) + + d_projectors = derham_h.dirichlet_projectors(kind='linop') + + nn = NormalVector('nn') + + for i in range(2): + print(f' - Test DBP{i}') + + u, v = elements_of(derham.spaces[i], names='u, v') + + if i == 0: + boundary_expr = u*v + expr = u*v + if (i == 1): + boundary_expr = cross(nn, u) * cross(nn, v) + expr = inner(u,v) + + a = BilinearForm((u, v), integral(domain, expr)) + ab = BilinearForm((u, v), integral(boundary, boundary_expr)) + + ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) + abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) + + I = IdentityOperator(derham_h.spaces[i].coeff_space) + DP = d_projectors[i] + + M = ah.assemble() + M_0 = DP @ M @ DP + (I - DP) + Mb = abh.assemble() + + f = funs[i] + fc = projectors[i](f).coeffs + + # 1. + # The coefficients of functions satisfying homogeneous Dirichlet + # boundary conditions should not change under application of the corresponding projector + fc2 = DP @ fc + diff = fc - fc2 + err = diff.inner(diff) + print(f' | f - P @ f |^2 = {err}') + assert err < tol**2 + + # 2.1 + # After applying a projector to a random vector, we want to verify that the + # corresponding boundary integral vanishes + rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() + print(' Random boundary integrals:') + for _ in range(3): + for patch in rdm_coeffs.blocks: + + if isinstance(patch.space, BlockVectorSpace): + for block in patch.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=patch._data.shape, dtype="float64", out=patch._data) + + rdm_coeffs2 = DP @ rdm_coeffs + boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 + boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 + print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') + assert boundary_int_proj_rdm < tol**2 + + # 2.2 + # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) + DP_arr = DP.toarray() + rdm_coeffs_arr = rdm_coeffs.toarray() + diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() + err = np.linalg.norm(diff_arr) + assert err < tol**2 + + # 3. + # We want to verify that applying a projector twice does not change the vector twice + fc3 = DP @ fc2 + diff = fc2 - fc3 + err = diff.inner(diff) + print(f' | P @ f - P @ P @ f |^2 = {err}') + assert err < tol**2 + + # 4. + # Finally, the modified mass matrix should still compute inner products correctly + l2_norm = M.dot_inner (fc, fc) + l2_norm2 = M_0.dot_inner(fc, fc) + diff = abs(l2_norm - l2_norm2) + print(f' || f ||^2 = {l2_norm} should be equal to') + print(f' || P @ f ||^2 = {l2_norm2}') + # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. + assert diff < tol + + print() + + +# =============================================================================== +# SCRIPT FUNCTIONALITY +#=============================================================================== + +if __name__ == "__main__": + import sys + pytest.main( sys.argv ) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 9316d37a6..da7305d0d 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -1,19 +1,7 @@ import numpy as np import pytest -from sympy import sin, pi, sqrt, Tuple - -from sympde.calculus import inner, cross -from sympde.expr import integral, LinearForm, BilinearForm -from sympde.topology import elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace -from sympde.topology.datatype import H1Space, HcurlSpace - -from psydac.api.discretization import discretize -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.fem.projectors import DirichletProjector -from psydac.linalg.basic import LinearOperator, IdentityOperator -from psydac.linalg.block import BlockVectorSpace from psydac.linalg.solvers import inverse from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector @@ -64,69 +52,6 @@ def define_data(n, p, matrix_data, dtype=float): xe[s:e + 1] = np.random.random(e + 1 - s) return(V, A, xe) -class SquareTorus(Mapping): - - _expressions = {'x': 'x1 * cos(x2)', - 'y': 'x1 * sin(x2)', - 'z': 'x3'} - - _ldim = 3 - _pdim = 3 - -class Annulus(Mapping): - - _expressions = {'x': 'x1 * cos(x2)', - 'y': 'x1 * sin(x2)'} - - _ldim = 2 - _pdim = 2 - -class SinMapping1D(Mapping): - - _expressions = {'x': 'sin((pi/2)*x1)'} - - _ldim = 1 - _pdim = 1 - -def _test_LO_equality_using_rng(A, B): - """ - A simple tool to check with almost certainty that two linear operators are identical, - by applying them repeatedly to random vectors. - - """ - - assert isinstance(A, LinearOperator) - assert isinstance(B, LinearOperator) - assert A.domain is B.domain - assert A.codomain is B.codomain - - rng = np.random.default_rng(42) - - x = A.domain.zeros() - y1 = A.codomain.zeros() - y2 = y1.copy() - - n = 10 - - for _ in range(n): - - x *= 0. - - if isinstance(A.domain, BlockVectorSpace): - for block in x.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=x._data.shape, dtype="float64", out=x._data) - - A.dot(x, out=y1) - B.dot(x, out=y2) - - diff = y1 - y2 - err = diff.inner(diff) / diff.space.dimension**2 - tol = 1e-15 - - assert err < tol**2 - #=============================================================================== @pytest.mark.parametrize('n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @@ -276,477 +201,6 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): assert errh_norm < tol assert solver == 'pcg' or errc_norm < tol -#=============================================================================== -@pytest.mark.parametrize('dim', [1, 2, 3]) - -def test_function_space_boundary_projector(dim): - - tol = 1e-15 - - ncells_3d = [8, 8, 8] - degree_3d = [2, 2, 2] - periodic_3d = [False, True, False] - - comm = None - backend = PSYDAC_BACKEND_GPYCCEL - - logical_domain_1d = Line ('L', bounds= (0, 1)) - logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) - logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) - logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] - - mapping_1d = SinMapping1D('LM') - mapping_2d = Annulus ('A' ) - mapping_3d = SquareTorus ('ST') - mappings = [mapping_1d, mapping_2d, mapping_3d] - - rng = np.random.default_rng(42) - - print() - print(f' ----- Test projectors in dimension {dim} -----') - print() - - domain = mappings[dim-1](logical_domains[dim-1]) - from sympde.utilities.utils import plot_domain - #plot_domain(domain, draw=True, isolines=True) - - # Obtain "true" boundary, i.e., remove periodic y-direction boundary - if dim == 1: - boundary = domain.boundary - elif dim == 2: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) - else: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), - domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) - - ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim] - degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim] - periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim] - - domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) - - nn = NormalVector('nn') - - for i in range(dim): - print(f' - Test DBP{i}') - - # The function defined here satisfy the corresponding homogeneous Dirichlet BCs - if dim == 1: - x = domain.coordinates - V = ScalarFunctionSpace('V', domain, kind='H1') # testing various kind arguments - f = sin(2*pi*x) - if dim == 2: - x, y = domain.coordinates - if i == 0: - V = ScalarFunctionSpace('V', domain, kind=H1Space) # testing various kind arguments - f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - else: - V = VectorFunctionSpace('V', domain, kind='hCuRl') # testing various kind arguments - f1 = x - f2 = y - f = Tuple(f1, f2) - if dim == 3: - x, y, z = domain.coordinates - if i == 0: - V = ScalarFunctionSpace('V', domain, kind='h1') # testing various kind arguments - f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1) - elif i == 1: - V = VectorFunctionSpace('V', domain, kind=HcurlSpace) # testing various kind arguments - f1 = z * (z - 1) * x - f2 = z * (z - 1) * y - f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f = Tuple(f1, f2, f3) - else: - V = VectorFunctionSpace('V', domain, kind='Hdiv') # testing various kind arguments - f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) - f3 = z * (z-1) * sin(x*y) - f = Tuple(f1, f2, f3) - - u, v = elements_of(V, names='u, v') - if i == 0: - boundary_expr = u*v - if (i == 1) and (dim == 2): - boundary_expr = cross(nn, u) * cross(nn, v) - if (i == 1) and (dim == 3): - boundary_expr = inner(cross(nn, u), cross(nn, v)) - if i == 2: - boundary_expr = inner(nn, u) * inner(nn, v) - - Vh = discretize(V, domain_h, degree=degree) - expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v - - a = BilinearForm((u, v), integral(domain, expr)) - ab = BilinearForm((u, v), integral(boundary, boundary_expr)) - - ah = discretize(a, domain_h, (Vh, Vh), backend=backend) - abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False) - - I = IdentityOperator(Vh.coeff_space) - DP = DirichletProjector(Vh) - - M = ah.assemble() - M_0 = DP @ M @ DP + (I - DP) - Mb = abh.assemble() - - # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc - lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f - l = LinearForm(v, integral(domain, lexpr)) - lh = discretize(l, domain_h, Vh, backend=backend) - rhs = lh.assemble() - A = M + 1e30*Mb - A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10) - fc = A_inv @ rhs - - # 1. - # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet - # boundary conditions should not change under application of the corresponding projector - fc2 = DP @ fc - diff = fc - fc2 - err = diff.inner(diff) - print(f' | f - P @ f | = {err}') - assert err < tol**2 - - # 2.1 - # After applying a projector to a random vector, we want to verify that the - # corresponding boundary integral vanishes - rdm_coeffs = Vh.coeff_space.zeros() - print(' Random boundary integrals:') - for _ in range(3): - if isinstance(rdm_coeffs.space, BlockVectorSpace): - for block in rdm_coeffs.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 - print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < tol**2 - - # 2.2 - # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) - DP_arr = DP.toarray() - rdm_coeffs_arr = rdm_coeffs.toarray() - diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() - err = np.linalg.norm(diff_arr) - assert err < tol**2 - - # 3. - # We want to verify that applying a projector twice does not change the vector twice - fc3 = DP @ fc2 - diff = fc2 - fc3 - err = diff.inner(diff) - print(f' | P @ f - P @ P @ f |^2 = {err}') - assert err < tol**2 - - # 4. - # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = M.dot_inner (fc, fc) - l2_norm2 = M_0.dot_inner(fc, fc) - diff = abs(l2_norm - l2_norm2) - # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. - assert diff < tol - - print() - -#=============================================================================== -@pytest.mark.parametrize('dim', [1, 2, 3]) - -def test_discrete_derham_boundary_projector(dim): - - tol = 1e-15 - - ncells = [8, 8, 8] - degree = [2, 2, 2] - periodic = [False, True, False] - - comm = None - backend = PSYDAC_BACKEND_GPYCCEL - - logical_domain_1d = Line ('L', bounds= (0, 1)) - logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) - logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) - logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d] - - mapping_1d = SinMapping1D('LM') - mapping_2d = Annulus ('A' ) - mapping_3d = SquareTorus ('ST') - mappings = [mapping_1d, mapping_2d, mapping_3d] - - rng = np.random.default_rng(42) - - # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs - - f11 = lambda x : np.sin(2*np.pi*x) - - r2 = lambda x, y : np.sqrt(x**2 + y**2) - f21 = lambda x, y : (r2(x, y) - 0.5) * (r2(x, y) - 1) - f22_1 = lambda x, y : x - f22_2 = lambda x, y : y - f22 = (f22_1, f22_2) - - f31 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) * z * (z - 1) - f32_1 = lambda x, y, z : z * (z - 1) * x - f32_2 = lambda x, y, z : z * (z - 1) * y - f32_3 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) - f32 = (f32_1, f32_2, f32_3) - f33_1 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) - f33_2 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) - f33_3 = lambda x, y, z : z * (z - 1) * np.sin(x*y) - f33 = (f33_1, f33_2, f33_3) - - funs = [[f11], [f21, f22], [f31, f32, f33]] - - print() - print(f' ----- Test projectors in dimension {dim} -----') - print() - - domain = mappings[dim-1](logical_domains[dim-1]) - from sympde.utilities.utils import plot_domain - #plot_domain(domain, draw=True, isolines=True) - - # Obtain "true" boundary, i.e., remove periodic y-direction boundary - if dim == 1: - boundary = domain.boundary - elif dim == 2: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1)) - else: - boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1), - domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1)) - - derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2']) - - ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim] - degree_dim = [degree[0], ] if dim == 1 else degree[0:dim] - periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim] - - domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm) - derham_h = discretize(derham, domain_h, degree=degree_dim) - - d_projectors = derham_h.dirichlet_projectors(kind='linop') - - if dim == 2: - conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True) - - nn = NormalVector('nn') - - for i in range(dim): - print(f' - Test DBP{i}') - - u, v = elements_of(derham.spaces[i], names='u, v') - - if i == 0: - boundary_expr = u*v - if (i == 1) and (dim == 2): - boundary_expr = cross(nn, u) * cross(nn, v) - if (i == 1) and (dim == 3): - boundary_expr = inner(cross(nn, u), cross(nn, v)) - if i == 2: - boundary_expr = inner(nn, u) * inner(nn, v) - - expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v - - a = BilinearForm((u, v), integral(domain, expr)) - ab = BilinearForm((u, v), integral(boundary, boundary_expr)) - - ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) - abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) - - I = IdentityOperator(derham_h.spaces[i].coeff_space) - DP = d_projectors[i] - - if dim == 2: - CP = conf_projectors[i] - _test_LO_equality_using_rng(DP, CP) - - M = ah.assemble() - M_0 = DP @ M @ DP + (I - DP) - Mb = abh.assemble() - - f = funs[dim-1][i] - fc = derham_h.projectors()[i](f).coeffs - - # 1. - # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet - # boundary conditions should not change under application of the corresponding projector - fc2 = DP @ fc - diff = fc - fc2 - err = diff.inner(diff) - print(f' | f - P @ f |^2 = {err}') - assert err < tol**2 - - # 2.1 - # After applying a projector to a random vector, we want to verify that the - # corresponding boundary integral vanishes - rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() - print(' Random boundary integrals:') - for _ in range(3): - if isinstance(rdm_coeffs.space, BlockVectorSpace): - for block in rdm_coeffs.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data) - rdm_coeffs2 = DP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 - print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < tol**2 - - # 2.2 - # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) - DP_arr = DP.toarray() - rdm_coeffs_arr = rdm_coeffs.toarray() - diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() - err = np.linalg.norm(diff_arr) - assert err < tol**2 - - # 3. - # We want to verify that applying a projector twice does not change the vector twice - fc3 = DP @ fc2 - diff = fc2 - fc3 - err = diff.inner(diff) - print(f' | P @ f - P @ P @ f |^2 = {err}') - assert err < tol**2 - - # 4. - # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = M.dot_inner (fc, fc) - l2_norm2 = M_0.dot_inner(fc, fc) - diff = abs(l2_norm - l2_norm2) - print(f' || f ||^2 = {l2_norm} should be equal to') - print(f' || P @ f ||^2 = {l2_norm2}') - # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. - assert diff < tol - - print() - -#=============================================================================== -def test_discrete_derham_boundary_projector_multipatch(): - - tol = 1e-15 - - ncells = [8, 8] - degree = [2, 2] - - comm = None - backend = PSYDAC_BACKEND_GPYCCEL - - from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain - domain = build_multipatch_domain(domain_name='annulus_3') - - rng = np.random.default_rng(42) - - # The following are functions satisfying homogeneous Dirichlet BCs - r = lambda x, y : np.sqrt(x**2 + y**2) - f1 = lambda x, y : (r(x, y) - 0.5) * (r(x, y) - 1) - f2_1 = lambda x, y : x - f2_2 = lambda x, y : y - f2 = (f2_1, f2_2) - funs = [f1, f2] - print() - - boundary = domain.boundary - - derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) - - ncells_h = {} - for D in domain.interior: - ncells_h[D.name] = ncells - - domain_h = discretize(domain, ncells=ncells_h, comm=comm) - derham_h = discretize(derham, domain_h, degree=degree) - - projectors = derham_h.projectors(nquads=[(d + 1) for d in degree]) - - d_projectors = derham_h.dirichlet_projectors(kind='linop') - - nn = NormalVector('nn') - - for i in range(2): - print(f' - Test DBP{i}') - - u, v = elements_of(derham.spaces[i], names='u, v') - - if i == 0: - boundary_expr = u*v - expr = u*v - if (i == 1): - boundary_expr = cross(nn, u) * cross(nn, v) - expr = inner(u,v) - - a = BilinearForm((u, v), integral(domain, expr)) - ab = BilinearForm((u, v), integral(boundary, boundary_expr)) - - ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend) - abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False) - - I = IdentityOperator(derham_h.spaces[i].coeff_space) - DP = d_projectors[i] - - M = ah.assemble() - M_0 = DP @ M @ DP + (I - DP) - Mb = abh.assemble() - - f = funs[i] - fc = projectors[i](f).coeffs - - # 1. - # The coefficients of functions satisfying homogeneous Dirichlet - # boundary conditions should not change under application of the corresponding projector - fc2 = DP @ fc - diff = fc - fc2 - err = diff.inner(diff) - print(f' | f - P @ f |^2 = {err}') - assert err < tol**2 - - # 2.1 - # After applying a projector to a random vector, we want to verify that the - # corresponding boundary integral vanishes - rdm_coeffs = derham_h.spaces[i].coeff_space.zeros() - print(' Random boundary integrals:') - for _ in range(3): - for patch in rdm_coeffs.blocks: - - if isinstance(patch.space, BlockVectorSpace): - for block in patch.blocks: - rng.random(size=block._data.shape, dtype="float64", out=block._data) - else: - rng.random(size=patch._data.shape, dtype="float64", out=patch._data) - - rdm_coeffs2 = DP @ rdm_coeffs - boundary_int_rdm = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2 - boundary_int_proj_rdm = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2 - print(f' rdm: {boundary_int_rdm} proj. rdm: {boundary_int_proj_rdm}') - assert boundary_int_proj_rdm < tol**2 - - # 2.2 - # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray()) - DP_arr = DP.toarray() - rdm_coeffs_arr = rdm_coeffs.toarray() - diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray() - err = np.linalg.norm(diff_arr) - assert err < tol**2 - - # 3. - # We want to verify that applying a projector twice does not change the vector twice - fc3 = DP @ fc2 - diff = fc2 - fc3 - err = diff.inner(diff) - print(f' | P @ f - P @ P @ f |^2 = {err}') - assert err < tol**2 - - # 4. - # Finally, the modified mass matrix should still compute inner products correctly - l2_norm = M.dot_inner (fc, fc) - l2_norm2 = M_0.dot_inner(fc, fc) - diff = abs(l2_norm - l2_norm2) - print(f' || f ||^2 = {l2_norm} should be equal to') - print(f' || P @ f ||^2 = {l2_norm2}') - # This test requires a higher tolerance. M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same up to order 1e-15. - assert diff < tol - - print() # =============================================================================== # SCRIPT FUNCTIONALITY From 9b133829bd4e1b4b54deb9c6f65c75f2ec0a3676 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 7 Nov 2025 11:56:19 +0100 Subject: [PATCH 20/25] Improve DirMultipatch constructor, fix typos, bcs always tuple --- psydac/fem/projectors.py | 42 +++++++++++----------------------------- 1 file changed, 11 insertions(+), 31 deletions(-) diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index 80dbd6b28..3aae76cb9 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -182,7 +182,7 @@ class DirichletProjector(LinearOperator): space_kind : str | SpaceType | None Necessary only if fem_space.kind.name is undefined and no bcs are passed. - Must not be passed if a space_kind argument is passed. + Must not be passed if a bcs argument is passed. Notes ----- @@ -194,10 +194,8 @@ def __init__(self, fem_space, *, bcs=None, space_kind=None): assert isinstance(fem_space, FemSpace) assert bcs is None or isinstance(bcs, Iterable) assert space_kind is None or isinstance(space_kind, (str, SpaceType)) - if bcs is not None: - assert space_kind is None - if space_kind is not None: - assert bcs is None + assert bcs is None or space_kind is None, \ + "Parameters bcs and space_kind are mutually exclusive" coeff_space = fem_space.coeff_space self._domain = coeff_space @@ -205,7 +203,7 @@ def __init__(self, fem_space, *, bcs=None, space_kind=None): if bcs is not None: assert all(isinstance(bc, Boundary) for bc in bcs) - self._bcs = bcs + self._bcs = tuple(bcs) else: self._bcs = self._get_bcs(fem_space, space_kind=space_kind) @@ -266,7 +264,7 @@ def bcs(self): def _get_bcs(self, fem_space, *, space_kind=None): """ - Returns a tuple of Boundries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + Returns a tuple of Boundaries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. Parameters ---------- @@ -368,7 +366,7 @@ class MultipatchDirichletProjector(DirichletProjector): space_kind : str | SpaceType | None Necessary only if fem_space.kind.name is undefined and no bcs are passed. - Must not be passed if a space_kind argument is passed. + Must not be passed if a bcs argument is passed. Notes ----- @@ -376,28 +374,10 @@ class MultipatchDirichletProjector(DirichletProjector): """ def __init__(self, fem_space, *, bcs=None, space_kind=None): - - assert isinstance(fem_space, FemSpace) - assert fem_space.is_multipatch - assert bcs is None or isinstance(bcs, Iterable) - assert space_kind is None or isinstance(space_kind, (str, SpaceType)) - if bcs is not None: - assert space_kind is None - if space_kind is not None: - assert bcs is None - - dim = fem_space.symbolic_space.domain.dim - assert dim == 2, 'The class MultipatchBoundaryProjector is implemented only in 2D.' - - coeff_space = fem_space.coeff_space - self._domain = coeff_space - self._codomain = coeff_space - - if bcs is not None: - assert all(isinstance(bc, Boundary) for bc in bcs) - self._bcs = bcs - else: - self._bcs = self._get_bcs(fem_space, space_kind=space_kind) + super().__init__(fem_space, bcs=bcs, space_kind=space_kind) + if fem_space.ldim != 2: + msg = f'The class {__class__.__name__} is implemented only in 2D.' + raise NotImplementedError(msg) #------------------------------------- # Abstract interface @@ -427,7 +407,7 @@ def dot(self, v, out=None): #-------------------------------------- def _get_bcs(self, fem_space, *, space_kind=None): """ - Returns a tuple of Boundries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. + Returns a tuple of Boundaries that allows to apply homogeneous Dirichlet BCs to functions belonging to fem_space. Parameters ---------- From 91fb8d5c8f215cb110969ae63fa5c238060d3569 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Nov 2025 11:46:17 +0100 Subject: [PATCH 21/25] minor fix in SumLinearOperator toarray, tosparse --- psydac/linalg/basic.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index ba8ca3663..e3996ffce 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -801,15 +801,14 @@ def dtype(self): return None def tosparse(self): - from scipy.sparse import csr_matrix - out = csr_matrix(self.shape, dtype=self.dtype) - for a in self._addends: + out = self.addends[0].tosparse() + for a in self.addends[1:]: out += a.tosparse() return out def toarray(self): - out = np.zeros(self.shape, dtype=self.dtype) - for a in self._addends: + out = self.addends[0].toarray() + for a in self.addends[1:]: out += a.toarray() return out From 34de3a2d40135f9e0cf54539c7381912f9eda4b3 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Nov 2025 11:46:53 +0100 Subject: [PATCH 22/25] allow for tol arg in _test_LO_equality_using_rng --- psydac/fem/tests/test_dirichlet_projectors.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/psydac/fem/tests/test_dirichlet_projectors.py b/psydac/fem/tests/test_dirichlet_projectors.py index 044affb48..1773425d1 100644 --- a/psydac/fem/tests/test_dirichlet_projectors.py +++ b/psydac/fem/tests/test_dirichlet_projectors.py @@ -40,7 +40,7 @@ class SinMapping1D(Mapping): _ldim = 1 _pdim = 1 -def _test_LO_equality_using_rng(A, B): +def _test_LO_equality_using_rng(A, B, tol=1e-15): """ A simple tool to check with almost certainty that two linear operators are identical, by applying them repeatedly to random vectors. @@ -75,7 +75,6 @@ def _test_LO_equality_using_rng(A, B): diff = y1 - y2 err = diff.inner(diff) / diff.space.dimension**2 - tol = 1e-15 assert err < tol**2 From 400454525916fde1e5e107c38d3a9b206cd105c8 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Nov 2025 11:51:44 +0100 Subject: [PATCH 23/25] fix KroneckerStencilMatrix diagonal in parallel --- psydac/linalg/kron.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index 77ddf91ba..83ce504e0 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -260,8 +260,8 @@ def diagonal(self, *, inverse=False, sqrt=False, out=None): # Obtain nested numpy array of diagonal entries (or their inverse (square root)) # by using the `diagonal` method of StencilMatrices diag = 1. - for mat in self.mats[::-1]: - diag = np.array([d*diag for d in mat.diagonal(inverse=inverse, sqrt=sqrt)._data], dtype='float64') + for mat, start, end in zip(self.mats[::-1], V.starts[::-1], V.ends[::-1]): + diag = np.array([d*diag for d in mat.diagonal(inverse=inverse, sqrt=sqrt)._data[start:end+1]], dtype='float64') if out is not None: np.copyto(diag, out._data) From d11d533a376c46fda9d180f90bd64d74f1ba3032 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Nov 2025 11:53:06 +0100 Subject: [PATCH 24/25] improve tests, add parallel test --- .../linalg/tests/test_kron_stencil_matrix.py | 70 +++++++++---------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/psydac/linalg/tests/test_kron_stencil_matrix.py b/psydac/linalg/tests/test_kron_stencil_matrix.py index d4b76d6da..5a3dd3076 100644 --- a/psydac/linalg/tests/test_kron_stencil_matrix.py +++ b/psydac/linalg/tests/test_kron_stencil_matrix.py @@ -1,20 +1,23 @@ -from functools import reduce - -import pytest -import numpy as np -from scipy.sparse import kron - -from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.linalg.stencil import StencilVectorSpace -from psydac.linalg.stencil import StencilVector -from psydac.linalg.stencil import StencilMatrix -from psydac.linalg.kron import KroneckerStencilMatrix - -from sympde.topology import Square, Line, Derham, elements_of -from sympde.expr import BilinearForm, integral -from sympde.calculus import inner -from psydac.linalg.block import BlockLinearOperator -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from functools import reduce + +import pytest +import numpy as np +from mpi4py import MPI +from scipy.sparse import kron + +from sympde.calculus import inner +from sympde.expr import integral, BilinearForm +from sympde.topology import elements_of, Square, Line, Derham + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from psydac.linalg.kron import KroneckerStencilMatrix +from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.stencil import StencilVectorSpace +from psydac.linalg.stencil import StencilVector +from psydac.linalg.stencil import StencilMatrix + #=============================================================================== def compute_global_starts_ends(domain_decomposition, npts): ndims = len(npts) @@ -120,11 +123,9 @@ def test_KroneckerStencilMatrix(dtype, npts, pads, periodic): assert np.array_equal(M_sp.dot(w.toarray()), M.dot(w).toarray()) #============================================================================== -def test_KroneckerStencilMatrix_diagonal(): +def test_KroneckerStencilMatrix_diagonal(comm=None): """We create three mass matrices (Stencil/Block and Kronecker) belonging to a 2D de Rham sequence, and compare their diagonals.""" - from psydac.api.discretization import discretize - ncells = [6, 7] degree = [3, 2] mult = [1, 2] @@ -137,7 +138,7 @@ def test_KroneckerStencilMatrix_diagonal(): domain = Square('S', bounds1=(0,1), bounds2=(0,2)) derham = Derham(domain, sequence=['h1', 'hcurl', 'l2']) - domain_h = discretize(domain, ncells=ncells, periodic=periodic) + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) derham_h = discretize(derham, domain_h, degree=degree, multiplicity=mult) V0, V1, V2 = derham.spaces @@ -162,9 +163,9 @@ def test_KroneckerStencilMatrix_diagonal(): # 2. Obtain KroneckerStencilMatrix / BlockLinearOperator (of KroneckerStencilMatrices) mass matrices - domain_1d_x = Line('L', bounds=(0,1)) - domain_1d_y = Line('L', bounds=(0,2)) - domains_1d = [domain_1d_x, domain_1d_y] + domain_1d_x = Line('L', bounds=domain.bounds1) + domain_1d_y = Line('L', bounds=domain.bounds2) + domains_1d = (domain_1d_x, domain_1d_y) M0s_1d = [] M1s_1d = [] @@ -195,22 +196,21 @@ def test_KroneckerStencilMatrix_diagonal(): [None, KroneckerStencilMatrix(V1cs[1], V1cs[1], M0s_1d[0], M1s_1d[1])]]) M2_kron = KroneckerStencilMatrix(V2cs, V2cs, *M1s_1d) - # 3. Test! - - for M, M_kron in zip([M0, M1, M2], [M0_kron, M1_kron, M2_kron]): - M_diag_arr = [] - M_kron_diag_arr = [] + # 3. Test whether M0/1/2.diagonal() is equal to M0/1/2_kron.diagonal() for all possible kwargs + for M, M_kron in zip((M0, M1, M2), (M0_kron, M1_kron, M2_kron)): options = [True, False] for inverse in options: for sqrt in options: M_diag = M.diagonal(inverse=inverse, sqrt=sqrt) M_kron_diag = M_kron.diagonal(inverse=inverse, sqrt=sqrt) - M_diag_arr.append(M_diag) - M_kron_diag_arr.append(M_kron_diag) - for M_diag, M_kron_diag in zip(M_diag_arr, M_kron_diag_arr): - diff = M_diag.toarray() - M_kron_diag.toarray() - err = np.linalg.norm(diff) - assert err < 1e-10 # arbitrary bound (largest occuring error = 6.7e-13) + from psydac.fem.tests.test_dirichlet_projectors import _test_LO_equality_using_rng + _test_LO_equality_using_rng(M_diag, M_kron_diag, tol=1e-13) + +#============================================================================== +@pytest.mark.parallel +def test_KroneckerStencilMatrix_diagonal_parallel(): + comm = MPI.COMM_WORLD + test_KroneckerStencilMatrix_diagonal(comm=comm) From f934ded3e6e57346c75d977ca02ba10e27d7d865 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Nov 2025 21:28:30 +0100 Subject: [PATCH 25/25] addressing comments --- psydac/fem/tests/test_dirichlet_projectors.py | 6 +++--- psydac/linalg/basic.py | 4 ++-- psydac/linalg/kron.py | 2 +- psydac/linalg/stencil.py | 14 ++++++-------- 4 files changed, 12 insertions(+), 14 deletions(-) diff --git a/psydac/fem/tests/test_dirichlet_projectors.py b/psydac/fem/tests/test_dirichlet_projectors.py index 1773425d1..4e9a9e677 100644 --- a/psydac/fem/tests/test_dirichlet_projectors.py +++ b/psydac/fem/tests/test_dirichlet_projectors.py @@ -73,10 +73,10 @@ def _test_LO_equality_using_rng(A, B, tol=1e-15): A.dot(x, out=y1) B.dot(x, out=y2) - diff = y1 - y2 - err = diff.inner(diff) / diff.space.dimension**2 + diff = y1 - y2 + scaled_err_sqr = diff.inner(diff) / diff.space.dimension**2 - assert err < tol**2 + assert scaled_err_sqr < tol**2 #=============================================================================== @pytest.mark.parametrize('dim', [1, 2]) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index e3996ffce..e721c5ff2 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -801,9 +801,9 @@ def dtype(self): return None def tosparse(self): - out = self.addends[0].tosparse() + out = self.addends[0].tosparse().tocsr() for a in self.addends[1:]: - out += a.tosparse() + out += a.tosparse().tocsr() return out def toarray(self): diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index 83ce504e0..b20645a29 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -261,7 +261,7 @@ def diagonal(self, *, inverse=False, sqrt=False, out=None): # by using the `diagonal` method of StencilMatrices diag = 1. for mat, start, end in zip(self.mats[::-1], V.starts[::-1], V.ends[::-1]): - diag = np.array([d*diag for d in mat.diagonal(inverse=inverse, sqrt=sqrt)._data[start:end+1]], dtype='float64') + diag = np.array([d*diag for d in mat.diagonal(inverse=inverse, sqrt=sqrt)._data[start:end+1]]) if out is not None: np.copyto(diag, out._data) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 3f5ec8d87..ac2dc283b 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -2097,15 +2097,13 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None): # Calculate entries, or set `out=self` in default case if inverse: data = np.divide(1, diag, out=data) - elif out: + if sqrt: + data = np.sqrt(data, out=data) + elif sqrt: + data = np.sqrt(diag, out=data) + elif out not in (None, self): np.copyto(data, diag) - - if sqrt: - if (not inverse) and (out is None): - data = diag.copy() - np.sqrt(data, out=data) - - if (not inverse) and (not sqrt) and (out is None): + else: out = self # If needed create a new StencilDiagonalMatrix object