Beginner’s cookbook

This page contains several examples of Capytaine’s main features, mostly about the computation of added mass, radiation damping and excitation force. The scripts can be downloaded individually as Python files from the Github folder with examples for version 3.0.dev. See also the template in Quickstart.

A1. Added mass of a single rigid body

This example generates the mesh of an horizontal cylinder, solves radiation problems for the six rigid-body degrees of freedom and then plots the added mass.

import numpy as np
import capytaine as cpt

cpt.set_logging()

# Generating a geometric mesh
mesh = cpt.mesh_horizontal_cylinder(
    length=10.0,
    radius=1.0,
    center=(0, 0, -2,),
    resolution=(10, 20, 30,)
)

# Define a rigid body using this mesh
cylinder = cpt.FloatingBody(
    mesh=mesh,
    dofs=cpt.rigid_body_dofs(),
    name="floating cylinder"
)

# Define the range of frequencies as a Numpy array
omega_range = np.linspace(0.1, 3.0, 10)

# Set up the problems: we will solve a radiation problem for each
# degree of freedom of the body and for each frequency in the
# frequency range.
problems = [
    cpt.RadiationProblem(body=cylinder, radiating_dof=dof, omega=omega)
    for dof in cylinder.dofs
    for omega in omega_range
]
# Water density, gravity and water depth have not been specified.
# Default values are used.

# Solve all radiation problems
solver = cpt.BEMSolver()
results = solver.solve_all(problems)

# Gather the computed added mass into a labelled array.
data = cpt.assemble_dataset(results)

# Plot the added mass of each dofs as a function of the frequency
import matplotlib.pyplot as plt
plt.figure()
for dof in cylinder.dofs:
    plt.plot(
        omega_range,
        data['added_mass'].sel(radiating_dof=dof, influenced_dof=dof),
        label=dof,
        marker='o',
    )
plt.xlabel('omega')
plt.ylabel('added mass')
plt.legend()
plt.tight_layout()
plt.show()

A2. Simulation with several bodies

import capytaine as cpt

# Define the first body
sphere = cpt.FloatingBody(
        cpt.mesh_sphere(radius=1.0, center=(0, 0, -2.0),
                        resolution=(20, 20)),
        name="sphere_1")
sphere.add_translation_dof(name="Surge")
sphere.add_translation_dof(name="Heave")

# Define the second body
other_sphere = cpt.FloatingBody(
        cpt.mesh_sphere(radius=0.5, center=(-2, -3, -1),
                        resolution=(20, 20)),
        name="sphere_2")
other_sphere.add_translation_dof(name="Surge")
other_sphere.add_translation_dof(name="Heave")

# Define the third body
cylinder = cpt.FloatingBody(
        cpt.mesh_horizontal_cylinder(
            length=5.0, radius=1.0,
            center=(1.5, 3.0, -3.0),
            resolution=(20, 20, 3)),
        name="cylinder")
cylinder.add_translation_dof(name="Surge")
cylinder.add_translation_dof(name="Heave")

# Combine the three individual bodies into a single body.
all_bodies = cylinder + sphere + other_sphere

print("Merged body name:", all_bodies.name)
print("Merged body dofs:", list(all_bodies.dofs.keys()))

# The merged body can be used to define the problems in the usual way
problems = [cpt.RadiationProblem(body=all_bodies, radiating_dof=dof, omega=1.0) for dof in all_bodies.dofs]
problems += [cpt.DiffractionProblem(body=all_bodies, wave_direction=0.0, omega=1.0)]

# Solves the problem
solver = cpt.BEMSolver()
results = solver.solve_all(problems)
data = cpt.assemble_dataset(results)

print(data)

A3. Influence of the water depth

This example runs the same simulation for several water depth and plot the results.

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

cpt.set_logging('INFO')

# Initialize floating body by generating a geometric mesh
mesh = cpt.mesh_horizontal_cylinder(
    length=10.0, radius=0.5,  # Dimensions
    center=(0, 0, 0),         # Position
    resolution=(5, 20, 40)    # Fineness of the mesh
    ).immersed_part()
body = cpt.FloatingBody(mesh, dofs=cpt.rigid_body_dofs())

