Issue with radiation wave elevation (probably related to phases)

[sdillenburg]

Dear all,

I have a question regarding the free surface elevation of radiated waves. I computed them based on the cookbook examples using the new compute_free_surface_elevation in version 2.0 (source: here) and get_free_surface_elevation in version 1.5 and 2.0 (source: here and here)

Each function yields different wave elevations (code to generate the plots is at the end of this issue):

2.

3.

4.

5.

6.

The results of get_free_surface_elevation (plots 1. and 2. as well as plots 4. and 5.) seem to differ in phase only. The results of compute_free_surface_elevation look very different. With motion RAOs for a wave direction equal pi/2, I would consider the results of the new function in 2.0 (plot 3.) to be correct, since the radiation should be in both directions. But edit: Since the radiated waves due to roll and sway cannot be symmetrical (opposite phases on port and starboard side) for a given time instant, the total radiated wave pattern for a certain time instant cannot be symmetrical either. Only heave radiated waves are symmetric to the x axis. This can be seen in the plots 7., 8. and 9. For 90°, roll and sway contributions are cancelled out by the heave contribution on the side where x is positive.

I would expect plot 6. (motion RAOs for 270° wave direction) to be symmetric to plot 3. (90°).

8.

9.

Additionally, results obtained with version 1.5 show a better agreement with tank tests (really small 1:144 model in a narrow tank, so I consider the agreement acceptable) compared to compute_free_surface_elevation results (plots consider incident, radiated and diffracted waves):

Here's the code to produce the first six plots in this issue

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

logging.basicConfig(level=logging.INFO, format='%(levelname)-8s: %(message)s')

version = float(cpt.__version__)

# Environment variables
directions = [90, 270]
omega = 1.0
water_depth = np.inf

# Free surface
xmin = -150
xmax = 150
ymin = -150
ymax = 150
nx = 60
ny = 60
x = np.linspace(xmin, xmax, nx, endpoint=True)
y = np.linspace(ymin, ymax, ny, endpoint=True)
fs = cpt.FreeSurface(x_range=(xmin, xmax), y_range=(ymin, ymax), nx=nx, ny=ny)

# Define box shaped body
length = 90  # m
beam = 30  # m
draft = 6 # m
zcg = 3

body = cpt.RectangularParallelepiped(size=(length, beam, draft), resolution=(int(length/2), int(beam/2), max(int(draft),5)),
                                           center=(0, 0, -0.5*draft), top=False, name='box')

body.add_all_rigid_body_dofs()
body.center_of_mass = (0, 0, zcg)
body.inertia_matrix = body.compute_rigid_body_inertia()
body.hydrostatic_stiffness = body.compute_hydrostatic_stiffness()

# Solver
bem_solver = cpt.BEMSolver()

