Post-processing and collection in dataset

Outputs stored in LinearPotentialFlowResult

Each LinearPotentialFlowResult object contains all the data defined for the corresponding problem:

problem = cpt.DiffractionProblem(body=body, period=1.0, wave_direction=3.14)
solver = cpt.BEMSolver()
result = solver.solve(problem, keep_details=True)

print(result.period)
# 1.0

They also contain a supplementary attribute named forces, containing the computed forces on the body:

print(result.forces)
# {'Surge': (23941.728129534735+9487.048661471363j), 'Sway': (-2.010438038269058e-09-2.8862814360763878e-09j), ...}

The force is stored as Python dictionary associating a complex value in SI units to the name of the influenced degree of freedom. In other words, in the example code above, the \(x\)-component of the force vector has magnitude 23941 Newton at \(t =0\), while the \(y\)-component of the force vector is negligible at all time.

For radiation problems, the result object also contain added_mass and radiation_damping attributes, with the same shape of a Python dictionary.

It the solver was called with keep_details=True, the result object also contains the potential and pressure fields on each face of the mesh used for the resolution (that is, including the lid faces, if a lid has been defined on the body to remove irregular frequencies):

print(result.potential)
# [-1.72534485e-03+0.14128629j -3.98932611e-03+0.33387497j
#  -6.54135830e-03+0.54208628j ... -2.09853806e+00-0.56106653j
#  -2.18441640e+00-0.58402701j -2.22777755e+00-0.59562008j]

print(result.pressure)
# [-1.72534485e-03+0.14128629j -3.98932611e-03+0.33387497j
#  -6.54135830e-03+0.54208628j ... -2.09853806e+00-0.56106653j
#  -2.18441640e+00-0.58402701j -2.22777755e+00-0.59562008j]

These magnitudes are stored in an one-dimensional array as long as the number of faces of the mesh, and stored in the same order as the faces in the mesh_including_lid object. In other words, result.pressure[3] contains the pressure on the face of center body.mesh_including_lid.faces_centers[3].

Since the pressure on the lid is not physically meaningful, the result object also provides result.pressure_on_hull, which is result.pressure restricted to the faces of body.mesh (that is, excluding the lid). When no lid has been defined, result.pressure_on_hull and result.pressure are the same.

Recall that the potential and the pressure are related by \(p = j \rho \omega \Phi\), where \(\rho\) is the fluid density and \(\omega\) is the angular frequency.

The potential, the pressure and the velocity in the rest of the fluid can be computed in post-processing as described below.

When using the indirect boundary integral equation (by default), the result object with details also contains the distribution of sources \(\sigma\) on the hull, under the same form as the potential above:

print(result.sources)
# [  0.0281344 -0.35719578j   0.06715056-1.34127138j
#    0.10891061-2.26921669j ... -13.58069557+2.13219904j
#  -14.13645749+2.21945488j -14.41706935+2.26351156j]

Building a dataframe or a dataset from LinearPotentialFlowResult

If you have a list of LinearPotentialFlowResult, you can assemble them in a pandas DataFrame or a xarray dataset for a more convenient post-processing. Use the following syntax:

dataframe = cpt.assemble_dataframe(list_of_results)

or:

dataset = cpt.assemble_dataset(list_of_results)

If you gave a test matrix to the BEMSolver.fill_dataset method, the output will directly be an xarray dataset.

Errored resolutions stored as a FailedDiffractionResult or a FailedRadiationResult will appear as NaN in the dataset.

Both assemble_dataset and fill_dataset accept some optional keyword arguments to store more information in the dataset:

• wavenumber, wavelength, period, omega, (default: all True): control whether which of the representations of the wave frequency are stored in the dataset. At least one should be included, by default they all are.

• hydrostatics (default: True): tries to compute the hydrostatic properties of the FloatingBody and add them to the dataset.