# Define the range of water depth
depth_range = np.concatenate([np.linspace(2, 100, 10), [np.inf]])

test_matrix = xr.Dataset(coords={
    "wavelength": 100.0,
    "water_depth": depth_range,
    "radiating_dof": ["Heave"]
    })

# Solve all radiation problems
solver = cpt.BEMSolver()
data = solver.fill_dataset(test_matrix, body)

# Note that the solver could not solve the most shallow case and returned NaN:
data.added_mass.sel(radiating_dof="Heave", influenced_dof="Heave", water_depth=2.0).values

# Plot the added mass of each dofs as a function of the water depth.
fig, ax = plt.subplots(layout="constrained")
ax.plot(
    data.water_depth * data.wavenumber,
    data['added_mass'].sel(radiating_dof="Heave", influenced_dof="Heave"),
    marker="s", label="Finite depth"
)
ax.hlines(
    data['added_mass'].sel(radiating_dof="Heave", influenced_dof="Heave", water_depth=np.inf),
    xmin=5, xmax=8, label="Infinite depth",
    linestyle="--"
)
ax.set(xlabel='wavenumber * water_depth', ylabel='Heave-heave added mass')
ax.legend()
plt.show()

A4. Custom degree of freedom

This example defines arbitrary degrees of freedom for a sphere and solves a diffraction problem.

import numpy as np
import capytaine as cpt

# Initialize floating body
sphere = cpt.FloatingBody(
        mesh=cpt.mesh_sphere(
            radius=1.0,          # Dimension
            center=(0, 0, -2),   # Position
            resolution=(20, 20), # Fineness of the mesh
            ))

# DEFINE THE DOFS

# Manually defined heave,
# that is a vertical unit vector for all faces.
sphere.dofs["Heave"] = np.array(
    [(0, 0, 1) for face in sphere.mesh.faces]
)

# Bulging of the sphere,
# that is a deformation vector normal to the body at all faces.
# We can simply point to the normal vectors of the mesh for that.
sphere.dofs["Bulge"] = sphere.mesh.faces_normals

# Shearing of the sphere in the x direction.
# The deformation vector on a face is computed from the position of the face.
sphere.dofs["x-shear"] = np.array(
    [(np.cos(np.pi*z/2), 0, 0) for x, y, z in sphere.mesh.faces_centers]
)

# SOLVE DIFFRACTION PROBLEMS
solver = cpt.BEMSolver()

# Solve the problem for β=0 (incoming wave in the x direction).
problem_1 = cpt.DiffractionProblem(body=sphere, wave_direction=0, omega=1.0)
result_1 = solver.solve(problem_1)

# Solve the problem for β=π/2 (incoming wave in the y direction).
problem_2 = cpt.DiffractionProblem(body=sphere, wave_direction=np.pi/2, omega=1.0)
result_2 = solver.solve(problem_2)

# Print the generalized diffraction forces
# for the three dofs we defined
# for both values of the wave_direction β.
for result in [result_1, result_2]:
    print(f"Angle: {result.wave_direction:.2f}")
    for dof in sphere.dofs:
        force = result.forces[dof]
        print(f"{dof}: {np.abs(force):.2f}·exp({np.angle(force):.2f}i) N")
    print()

The diffraction force on the “Heave” and “Bulge” dofs should be the same for both incoming wave directions. The diffraction force on the “x-shear” dof is zero when the wave comes from the y direction.

A5. Convergence study

This example runs a mesh convergence study for a submerged cylinder.

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

cpt.set_logging('INFO')

def make_cylinder(resolution):
    """Make cylinder with a mesh of a given resolution in panels/meter."""
    radius = 1.0
    length = 5.0
    mesh = cpt.mesh_horizontal_cylinder(
        length=length, radius=radius,
        center=(0, 0, -1.5*radius),
        resolution=(int(resolution*radius),
                    int(2*pi*length*resolution),
                    int(length*resolution))
    )
    body = cpt.FloatingBody(mesh, name=f"cylinder_{mesh.nb_faces:04d}")
    # Store the number of panels in the name of the body
    body.add_translation_dof(name="Heave")
    return body