for direction in directions:
    # Problem setup and solve
    if version < 2:
        problems = [cpt.RadiationProblem(omega=omega, body=body, sea_bottom=-water_depth, radiating_dof=dof) for dof in body.dofs]
        problems.append(cpt.DiffractionProblem(omega=omega, body=body, sea_bottom=-water_depth, wave_direction=np.radians(direction)))
    elif version >= 2:
        problems = [cpt.RadiationProblem(omega=omega, body=body, water_depth=water_depth, radiating_dof=dof) for dof in body.dofs]
        problems.append(cpt.DiffractionProblem(omega=omega, body=body, water_depth=water_depth, wave_direction=np.radians(direction)))
    
    results = bem_solver.solve_all(problems)
    *radiation_results, diffraction_result = results
    
    dataset = cpt.assemble_dataset(results) 
    
    # Postprocessing   
    rao = cpt.post_pro.rao(dataset, wave_direction=np.radians(direction))

    print(direction, rao)
            
    radiation_elevations_per_dof = {res.radiating_dof: (-1j*omega)*bem_solver.get_free_surface_elevation(res, fs) 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)
    
    def plot_free_surface_elevation(free_surface_elevation, filename):
        fig, ax = plt.subplots(figsize=(12, 9))
    
        ax.set_title(filename)
        ax.set_aspect('equal', 'box')
        ax.set_xlabel('x [m]')
        ax.set_ylabel('y [m]')
        
        
        fs_elevation = np.abs(free_surface_elevation)*np.sin(np.angle(free_surface_elevation))
        fs_elevation = fs_elevation.reshape((nx, ny)).T
        
        h = ax.contourf(x, y, fs_elevation, levels=np.linspace(-1.,1.,9),extend='both')
        
        box = patches.Rectangle((-length/2, -beam/2), length, beam, linewidth=1, edgecolor='k', facecolor='gray')
        ax.add_patch(box)
        
        plt.colorbar(h, ax=ax, aspect=30, orientation='horizontal', label='Wave elevation [m]')
    
        plt.savefig(filename)
    
        
    plot_free_surface_elevation(radiation_elevation, 'free_surface_dir_{}__version_{}.svg'.format(direction, version))
    
    if version >= 2:
        # Postprocessing using new free surface elevation functions
        radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs.mesh, diffraction_result) 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)
        
        plot_free_surface_elevation(radiation_elevation, 'free_surface_version_dir_{}_{}_new_functions.svg'.format(direction, version))

Summarizing my questions:

1. Which approach is the correct one?
2. Although get_free_surface_elevation is deprecated in 2.0, the results should be that same in 1.5 and 2.0, shouldn't they?
3. Why do the radiated waves computed with compute_free_surface_elevation differ for 90° and 270°? This is not due to the instant of time, see following plots which show the absolute values:

Any hint for further investigation is appreciated!

[mancellin]

Thank you for reporting this.

I don't have the time to look in depth into it right now, but here are some comments:

• the result is expected to be different by a factor of $-j \omega$ between v1.5 and v2.0. Some background and motivation can be found in Unit amplitude or unit velocity in solver #173. Radiation problems (and not diffraction problems) used to need a manual correction of their output. That is the (-1j*omega) term in the following line, that is only required in v1.5 and should be removed in v2.0¹.

radiation_elevations_per_dof = {res.radiating_dof: (-1j*omega)*bem_solver.get_free_surface_elevation(res, fs) for res in radiation_results}

• in v2.0, compute_free_surface_elevation and get_free_surface_elevation should give the same result. The former just has a different call signature that is meant to be the same in all of Capytaine (e.g. the same as airy_wave_free_surface_elevation).

• in one of the last lines of your code, you are iterating other the radiation results but use the diffraction result to compute the free surface elevation:

radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs.mesh, diffraction_result) for res in radiation_results}

But there has been some reports of suspicious results (such as #191 and #387) so there might be a bug somewhere in there.

Footnotes

1. In retrospect, I could have made the change only in compute_free_surface_elevation and kept the old convention in the legacy get_free_surface_elevation. These were two different changes and I did not think of that in time. ↩

[sdillenburg]

@mancellin

Your comments solved the entire issue ;)
Many thanks!!

Your last comment is a copy & paste error from here: https://github.com/capytaine/capytaine/blob/v2.0/docs/user_manual/examples/boat_animation.py#L44

So the example should be corrected.

For the sake of completeness, here's the adjusted code:

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

logging.basicConfig(level=logging.INFO, format='%(levelname)-8s: %(message)s')

version = float(cpt.__version__)

# Environment variables
directions = [90, 270]
omega = 1.0
water_depth = np.inf

# Free surface
xmin = -150
xmax = 150
ymin = -150
ymax = 150
nx = 60
ny = 60
x = np.linspace(xmin, xmax, nx, endpoint=True)
y = np.linspace(ymin, ymax, ny, endpoint=True)
fs = cpt.FreeSurface(x_range=(xmin, xmax), y_range=(ymin, ymax), nx=nx, ny=ny)

