Advanced cookbook

This page contains several examples of Capytaine’s advanced features, about the fine control of the internal of the solver. The scripts can be downloaded individually as Python files from the Github folder with examples for version 3.0.dev.

C1. Parallel computation with Joblib

TODO

C2. Direct or indirect boundary integral equation

TODO

C3. Single precision computations

TODO

C4. Using an iterative solver

TODO

C5. Plot the influence matrix

This example plots the influence matrix for an horizontal cylinder.

import numpy as np
import capytaine as cpt

# Generate the mesh of a cylinder
cylinder = cpt.FloatingBody(
        mesh=cpt.mesh_horizontal_cylinder(
            length=10.0, radius=1.0,
            center=(0, 0, -2),
            resolution=(1, 6, 8)
            ))

engine = cpt.BasicMatrixEngine()
green_function = cpt.Delhommeau()

S, K = engine.build_matrices(
    cylinder.mesh, cylinder.mesh,
    free_surface=0.0, water_depth=np.inf,
    wavenumber=1.0,
    green_function=green_function,
)

# Plot the absolute value of the matrix S
import matplotlib.pyplot as plt
plt.imshow(abs(np.array(S)))
plt.colorbar()
plt.title("$|S|$")
plt.tight_layout()
plt.show()

C6. Toeplitz matrix of an axisymmetric buoy

import numpy as np
import capytaine as cpt
import capytaine.io.xarray

cpt.set_logging('INFO')

# Profile of the axisymmetric body
def shape(z):
    return 0.1*(-(z+1)**2 + 16)

# Generate the mesh and display it with VTK.
buoy = cpt.FloatingBody(
    mesh=cpt.AxialSymmetricMesh.from_profile(shape, z_range=np.linspace(-5, 0, 30), nphi=40)
)
buoy.add_translation_dof(name="Heave")
buoy.show()

# Set up problems
omega_range = np.linspace(0.1, 5.0, 60)
problems = [cpt.RadiationProblem(body=buoy, radiating_dof='Heave', omega=omega)
            for omega in omega_range]

# Solve the problems using the axial symmetry
solver = cpt.BEMSolver(engine=cpt.HierarchicalToeplitzMatrixEngine())
results = solver.solve_all(problems)
dataset = capytaine.io.xarray.assemble_dataset(results)

# Plot results
import matplotlib.pyplot as plt
plt.figure()
plt.plot(
    omega_range,
    dataset['added_mass'].sel(radiating_dof='Heave',
                              influenced_dof='Heave'),
    label="Added mass",
)
plt.plot(
    omega_range,
    dataset['radiation_damping'].sel(radiating_dof='Heave',
                                     influenced_dof='Heave'),
    label="Radiation damping",
)
plt.xlabel('omega')
plt.grid()
plt.legend()
plt.show()

C7. Hierarchical matrices with preconditionning

import numpy as np
import capytaine as cpt
from time import time

radius = 2.
spacing = 5.
nb_cols = 14
nb_rows = 14

omega = 1.

x, y = np.meshgrid(np.arange(nb_cols)*spacing, np.arange(nb_rows)*spacing)
x = x.flatten()
y = y.flatten()

bodies = [cpt.FloatingBody(mesh=cpt.mesh_sphere(radius=radius,
                                                center=(x[ii], y[ii], 0),
                                                resolution=(10,10)))
            for ii in range(len(x))]

for body in bodies:
    body.keep_immersed_part()
    body.add_translation_dof(name='Heave')

# Dense problem setup and solution
print('Solving dense problem')
all_bodies = cpt.FloatingBody.join_bodies(*bodies, name="joined bodies")
problems = [cpt.RadiationProblem(body=all_bodies, radiating_dof=dof, omega=omega)
        for dof in all_bodies.dofs]
solver = cpt.BEMSolver()
dresults = solver.solve_all(problems)
ddata = cpt.assemble_dataset(dresults)


# Hierarchical problem setup
all_bodies = cpt.FloatingBody.cluster_bodies(*bodies, name="clustered bodies")
problems = [cpt.RadiationProblem(body=all_bodies, radiating_dof=dof, omega=omega)
            for dof in all_bodies.dofs]


cpt.set_logging(level='INFO')  # prints the number of iterations
# Hierarchical solve, with preconditioner
print('Solving hierarchical problem with preconditioner')
t0 = time()
sparse_engine = cpt.HierarchicalPrecondMatrixEngine(ACA_distance=2.0, ACA_tol=1e-2)
solver = cpt.BEMSolver(engine=sparse_engine)
presults = solver.solve_all(problems)
tP = time() - t0

data = cpt.assemble_dataset(presults)

# Hierarchical solve, no preconditioner
print('Solving hierarchical problem without preconditioner')
t0 = time()
sparse_engine = cpt.HierarchicalToeplitzMatrixEngine(ACA_distance=2.0, ACA_tol=1e-2)
solver = cpt.BEMSolver(engine=sparse_engine)
hresults = solver.solve_all(problems)
tNP = time() - t0

