Skip to content

Commit

Permalink
Expose inter-/histopolation matrices in GlobalProjector, add unit tes…
Browse files Browse the repository at this point in the history 
…ts (#347)

Fixes #346.

The inter-/histopolation matrices are now exposed in
`GlobalProjector.imat_kronecker` as a `LinearOperator` object of type
`KroneckerStencilMatrix` (in the case where both domain and codomain are
scalar spaces) or `BlockLinearOperator` (in the case where at least one
space is vector valued).

The unit tests in `feec/tests/test_commuting_projections.py` have been
enhanced: they now also test how close `imat_kronecker` is to being the
right-inverse of the `solver`attribute (which uses a Kronecker product
of 1D linear solvers).

---------

Co-authored-by: Yaman Güçlü <[email protected]>
Co-authored-by: Max Lindqvist <[email protected]>
Co-authored-by: Martin Campos Pinto <[email protected]>
  • Loading branch information
@spossann @yguclu
@max-models @campospinto
4 people authored Feb 12, 2025
1 parent 66205a5 commit 5669252
Show file tree
Hide file tree
Showing 4 changed files with 228 additions and 18 deletions.
82 changes: 78 additions & 4 deletions psydac/feec/global_projectors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,16 +2,17 @@

import numpy as np

from psydac.linalg.kron import KroneckerLinearSolver
from psydac.linalg.stencil import StencilVector
from psydac.linalg.block import BlockLinearOperator, BlockVector
from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix
from psydac.linalg.stencil import StencilMatrix, StencilVectorSpace
from psydac.linalg.block import BlockLinearOperator
from psydac.core.bsplines import quadrature_grid
from psydac.utilities.quadratures import gauss_legendre
from psydac.fem.basic import FemField

from psydac.fem.tensor import TensorFemSpace
from psydac.fem.vector import VectorFemSpace

from psydac.ddm.cart import DomainDecomposition, CartDecomposition
from psydac.utilities.utils import roll_edges

from abc import ABCMeta, abstractmethod
Expand Down Expand Up @@ -120,18 +121,34 @@ def __init__(self, space, nquads = None):
self._grid_x = []
self._grid_w = []
solverblocks = []
matrixblocks = []
blocks = [] # for BlockLinearOperator
for i,block in enumerate(structure):
# do for each block (i.e. each TensorFemSpace):

block_x = []
block_w = []
solvercells = []
matrixcells = []
blocks += [[]]
for j, cell in enumerate(block):
# for each direction in the tensor space (i.e. each SplineSpace):

V = tensorspaces[i].spaces[j]
s = tensorspaces[i].vector_space.starts[j]
e = tensorspaces[i].vector_space.ends[j]
p = tensorspaces[i].vector_space.pads[j]
n = tensorspaces[i].vector_space.npts[j]
m = tensorspaces[i].multiplicity[j]
periodic = tensorspaces[i].vector_space.periods[j]
ncells = tensorspaces[i].ncells[j]
blocks[-1] += [None] # fill blocks with None, fill the diagonals later

# create a distributed matrix for the current 1d SplineSpace
domain_decomp = DomainDecomposition([ncells], [periodic])
cart_decomp = CartDecomposition(domain_decomp, [n], [[s]], [[e]], [p], [m])
V_cart = StencilVectorSpace(cart_decomp)
M = StencilMatrix(V_cart, V_cart)

if cell == 'I':
# interpolation case
Expand All @@ -143,6 +160,23 @@ def __init__(self, space, nquads = None):
local_x = local_intp_x[:, np.newaxis]
local_w = np.ones_like(local_x)
solvercells += [V._interpolator]

# make 1D collocation matrix in stencil format
row_indices, col_indices = np.nonzero(V.imat)

for row_i, col_i in zip(row_indices, col_indices):

# only consider row indices on process
if row_i in range(V_cart.starts[0], V_cart.ends[0] + 1):
row_i_loc = row_i - s

M._data[row_i_loc + m*p, (col_i + p - row_i)%V.imat.shape[1]] = V.imat[row_i, col_i]

# check if stencil matrix was built correctly
assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1])
# TODO Fix toarray() for multiplicity m > 1
matrixcells += [M.copy()]

elif cell == 'H':
# histopolation case
if quad_x[j] is None:
Expand All @@ -159,11 +193,37 @@ def __init__(self, space, nquads = None):

local_x, local_w = quad_x[j], quad_w[j]
solvercells += [V._histopolator]

# make 1D collocation matrix in stencil format
row_indices, col_indices = np.nonzero(V.hmat)

for row_i, col_i in zip(row_indices, col_indices):

