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. Single rigid body hydrodynamics

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

cpt.set_logging("INFO")

# Generating a geometric mesh.
mesh = cpt.mesh_horizontal_cylinder(
    length=10.0,
    radius=1.0,
    center=(0, 0, 0,),
    resolution=(10, 20, 30,),
    name="cylinder_mesh",
)
# Consider also using `cpt.load_mesh` to open an existing file.

# Define a rigid body using this mesh
body = cpt.FloatingBody(
    mesh=mesh,
    lid_mesh=mesh.generate_lid(),  # Generate a lid for irregular frequency removal
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)),
    name="floating cylinder",
)

# Set up the problems: we will solve a radiation problem for each
# degree of freedom and a diffraction problem for each wave direction,
# both for each frequency in the frequency range.
test_matrix = xr.Dataset(
    {
        "omega": np.linspace(0.1, 3.0, 20),
        "wave_direction": [0.0, np.pi/2],
        "radiating_dof": list(body.dofs),
    }
)
# Water density, gravity and water depth have not been specified.
# Default values are used.

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


# Export data in various formats
from pathlib import Path
output_dir = Path("A1_example_outputs")
output_dir.mkdir(exist_ok=True)
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")
cpt.export_dataset(output_dir / "nemoh_data", dataset, format="nemoh")


