I have problems in the coupling of OpenFOAM and FEniCSx with respect to the RBF interpolation in the FEniCSx adapter in python.
The interpolation from the OpenFOAM to the FEniCSx grid works well, however in the internal interpolation in the adapter I get the following warning
C:\Users\enders-seidlitz\Anaconda3\lib\site-packages\scipy\interpolate\rbf.py:257: LinAlgWarning: Ill-conditioned matrix (rcond=9.16567e-22): result may not be accurate.
self.nodes = linalg.solve(self.A, self.di)
and the values written to my coupling_expression contain large errors.
I could find the following problem in the RBF interpolation on my coordinates:
from scipy.interpolate import Rbf
import numpy as np
coords= np.loadtxt("interpolation_coords.txt")
data = np.random.rand(coords.shape[1])
rbf = Rbf(coords[0, :], coords[1, :], data)
print("Max Rbf interpolation error:", max(abs((rbf(coords[0, :], coords[1, :]) - data)/data)))
This code yields, e.g., “Max Rbf interpolation error: 46208.524712400205”
I first thought there would be duplicate nodes in interpolation_coords.txt, but this looks fine (I checked this by comparing np.unique(coords,axis=1).shape to coords.shape). Can you try to make the mesh a bit more “friendly” - e.g. by removing some nodes at the left part of the domain? Maybe this resolves the conditioning problems indicated by the warning above.
Changing the mesh is obviously only a work-around and preCICE’s RBF interpolation seems to be able to handle the mesh. So there should also be a way to get the RBFs in the FEniCS adapter working as expected. As far as I know, preCICE uses some advanced techniques for RBF interpolation (I cannot give details here, because I’m not really a RBF person). In the FEniCS adapter I only implemented the bare minimum (so far only very “friendly” meshes tested).
Maybe somebody with a better RBF background can help out here? Is the observed error something one would expect for RBFs without special treatment and the given mesh?
It appears that the matrix is badly conditioned with the default multiquadratics basis function.
Thin-plate: Max Rbf interpolation error: 7.032982565798581e-07
Default: Ill-conditioned matrix (rcond=1.4045e-21): result may not be accurate.
self.nodes = linalg.solve(self.A, self.di)
Max Rbf interpolation error: 39347.20266935741
Great to hear that you got it working and thanks @KyleDavisSA for lending us your expertise
Do you see any possibility to integrate that in the adapter?
@arvedes: Absolutely! I did not spend much of thought on how to perform the interpolation the best way. I just picked the standard configuration from scipy when developing the adapter and tested it with few test cases (with friendly meshes), where it worked well. Feel free to open an issue in the adapter. We can discuss using function=thin-plate as default there (or make it configurable, at least). With your case there is now at least one good use-case for this.
Note: The adapter allows to define different interpolation schemes and is extensible in this direction, because we already expected that a single interpolation scheme might not be sufficient at some point in time.
@arvedes previous work on RBF’s have shown that thin-plate works quite well, however it is a global basis function and not really suitable for large meshes (probably somewhere less than 10^4 vertices). I think the default basis function is not as stable. Another option is Compact Polynomials of Gaussian basis functions, but then this introduces the shape parameter as another input variable, hence me choosing the Thin-Plate basis function instead.
