ZeppSav / BEMUse

A lightweight boundary element method library for hydrodynamics and aerodynamics.
GNU General Public License v3.0
9 stars 3 forks source link

Possible bug in Hydrodynamic_Solver.cpp? #1

Closed jnvn7 closed 10 months ago

jnvn7 commented 1 year ago

Thank you for this great tool!

One small thing pops out when looking at the code, should Omega = 2*PI*Frequency on the following line and not Omega = Frequency/2/PI?

https://github.com/ZeppSav/BEMUse/blob/ecab031dfcf0bf0526beede962557668968654d9/src/Solver/Hydrodynamic_Solver.cpp#L374

ZeppSav commented 1 year ago

Hi there, Thanks for the kudos, I hope you are enjoying using BEMUse. If you look a few lines below this point you will see where I calculate the wave number Kappa. A quick dimension check indicates that indeed the variable Omega represents the frequency. The code was verified against a couple of other solvers, some of which use (f) as input and some use (Omega) as input. While checking this I probably starting using the wrong naming convention.

This should not influence the output of the solver however, besides your requirement to scale the input / output by a factor of either 2PI or 1/2PI. If you could report back to me on this I will try to correct the nomenclature in the code ASAP. Thanks!

jnvn7 commented 1 year ago

@ZeppSav Thank you for your quick reply. I also thought that it was a naming convention but then seeing the wavenumber (Kappa) is calculated as omega*omega/g for deep (infinite) water, omega should be (rad/s) and not (1/s or Hz), which is likely represented by f in this case since it is used to calculate the Period (1/Frequency) a few lines down. With this convention, omega should be 2*PI*Frequency.

With that said, let me know what you think. I haven't look too deep into the code and definitely are not familiar with the program's structure as much as you do. Thanks!

ZeppSav commented 1 year ago

Thank you for pointing out the issue with the naming, I will try to correct this ASAP. This however is quite a superficial issue and is simply a matter of modifying the names and double cheking that the outputs correspond correctly.

Is there actually an issue with the functioning of the software besides this? I am guessing you were going through the code to try to debug something right? As mentioned above the software should still calculate correctly, besides the necessity to scale the inputs / output frequencies.

jnvn7 commented 1 year ago

Sorry for the late response. I wanted to do some more benchmarking tests before reaching out.

I went through the code to check the unit of the input frequencies and whether that should be in Hz or rad/s, or something else. I agree it does not affect the function of the program but it is good to know the software is performing calculations on the desired frequencies.

On a different note, I also performed some benchmark tests comparing BEMUse, Wamit, HAMS, and Capytaine. I did two cases including a truncated cylinder and a hemisphere. While the results show good correlation among the programs for translational modes (surge, sway, heave), the outputs from BEMUse are different compared to the rest in rotational. I'm not sure if there is an issue with the setup or due to something else.

I've attached the plots and the case setup files for your reference. I only include BEMUse and HAMS. Setup for Wamit and Capytaine are similarly done and they all share the same mesh.

Thank you!

Truncated Cylinder:

CylinderTrans

CylinderRot

Hemisphere:

SphereTrans

SphereRot

CaseSetup.zip

ZeppSav commented 1 year ago

Hi nhu-vnguyen,

No problems regarding taking time, I understand how time consuming testing can be. I really appreciate you taking the time to go through the results and report them. So let's go through the above plots/issues step by step.

i) Re. Frequency input / output: As we have previously written, you may have had to scale the input/output frequency to align with the other softwares. If indeed you did do this, what exactly were the modifications you had to make? I can then correct this.

ii) Re. Translational DoFs. I observe also the agreement, which is obviously pleasing and aligns well with the validaitons I did myself.

iii) Re. Rotational DoFs- Hemisphere: For the Hemipshere due to the nature of potential flow and the symmetry of the body in all axes of rotation, these values are theoretically zero. All that you are seeing here (in all codes) is the result of discretising the surface and using floating point operations- there is inevitably some "numerical noise" which comes through to the results. This is observable as the magnitude of these values are about five orders of magnitude below the corresponding results for translation. In all cases the actual magnitude of the loads acting in these DoFs are vanishingly small, and for all codes I would expect it to get even smaller with increased mesh resolution. The fact that BEMUse is greater than the values of the others will have practically no influence on calculations of dynamics. It is possibly a result of the fact that BEMUse is compiled with single floating point precision in order to reduce memory overhead and slightly improve performance. I am unsure if this is also the case for the other codes, which may have been compiled with double floating point precision. Either way you can test this yourself by compiling BEMUse with the compilation flag DoublePrec instead of the defaultSinglePrec.

iv) Re. Rotational DoFs- Cylinder- Yaw: The previous argument applies equally to the yaw DoF for the cylinder - symmetry means it should practically be zero besides numerical noise.

v) Re. Rotational DoFs- Cylinder- Pitch and Roll: Here indeed a difference can be observed which is somewhat harder to explain. For the calculation of rotational DoFs A_ij and B_ij , a cross product is taken between the surface node position and the specified "Reference position" of the geometry. Within BEMuse this is currently always specified to be the origin [0,0,0]. If another position was chosen for this, then the integrated value over the geometry will be different. I observe in your HAMS inputs that a variable _Reference_bodycentre has been defined at the (I assume) CoG of the body. If this is being used for this calculation rather than the origin, this very well might explain the discrepancy. You can check this by re-running the simulation (in HAMS or any other solver) with this position set to [0,0,0]. This actually might also explain the differences in point iii) above if those codes are also compiled for single floating point precision.

If you could try that modification it would be greatly appreciated, otherwise we can perhaps also have a telcon. My email address is in the source files. Thanks again and looking forward to your response.

ZeppSav commented 1 year ago

Hi nhu-vnguyen, Were you able to try the changes I suggested above and resolve the problem? I would like to be able to close this issue ;)

jnvn7 commented 1 year ago

Hi @ZeppSav, sorry for the late reply, and thank you for the suggestions. I have rerun the simulations with CG at [0,0,0] and the issues were resolved. The comparisons with other solvers are great!