# Plot the added mass of each dofs as a function of the frequency
import matplotlib.pyplot as plt
plt.figure()
for dof in body.dofs:
    plt.plot(
        dataset.coords["omega"],
        dataset["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()

language:

python

A2. Simulation with several rigid bodies

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

cpt.set_logging("INFO")

# Define the first body
mesh_1 = cpt.mesh_sphere(radius=1.0, center=(0, 0, 0), resolution=(20, 20))
sphere = cpt.FloatingBody(
    mesh=mesh_1,
    lid_mesh=mesh_1.generate_lid(),
    dofs=cpt.rigid_body_dofs(only=["Surge", "Heave"], rotation_center=(0, 0, 0)),
    center_of_mass=(0, 0, 0),
    name="sphere_1",
)

# Define the second body
mesh_2 = cpt.mesh_sphere(radius=0.5, center=(-2, -3, 0), resolution=(20, 20))
other_sphere = cpt.FloatingBody(
    mesh=mesh_2,
    lid_mesh=mesh_2.generate_lid(),
    dofs=cpt.rigid_body_dofs(only=["Surge", "Heave"], rotation_center=(-2, -3, 0)),
    center_of_mass=(-2, -3, 0),
    name="sphere_2",
)

# Define the third body
mesh_3 = cpt.mesh_horizontal_cylinder(
    length=5.0, radius=1.0, center=(0.5, 3.0, -0.5), resolution=(3, 20, 20)
)
cylinder = cpt.FloatingBody(
    mesh=mesh_3,
    lid_mesh=mesh_3.generate_lid(),
    dofs=cpt.rigid_body_dofs(only=["Surge", "Heave"], rotation_center=(0.5, 3.0, -0.5)),
    center_of_mass=(0.5, 3.0, -0.5),
    name="cylinder",
)

# Combine the three individual bodies into a single body.
all_bodies = cpt.Multibody([sphere, other_sphere, cylinder])

# all_bodies.mesh.immersed_part().show()

print("Merged body name:", all_bodies.name)
print("Merged body dofs:", list(all_bodies.dofs.keys()))
print("Multibody centers of mass:", all_bodies.center_of_mass)

# Compute multibody hydrostatics
M = all_bodies.compute_rigid_body_inertia()
K = all_bodies.compute_hydrostatic_stiffness()

with np.printoptions(precision=2, suppress=True):
    print(M)
    print(K)

# Compute multibody hydrodynamics
test_matrix = xr.Dataset(
    {
        "omega": np.linspace(0.5, 2.0, 2),
        "wave_direction": [0.0],
        "radiating_dof": list(all_bodies.dofs),
        "water_depth": [np.inf],
    }
)
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, all_bodies.immersed_part())

print(dataset)

A3. Finite depth flap

This example illustrates the following features: finite depth, rotation around an arbitrary axis (flap’s hinge) and fixed obstacle without degrees of freedom.

"""Single body with a rotation center different from the center of mass.
Also a second body that is just a fixed obstacle.
"""

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

cpt.set_logging("INFO")

water_depth = 10.0
flap_height = 9.0
flap_draft = 8.0
flap_width = 6.0
base_height = water_depth - flap_draft - 0.1

flap_mesh = cpt.mesh_parallelepiped(
    size=(flap_width, 0.5, flap_height),
    center=(0.0, 0.0, -flap_draft + flap_height / 2),
    faces_max_radius=0.2,
    missing_sides={"top"},
)
flap_body = cpt.FloatingBody(
    mesh=flap_mesh,
    lid_mesh=flap_mesh.generate_lid(),
    dofs=cpt.rigid_body_dofs(only=["Pitch"], rotation_center=(0, 0, -flap_draft)),
    # The Pitch dof is defined around the hinge at z=-flap_draft.
    center_of_mass=(0, 0, -flap_draft + flap_height / 2),
    name="flap",
)

base_mesh = cpt.mesh_parallelepiped(
    size=(1.1 * flap_width, 0.8, base_height),
    center=(0.0, 0.0, -water_depth + base_height / 2),
    faces_max_radius=0.2,
    missing_sides={"bottom"},
)
base_body = cpt.FloatingBody(
    mesh=base_mesh,
    # no lid for fully immersed base
    # no dofs for fixed body
    center_of_mass=(0, 0, -10.90),
    name="base",
)

full_flap = flap_body + base_body

test_matrix = xr.Dataset(
    coords={
        "wavelength": np.linspace(5.0, 30.0, 10),
        "radiating_dof": ["flap__Pitch"],
        "wave_direction": [0.0],
        "water_depth": [water_depth],
        "rho": [1000.0],
    }
)

solver = cpt.BEMSolver()
dataset = solver.fill_dataset(
    test_matrix, full_flap.immersed_part(water_depth=water_depth)
)

dataset.added_mass.sel(radiating_dof="flap__Pitch", influenced_dof="flap__Pitch").plot(
    x="wavelength"
)
plt.show()

A4. Forward speed

"""
Example file based on Malenica 1995:
Š. Malenica, P.J. Clark, B. Molin,
Wave and current forces on a vertical cylinder free to surge and sway,
Applied Ocean Research, Volume 17, Issue 2, 1995, Pages 79-90,
https://doi.org/10.1016/0141-1187(95)00002-I.

This script is based on the work of the Maritime Technology Division (MTD) of Ghent University, Belgium.
Please refer to following papers:

I. Herdayanditya, L. Donatini, G. Verao Fernandez, A. B. K. Pribadi, and P. Rauwoens,
“Waves-current effect investigation on monopile excitation force employing approximate forward speed approach,”
in 6th MASHCON : international conference on ship manoeuvring in shallow and confined water
with special focus on port manoeuvres, Glasgow, UK, 2022, pp. 73–80.
http://hdl.handle.net/1854/LU-8756871

L. Donatini, I. Herdayanditya, G. Verao Fernandez, A. B. K. Pribadi, E. Lataire, and G. Delefortrie,
“Implementation of forward speed effects on an open source seakeeping solver,”
in 6th MASHCON : international conference on ship manoeuvring in shallow and confined water
with special focus on port manoeuvres, Glasgow, UK, 2022, pp. 20–33.
http://hdl.handle.net/1854/LU-8756868

"""

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

g = 9.81
rho = 1025.0

water_depth = 1.0
radius = 1.0
center = (0, 0, -0.5)
resolution = (2, 20, 20)

mesh = cpt.mesh_vertical_cylinder(
    length=1.001 * water_depth, radius=radius, center=center, resolution=resolution
)
cylinder = cpt.FloatingBody(
    mesh=mesh, dofs=cpt.rigid_body_dofs(rotation_center=center), center_of_mass=center
)
cylinder = cylinder.immersed_part(water_depth=water_depth)

wavenumber_range = np.linspace(0.2, 2, 25)
froude_number_range = np.array([0, 0.05, -0.05])
forward_speed_range = froude_number_range * np.sqrt(g * radius)

test_matrix = xr.Dataset(
    coords={
        "wavenumber": wavenumber_range,
        "forward_speed": forward_speed_range,
        "wave_direction": [0.0],
        "radiating_dof": ["Surge"],
        "water_depth": [water_depth],
        "rho": [rho],
    }
)
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, cylinder)
dataset["excitation_force"] = (
    dataset["Froude_Krylov_force"] + dataset["diffraction_force"]
)

