Wall tests

[rykerfish]

Adding some tests for DPStokes (and maybe other solvers later) that compare to some reference results.

[rykerfish]

@RaulPPelaez This plot is to illustrate an issue I've been running into with DPStokes in libMobility. Here, the square points and the lines are reference data from the DPStokes repository, and the triangular points are from libMobility generated from the bottom wall test on this branch. In general, libMobility lines up well, but the inset shows a small error between libMobility and the reference that is consistently of order 1e-3. The error goes to zero as $z/R_h \to 0$. This same mismatch occurs in slit channel examples as well.

Brennan suspects this mismatch may be from suboptimal grid/kernel parameters. Do you agree? I'm going to attempt to integrate the parameter selection from DPStokes to see if that resolves the mismatch.

Let me know if have a different idea of what could be causing this systematic underestimation so I can check that out as well!

[RaulPPelaez]

I agree with Brennan, it is most probably suboptimal grid parameters.

[RaulPPelaez]

We should be merging this before #19

[rykerfish]

@RaulPPelaez I haven't merged this yet as there are still problems with the DPStokes results not agreeing perfectly with reference results and Brennan and I think it should be possible to get the results to match to more than a tolerance of ~1e-3. Even when using the parameter selection code, there remains a persistent underestimation on the order of 1e-4.

I have the reference code running from the DPStokes repository. As far as I can tell, all the grid and kernel parameters are the same between the two simulations. Do you know if there is anything that got changed between the libMobility version of DPStokes and the original version that could be causing this mismatch? We aren't sure where to look to try and resolve this as I've (mostly) ruled out the difference in parameters.

[RaulPPelaez]

1. Does this happen in all modes (bottom wall etc) consistently or is just one?
2. Try to run it in double just to discard that...
3. Can you post the error, (ref-libmo) or smth, vs z/R_h?

I cannot think of anything that could cause this given identical parameters.

[rykerfish]

Here's a relative error plot for the xx and zz trans-trans components of the mobility matrix, computed as $\displaystyle \big\vert \frac{ref-libm}{ref} \big\vert$, versus z/R_h with R_h=1. The dashed vertical line is z=1. The error is 0 at z=0 for both bottom and and slit channel geometries but small but non-zero elsewhere. There's a strange local minimum at z=1.27 for the slit channel geometry where the error goes down to $O(1e^{-6})$ that I can't explain.

The plot is the same for double and single precision, with negligible changes in the exact values of the computed mobilities.

Here, the values for $\beta$, $w$, $\alpha$, $N_x, N_y, N_z$, and $h$ are the same as given by the UAMMD printout for libMobility and the parameter selection printout for the original DPStokes code.

[RaulPPelaez]

Between the version from the reference and the current one used in libmobility there is this big commit in UAMMD https://github.com/RaulPPelaez/UAMMD/commit/c662443254bffefc4d3f23e74c9475fa79c07d80 You could try to undo this commit and try, maybe I messed up something subtle. The version in the reference (which I understood you have ran and confirmed the results are consistent with the paper) is this https://github.com/RaulPPelaez/UAMMD/commits/7f6119fa03986d3ef03c46533f75e5a69a81340a/ . You could also try taking UAMMD to the latest version in the reference repo and see if it still reproduces the results. Just to confirm/discard this difference is rooted in some change to UAMMD since the paper.

The reference uses w=6 or 12, but this PR has w=4, I reckon that is taken into account?

[rykerfish]

@RaulPPelaez finally think I figured out the rest of the issues. Turns out that all the reference data from the DPStokes repository was generated using the CPU implementations. Additionally, although I was using the reference data for the w=4 kernel, their script for the pair mobility had a bug that could force the wrong kernel to be used. I regenerated all the data in that repository using the GPU implementations and with correct kernels and ended up with the relative error plots below, which Brennan and I are satisfied with.

The pair mobility plot (2nd plot) isn't great since I didn't mess with the formatting enough, so a few notes on that one: The green curve with the error spikes is for two particles spaced 8R_h apart. The errors in that curve are consistently higher than the rest, but mostly small for z/R_h > 1. The large spike in the bottom right plot is because that height is at the exact middle of the slit channel, so the mobility gets small quickly. The actual error there is on the scale of 1e-9, but the relative error appears larger.

To get these small relative errors, I had to use the exact same parameters as the DPStokes repository code. To avoid doing that in the future, I regenerated all the reference data using our implementations with our current default parameters in single precision to use as the reference going forwards.

If all this looks good, I think this branch can finally get merged. Let me know what you think!

[RaulPPelaez]

I am not sure I am following. In these plots, you are comparing the GPU impl from the article with the GPU impl in libmobility?. What this tells us is that the GPU code has not changed behavior since the publication of the article. Which is nice, but it might still have an issue, right? Or you say you do not trust the CPU implementation? If you are comparing GPU with GPU then any error you see can probably be attributed to non deterministic execution order.