# The clustering process changes the order of the bodies with respect to
# the one in list 'bodies'. The correspondence with keys is maintained correct.
# However, if one wants to go back to the initial ordering, the following
# lines can be used.

ordered_dofs_names = []
for body in bodies:
    ordered_dofs_names += [body.name + '__' + dofkey for dofkey in list(body.dofs.keys())  ]
reordered_dofs = {'radiating_dof': ordered_dofs_names,
                  'influenced_dof': ordered_dofs_names}

reordered_data = data.reindex(reordered_dofs)

print('error = ', np.linalg.norm(presults[0].sources - hresults[0].sources))
print('time with preconditioner = ', tP)
print('time without preconditioner = ', tNP)

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 3)
ax[0].matshow(ddata['added_mass'].data[0,:,:])
ax[1].matshow(reordered_data['added_mass'].data[0,:,:])
ax[2].matshow(data['added_mass'].data[0,:,:])
ax[0].set_title('Dense')
ax[1].set_title('HP, reordered')
ax[2].set_title('HP, before reordering')
fig.suptitle('Added mass matrix')

plt.show()

C8. Compare two implementations of the Green function

This is an example of comparison of two implementations of the Green function.

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import capytaine as cpt

cpt.set_logging('INFO')

# Generate body
mesh = cpt.mesh_horizontal_cylinder(
    length=3.0, radius=1.0, center=(0, 0, -1.01), resolution=(5, 30, 15)
)
body = cpt.FloatingBody(mesh)
body.add_translation_dof(name="Heave")

test_matrix = xr.Dataset(coords={
    'omega': np.linspace(0.5, 4, 20),
    'radiating_dof': list(body.dofs.keys()),
    'water_depth': [np.inf, 4.0]
})

green_functions = [
        cpt.Delhommeau(gf_singularities="low_freq"),
        # cpt.Delhommeau(gf_singularities="high_freq"),  # For this problem, more difficult to converge
        cpt.HAMS_GF(),
        ]

data = []
for gf in green_functions:
    data.append(cpt.BEMSolver(green_function=gf).fill_dataset(test_matrix, body))

fig, axs = plt.subplots(2, 1, sharex=True, layout="constrained")
for gf, ds in zip(green_functions, data):
    ds['added_mass'].sel(water_depth=np.inf).plot(ax=axs[0], x='omega', linestyle="--", label=f"Infinite depth {gf}")
    ds['added_mass'].sel(water_depth=4.0).plot(ax=axs[0], x='omega', label=f"Finite depth {gf}")
    ds['radiation_damping'].sel(water_depth=np.inf).plot(ax=axs[1], x='omega', linestyle="--", label=f"Infinite depth {gf}")
    ds['radiation_damping'].sel(water_depth=4.0).plot(ax=axs[1], x='omega', label=f"Finite depth {gf}")
axs[0].set_title("Added mass")
axs[0].legend()
axs[1].set_title("Radiation damping")

plt.show()

C9. Use a custom Green function

This is an example of how to implement a custom Green function.

import numpy as np
from numpy import dot, pi
from numpy.linalg import norm

import capytaine as cpt
from capytaine.green_functions.abstract_green_function import AbstractGreenFunction


class MyGreenFunction(AbstractGreenFunction):
    """An example of a custom routine to evaluate the Green function."""

    floating_point_precision = "float64"  # Optional, could allow for RAM usage optimisations

    def evaluate(
        self,
        mesh1,
        mesh2,
        free_surface,
        water_depth,
        wavenumber,
        adjoint_double_layer=True,
        early_dot_product=True,
    ):
        """The main method that needs to be implemented in the class."""

        if not adjoint_double_layer:
            raise NotImplementedError(
                """MyGreenFunction only implements the adjoint double layer"""
            )

        if free_surface != np.inf and water_depth != np.inf:
            raise NotImplementedError(
                """MyGreenFunction only implements the case without a free surface"""
            )

        # Initialize the matrices
        S = np.zeros((mesh1.nb_faces, mesh2.nb_faces))
        if early_dot_product:
            gradient_shape = (mesh1.nb_faces, mesh2.nb_faces)
        else:
            gradient_shape = (mesh1.nb_faces, mesh2.nb_faces, 3)
        K = np.zeros(gradient_shape)

        for i in range(mesh1.nb_faces):
            for j in range(mesh2.nb_faces):
                area = mesh1.faces_areas[j]
                if i == j:
                    # Approximation of the integral of the 1/r Green
                    # function over the panel by taking the integral over a
                    # circle of same center and same area.
                    S[i, i] = -1 / (4 * pi) * 2 * np.sqrt(area * pi)
                    if mesh1 is mesh2:
                        if early_dot_product:
                            K[i, i] = 1 / 2
                        else:
                            K[i, i, :] = 1 / 2 * mesh1.faces_normals[i, :]
                    else:
                        if early_dot_product:
                            K[i, i] = 0.0
                        else:
                            K[i, i, :] = 0.0

                else:
                    # The integral on the face is computed using the value
                    # at the center of the face.
                    r = mesh1.faces_centers[i, :] - mesh2.faces_centers[j, :]

                    S[i, j] = -area / (4 * pi * norm(r))
                    gradient_rankine = area * r / (4 * pi * norm(r) ** 3)
                    if early_dot_product:
                        K[i, j] = dot(gradient_rankine, mesh1.faces_normals[i, :])
                    else:
                        K[i, j, :] = gradient_rankine

        return S, K


