Mismatch between Nemoh and Capytaine when computing potential in domain

[ffrancesco]

I m trying to estimate the velocity potential at a point different then the body mesh, using the get_potential_on_mesh function of the BEMSolver class. During the solution of the problem I enforce the keep_details=True so I can reuse the sources.
I compared the results with Nemoh and WAMIT, but the results are not similar.

Did I interpreted the function correctly by assuming that it can be used to estimate the velocity potential at a different location than the body mesh?

[mancellin]

Yes, that should be the way to do it.

Could you share a bit more about what you tried? How different are the results?

[ffrancesco]

The test case is an Heaving cylinder with radius 5m and draft 10m (supporting file/Cylinder.dat)
Two frequencies are used 0.1 and 2 rad/s with a single wave direction of 0deg.

The "correct" potential at the inscribing cylinder can be found in the supporting file.
cylsurface 1.dat -> diffraction at frequency 0.1
cylsurface 2.dat -> radiation at 0.1 frequency
etc.

I m using this fork of Nemoh (https://github.com/DTOcean/NEMOH) to estimate the potential.
This Nemoh branch was verified against WAMIT.

Running the python script issue_191.py it is possible to see that for the low frequency (0.1rad/s) there is a good agreement between Nemoh and Capytaine, but at the frequency 2.0rad/s, both radiation and diffraction results differ significantly.

If you need more information or additional test please let me know.

supporting files

[mancellin]

Interesting, thank you. I can reproduce the discrepancy, but I've no idea where it could come from.

@H0R5E Since the only commit differing between mainline Nemoh and https://github.com/DTOcean/NEMOH is yours: can you confirm that the change is only meant to add a new output format but does not change the way the potential is computed internally?

[ffrancesco]

@mancellin the work on the cylinder potential was done originally by Cameron McNat and verified against WAMIT, the fork on the DTOcean is just there for installation purposes, so there should be no modification to the Fortram files.

[H0R5E]

@ffrancesco, is that you, Francesco? It's fair to say that the DTOcean version isn't vanilla NEMOH. There is an additional cylinder added in the computation (for calculation of interactions between devices). I think this is a relevant reference.

[ffrancesco]

@H0R5E yes it is Francesco (AAU). Yes, the version of Nemoh is a relatively old branch of Nemoh, in which the major modification to the code (in term of solver) can be found at line 117 of the SUBROUTINE SOLVE_BVP (SOLVE_BEM.f90)

Cylindrical surface computation
IF (Switch_Cyl.EQ.1) THEN
DO j=1,MeshCyl%Npoints
CALL COMPUTE_POTENTIAL_DOMAIN(ID,PhiC(j),MeshCyl%X(1,j),MeshCyl%X(2,j),MeshCyl%X(3,j),kwave,AMH,NEXP)
EtaC(j)=IIW/GPhiC(j)
END DO
CALL WRITE_SURF(ID,ProblemNumber,EtaC,MeshCyl,'cylsurface.',0)
END IF

what I wanted to say is that the branch DTOcean/NEMOH is the same copy of the code that has been verified against WAMIT.

@mancellin I can try to see if I can find some of the old WAMIT verification and post it here. Can that be helpful?
if not what else can I do to help?

[mancellin]

So, there is one more output format, but the result itself should be the same as in the reference stable version of Nemoh.

You can look for the WAMIT verification results, but whatever it says, an inconstancy of Capytaine with NEMOH is already a bug that deserves some investigation.

[ffrancesco]

@mancellin if there is a way I can help out please let me know.

[mancellin]

@ffrancesco Could you run the Nemoh/Capytaine comparison with some other meshes and parameters, to see if it affects the issue? For instance, what about infinite depth?

[ffrancesco]

@mancellin I cannot really run Nemoh with infinity water depth and estimate the potential on the outer cylinder, since the outer cylinder is automatically meshed and it requires the water depth to be finite. I run a case with 500m water depth and then compared with Capytaine with inf water depth. I also run a case with smaller outer cylinder (5.5m) and a case with a larger outer cylinder (10m). Additionally, I run a case with a surging cylinder. The low frequency solutions seem correct, but the high frequency solutions seem to differ significantly from Nemoh.
You can find all the material to run the different case inside the zipped file issue_191_followup.zip

[mancellin]

Other outputs seem to match between Nemoh and Capytaine. However the match seems better at low frequency (3 or 4 significant digits) than high frequencies (2 significant digits). Maybe this small mismatch is blowing up when computing the potential far from the body...

[ffrancesco]

@mancellin i ll try to estimate the potential at the mesh of the original body and see if there is already any error there. Cos if the original potential is slightly off at the mesh moving it out will not help. I ll should be able to add the material later today.

[ffrancesco]

@mancellin I ve run the Heaving cylinder case again, and compared the potential at the body mesh.
The results show some difference between the potential at the body, if I did not do any major mistake, this can justify the differences in the estimation of the potential at another mesh.
The attached files should show the difference between real and imaginary part of the potential for different cases.
I m not sure the presentation of the results is the best, but it should give you a quick glance at the differences.

PotentialATbody.zip

[mancellin]

@mancellin I ve run the Heaving cylinder case again, and compared the potential at the body mesh. The results show some difference between the potential at the body, if I did not do any major mistake, this can justify the differences in the estimation of the potential at another mesh. The attached files should show the difference between real and imaginary part of the potential for different cases. I m not sure the presentation of the results is the best, but it should give you a quick glance at the differences.

Thank you for this.
I tried to get some quantitative estimates of the mismatch using your last test.

for i in range(len(problems)):
    ...
    print(problems[i])
    def norm(x):
        return np.sum(np.abs(x))
    print(norm(PHI_faces.real - phi_faces.real)/norm(PHI_faces.real))
    print(norm(PHI_faces.imag - phi_faces.imag)/norm(PHI_faces.imag))

But I'm confused because the mismatch seems higher for low frequencies there.

RadiationProblem(body=None, omega=0.100, depth=inf, radiating_dof=Heave)
0.10436866652546821
9.16388154725368
[...]
RadiationProblem(body=None, omega=2.000, depth=inf, radiating_dof=Heave)
2.9106719430155306e-07
5.667434431347697e-06

[mancellin]

I was wondering if the discrepancy could come from floating point precision (Nemoh is single precision whereas Capytaine uses double precision by default).
But after a couple of tries with the new #224, the floating point precision does not seem to change the result of Capytaine on this problem.
That is probably not the issue.

[ffrancesco]

@mancellin, can it be a difference associated with the solver?
Capytaine uses GMRES while I m not sure if Nemoh has implemented it.
I tried to run capytaine with direct solver, but it was a quick check so I might had overlook at the results. I can try it again in the coming weeks.

[mancellin]

@ffrancesco The latest version of Capytaine uses a direct solver by default, so that's probably what you've already been using.

[mancellin]

@nhu-vnguyen Do you think you could look for an analytical solution for this problem (potential away from a floating vertical cylinder) and compare it with the same computation in Capytaine? It would help us figure out what is going on.

[jnvn7]

@mancellin yes, I have a few references on this problem. I can take a look at this and build a case for comparison.