test_matrix = xr.Dataset(coords={
    'omega': [1.0],
    'radiating_dof': ['Heave'],
})

bodies = [make_cylinder(n) for n in np.linspace(1, 5, 10)]
ds1 = cpt.BEMSolver().fill_dataset(test_matrix, bodies)

def read_nb_faces_in_mesh_name(ds):
    """Read the name of the body to guess the resolution of the mesh."""
    ds.coords['nb_faces'] = xr.DataArray([int(name[9:]) for name in ds['body_name'].values], coords=[ds['body_name']])
    ds = ds.swap_dims({'body_name': 'nb_faces'})
    return ds
ds1 = read_nb_faces_in_mesh_name(ds1)

ds1['added_mass'].plot(x='nb_faces')

plt.show()

A5. Forward speed

TODO

A6. Irregular frequencies removal

import numpy as np
import capytaine as cpt

cpt.set_logging('INFO')

radius = 12.5
draught = 37.5
n_radius = 15
n_theta = 25
n_z = 30

# Define the range of frequencies to solve
omega_range = np.linspace(1.5, 2.25, 50)

solver = cpt.BEMSolver()

def compute_excitation_force(body, omega_range):
    problems = [
        cpt.DiffractionProblem(body=body, omega=omega)
        for omega in omega_range
    ]
    results = solver.solve_all(problems)
    data = cpt.assemble_dataset(results)
    return data['excitation_force']


# Initialize floating body by generating a geometric mesh
cylinder_mesh = cpt.mesh_vertical_cylinder(
    length=draught*2, radius=radius, center=(0, 0, 0), faces_max_radius=1.0,
    ).immersed_part()

### Body Without Lid (Body A)
body_1 = cpt.FloatingBody(mesh=cylinder_mesh)
body_1.add_translation_dof(name="Surge")

### Body With Lid at z = -3.699 (Body B)
lid_mesh = cylinder_mesh.generate_lid(z=-3.699)
body_2 = cpt.FloatingBody(mesh=cylinder_mesh, lid_mesh=lid_mesh)
body_2.add_translation_dof(name="Surge")

### Body With Lid at z = auto (Body C)
lid_mesh = cylinder_mesh.generate_lid(z=cylinder_mesh.lowest_lid_position(omega_max=omega_range.max()))
body_3 = cpt.FloatingBody(mesh=cylinder_mesh, lid_mesh=lid_mesh)
body_3.add_translation_dof(name="Surge")

force_1 = compute_excitation_force(body_1, omega_range)
force_2 = compute_excitation_force(body_2, omega_range)
force_3 = compute_excitation_force(body_3, omega_range)

# Plot the added mass of each dofs as a function of the frequency
import matplotlib.pyplot as plt
plt.figure()
plt.plot(
    omega_range,
    np.abs(force_1.sel(influenced_dof='Surge')),
    marker='+',
    label='no lid'
)
plt.plot(
    omega_range,
    np.abs(force_2.sel(influenced_dof='Surge')),
    marker='*',
    label=r'$z = -3.699$m'
)
plt.plot(
    omega_range,
    np.abs(force_3.sel(influenced_dof='Surge')),
    marker='x',
    label='auto'
)

plt.xlabel(r'$\omega$')
plt.ylabel(r'$F_1 (abs)$')
plt.legend()
plt.tight_layout()
plt.show()

A7. Elastic beam

import numpy as np
import xarray as xr
import capytaine as cpt
import matplotlib.pyplot as plt
from scipy.optimize import newton

# We consider a vertical cylinder clamped at the sea bottom and freely moving at the top above the free surface.

water_depth = 10

cylinder_length = 12.0
cylinder_density = 1000.0
cylinder_diameter = 1.0
young_modulus = 1e9
inertia_moment = np.pi * cylinder_diameter**4 / 64  # Cylinder with circular section

# The following code computes the shapes of the deformation modes of the body and the associated eigenfrequencies.

def deformation_mode(i_mode, z):
    alpha = newton(lambda a: np.cos(a)*np.cosh(a) + 1, (2*i_mode + 1)*np.pi/2)
    coeff = (-np.cos(alpha) - np.cosh(alpha))/(np.sin(alpha) + np.sinh(alpha))
    k = alpha/cylinder_length
    z_ = z + water_depth
    return 0.5 * (np.cos(k*z_) - np.cosh(k*z_) + coeff*(np.sin(k*z_) - np.sinh(k*z_)))