dataset.coords["encounter_omega"] = dataset["omega"] - dataset["wavenumber"] * dataset[
    "forward_speed"
] * np.cos(dataset["wave_direction"])
dataset.coords["Froude_number"] = dataset.coords["forward_speed"] / np.sqrt(g * radius)

plt.figure()
dataset["norm_excitation"] = np.abs(dataset["excitation_force"] / (g * rho * radius**2))
dataset["norm_excitation"].attrs["long_name"] = "Normalized excitation force F/ρgr²"
dataset["norm_excitation"].sel(wave_direction=0.0, influenced_dof="Surge").plot(
    x="wavenumber", hue="Froude_number"
)

plt.figure()
dataset["norm_added_mass"] = dataset["added_mass"] / (rho * radius**3)
dataset["norm_added_mass"].attrs["long_name"] = "Normalized added mass A/ρr³"
dataset["norm_added_mass"].sel(influenced_dof="Surge", radiating_dof="Surge").plot(
    x="wavenumber", hue="Froude_number"
)

plt.figure()
dataset["norm_rad_damping"] = dataset["radiation_damping"] / (
    rho * dataset.encounter_omega * radius**3
)
dataset["norm_rad_damping"].attrs["long_name"] = "Normalized radiation damping B/ρωr³"
dataset["norm_rad_damping"].sel(influenced_dof="Surge", radiating_dof="Surge").plot(
    x="wavenumber", hue="Froude_number"
)

print(dataset)
plt.show()

A5. Benchmark plane symmetry

"""
This example benchmarks the performance of the mesh plane symmetry.
The same problems are solved, first on the full reference mesh, then using one symmetry plane, then using two symmetry planes.
"""

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

cpt.set_logging("WARNING")

reference_mesh = cpt.mesh_parallelepiped(resolution=(30, 30, 30), name="reference_mesh")
print("Number of faces in full mesh:", reference_mesh.nb_faces)

# Cut the mesh to build symmetric mesh object
half_mesh = reference_mesh.clipped(origin=(0, 0, 0), normal=(1, 0, 0))
simple_symmetric_mesh = cpt.ReflectionSymmetricMesh(
    half=half_mesh, plane="yOz", name="simple_symmetric_mesh"
)

quarter_mesh = half_mesh.clipped(origin=(0, 0, 0), normal=(0, 1, 0))
nested_symmetric_mesh = cpt.ReflectionSymmetricMesh(
    half=cpt.ReflectionSymmetricMesh(half=quarter_mesh, plane="xOz"),
    plane="yOz",
    name="nested_symmetric_mesh",
)

for mesh in [reference_mesh, simple_symmetric_mesh, nested_symmetric_mesh]:
    body = cpt.FloatingBody(
        mesh=mesh,
        lid_mesh=mesh.generate_lid(),
        dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, -0.25)),
        name=mesh.name,
    )

    test_matrix = xr.Dataset(
        {
            "omega": np.linspace(0.5, 2.0, 5),
            "wave_direction": np.linspace(0, np.pi, 3),
            "radiating_dof": list(body.dofs),
        }
    )

    solver = cpt.BEMSolver()
    dataset = solver.fill_dataset(
        test_matrix,
        body.immersed_part(),
        n_threads=1,  # No parallel resolution for cleaner benchmark
    )
    print()
    print("--->", body.name)
    print(
        "Added mass",
        dataset.added_mass.sel(radiating_dof="Heave", influenced_dof="Heave").values,
    )
    print("Computation time", solver.timer_summary())
    print()

