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

plt.pcolormesh(grid[0], grid[1], np.real(fse + incoming_fse))
plt.xlabel("x")
plt.ylabel("y")
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']
     .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.

import capytaine as cpt
from capytaine.bem.airy_waves import airy_waves_free_surface_elevation
from capytaine.ui.vtk.animation import Animation

cpt.set_logging('INFO')

# Generate the mesh of a sphere

full_mesh = cpt.mesh_sphere(radius=3, center=(0, 0, 0), resolution=(20, 20))
full_sphere = cpt.FloatingBody(mesh=full_mesh)
full_sphere.add_translation_dof(name="Heave")

# Keep only the immersed part of the mesh
sphere = full_sphere.immersed_part()

# Set up and solve problem
solver = cpt.BEMSolver()

diffraction_problem = cpt.DiffractionProblem(body=sphere, wave_direction=0.0, omega=2.0)
diffraction_result = solver.solve(diffraction_problem)

radiation_problem = cpt.RadiationProblem(body=sphere, radiating_dof="Heave", omega=2.0)
radiation_result = solver.solve(radiation_problem)

# Define a mesh of the free surface and compute the free surface elevation
free_surface = cpt.FreeSurface(x_range=(-50, 50), y_range=(-50, 50), nx=150, ny=150)
diffraction_elevation_at_faces = solver.compute_free_surface_elevation(free_surface, diffraction_result)
radiation_elevation_at_faces = solver.compute_free_surface_elevation(free_surface, radiation_result)

# Add incoming waves
diffraction_elevation_at_faces = diffraction_elevation_at_faces + airy_waves_free_surface_elevation(free_surface, diffraction_problem)

# Run the animations
animation = Animation(loop_duration=diffraction_result.period)
animation.add_body(full_sphere, faces_motion=None)
animation.add_free_surface(free_surface, faces_elevation=0.5*diffraction_elevation_at_faces)
animation.run(camera_position=(-30, -30, 30))  # The camera is oriented towards (0, 0, 0) by default.
# animation.save("path/to/the/video/file.ogv", camera_position=(-30, -30, 30))

animation = Animation(loop_duration=radiation_result.period)
animation.add_body(full_sphere, faces_motion=full_sphere.dofs["Heave"])
animation.add_free_surface(free_surface, faces_elevation=3.0*radiation_elevation_at_faces)
animation.run(camera_position=(-30, -30, 30))
# animation.save("path/to/the/video/file.ogv", camera_position=(-30, -30, 30))

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.vtk import Animation

cpt.set_logging('INFO')

bem_solver = cpt.BEMSolver()


def generate_boat():
    boat_mesh = cpt.load_mesh("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])
    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, diffraction_result)
    diffraction_elevation = bem_solver.compute_free_surface_elevation(fs, diffraction_result)

    # Compute the wave pattern radiated by the RAO
    radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs, 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
    rao_faces_motion = sum(rao.sel(omega=omega, radiating_dof=dof).data * body.dofs[dof] for dof in body.dofs)

    # Set up scene
    animation = Animation(loop_duration=2*pi/omega)
    animation.add_body(body, faces_motion=wave_amplitude*rao_faces_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.FreeSurface(x_range=(-100, 75), y_range=(-100, 75), nx=100, ny=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.ogv", camera_position=(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}")