Intermediate cookbook

This page contains several examples of Capytaine’s features, mostly focused on post-processing features such as reconstruction of the free surface elevation or computation of the Kochin functions. The scripts can be downloaded individually as Python files from the Github folder with examples for version 3.0.dev.

B1. Plot pressure on hull

This example solves a BEM problem and display the computed pressure on the hull on a 3D view of the mesh.

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

mesh = cpt.mesh_sphere(faces_max_radius=0.1).immersed_part()
body = cpt.FloatingBody(
    mesh=mesh,
    dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0))
)

pb = cpt.RadiationProblem(wavelength=1.0, body=body, radiating_dof="Surge", water_depth=10.0)
solver = cpt.BEMSolver()
res = solver.solve(pb, keep_details=True)

body.show_matplotlib(
        color_field=np.real(res.pressure),
        cmap=plt.get_cmap("viridis"),  # Colormap
        )

B2. Haskind’s relation

This example computes the excitation force from the radiation potential using Haskind’s relation. The result is compared with the one obtained by direct integration of the potentials from incident waves and from the diffraction problem.

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_potential, airy_waves_velocity, froude_krylov_force

r = 1.
draft = 0.5
depth = 10.
omega = 1.
rho = 1000

# Define geometry and heave degree of freedom
body = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            length=2*draft, radius=r, center=(0.,0.,0.),
            resolution=(10, 50, 10)
            ))
body.add_translation_dof(name='Heave')
body = body.immersed_part()

solver = cpt.BEMSolver()

# Define and solve the diffraction and radiation problems
diff_problem = cpt.DiffractionProblem(body=body, water_depth=depth,
                                      omega=omega, wave_direction=0.)
rad_problem = cpt.RadiationProblem(body=body, water_depth=depth,
                                   omega=omega, radiating_dof='Heave')

diff_solution = solver.solve(diff_problem)
rad_solution = solver.solve(rad_problem)


# Read mesh properties
faces_centers = body.mesh.faces_centers
faces_normals = body.mesh.faces_normals
faces_areas = body.mesh.faces_areas

# Computation from the diffraction solution (Capytaine)
FK = froude_krylov_force(diff_problem)['Heave']
diff_force = diff_solution.forces['Heave']

# Get potentials
phi_inc = airy_waves_potential(faces_centers, diff_problem)
v_inc = airy_waves_velocity(faces_centers, diff_problem)
phi_rad = rad_solution.potential

# Indirect computation from the radiation solution, via the Haskind relation
integrand = - (phi_inc * faces_normals[:,2]
              - 1j/omega * phi_rad * np.diag(v_inc@faces_normals.T))
F_hask = 1j * omega * rho * np.sum(integrand*faces_areas)

# Direct integration of the potentials
diff_force_recomputed = - 1j * omega * rho * np.sum(
                        diff_solution.potential*faces_areas*faces_normals[:,2])
FK_recomputed = - 1j * omega * rho * np.sum(
                        phi_inc*faces_areas*faces_normals[:,2])

print('Result from Capytaine = ', FK + diff_force)
print('Result from recomputed direct integration',
        FK_recomputed + diff_force_recomputed)
print('Result from Haskind relation = ', F_hask)

B3. Free surface elevation

This example computes the free surface elevation as a post-processing step of a diffraction problem.

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
import matplotlib.pyplot as plt

mesh = cpt.mesh_sphere(radius=1, resolution=(10, 10)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs())

pb = cpt.DiffractionProblem(body=body, wave_direction=np.pi/2, wavelength=3.0)
solver = cpt.BEMSolver()
res = solver.solve(pb)

grid = np.meshgrid(np.linspace(-10.0, 10.0, 100), np.linspace(-10.0, 10.0, 100))
fse = solver.compute_free_surface_elevation(grid, res)
incoming_fse = airy_waves_free_surface_elevation(grid, res)

fig, ax = plt.subplots(1, 1)
ax.pcolormesh(grid[0], grid[1], np.real(fse + incoming_fse))
ax.add_patch(plt.Circle((0, 0), radius=1, color="black"))
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_aspect("equal")
ax.set_title("Wave crests and troughs at t=0 seen from above")
plt.show()

B4. Kochin function

This example computes the Kochin function for a surging buoy and plot the results.

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')

# SET UP THE MESH
buoy = cpt.FloatingBody(
        mesh=cpt.mesh_vertical_cylinder(
            radius=1.0, length=3.0, center=(0, 0, -1.0),
            resolution=(15, 15, 3)
            ))
buoy.add_translation_dof(name="Surge")

# SOLVE THE BEM PROBLEMS AND COMPUTE THE KOCHIN FUNCTIONS
theta_range = np.linspace(0, 2*pi, 40)
omega_range = np.linspace(1.0, 5.0, 3)

test_matrix = xr.Dataset(coords={
    'omega': omega_range, 'theta': theta_range, 'radiating_dof': ["Surge"],
})
solver = cpt.BEMSolver()
dataset = solver.fill_dataset(test_matrix, buoy.immersed_part())