• mesh (default: False): stores in the dataset the full mesh for the hull and the lid, as well as the dofs. - The hull mesh is stored using the hull_face dimension in the arrays mesh_vertices (shape (hull_face, 4, 3)` and mesh_faces_center (shape: (hull_face, 3)). - If the body has a lid, it is stored similarly using the lid_face dimension in the arrays lid_mesh_vertices (shape (lid_face, 4, 3) and lid_mesh_faces_center (shape: (lid_face, 3)). - If the body has dofs, the array dof_motions (shape: (influenced_dof, hull_face, 3)) stores the displacement at each face center due to each dof. - If the dofs have a Jacobian defined (as do rigid body dofs or CustomDof when setting the gradient_of_motion parameter), it is stored also in an array of shape (influenced_dof, hull_face, 3, 3). - The name of the quadrature scheme used for the mesh is stored as a string in the quadrature_method variable.

Note

If the body’s mesh takes advantage of a planar or rotational symmetry (e.g. a ReflectionSymmetricMesh), the symmetry itself is not preserved by mesh=True, as of the current version of Capytaine: mesh_vertices, mesh_faces_center and dof_motions always contain the full, explicit list of faces (as if merged() had been called first), with nothing in the dataset indicating that the original mesh was symmetric. This only affects the exported mesh/dofs data itself; the hydrodynamic coefficients are unaffected and are still computed by taking advantage of the symmetry internally.

In addition, fill_dataset accepts a keep_details keyword argument (default: False), forwarded to solve_all(). When set to True, the potential, pressure and sources computed for each problem are kept and the resulting dataset contains the pressure distribution on the hull as the data variables pressure_diffraction, pressure_Froude_Krylov and pressure_radiation (with a hull_face dimension, in the same order as the faces of body.mesh):

dataset = solver.fill_dataset(test_matrix, body, keep_details=True)
print(dataset["pressure_diffraction"])

Keeping these details increases the memory usage of the dataset, so this option is disabled by default.

Note

The code does its best to keep the degrees of freedom in the same order as they have been provided by the user, but there is no guarantee that this will always be the case. Nonetheless, in all cases the labels will always match the data. So selecting a value by the name of the dof will always return the right one:

data.sel(radiating_dof="Heave", influenced_dof="Heave")

You can also manually reorder the dofs with the following syntax:

sorted_dofs = ["Surge", "Sway", "Heave", "Roll", "Pitch", "Yaw"]
print(data.sel(radiating_dof=sorted_dofs, influenced_dof=sorted_dofs))

Note

Datasets created with assemble_dataset only include data on cases with a free surface. Cases without a free surface (free_surface=inf) are ignored.

The results can also be collected by assemble_matrices(), which returns some matrices of assemble_dataset() as numpy arrays stripped of their metadata. This function is meant to be used for teaching, to assemble the matrices without getting the students in contact with xarray.

Pressure, velocity, free surface elevation

Once the problem has been solved, several fields of interest can be computed at post-processing:

Code	Description
solver.compute_potential(points, result)	The velocity potential \(\phi(x, y, z)\)
solver.compute_pressure(points, result)	The pressure in the fluid \(p(x, y, z)\)
solver.compute_velocity(points, result)	The velocity of the fluid \(u(x, y, z)\)
solver.compute_free_surface_elevation(points, result)	The elevation of the free surface \(\eta(x, y)\)

All the methods listed above work in the same way: they require the LinearPotentialFlowResult object containing the required data about the solved problem and some points at which the field should be evaluated.

The result object should have been computed with the indirect method (on by default) and the option keep_details=True (on by default for solve(), off by default for solve_all()). The solver does not need to be the one that computed the result object.

Note

The functions in the airy_waves module, used to compute the same magnitudes for an undisturbed incoming wave field, have the same structure.

The point(s) can be given in several ways:

• Either a single point, given as a list, a tuple, or an 1d-array:

  solver.compute_potential([3.0, -2.0, -5.0], result)
  

• or a list of points, given as a list of lists, or a list of tuples, or a 2d-array:

  solver.compute_potential([[3.0, -2.0, -5.0], [4.0, 5.0, -2.0]], result)
  

• or the return of a call to meshgrid:

  points = np.meshgrid(np.linspace(-2.0, 2.0, 10), np.linspace(-3.0, 3.0, 20), np.linspace(-2.0, 0.0, 30))
  solver.compute_potential(points, result)
  

• or a mesh, in which case the centers of the faces of the mesh are used:

  solver.compute_potential(mesh, result)
  

in particular, the following can be used to compute the value on the mesh used for the resolution:

solver.compute_potential(result.body.mesh, result)

• or a floating body, in which case the corresponding mesh will be used:

  solver.compute_potential(body, result)
  

The returned values is an array of shape matching the shape of the input points.

Warning

There is a single case in which passing a mesh is not equivalent to a list of point: if you want the compute the velocity on the hull of the floating body. In this case, you should give the same mesh object that has been used for the resolution:

solver.compute_velocity(result.body.mesh, result)

Other Python objects might return incorrect values or errors.

For potential, pressure and velocity, 3 coordinates \((x, y, z)\) are expected for each points. For the free surface elevation, 2 coordinates \((x, y)\) are sufficient.

Note

In the limit cases of zero and infinite frequencies, these magnitudes can also be computed. Strictly speaking, their value is 0 (resp. infinity) at these frequencies. Nonetheless, the following magnitudes \(\phi/\omega\), \(u/\omega\), \(\eta/\omega^2\) and \(p/\omega^2\) are well defined and have a finite value. They can be accessed in Capytaine as in the following example:

[...]
pb = cpt.RadiationProblem(body=body, omega=0, radiating_dof="Heave")
res = solver.solve(pb, keep_details=True)
pressure_over_omega_2 = solver.compute_pressure(points, res)/(pb.omega*pb.omega)

Impedance and RAO

The intrinsic impedance can be computed based on the hydrodynamics, hydrostatics, and inertial properties:

import numpy as np
import xarray as xr
from capytaine import BEMSolver
from capytaine.bodies.predefined.spheres import Sphere
from capytaine.post_pro import impedance

f = np.linspace(0.1, 2.0)
omega = 2*np.pi*f
rho_water = 1e3
r = 1

sphere = Sphere(radius=r, ntheta=3, nphi=12, clip_free_surface=True)
sphere.center_of_mass = np.array([0, 0, 0])
sphere.add_all_rigid_body_dofs()

sphere.inertia_matrix = sphere.compute_rigid_body_inertia(rho=rho_water)
sphere.hydrostatic_stiffness = sphere.compute_hydrostatic_stiffness(rho=rho_water)

solver = BEMSolver()
test_matrix = xr.Dataset(coords={
    'rho': rho_water,
    'water_depth': [np.inf],
    'omega': omega,
    'wave_direction': 0,
    'radiating_dof': list(sphere.dofs.keys()),
    })

data = solver.fill_dataset(test_matrix, sphere_fb,
                           hydrostatics=True,
                           mesh=True,
                           wavelength=True,
                           wavenumber=True)

Zi = impedance(data)

Note that we assigned the inertia and stiffness to attributes of body called inertia_matrix and hydrostatic_stiffness. These are the names expected by the fill_dataset and impedance functions to compute the impedance matrix.

By simple extension of incorporating the excitation transfer function response amplitude operator (RAO):

from capytaine.post_pro import rao
rao = rao(data)

Kochin functions

The function compute_kochin() can compute the Kochin function of the waves for a LinearPotentialFlowResult containing the source distribution (that is, if it has been solved with the indirect method):

from capytaine.post_pro.kochin import compute_kochin
thetas = np.linspace(0, 2*np.pi, 20)
kochin_values = compute_kochin(result, thetas)

Alternatively, if the test_matrix of fill_dataset() contains a coordinate called theta, it is used to compute the Kochin function of all solved problems (see cookbook example).

Mean drift forces

The horizontal mean drift forces can be computed using the function far_field_mean_drift_force(), here is an example:

import numpy as np
import xarray as xr
import capytaine as cpt
from capytaine.post_pro.mean_drift_force import far_field_mean_drift_force

mesh = cpt.mesh_sphere().immersed_part()
body = cpt.FloatingBody(mesh=mesh, dofs=cpt.rigid_body_dofs(), center_of_mass=(0,0,0))
body.inertia_matrix = body.compute_rigid_body_inertia()
body.hydrostatic_stiffness = body.compute_hydrostatic_stiffness()
solver = cpt.BEMSolver()
k = [1.16, 1.24, 1.39]
theta = np.linspace(0, 2*np.pi, 10)
test_matrix = xr.Dataset(coords={
            'wavenumber': k ,
            'wave_direction' : 0,
            'theta': theta,
            'radiating_dof': list(body.dofs.keys())
        })
dataset = solver.fill_dataset(test_matrix, body, hydrostatics=True)
X = cpt.post_pro.rao(dataset)
mean_drift_forces = far_field_mean_drift_force(X, dataset)