def eigenfrequency(i_mode):
    alpha = newton(lambda a: np.cos(a)*np.cosh(a) + 1, (2*i_mode + 1)*np.pi/2)
    return alpha/cylinder_length**2 * np.sqrt(young_modulus*inertia_moment/(cylinder_density*np.pi*(cylinder_diameter/2)**2))


# Plotting the x-deformation of the first 4 modes
z_range = np.linspace(-water_depth, -water_depth + cylinder_length, 100)

# for i_mode in range(4):
#     plt.plot(deformation_mode(i_mode, z_range), z_range, label=f"mode_{i_mode}")
# plt.legend()


mesh = cpt.mesh_vertical_cylinder(
        center=(0.0, 0.0, cylinder_length/2 - water_depth - 0.1),
        length=cylinder_length,
        radius=cylinder_diameter/2,
        resolution=(4, 15, 36),
        )

nb_modes = 10
dofs = {}
for i in range(nb_modes):
    dofs[f"mode_{i}"] = np.array([(deformation_mode(i, z), 0.0, 0.0) for x, y, z in mesh.faces_centers])

full_body = cpt.FloatingBody(mesh=mesh, dofs=dofs, center_of_mass=(0, 0, cylinder_length/2 - water_depth))
body = full_body.immersed_part(water_depth=water_depth)

body.inertia_matrix = body.add_dofs_labels_to_matrix(np.eye(nb_modes))
body.internal_stiffness_matrix = body.add_dofs_labels_to_matrix(np.diag([eigenfrequency(i_mode) for i_mode in range(nb_modes)]))
body.hydrostatic_stiffness = body.compute_hydrostatic_stiffness()

test_matrix = xr.Dataset(coords={
    "omega": [1.0],
    "radiating_dof": list(body.dofs),
    "wave_direction": [0.0],
    "water_depth": [water_depth],
    })
ds = cpt.BEMSolver().fill_dataset(test_matrix, body)
rao = cpt.post_pro.rao(ds, stiffness=body.internal_stiffness_matrix)


deformation_profile = sum(
    rao.isel(radiating_dof=i_mode).values[0] * deformation_mode(i_mode, z_range)
    for i_mode, dof in enumerate(body.dofs)
    )
plt.plot(np.real(deformation_profile), z_range)

plt.show()

A8. Export dataset

import numpy as np
from pathlib import Path
import xarray as xr

import capytaine as cpt

# --- Parameters ---
output_dir = Path("outputs")
output_dir.mkdir(exist_ok=True)

# --- Create mesh and set up body ---
mesh = cpt.mesh_sphere()
dofs = cpt.rigid_body_dofs(rotation_center=(0, 0, 0))
full_body = cpt.FloatingBody(mesh, dofs)
full_body.center_of_mass = np.copy(full_body.mesh.center_of_mass_of_nodes)
immersed_body = full_body.immersed_part()
immersed_body.compute_hydrostatics()

# --- Define parameter grid ---
omegas = np.concatenate(([0], np.linspace(0.1, 1.0, 10), [np.inf]))  # 0, 5 values, inf
wave_directions = np.linspace(0, np.pi, 3)  # 3 directions: 0, pi/2, pi

test_matrix = xr.Dataset(
    {
        "omega": omegas,
        "wave_direction": wave_directions,
        "radiating_dof": list(immersed_body.dofs),
        "water_depth": [np.inf],
        "rho": [1025],
    }
)

# --- Run simulation ---
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, immersed_body)


cpt.export_dataset(output_dir / "test_boat200.nc", dataset)  # Format is deduced from filename suffix.
cpt.export_dataset(output_dir / "wamit_data", dataset, format="wamit")
(output_dir / 'nemoh_data').mkdir(exist_ok=True)
cpt.export_dataset(output_dir / "nemoh_data", dataset, format="nemoh")

print("Exported files:")
for file in output_dir.glob("**/*"):
    print(file)