# only consider row indices on process
if row_i in range(V_cart.starts[0], V_cart.ends[0] + 1):
row_i_loc = row_i - s

M._data[row_i_loc + m*p, (col_i + p - row_i)%V.hmat.shape[1]] = V.hmat[row_i, col_i]

# check if stencil matrix was built correctly
assert np.allclose(M.toarray()[s:e + 1], V.hmat[s:e + 1])

matrixcells += [M.copy()]

else:
raise NotImplementedError('Invalid entry in structure array.')

block_x += [local_x]
block_w += [local_w]

# build Kronecker out of single directions
if isinstance(self.space, TensorFemSpace):
matrixblocks += [KroneckerStencilMatrix(self.space.vector_space, self.space.vector_space, *matrixcells)]
else:
matrixblocks += [KroneckerStencilMatrix(self.space.vector_space[i], self.space.vector_space[i], *matrixcells)]

# fill the diagonals for BlockLinearOperator
blocks[i][i] = matrixblocks[-1]

# finish block, build solvers, get dataslice to project to
self._grid_x += [block_x]
Expand All @@ -173,6 +233,13 @@ def __init__(self, space, nquads = None):

dataslice = tuple(slice(p*m, -p*m) for p, m in zip(tensorspaces[i].vector_space.pads,tensorspaces[i].vector_space.shifts))
dofs[i] = rhsblocks[i]._data[dataslice]

# build final Inter-/Histopolation matrix (distributed)
if isinstance(self.space, TensorFemSpace):
self._imat_kronecker = matrixblocks[0]
else:
self._imat_kronecker = BlockLinearOperator(self.space.vector_space, self.space.vector_space,
blocks=blocks)

# finish arguments and create a lambda
args = (*intp_x, *quad_x, *quad_w, *dofs)
Expand Down Expand Up @@ -238,6 +305,13 @@ def solver(self):
"""
return self._solver

@property
def imat_kronecker(self):
"""
Inter-/Histopolation matrix in distributed format.
"""
return self._imat_kronecker

@abstractmethod
def _structure(self, dim):
"""
Expand Down
154 changes: 145 additions & 9 deletions psydac/feec/tests/test_commuting_projections.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,19 +3,20 @@
from psydac.feec.global_projectors import Projector_H1, Projector_L2, Projector_Hcurl, Projector_Hdiv
from psydac.fem.tensor import TensorFemSpace, SplineSpace
from psydac.fem.vector import VectorFemSpace
from psydac.linalg.block import BlockVector
from psydac.core.bsplines import make_knots
from psydac.feec.derivatives import Derivative_1D, Gradient_2D, Gradient_3D
from psydac.feec.derivatives import ScalarCurl_2D, VectorCurl_2D, Curl_3D
from psydac.feec.derivatives import Divergence_2D, Divergence_3D
from psydac.ddm.cart import DomainDecomposition
from psydac.linalg.solvers import inverse
from psydac.linalg.basic import IdentityOperator

from mpi4py import MPI
import numpy as np
import pytest

#==============================================================================
# Run test
# 3D tests
#==============================================================================
@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [5])
Expand Down Expand Up @@ -76,7 +77,21 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#==============================================================================
#--------------------------
# check BlockLinearOperator
#--------------------------
Id_0 = IdentityOperator(H1.vector_space)
Err_0 = P0.solver @ P0.imat_kronecker - Id_0
e0 = Err_0 @ u0.coeffs # random vector could be used as well
norm2_e0 = np.sqrt(e0.dot(e0))
assert norm2_e0 < 1e-12

Id_1 = IdentityOperator(Hcurl.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs # random vector could be used as well
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e1 < 1e-12

@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [8])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -153,7 +168,22 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#==============================================================================
#--------------------------
# check BlockLinearOperator
#--------------------------

Id_1 = IdentityOperator(Hcurl.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs # random vector could be used as well
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e1 < 1e-12

Id_2 = IdentityOperator(Hdiv.vector_space)
Err_2 = P2.solver @ P2.imat_kronecker - Id_2
e2 = Err_2 @ u2.coeffs # random vector could be used as well
norm2_e2 = np.sqrt(e2.dot(e2))
assert norm2_e2 < 1e-12

@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [8])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -221,6 +251,26 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_2 = IdentityOperator(Hdiv.vector_space)
Err_2 = P2.solver @ P2.imat_kronecker - Id_2
e2 = Err_2 @ u2.coeffs # random vector could be used as well
norm2_e2 = np.sqrt(e2.dot(e2))
assert norm2_e2 < 1e-12

