Added Mass different between v2.0 and v2.1

[jtgrasb]

I am using Capytaine to run the boundary element method for a 2-body attenuator type WEC. Recently, I updated from v2.0 to v2.1 and noticed that I am now getting different results, specifically for the added mass at low frequencies. You can see the added mass plotted for one of the bodies in the surge, heave, and pitch dofs below. v2.0: v2.1:

For context, the body shown is a small disk-like cylinder with the radius of about 0.15 m and height of about 0.075 m. and a relatively small draft (~.03 m). What is causing this change in the added mass between the two Capytaine versions? Should the results from v2.1 be trusted over v2.0?

[mancellin]

Thank you for reporting this.

My main guess is that it is due to the increased resolution of the tabulation, meaning that v2.1 is more reliable.

You could check that by solving in v2.1 with

solver = cpt.BEMSolver(green_function=cpt.Delhommeau(
            tabulation_nr=400, tabulation_nz=80,
            tabulation_nb_integration_points=251,
            tabulation_method="legacy"))

which should give the same result as v2.0 or by solving with

solver = cpt.BEMSolver(green_function=cpt.Delhommeau(tabulation_nr=0))

which should give the same result in both versions, and be closer to v2.1 than v2.0. (Beware that the latter setting is up to 100 times slower, so maybe reduce the number of frequencies.)

Otherwise, could you share a code snippet for me to investigate?

[mancellin]

it is due to the increased resolution of the tabulation

Which is surprisingly missing from the changelog... I need to fix that.

[jtgrasb]

Thank you for the quick response. I checked the two cases you suggested and found the expected results

solver = cpt.BEMSolver(green_function=cpt.Delhommeau(
            tabulation_nr=400, tabulation_nz=80,
            tabulation_nb_integration_points=251,
            tabulation_method="legacy"))

gives me the same results as v2.0. And

solver = cpt.BEMSolver(green_function=cpt.Delhommeau(tabulation_nr=0))

gives me the same results for both. I've included a code snippet for a similar thin cylinder.

And here is the resultant added mass plots using v2.1. I believe the large negative added mass terms at low frequencies are leading to instabilities in our model. I am working on comparing these results to WAMIT for this same geometry.

Python code

```python import capytaine as cpy import matplotlib.pyplot as plt import pygmsh import gmsh import os import xarray as xr import numpy as np print(cpy.__version__) # create mesh r = .35 # cylinder radius h = .1 # cylinder height freeboard = h/2 mesh_size_factor = 0.5 with pygmsh.occ.Geometry() as geom: gmsh.option.setNumber('Mesh.MeshSizeFactor', mesh_size_factor) cyl = geom.add_cylinder([0, 0, 0],[0, 0, -h],r) geom.translate(cyl, [0, 0, freeboard]) mesh = geom.generate_mesh() fb = cpy.FloatingBody.from_meshio(mesh, name="Pioneer") fb.show_matplotlib() fb.keep_immersed_part() fb.add_all_rigid_body_dofs() fb.show_matplotlib() # set frequencies wCapy = np.linspace(.1, 10, 100) # wave frequencies headings = np.linspace(0,np.pi/2,1) depth = np.infty density = 1025.0 # call Capytaine solver print(f'\n-------------------------------\n' f'Calling Capytaine BEM solver...\n' f'-------------------------------\n' f'mesh = {fb.name}\n' f'w range = {wCapy[0]:.3f} - {wCapy[-1]:.3f} rad/s\n' f'dw = {(wCapy[1]-wCapy[0]):.3f} rad/s\n' f'no of headings = {len(headings)}\n' f'no of radiation & diffraction problems = {len(wCapy)*(len(headings) + len(fb.dofs))}\n' f'-------------------------------\n') problems = xr.Dataset(coords={ 'omega': wCapy, 'wave_direction': headings, 'radiating_dof': list(fb.dofs), 'water_depth': [depth], 'rho': [density], }) solver = cpy.BEMSolver() capyData = solver.fill_dataset(problems, [fb], n_jobs = 2, hydrostatics=False) def plot_hydrodynamic_coefficients(bem_data, wave_dir, ): """Plots hydrodynamic coefficients (added mass, radiation damping, and wave excitation) based on BEM data. Parameters ---------- bem_data Linear hydrodynamic coefficients obtained using the boundary element method (BEM) code Capytaine, with sign convention corrected. wave_dir Wave direction(s) to plot the """ bem_data = bem_data.sel(wave_direction = wave_dir, method='nearest') radiating_dofs = bem_data.radiating_dof.values radiating_dofs = ['Surge','Heave','Pitch'] influenced_dofs = bem_data.influenced_dof.values influenced_dofs = ['Surge','Heave','Pitch'] # plots fig_am, ax_am = plt.subplots( 1, len(radiating_dofs), tight_layout=True, sharex=True, figsize=(3*len(radiating_dofs), 3), squeeze=False ) fig_rd, ax_rd = plt.subplots( 1, len(radiating_dofs), tight_layout=True, sharex=True, figsize=(3*len(radiating_dofs), 3), squeeze=False ) fig_ex, ax_ex = plt.subplots( 1, len(radiating_dofs), tight_layout=True, sharex=True, figsize=(3*len(radiating_dofs), 3), squeeze=False ) [ax.grid(True) for axs in (ax_am, ax_rd, ax_ex) for ax in axs.flatten()] # plot titles fig_am.suptitle('Added Mass Coefficients', fontweight='bold') fig_rd.suptitle('Radiation Damping Coefficients', fontweight='bold') fig_ex.suptitle('Wave Excitation Coefficients', fontweight='bold') sp_idx = 0 for i, rdof in enumerate(radiating_dofs): sp_idx += 1 np.abs(bem_data.diffraction_force.sel(influenced_dof=rdof)).plot( ax=ax_ex[0,i], linestyle='dashed', label='Diffraction') np.abs(bem_data.Froude_Krylov_force.sel(influenced_dof=rdof)).plot( ax=ax_ex[0,i], linestyle='dashdot', label='Froude-Krylov') ex_handles, ex_labels = ax_ex[0,i].get_legend_handles_labels() ax_ex[0,i].set_title(f'{rdof}') ax_ex[0,i].set_xlabel('') ax_ex[0,i].set_ylabel('') bem_data.added_mass.sel( radiating_dof=rdof, influenced_dof=rdof).plot(ax=ax_am[0,i]) bem_data.radiation_damping.sel( radiating_dof=rdof, influenced_dof=rdof).plot(ax=ax_rd[0, i]) ax_am[0, i].set_xlabel(f'$\omega$', fontsize=10) ax_rd[0, i].set_xlabel(f'$\omega$', fontsize=10) ax_ex[0, i].set_xlabel(f'$\omega$', fontsize=10) fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False) return #[(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)] plot_hydrodynamic_coefficients(capyData,headings) ```

