Since then, I have tried to validate my attempts on the proposed formulation without any success. The problem now strikes me again as I am finishing my acoustical engineering undergrad work. I have seen the example notebook on BEM-BEM coupling was recently added to github, so I figured I should give it a try again.

I am trying to model a porous sample with complex speed of sound and density and the boundary data of an incident wave. I understand that the speed of sound is included in the wavenumber that goes in the multitrace_operator. But I’m confused of where the interior’s medium density should go. According to the cited paper, the density ratio shall scale the interior operator, but all my attempts to do so have resulted in incorrect results.

Does anyone have any idea of how I should correctly scale the operator and proceed to solve this problem? I don’t think posting my code will be useful since I’m basing it all on the provided notebook example.

I hope my problem was made clear, and I figured it was best to move the discussion to the new platform.

The speed of sound goes into the wavenumber as k=2pi*f/c with f the frequency and c the speed of sound. The density goes into the Neumann boundary condition. This condition models continuity of particle velocity and typically is given by the normal derivative of the pressure divided by the density.
You need to scale the interior Calderón operators with density ratios as in Eq. (7) of the Arxiv paper. Specifically, scale the interior hypersingular operator with rho_ext/rho_int and the interior single-layer operator with rho_int/rho_ext.

I have attempted said solution but my results are still not matching my reference, which have been validated in two different softwares, I reckon I’m still doing something wrong.

I am quite positive that the Neumann boundary condition is correct, as the only step taken is dividing it by the density. So I assume my scaling procedure is incorrect. I am using the following lines to scale:

Please remember that there is a scaling in the surface potentials as well. See Eq. (6b) in the Arxiv paper. The solution vector consists of the Dirichlet and Neumann trace of the exterior field. However, the interior field is given by the potential operators over the traces of the interior field, see Eq. (A.1) for the representation formula.
Hence, when calculating the potential operators for visualisation of the interior field, you need to use the interior Neumann trace of the field. With the interface condition, this is (rho_int/rho_ext)*sol_neu.
Please check if you have this scaling correctly implemented.

Thanks for the response, I indeed was not doing this scaling related to the surface potentials for the interior field.

Nonetheless, the problem I am attempting to solve is only concerned with the exterior field, and therefore my validation reference only takes into account field points in the exterior field. So I assume for the present case this shouldn’t be a problem.

Well, I will continue to seek for possible mistakes in my code, if you think there are any possible areas in the code that are prone to mistakes and I should pay special attention, let me know.

A question just occurred, when setting the Neumann boundary condition, should I divide it by the density ratio, or just by rho_ext? When comparing to my reference validation, I have noticed that the results matched way better when considering a division by rho_ext. While, when considering the density ratio, the mismatch was big. I am confused if Eq 6b refers to this Neumman boundary condition scaling process too, or only the surface potential scaling you mentioned in your last reply.

Furthermore, are there any known wavenumber dependent limitations to the method, apart from the slower convergence at higher wavenumbers? I am dealing with problems in the human hearing frequency range.

For this BEM formulation you need to pick either the interior or exterior trace as unknown surface potential. In the Arxiv paper, the exterior trace of the total field was chosen, hence you need to scale with the density ratio to calculate the interior trace.

Based on the information provided, it is difficult to assess what is going wrong in your implementation. I’d recommend checking all parts of your implementation. Check the number of elements per wavelength and increase the mesh density if possible. Check the convergence of the linear solver and reduce the tolerance if necessary. Check what happens at different frequencies and wavespeed/density ratios. Check if your mesh has high quality and all normals pointing outwards.