# PLOT THE KOCHIN FUNCTION
plt.figure()
dataset = dataset.swap_dims({'omega': 'wavenumber'})  # Use wavenumber to index data
for k in dataset.coords['wavenumber']:
    (dataset['kochin_radiation']
     .sel(radiating_dof="Surge", wavenumber=k)
     .real
     .plot(x='theta', label=f"k={float(k):.2f} m¯¹", ax=plt.gca()))
plt.legend()
plt.title("Real part of the Kochin function as a function of theta")
plt.tight_layout()
plt.show()

B5. Plot velocity in domain

This example reconstruct the velocity field in the domain and display it using matplotlib.

import numpy as np
import capytaine as cpt

# Solve the potential flow problems
radius = 1.0
mesh = cpt.mesh_horizontal_cylinder(radius=radius, length=1.0, center=(0, 0, 0), resolution=(5, 20, 5)).immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs())
solver = cpt.BEMSolver()
results = solver.solve_all([cpt.RadiationProblem(body=body, omega=1.0, radiating_dof=dof) for dof in body.dofs])

# Define a set of points in the fluid on which the fluid velocities will be computed
Y, Z = np.meshgrid(np.linspace(-3.0, 3.0, 40), np.linspace(-3.0, -0.05, 20))
points = np.vstack([np.zeros(Y.size), Y.ravel(), Z.ravel()]).T
# Filter out points that are inside the floating body
points = np.array([[x, y, z] for (x, y, z) in points if y**2 + z**2 > radius**2])

# Compute the velocity field for each of the result object, that is each of the radiating dof
velocities = {}
for res in results:
    velocities[res.radiating_dof] = solver.compute_velocity(points, res)

import matplotlib.pyplot as plt
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
for ax, dof in zip(axs.ravel(), ["Sway", "Heave"]):
    # Display the arrows of the velocity
    y_velocities = -velocities[dof][:, 1].imag
    z_velocities = -velocities[dof][:, 2].imag
    ax.quiver(points[:, 1], points[:, 2], y_velocities, z_velocities, angles="xy")
    ax.axis('equal')

    # Display the boundary of the floating body
    theta = np.linspace(-np.pi, 0.0, 100)
    ax.plot(radius*np.cos(theta), radius*np.sin(theta), color="red")
    ax.set_title(dof)
    ax.set_xlabel("y")
    ax.set_ylabel("z")

plt.show()

B6. Animated free surface elevation

This example solves a diffraction problem, it computes the free surface elevation and shows it as a 3D animation.

"""This example shows how to animate the free surface elevation and pressure on a floating body using vedo."""

import numpy as np
import capytaine as cpt
from capytaine.bem.airy_waves import (
    airy_waves_free_surface_elevation,
    airy_waves_pressure,
)
from capytaine.ui.vedo_animations import Animation


# Parameters
omega = 3.0
T = 2 * np.pi / omega


## RESOLUTION WITH CAPYTAINE
positions = [
    (np.cos(angle) * 5.0, np.sin(angle) * 5.0, 0.0)
    for angle in [0.0, 2 * np.pi / 3, -2 * np.pi / 3]
]
meshes = [
    cpt.mesh_vertical_cylinder(length=4.0, radius=2.0, center=p, faces_max_radius=0.3)
    for p in positions
]
mesh = cpt.Mesh.join_meshes(*meshes)
sphere = cpt.FloatingBody(mesh=mesh.immersed_part(), lid_mesh=mesh.generate_lid())
solver = cpt.BEMSolver()
diffraction_problem = cpt.DiffractionProblem(body=sphere, wave_direction=0.0, omega=omega)
print("Solving the diffraction problem...")
diffraction_result = solver.solve(diffraction_problem)
print("Diffraction problem solved.")


## POST-PROCESSING THE PRESSURE AND FREE SURFACE ELEVATION
free_surface_mesh = cpt.mesh_rectangle(
    size=(100.0, 100.0), center=(0.0, 0.0, 0.0), resolution=(300, 300)
)
amplitude = 0.2
print("Computing the free surface elevation and pressure...")
diffraction_elevation_at_faces = amplitude * (
    solver.compute_free_surface_elevation(free_surface_mesh.vertices, diffraction_result)
    + airy_waves_free_surface_elevation(free_surface_mesh.vertices, diffraction_problem)
)
pressure_on_sphere = amplitude * (
    solver.compute_pressure(sphere.mesh.vertices, diffraction_result)
    + airy_waves_pressure(sphere.mesh.vertices, diffraction_problem)
)
print("Free surface elevation and pressure computed.")


## ANIMATION WITH VEDO
anim = Animation(loop_duration=T, fps=24)
anim.add_body(sphere, vertices_pressure=pressure_on_sphere)
anim.add_free_surface(free_surface_mesh, vertices_elevation=diffraction_elevation_at_faces)
anim.save("diffraction.mp4", size=(1024, 512))