# Define box shaped body
length = 90  # m
beam = 30  # m
draft = 6 # m
zcg = 3

body = cpt.RectangularParallelepiped(size=(length, beam, draft), resolution=(int(length/2), int(beam/2), max(int(draft),5)),
                                           center=(0, 0, -0.5*draft), top=False, name='box')

body.add_all_rigid_body_dofs()
body.center_of_mass = (0, 0, zcg)
body.inertia_matrix = body.compute_rigid_body_inertia()
body.hydrostatic_stiffness = body.compute_hydrostatic_stiffness()

# Solver
bem_solver = cpt.BEMSolver()

for direction in directions:
    # Problem setup and solve
    if version < 2:
        radiation_problems = [cpt.RadiationProblem(omega=omega, body=body, sea_bottom=-water_depth, radiating_dof=dof) for dof in body.dofs]
        diffraction_problem  = cpt.DiffractionProblem(omega=omega, body=body, sea_bottom=-water_depth, wave_direction=np.radians(direction))
    elif version >= 2:
        radiation_problems = [cpt.RadiationProblem(omega=omega, body=body, water_depth=water_depth, radiating_dof=dof) for dof in body.dofs]
        diffraction_problem  = cpt.DiffractionProblem(omega=omega, body=body, water_depth=water_depth, wave_direction=np.radians(direction))
    
    radiation_results = bem_solver.solve_all(radiation_problems)
    diffraction_result = bem_solver.solve(diffraction_problem)
    
    results = radiation_results + [diffraction_result]
    
    dataset = cpt.assemble_dataset(results) 
    
    # Postprocessing   
    rao = cpt.post_pro.rao(dataset, wave_direction=np.radians(direction))
    
    if version < 2:    
        radiation_elevations_per_dof = {res.radiating_dof: (-1j*omega)*bem_solver.get_free_surface_elevation(res, fs) for res in radiation_results}
    elif version >= 2:
        radiation_elevations_per_dof = {res.radiating_dof: bem_solver.get_free_surface_elevation(res, fs) 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)
    
    def plot_free_surface_elevation(free_surface_elevation, filename):
        fig, ax = plt.subplots(figsize=(12, 9))
    
        ax.set_title(filename)
        ax.set_aspect('equal', 'box')
        ax.set_xlabel('x [m]')
        ax.set_ylabel('y [m]')
        
        
        fs_elevation = np.abs(free_surface_elevation)*np.sin(np.angle(free_surface_elevation))
        fs_elevation = fs_elevation.reshape((nx, ny)).T
        
        h = ax.contourf(x, y, fs_elevation, levels=np.linspace(-1.,1.,9),extend='both')
        
        box = patches.Rectangle((-length/2, -beam/2), length, beam, linewidth=1, edgecolor='k', facecolor='gray')
        ax.add_patch(box)
        
        plt.colorbar(h, ax=ax, aspect=30, orientation='horizontal', label='Wave elevation [m]')
    
        plt.savefig(filename)
    
        
    plot_free_surface_elevation(radiation_elevation, 'free_surface_dir_{}__version_{}_w_get_free_surface_elevation.svg'.format(direction, version))
    
    if version >= 2:
        # Postprocessing using new free surface elevation functions
        radiation_elevations_per_dof = {res.radiating_dof: bem_solver.compute_free_surface_elevation(fs.mesh, 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)
        
        plot_free_surface_elevation(radiation_elevation, 'free_surface_dir_{}__version_{}_w_compute_free_surface_elevation.svg'.format(direction, version))

[mancellin]

Your last comment is a copy & paste error from here: https://github.com/capytaine/capytaine/blob/v2.0/docs/user_manual/examples/boat_animation.py#L44

Ah, yes indeed. Would you like to submit a PR with the correction?

[sdillenburg]

Yes, will do (I have not done this befoe, so I probably need some time to dig in the github details ;))