# Copyright 2026 Capytaine developers
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Definition of the methods to build influence matrices, using possibly some sparse structures."""
import logging
from abc import ABC, abstractmethod
from typing import Tuple, Union, Optional, Callable
import numpy as np
import scipy.sparse.linalg as ssl
from capytaine.meshes.symmetric_meshes import ReflectionSymmetricMesh, RotationSymmetricMesh
from capytaine.green_functions.abstract_green_function import AbstractGreenFunction, GreenFunctionEvaluationError
from capytaine.green_functions.delhommeau import Delhommeau
from capytaine.tools.lazy_matrices import LazyMatrix
from capytaine.tools.lists_of_points import _normalize_points
from capytaine.tools.block_circulant_matrices import (
BlockCirculantMatrix,
NestedBlockCirculantMatrix,
lu_decompose,
has_been_lu_decomposed,
MatrixLike,
LUDecomposedMatrixLike
)
LOG = logging.getLogger(__name__)
####################
# ABSTRACT CLASS #
####################
[docs]
class MatrixEngine(ABC):
"""Abstract method to build a matrix."""
[docs]
@abstractmethod
def build_matrices(self, mesh1, mesh2, free_surface, water_depth, wavenumber, adjoint_double_layer):
pass
[docs]
@abstractmethod
def build_S_matrix(self, mesh1, mesh2, free_surface, water_depth, wavenumber):
pass
[docs]
@abstractmethod
def build_fullK_matrix(self, mesh1, mesh2, free_surface, water_depth, wavenumber):
pass
##################
# BASIC ENGINE #
##################
[docs]
class Counter:
def __init__(self):
self.nb_iter = 0
def __call__(self, *args, **kwargs):
self.nb_iter += 1
[docs]
def solve_gmres(A, b):
LOG.debug(f"Solve with GMRES for {A}.")
if LOG.isEnabledFor(logging.INFO):
counter = Counter()
x, info = ssl.gmres(A, b, atol=1e-6, callback=counter)
LOG.info(f"End of GMRES after {counter.nb_iter} iterations.")
else:
x, info = ssl.gmres(A, b, atol=1e-6)
if info > 0:
raise RuntimeError(f"No convergence of the GMRES after {info} iterations.\n"
"This can be due to overlapping panels or irregular frequencies.")
return x
LUDecomposedMatrixOrNot = Union[MatrixLike, LUDecomposedMatrixLike]
[docs]
class DefaultMatrixEngine(MatrixEngine):
"""
Features:
- Caching of the last computed matrices.
- Supports plane symmetries and nested plane symmetries.
- Linear solver can be customized. Default is `lu_decomposition` with caching of the LU decomposition.
Parameters
----------
green_function: AbstractGreenFunction
the low level implementation used to compute the coefficients of the matrices.
linear_solver: str or function, optional
Setting of the numerical solver for linear problems Ax = b.
It can be set with the name of a preexisting solver
(available: "lu_decomposition", "lu_decompositon_with_overwrite" and "gmres", the former is the default choice)
or by passing directly a solver function.
"""
green_function: AbstractGreenFunction
_linear_solver: Union[str, Callable]
last_computed_matrices: Optional[Tuple[MatrixLike, LUDecomposedMatrixOrNot]]
def __init__(self, *, green_function=None, linear_solver='lu_decomposition'):
self.green_function = Delhommeau() if green_function is None else green_function
self._linear_solver = linear_solver
self.last_computed_inputs = None
self.last_computed_matrices = None
self.exportable_settings = {
'engine': 'DefaultMatrixEngine',
'linear_solver': str(linear_solver),
**self.green_function.exportable_settings,
}
def __str__(self):
params= [f"green_function={self.green_function}", f"linear_solver={repr(self._linear_solver)}"]
return f"DefaultMatrixEngine({', '.join(params)})"
def __repr__(self):
return self.__str__()
def _repr_pretty_(self, p, cycle):
p.text(self.__str__())
[docs]
def build_S_matrix(self, points, mesh, **gf_params) -> np.ndarray:
"""Similar to :code:`build_matrices`, but returning only :math:`S`
points: np.ndarray of shape (nb_collocations_points, 3)
List of collocation points. Usually went through `_normalize_points`.
"""
full_mesh = mesh.merged()
def build_S_rows(slice_):
S, _ = self.green_function.evaluate(points[slice_, :], full_mesh, **gf_params)
check_if_nan_in_matrix([S])
return S
nb_points = points.shape[0]
if nb_points < 500:
return build_S_rows(slice(0, nb_points)) # Just the full matrix
else:
return LazyMatrix(
build_S_rows,
shape=(points.shape[0], mesh.nb_faces),
dtype=complex,
chunk_size=100
)
[docs]
def build_fullK_matrix(self, mesh1, mesh2, **gf_params) -> np.ndarray:
"""Similar to :code:`build_matrices`, but returning only full :math:`K`
(that is the three components of the gradient, not just the normal one)"""
# TODO: could use symmetries. In particular for forward, we compute the
# full velocity on the same mesh so symmetries could be used.
gf_params.setdefault("diagonal_term_in_double_layer", True) # Unclear if this is a good default
gf_params.setdefault("adjoint_double_layer", True)
gf_params.setdefault("early_dot_product", False)
_, fullK = self.green_function.evaluate(mesh1, mesh2.merged(), **gf_params)
check_if_nan_in_matrix([fullK])
return fullK
def _build_matrices_with_symmetries(self, mesh1, mesh2, *, diagonal_term_in_double_layer=True, **gf_params) -> Tuple[MatrixLike, MatrixLike]:
gf_params.setdefault("early_dot_product", True)
if (isinstance(mesh1, ReflectionSymmetricMesh)
and isinstance(mesh2, ReflectionSymmetricMesh)
and mesh1.plane == mesh2.plane):
if (isinstance(mesh1.half, ReflectionSymmetricMesh)
and isinstance(mesh2.half, ReflectionSymmetricMesh)
and mesh1.half.plane == mesh2.half.plane):
# Nested plane symmetries
S, K = self.green_function.evaluate(
mesh1.merged(), mesh2.half.half,
diagonal_term_in_double_layer=diagonal_term_in_double_layer,
**gf_params
)
check_if_nan_in_matrix([S, K])
block_shape = (mesh2.half.half.nb_faces, mesh2.half.half.nb_faces)
return NestedBlockCirculantMatrix(S.reshape((4, *block_shape))), NestedBlockCirculantMatrix(K.reshape((4, *block_shape)))
# elif (isinstance(mesh1.half, RotationSymmetricMesh)
# and isinstance(mesh2.half, RotationSymmetricMesh)
# and mesh1.half.n == mesh2.half.n):
else:
# Single plane symmetries
S, K = self.green_function.evaluate(
mesh1.merged(), mesh2.half,
diagonal_term_in_double_layer=diagonal_term_in_double_layer,
**gf_params
)
check_if_nan_in_matrix([S, K])
block_shape = (mesh2.half.nb_faces, mesh2.half.nb_faces)
return BlockCirculantMatrix(S.reshape((2, *block_shape))), BlockCirculantMatrix(K.reshape((2, *block_shape)))
elif (isinstance(mesh1, RotationSymmetricMesh)
and isinstance(mesh2, RotationSymmetricMesh)
and mesh1.n == mesh2.n):
if (isinstance(mesh1.wedge, ReflectionSymmetricMesh)
and isinstance(mesh2.wedge, ReflectionSymmetricMesh)
and mesh1.wedge.plane == mesh2.wedge.plane):
# Dihedral symmetry
S, K = self.green_function.evaluate(
mesh1.merged(), mesh2.wedge.half,
diagonal_term_in_double_layer=diagonal_term_in_double_layer,
**gf_params
)
check_if_nan_in_matrix([S, K])
n_blocks = 2*mesh1.n # == mesh2.n
block_shape = (mesh2.wedge.half.nb_faces, mesh2.wedge.half.nb_faces)
return NestedBlockCirculantMatrix(S.reshape((n_blocks, *block_shape))), NestedBlockCirculantMatrix(K.reshape((n_blocks, *block_shape)))
else:
# Rotation symmetry
S_cols, K_cols = self.green_function.evaluate(
mesh1.merged(), mesh2.wedge,
diagonal_term_in_double_layer=diagonal_term_in_double_layer,
**gf_params,
)
# Building the first column of blocks, that is the interactions of all of mesh1 with the reference wedge of mesh2.
check_if_nan_in_matrix([S_cols, K_cols])
n_blocks = mesh1.n # == mesh2.n
block_shape = (mesh2.wedge.nb_faces, mesh2.wedge.nb_faces)
return (
BlockCirculantMatrix(S_cols.reshape((n_blocks, *block_shape))),
BlockCirculantMatrix(K_cols.reshape((n_blocks, *block_shape))),
)
else:
S, K = self.green_function.evaluate(
mesh1, mesh2,
diagonal_term_in_double_layer=diagonal_term_in_double_layer,
**gf_params
)
check_if_nan_in_matrix([S, K])
return S, K
def _build_and_cache_matrices_with_symmetries(
self, mesh1, mesh2, **gf_params
) -> Tuple[MatrixLike, LUDecomposedMatrixOrNot]:
if (mesh1, mesh2, gf_params) == self.last_computed_inputs:
LOG.debug("%s: reading cache.", self.__class__.__name__)
return self.last_computed_matrices
else:
LOG.debug("%s: computing new matrices.", self.__class__.__name__)
self.last_computed_matrices = None # Unlink former cached values, so the memory can be freed to compute new matrices.
S, K = self._build_matrices_with_symmetries(mesh1, mesh2, **gf_params)
self.last_computed_inputs = (mesh1, mesh2, gf_params)
self.last_computed_matrices = (S, K)
return self.last_computed_matrices
# Main interface for compliance with AbstractGreenFunction interface
[docs]
def build_matrices(self, mesh1, mesh2, **gf_params):
r"""Build the influence matrices between mesh1 and mesh2.
Parameters
----------
mesh1: MeshLike or list of points
mesh of the receiving body (where the potential is measured)
mesh2: MeshLike
mesh of the source body (over which the source distribution is integrated)
free_surface: float
position of the free surface (default: :math:`z = 0`)
water_depth: float
position of the sea bottom (default: :math:`z = -\infty`)
wavenumber: float
wavenumber (default: 1.0)
adjoint_double_layer: bool, optional
compute double layer for direct method (F) or adjoint double layer for indirect method (T) matrices (default: True)
Returns
-------
tuple of matrix-like (Numpy arrays or BlockCirculantMatrix)
the matrices :math:`S` and :math:`K`
"""
return self._build_and_cache_matrices_with_symmetries(
mesh1, mesh2, **gf_params
)
[docs]
def linear_solver(self, A: LUDecomposedMatrixOrNot, b: np.ndarray) -> np.ndarray:
"""Solve a linear system with left-hand side A and right-hand-side b
Parameters
----------
A: matrix-like
Expected to be the second output of `build_matrices`
b: np.ndarray
Vector of the correct length
Returns
-------
x: np.ndarray
Vector such that A@x = b
"""
if not isinstance(self._linear_solver, str):
# If not a string, it is expected to be a custom function that can
# be called to solve the system
x = self._linear_solver(A, b)
if not x.shape == b.shape:
raise ValueError(f"Error in linear solver of {self}: the shape of the output ({x.shape}) "
f"does not match the expected shape ({b.shape})")
return x
elif self._linear_solver in ("lu_decomposition", "lu_decomposition_with_overwrite") :
overwrite_a = (self._linear_solver == "lu_decomposition_with_overwrite")
if not has_been_lu_decomposed(A):
luA = lu_decompose(A, overwrite_a=overwrite_a)
if A is self.last_computed_matrices[1]:
# In normal operation of Capytaine, `A` is always the $D$
# or $K$ matrix stored in the cache of the solver.
# Here we replace the matrix by its LU decomposition in the
# cache to avoid doing the decomposition again.
self.last_computed_matrices = (self.last_computed_matrices[0], luA)
else:
luA: LUDecomposedMatrixLike = A
return luA.solve(b)
elif self._linear_solver == "gmres":
return solve_gmres(A, b)
else:
raise NotImplementedError(
f"Unknown `linear_solver` in DefaultMatrixEngine: {self._linear_solver}"
)
[docs]
def compute_ram_estimation(self, problem):
nb_faces = problem.body.mesh.nb_faces
nb_matrices = 2
nb_bytes = 16
if self._linear_solver == "lu_decomposition":
nb_matrices += 1
if self.green_function.floating_point_precision == "float32":
nb_bytes = 8
# In theory a simple symmetry is a gain of factor 1/2
# and a nested symmetry is a gain of factor 1/4.
# For the solvers that use LU decomposition the gain is a bit less.
solver_factors = {
# Formula to compute the factor of gain:
# (2 matrices * theoretical symmetry factor + LU decomposition + intermediate_step) / nb matrices without symmetry
"lu_decomposition": {
"simple": 2 / 3, # (2 * 1/2 + 1/2 + 1/2) / 3
"nested": 5 / 12, # (2 * 1/4 + 1/4 + 1/2) / 3
"rotation": 4 / 3,
},
# Formula to compute the factor of gain:
# (2 matrices * theoretical symmetry factor + intermediate step) / nb matrices without symmetry
"lu_decomposition_with_overwrite": {
"simple": 3 / 4, # (2 * 1/2 + 1/2) / 2
"nested": 1 / 2, # (2 * 1/4 + 1/2) / 2
"rotation": 3 / 2,
},
"gmres": {
"simple": 1 / 2,
"nested": 1 / 4,
"rotation": 1,
},
}
if isinstance(problem.body.mesh, ReflectionSymmetricMesh):
if isinstance(problem.body.mesh.half, ReflectionSymmetricMesh):
# Should not go deeper than that, there is currently only two
# symmetries available
symmetry_type = "nested"
else:
symmetry_type = "simple"
symmetry_factor = solver_factors[self._linear_solver][symmetry_type]
elif isinstance(problem.body.mesh, RotationSymmetricMesh):
symmetry_factor = solver_factors[self._linear_solver]["rotation"] / problem.body.mesh.n
else:
symmetry_factor = 1.0
memory_peak = symmetry_factor * nb_faces**2 * nb_matrices * nb_bytes/1e9
return memory_peak
[docs]
def check_if_nan_in_matrix(matrices):
for matrix in matrices:
if np.any(np.isnan(matrix)):
raise GreenFunctionEvaluationError(
"Green function returned a NaN in the interaction matrix.\n"
"It could be due to overlapping panels.")