BTW one difference between CPU and GPU impls is that in the CPU particlesare spread to the cell corners in XY (pos_cell(i) = hi), but in GPU I spread to cell centers (pos_cell(i) = (i+0.5)h). This will give you differences in the results due to the variation of the hydrodynamic radius inside a cell. You should be able to accommodate for this by comparing CPU(particle_x, particle_y) with GPU(particle_x-0.5h, particle_y-0.5h).

[RaulPPelaez]

Also you should run all these correctness tests in double, numerical accuracy could be hiding some effects. For instance, the correction expressions are really bad for float as is and are written in a really specific way to reduce numerical issues. I bet the increased error you see near the wall is due to this.

[rykerfish]

@RaulPPelaez sorry about the delay in getting this finished- I was traveling for a bit, got busy with another project, and Brennan and I had spent some extra time talking about the DPStokes parameters til we felt comfortable with it.

What we landed on was leaving $h$ and $\beta$ at the values determined from the DPStokes paper as Brennan wanted to best preserve the hydrodynamic radius due to the high level of accuracy with which a particle's radius can be measured in experiments. A fixed value for $h$ requires altering the periodic length $L$ to make $N=L/h$ an integer for the FFT grid, which is now what is done in the libMobility parameter selection.

One consequence of this is that the reference results from the DPStokes repository could no longer be used exactly since the hydrodynamics are a bit sensitive to $L$, so I regenerated reference results from this repository with the new parameters.

If this method of selection seems good, I think this branch is finally good to go.

[RaulPPelaez]

All the contrary man, thank you so much for working on this.

Have you generated the reference using libMobility now then? That is dangerous bc it assumes what you are trying to test is correct beforehand.

[RaulPPelaez]

What we settled in the past with Aleks on this "you have to change r_h or L" is to have a parameter that allows you to select if L can be changed or not.

Note that by changing L you also change r_h due to periodic corrections.

[rykerfish]

I did regenerate the DPStokes references using libMobility, although I was able to use the old NBody references.

What I did to have confidence when regenerating the reference was to first match the parameters used to generate the original reference and see that I matched the original reference to around $10^{-6}$ or so. Then, I'd change the parameters back to the ones I implemented and regenerate. Obviously not 100% bulletproof, though.

If we want to have both a parameter selection that preserves $L$ and one that is allowed to modify $L$, I could finish the implementation in #19 that emulates the fancy parameter selection. That should allow us to use the old reference data too by using the same parameter selection code that they used to generate the original refs. Do you think that's a good route to go?

Good point about modifying $L$ still changing $r_h$. Maybe our hope was that the effects from the periodic corrections could be minimized if necessary by making $L$ larger? Or maybe that's just the sign that we'll have to settle for changing something that effects $r_h$ in the end, like you said.

[RaulPPelaez]

What you did makes sense, I was not getting the full picture.

My advise here is to keep it simple for now. Add a parameter called "allowChangingBoxSize" or smth like that. If set to true (you can make this the default if you decide) change L so you can fix h, otherwise keep h fixed to the nearest possible value (changing r_h).

If the latter is too troublesome then just make the module change L and document accordingly. We can add that functionality when needed in the future.

[rykerfish]

@RaulPPelaez I've followed your suggestion and implemented a mode that changes $h$ while keeping $L$ fixed (as the new default), but I did it with the inverse-polynomial interpolation used in the DPStokes repository to re-select $\beta$ according to the new $h$. Theoretically, this should better preserve the relationship between $r_h$, $h$, and $\beta$.

I added the parameter to select which mode to use- I did it in setParameters() since I thought it made sense to go with where you set $L$, but it could get moved to initialize() if you think that's a better place.

Finally, I regenerated the ref data for hopefully the last time with the new default. Once again, I compared against the data from the original DPStokes repository to verify a decent match before regenerating.

[rykerfish]

@RaulPPelaez I've followed your suggestion and implemented a mode that changes $h$ while keeping $L$ fixed (as the new default), but I did it with the inverse-polynomial interpolation used in the DPStokes repository to re-select $\beta$ according to the new $h$. Theoretically, this should better preserve the relationship between $r_h$, $h$, and $\beta$.

I added the parameter to select which mode to use- I did it in setParameters() since I thought it made sense to go with where you set $L$, but it could get moved to initialize() if you think that's a better place.

Finally, I regenerated the ref data for hopefully the last time with the new default. Once again, I compared against the data from the original DPStokes repository to verify a decent match before regenerating.

Edit: I also changed the threshold for the mobility matrix symmetric test to 1e-5 to make them all pass. I think this is okay since if I use the default parameters from the original DPStokes repository, they only pass this test at the 1e-5 level as well.

[RaulPPelaez]

Sounds good Ryker. Is this ready to merge?

[rykerfish]

@RaulPPelaez Ready to merge unless you find anything to change on review!

[rykerfish]

@RaulPPelaez This PR should be ready to merge- I think I addressed all your comments but let me know if there is anything else you notice!

I talked with Brennan and we're deferring implementing rectangular boxes (Lx != Ly) for DPStokes to a later date, so I put an error in there during initialization.

[RaulPPelaez]

I talked with Brennan and we're deferring implementing rectangular boxes (Lx != Ly) for DPStokes to a later date, so I put an error in there during initialization.

That's fine, we will get around it.