# Define a solver using the Green function above.
custom_solver = cpt.BEMSolver(green_function=MyGreenFunction())
default_solver = cpt.BEMSolver()  # for comparison

# Example of a problem
radius = 1.4
body = cpt.FloatingBody(
    mesh=cpt.mesh_sphere(radius=radius, resolution=(20, 20)),
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)),
)
problem = cpt.RadiationProblem(
    body=body, free_surface=np.inf, water_depth=np.inf, radiating_dof="Surge"
)

custom_result = custom_solver.solve(problem)
default_result = default_solver.solve(problem)

print("Added mass:")
print(f"    Analytical value: {2/3*pi*problem.rho*radius**3:.2e}")
print(f"    Custom solver:    {custom_result.added_masses['Surge']:.2e}")
print(f"    Default solver:   {default_result.added_masses['Surge']:.2e}")
# The analytical solution is known for this particular test case.
# The custom solver might be more accurate. That is probably a coincidence.
# Both should converge to the same solution when the mesh is refined.

print()
print("Timer:")
print(f"    Custom solver:  {custom_solver.timer['  Green function'].mean:.1e} s")
print(f"    Default solver: {default_solver.timer['  Green function'].mean:.1e} s")

C10. Implementation of a GPU Custom Linear Solver

This is an example of how to implement a custom solver for linear systems that leverages your computers’s GPU, here with PyTorch. Your mileage may vary as performance will vary wildly across different hardware. Please also note that you are expected to already be familiar with GPU frameworks and have a working setup to use this; Capytaine’s developers can offer no support for GPU issues unrelated to Capytaine.

import time
import numpy as np
import capytaine as cpt
import xarray as xr
import torch

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

reference_cpu_solver = cpt.BEMSolver()


# A simple implementation of a linear solver on GPU using PyTorch
def linear_solver_on_GPU(A, b):
    A_on_GPU = torch.tensor(np.array(A), dtype=torch.complex64, device=device)
    b_on_GPU = torch.tensor(b, dtype=torch.complex64, device=device)
    res = torch.linalg.solve(A_on_GPU, b_on_GPU)
    return res.cpu().numpy()

naive_gpu_solver = cpt.BEMSolver(
    engine=cpt.BasicMatrixEngine(linear_solver=linear_solver_on_GPU)
)


# Alternatively, an optimization of the above to avoid doing the same GPU
# transfer and LU factorization twice
latest_A = None
latest_LU_decomp = None
def lu_linear_solver_with_cache_on_GPU(A, b):
    global latest_A, latest_LU_decomp
    if A is not latest_A:
        latest_A = A
        A_on_GPU = torch.tensor(np.array(A), dtype=torch.complex64, device=device)
        latest_LU_decomp = torch.linalg.lu_factor(A_on_GPU)
    b_on_GPU = torch.tensor(np.array(b), dtype=torch.complex64, device=device).reshape(-1, 1)
    # Reshaping because `torch.linalg.lu_solve` wants a column vector on the
    # right-hand-side (unlike `torch.linalg.solve`).
    res = torch.linalg.lu_solve(*latest_LU_decomp, b_on_GPU).reshape(-1)
    return res.cpu().numpy()

gpu_solver = cpt.BEMSolver(
    engine=cpt.BasicMatrixEngine(linear_solver=lu_linear_solver_with_cache_on_GPU)
)


###############
#  Benchmark  #
###############

mesh = cpt.mesh_sphere(resolution=(50, 50)).immersed_part()
print(mesh.nb_faces, "faces")
body = cpt.FloatingBody(mesh, dofs=cpt.rigid_body_dofs())

test_matrix = xr.Dataset(
    coords={
        "omega": np.linspace(0.0, 4.0, 20),
        "radiating_dof": list(body.dofs),
        "water_depth": [np.inf],
    }
)

start = time.perf_counter()
ds1 = reference_cpu_solver.fill_dataset(test_matrix, body)
print("CPU:", time.perf_counter() - start)

start = time.perf_counter()
ds2 = naive_gpu_solver.fill_dataset(test_matrix, body)
print("GPU:", time.perf_counter() - start)

start = time.perf_counter()
ds3 = gpu_solver.fill_dataset(test_matrix, body)
print("Slightly faster GPU:", time.perf_counter() - start)


# Check outputs
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(
    ds1.omega.values,
    ds1.added_mass.sel(radiating_dof="Roll", influenced_dof="Roll").values,
)
ax.plot(
    ds2.omega.values,
    ds2.added_mass.sel(radiating_dof="Roll", influenced_dof="Roll").values,
)
ax.plot(
    ds3.omega.values,
    ds3.added_mass.sel(radiating_dof="Roll", influenced_dof="Roll").values,
)
plt.show()