A6. Benchmark axial symmetry

"""
This example benchmarks the performance of the axial symmetry.
The same problems are solved, first on the full reference mesh, then using the rotation-symmetric mesh.
"""

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

cpt.set_logging("WARNING")

meridian_points = np.array(
    [(np.sqrt(1 - z**2), 0.0, z) for z in np.linspace(-1.0, 1.0, 50)]
)
sphere_mesh = cpt.RotationSymmetricMesh.from_profile_points(
    meridian_points, n=50
).immersed_part()


symmetric_body = cpt.FloatingBody(
    mesh=sphere_mesh,
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)),
    name="symmetric",
)

full_body = cpt.FloatingBody(
    mesh=sphere_mesh.merged(),  # Discard the symmetry
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)),
    name="full",
)

for body in [full_body, symmetric_body]:
    test_matrix = xr.Dataset(
        {
            "omega": np.linspace(0.5, 2.0, 5),
            "wave_direction": np.linspace(0, np.pi, 3),
            "radiating_dof": list(body.dofs),
            "water_depth": [8.0],
            "rho": [1025],
        }
    )

    solver = cpt.BEMSolver()
    dataset = solver.fill_dataset(
        test_matrix,
        body,
        n_threads=1,  # No parallel resolution for cleaner benchmark
        progress_bar=False,
    )

    print()
    print("--->", body.name)
    print(
        "Added mass",
        dataset.added_mass.sel(radiating_dof="Surge", influenced_dof="Surge").values,
    )
    print(
        "Diffraction",
        dataset.diffraction_force.real.sel(
            wave_direction=0.0, influenced_dof="Surge"
        ).values,
    )
    print("Computation time", solver.timer_summary())
    print()

A7. 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

mesh = cpt.mesh_sphere(radius=1.0, center=(0, 0, 0), faces_max_radius=0.2)

# Initialize floating body
body = cpt.FloatingBody(mesh, lid_mesh=mesh.generate_lid())

# DEFINE THE DOFS

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

# Bulging of the body,
# 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.
body.dofs["Bulge"] = body.mesh.faces_normals

# Shearing of the body in the x direction.
# The deformation vector on a face is computed from the position of the face.
body.dofs["x-shear"] = np.array(
    [(np.cos(np.pi * z / 2), 0, 0) for x, y, z in body.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=body, 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=body, 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 body.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.

A8. 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, 0),
        resolution=(
            int(resolution * radius),
            int(2 * pi * radius * resolution),
            int(length * resolution),
        ),
    ).immersed_part()
    body = cpt.FloatingBody(
        mesh,
        lid_mesh=mesh.generate_lid(),
        dofs=cpt.rigid_body_dofs(only=["Heave"]),
        name=f"cylinder_{mesh.nb_faces:04d}"
    )
    # Store the number of panels in the name of the body
    return body


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

bodies = [make_cylinder(n) for n in np.linspace(1.0, 9.0, 8)]
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()

A9. Test irregular frequencies removal

This example compare the output without a lid and with a lid at different position.

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, name="body_without_lid")
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, name="body_with_immersed_lid"
)
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, name="body_with_lid_on_free_surface"
)
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()

A10. 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()

A11. Parametric study: water depth

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

"""Example of parametric study by giving several water depth in the test matrix"""

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,
    center=(0, 0, 0),
    faces_max_radius=0.2
    ).immersed_part()

body = cpt.FloatingBody(
    mesh=mesh,
    lid_mesh=mesh.generate_lid(),
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0))
)

# 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()