Id_3 = IdentityOperator(L2.vector_space)
Err_3 = P3.solver @ P3.imat_kronecker - Id_3
e3 = Err_3 @ u3.coeffs # random vector could be used as well
norm2_e3 = np.sqrt(e3.dot(e3))
assert norm2_e3 < 1e-12

#==============================================================================
# 2D tests
#==============================================================================
@pytest.mark.parallel
@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [5])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -277,6 +327,23 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_0 = IdentityOperator(H1.vector_space)
Err_0 = P0.solver @ P0.imat_kronecker - Id_0
e0 = Err_0 @ u0.coeffs # random vector could be used as well
norm2_e0 = np.sqrt(e0.dot(e0))
assert norm2_e0 < 1e-12

Id_1 = IdentityOperator(Hcurl.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs # random vector could be used as well
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e1 < 1e-12

@pytest.mark.parallel
@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [5])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -333,6 +400,23 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_0 = IdentityOperator(H1.vector_space)
Err_0 = P0.solver @ P0.imat_kronecker - Id_0
e0 = Err_0 @ u0.coeffs # random vector could be used as well
norm2_e0 = np.sqrt(e0.dot(e0))
assert norm2_e0 < 1e-12

Id_1 = IdentityOperator(Hdiv.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs # random vector could be used as well
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e0 < 1e-12

@pytest.mark.parallel
@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [8])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -396,6 +480,23 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc, m):
error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_2 = IdentityOperator(Hdiv.vector_space)
Err_2 = P2.solver @ P2.imat_kronecker - Id_2
e2 = Err_2 @ u2.coeffs
norm2_e2 = np.sqrt(e2.dot(e2))
assert norm2_e2 < 1e-12

Id_3 = IdentityOperator(L2.vector_space)
Err_3 = P3.solver @ P3.imat_kronecker - Id_3
e3 = Err_3 @ u3.coeffs
norm2_e3 = np.sqrt(e3.dot(e3))
assert norm2_e3 < 1e-12

@pytest.mark.parallel
@pytest.mark.parametrize('Nel', [8, 12])
@pytest.mark.parametrize('Nq', [8])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -451,14 +552,33 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc, m):
# Projections and discrete derivatives
#-------------------------------------

u2 = P1((fun1, fun2))
u3 = P2(difun)
Dfun_h = curl(u2)
Dfun_proj = u3
u1 = P1((fun1, fun2))
u2 = P2(difun)
Dfun_h = curl(u1)
Dfun_proj = u2

error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_1 = IdentityOperator(Hcurl.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e1 < 1e-12

Id_2 = IdentityOperator(L2.vector_space)
Err_2 = P2.solver @ P2.imat_kronecker - Id_2
e2 = Err_2 @ u2.coeffs
norm2_e2 = np.sqrt(e2.dot(e2))
assert norm2_e2 < 1e-12

#==============================================================================
# 1D tests
#==============================================================================
@pytest.mark.parametrize('Nel', [16, 20])
@pytest.mark.parametrize('Nq', [5])
@pytest.mark.parametrize('p', [2,3])
Expand Down Expand Up @@ -509,6 +629,23 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc, m):

error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
assert error < 1e-9

#--------------------------
# check BlockLinearOperator
#--------------------------

Id_0 = IdentityOperator(H1.vector_space)
Err_0 = P0.solver @ P0.imat_kronecker - Id_0
e0 = Err_0 @ u0.coeffs
norm2_e0 = np.sqrt(e0.dot(e0))
assert norm2_e0 < 1e-12

Id_1 = IdentityOperator(L2.vector_space)
Err_1 = P1.solver @ P1.imat_kronecker - Id_1
e1 = Err_1 @ u1.coeffs
norm2_e1 = np.sqrt(e1.dot(e1))
assert norm2_e1 < 1e-12

#==============================================================================
if __name__ == '__main__':

Expand All @@ -526,4 +663,3 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc, m):
test_2d_commuting_pro_3(Nel, Nq, p, bc, m)
test_2d_commuting_pro_4(Nel, Nq, p, bc, m)
test_1d_commuting_pro_1(Nel, Nq, p, bc, m)

2 changes: 1 addition & 1 deletion psydac/fem/splines.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ def init_interpolation( self, dtype=float ):
for i,j in zip( *cmat.nonzero() ):
bmat[u+l+i-j,j] = cmat[i,j]
self._interpolator = BandedSolver( u, l, bmat )
self.imat = imat
self.imat = imat

# Store flag
self._interpolation_ready = True
Expand Down
Loading

0 comments on commit 5669252

Please sign in to comment.