NBody wall kernels lose SPD property close to a wall

[rykerfish]

For this plot, I took a single particle at different heights approaching a bottom wall with DPStokes and NBody and applied only forces, then used the same procedure as in the fluctuation dissipation test to compute the full mobility matrix. This plot is the smallest eigenvalue of that 3x3 matrix. The plot shows that the NBody eigenvalues seem to become unbounded as the particle goes toward the wall.

I noticed this because there are times when the NBody wall kernels produce NaNs in the guesses while iterating in the Lancoz algorithm. This is configuration dependent and seemingly only depends on the distance between the wall and the closest particle.

I squinted at the NBody wall kernel for UF and it looked okay compared to the formulas in the Brady/Swan paper, but haven't investigated much further to see what's causing it.

[RaulPPelaez]

I wonder if this is a bug in the code or a limitation of the theory that's being used in NBody. What is the scale in X here? units of "a" I reckon?

[rykerfish]

Yup, the x-axis is in units of a!

I had the same question. The Swan/Brady paper declares that the corrections are SPD by construction but doesn't provide explicit proof. I wouldn't be shocked if this is a limitation of the theory since the Green's function for the wall is from Blake and uses the image system in the construction, but I don't know enough to say for certain.

If it's a bug in the code it's been around for a while since I believe these are tested against the old Python implementations from RigidMultiBlobs

[RaulPPelaez]

I am inclined to trust Swan/Brady on this... Might it be a precision issue? Does it happen also in double? Perhaps the implementation is mathematically correct, but accumulates too much error at some operation.

[RaulPPelaez]

BTW this is a cool test that could be included in the test suite. i.e. checking the smallest eigen value for particles throughout the domain.

[rykerfish]

This plot was made using double precision, and the results were visually identical between precisions. Accumulation of error is possible, but it starts happening for values of h that seemingly aren't that small numerically. Do you have any suggestions to look into? I'll keep thinking about it as well!

[RaulPPelaez]

I do not think it is accumulation of error then. Assuming there is an error, given that there is only one particle, I reckon it should be in this code: https://github.com/stochasticHydroTools/libMobility/blob/4dafab7f863842b5332e45cd3fdc78ce4cd4492a/solvers/NBody/extra/hydrodynamicKernels.cuh#L164-L173

[RaulPPelaez]

That looks identical to the expression in the article to me.

[RaulPPelaez]

Try to run it with a=1 and viscosity = 1/(6pi). The article states that correction formulas are normalized by 6pietaa^n, maybe the a^n part is being introduced incorrectly in this case?