B7. Animation of the RAO

This script generates the animation of the RAO motion for a wave incoming in front of a ship, such as the one used on the main page of this documentation. This script requires the mesh of the ship boat_200.mar. It can be downloaded from: https://raw.githubusercontent.com/capytaine/capytaine/master/docs/examples/src/boat_200.mar

from numpy import pi

import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
from capytaine.ui.vedo_animations import Animation


bem_solver = cpt.BEMSolver()


def generate_boat():
    boat_mesh = cpt.load_mesh("docs/examples/src/boat_200.mar", file_format="mar")
    boat = cpt.FloatingBody(
            mesh=boat_mesh,
            dofs=cpt.rigid_body_dofs(rotation_center=boat_mesh.center_of_buoyancy),
            center_of_mass = boat_mesh.center_of_buoyancy,
            name="pirate ship"
            )
    boat.inertia_matrix = boat.compute_rigid_body_inertia() / 10 # Artificially lower to have a more appealing animation
    boat.hydrostatic_stiffness = boat.immersed_part().compute_hydrostatic_stiffness()
    return boat


def setup_animation(body, fs, omega, wave_amplitude, wave_direction):
    # SOLVE BEM PROBLEMS
    radiation_problems = [cpt.RadiationProblem(omega=omega, body=body.immersed_part(), radiating_dof=dof) for dof in body.dofs]
    radiation_results = bem_solver.solve_all(radiation_problems)
    diffraction_problem = cpt.DiffractionProblem(omega=omega, body=body.immersed_part(), wave_direction=wave_direction)
    diffraction_result = bem_solver.solve(diffraction_problem)

    dataset = cpt.assemble_dataset(radiation_results + [diffraction_result])
    dataset["inertia_matrix"] = body.inertia_matrix
    dataset["hydrostatic_stiffness"] = body.hydrostatic_stiffness
    rao = cpt.post_pro.rao(dataset, wave_direction=wave_direction)

    # COMPUTE FREE SURFACE ELEVATION
    # Compute the diffracted wave pattern
    incoming_waves_elevation = airy_waves_free_surface_elevation(fs.vertices, diffraction_result)
    diffraction_elevation = bem_solver.compute_free_surface_elevation(fs.vertices, diffraction_result)

    # Compute the wave pattern radiated by the RAO
    radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs.vertices, res) for res in radiation_results}
    radiation_elevation = sum(rao.sel(omega=omega, radiating_dof=dof).data * radiation_elevations_per_dof[dof] for dof in body.dofs)

    # SET UP ANIMATION
    # Compute the motion of each face of the mesh for the animation
    dofs_vertices_motions = {dof: body.dofs[dof].evaluate_motion_at_points(body.mesh.vertices) for dof in body.dofs}
    rao_vertices_motion = sum(rao.sel(omega=omega, radiating_dof=dof).data * dofs_vertices_motions[dof] for dof in body.dofs)

    # Set up scene
    animation = Animation(loop_duration=2*pi/omega)
    animation.add_body(body, vertices_motion=wave_amplitude*rao_vertices_motion)
    animation.add_free_surface(fs, wave_amplitude * (incoming_waves_elevation + diffraction_elevation + radiation_elevation))
    return animation


if __name__ == '__main__':
    body = generate_boat()
    fs = cpt.mesh_rectangle(size=(200.0, 200.0), resolution=(100, 100))

    anim = setup_animation(body, fs, omega=1.5, wave_amplitude=0.5, wave_direction=pi)
    # anim.run(camera_position=(70, 70, 100), resolution=(800, 600))
    anim.save("animated_boat.mp4", camera={"pos": (70, 70, 100)}, resolution=(800, 600))

B8. Pressure field for zero or infinite frequency

The computations for zero and infinite frequencies are tricky, because several magnitudes such as the pressure are zero or infinity, but they can be expressed as \(\omega\) or \(\omega^2\) times some finite value. The example below shows an example of extracting this value.

import numpy as np
import capytaine as cpt

# Initialize an example mesh and test case.
mesh = cpt.mesh_sphere().immersed_part()
body = cpt.FloatingBody(mesh, dofs=cpt.rigid_body_dofs(rotation_center=(0, 0, 0)))

pb = cpt.RadiationProblem(body=body, omega=np.inf, radiating_dof="Heave")
# Zero and infinity are accepted values for "omega", as well as for "period", "wavelength" and "wavenumber"

solver = cpt.BEMSolver()
res = solver.solve(pb)

print(f"Heave-heave added mass: {res.added_mass['Heave']:.2e}")

p = np.real(res.pressure)
# p is a custom Python object storing the pressure on each panel

i_panel = 0  # Some panel for testing

# All pressures are actually infinite
print(f"Pressure on panel {i_panel}:", float(p[i_panel]))

# However p/omega^2 is finite (hence the finite added mass)
p_over_omega2 = p/pb.omega**2
print(f"Pressure/omega^2 on panel {i_panel}: {p_over_omega2[i_panel]:.2e}")