[jtgrasb]

Confirming that when I run the same geometry in WAMIT, I do not get the low frequency peaks in added mass as shown below (different units cause the difference in magnitude):

[mancellin]

Ok, I can confirm that the default result of 2.1 is wrong here. By comparing with the output for $\omega=0$, the asymptotic is clearly wrong.

[mancellin]

But for some reasons, the current development version on the master branch gives good results :confused:

[mancellin]

Works also well with cpt.BEMSolver(green_function=cpt.XieDelhommeau()), so this might be the low-frequency counterparts of #298

[mancellin]

I can also observe that the low-frequency asymptotics was wrong for the Green function of 2.0 (diverging for $\omega < 0.01$ rad/s). That is expected and I was intending to fix that it in the future at the same time as #298.

It is completely unexpected to me that the change of version 2.1 would make it much worse.

Details

```python import numpy as np import xarray as xr import matplotlib.pyplot as plt import pygmsh import gmsh import capytaine as cpt print(cpt.__version__) # create mesh r = .35 # cylinder radius h = .1 # cylinder height freeboard = h/2 mesh_size_factor = 0.5 with pygmsh.occ.Geometry() as geom: gmsh.option.setNumber('Mesh.MeshSizeFactor', mesh_size_factor) cyl = geom.add_cylinder([0, 0, 0],[0, 0, -h],r) geom.translate(cyl, [0, 0, freeboard]) mesh = geom.generate_mesh() fb = cpt.FloatingBody.from_meshio(mesh, name="Pioneer") # fb.show_matplotlib() fb.keep_immersed_part() fb.add_all_rigid_body_dofs() # fb.show_matplotlib() # set frequencies wCapy = np.linspace(0.0, 0.1, 100) depth = np.infty density = 1025.0 # call Capytaine solver problems = xr.Dataset(coords={ 'omega': wCapy, 'radiating_dof': list(fb.dofs), 'water_depth': [depth], 'rho': [density], }) for solver in [ cpt.BEMSolver(green_function=cpt.Delhommeau(tabulation_nr=0)), cpt.BEMSolver(green_function=cpt.XieDelhommeau(tabulation_nr=0)), cpt.BEMSolver(green_function=cpt.Delhommeau()), cpt.BEMSolver(green_function=cpt.XieDelhommeau()), cpt.BEMSolver(green_function=cpt.Delhommeau(tabulation_nr=400, tabulation_nz=80, tabulation_method="legacy")), ]: ds = solver.fill_dataset(problems, [fb], n_jobs=2, hydrostatics=False) ds.added_mass.sel(radiating_dof="Heave", influenced_dof="Heave").plot(x="omega", label=str(solver)) plt.legend() plt.show() ```

[mancellin]

The low frequency asymptotics should be better in the next version (2.2)