Support for non square boxes in DPStokes

[RaulPPelaez]

I addressed a warning about non square boxes in DPStokes. I added a test to check all is fine.
I also added similar tests to UAMMD regarding this.

[github-actions]

✅ Linter reported no issues

All Python files are correctly formatted with Black.

[github-actions]

✅ Linter reported no issues

All C/C++ files are correctly formatted with clang-format.

[RaulPPelaez]

@rykerfish do you have the bandwidth to take a look at this to add some additional tests?

[rykerfish]

I made some time to look at it today. I tinkered around with using the reference datasets from the wall tests as another set of tests, but I don't think it's possible/a good idea since the reference domain size is roughly $80a$ and so periodic corrections still play enough of a roll to make it unreliable when modifying the domain to be non-square. Instead, I have the essence of some tests for isotropy that I'll finish up tomorrow.

@RaulPPelaez something of note- your test for test_non_square_box fails for me. I think it's because the spread/interp kernel in UAMMD is initialized only based on the x grid spacing, in the code here. The test fails when the domain is made non-square by extending in x. If I change the UAMMD code to use h as grid.cellSize.y inside the kernel in UAMMD, the test instead fails when the domain is made non-square by extending in y.

I tried to fix it here, but if anything I think I made it worse since now the test seems to randomly fail in both x and y. I may have set the z-radius incorrectly, though. Do you think it's correct to somehow modify the kernel parameters?

[rykerfish]

Two other things that might've caused my fix to go wrong:

1. Do we need to set a different beta for x/y based on the different h values? That would affect the radius
2. Similar question for alpha- I think I remember you saying a long time ago that alpha is in units of h (since we only set it as 0.5*w in libMobility). I couldn't figure out where this implicit multiplication by h was happening inside UAMMD to check if it was using the correct h values for each dimension.

[RaulPPelaez]

Two other things that might've caused my fix to go wrong:

1. Do we need to set a different beta for x/y based on the different h values? That would affect the radius

The way to do this is to choose beta, etc, using the smallest h. You get slightly more accuracy in one direction, but the kernel is not asymmetric.
This is how I imagine it in my mind (highly exaggerated hx/hy difference):

Note that the small h gets larger support in number of cells.

2. Similar question for alpha- I think I remember you saying a long time ago that alpha is in units of `h` (since we only set it as 0.5*w in libMobility). I couldn't figure out where this implicit multiplication by `h` was happening inside UAMMD to check if it was using the correct `h` values for each dimension.

I cannot remember -.-, the BM kernel is defined here: https://github.com/RaulPPelaez/UAMMD/blob/01cd9d175aadf296603989a1f4a3bbd65ccc0329/src/misc/IBM_kernels.cuh#L83-L113

your test for test_non_square_box fails for me. I think it's because the spread/interp kernel in UAMMD is initialized only based on the x grid spacing

That is strange that it fails for you and not for me. I have not run the test since I wrote it in this repo, though. But it was passing then. I can see things being weird when hx!=hy the way I left things.

[rykerfish]

I added a test to check that a particle responds the same when the domain is extended in x and y. It fails with the BM kernel as it currently is, but passes if you modify h in the BM kernel for x and y separately. I've done so here, and you can use this in the CMakeLists to test it:

FetchContent_Declare(
  uammd
  GIT_REPOSITORY https://github.com/rykerfish/UAMMD
  GIT_TAG        origin/dpstokes_kernel_params
  EXCLUDE_FROM_ALL
)

I'll need to modify and test the changes in UAMMD for torques as well if this is the change we want.

[RaulPPelaez]

That makes total sense, and your implementation looks correct. Let's merge RaulPPelaez/UAMMD#62
I will make a release, which you can then use in this PR. Thanks @rykerfish.

[rykerfish]

I have a somewhat insane-sounding bug.

With the test how you left it, I will get random test failures. More particles makes it more likely to fail, but it will pass sometimes. The max absolute error is usually between 1e-2 and 1e-3.

I've tracked it down to it being related to the number of grid points. Since we double the domain we'd expect double the number of grid points, but this is not always the case. With the parameters as-is (which uses h=a/1.205), we get the long edge having 38 grid points and the short edge having 20 grid points. If I fiddle with the initialization and set h=a/1.25 then we get exactly (40,20) for the number of grid points along the long and short sides. Running the test with this causes max absolute error to be substantially better, between 1e-6 and 1e-9. This behavior seems to stay consistent even when the number of grid points is doubled.

Is this something you're concerned by? It may just be what we get from the FFT grids and, in the bad case, is still around the error threshold that DPStokes is designed to provide. This explains the random failures I was experiencing.

I'm running the UAMMD tests as we speak and will let you know on the other PR once they finish!

[RaulPPelaez]

When DPStokes chooses the number of grid points, it takes into account FFT-friendlyness. That means that doubling the size in one direction can result in a different cell size, h.

You get different results perhaps because h is related to the width of the kernel and then to the hydrodynamic radius, so you are comparing systems with slightly differently sized particles.

Now, is that change in hydrodynamic radius ok? I just cannot give you an answer -.-
Take into account this is just my guess of what is going on. I reckon the "large" 1e-2 errors you are getting are related to the guarantee of these funny beta parameters on accuracy.

If I am correct, if you increase the accuracy (by using a larger support, etc) the effect would vanish.

[rykerfish]

[From slack] I think we came up with an okay test to answer the question you posed on git (since we were thinking about as well) for how much the different grid spacing affects the hydrodynamic radius. I set up a square domain and manually tinkered with the number of grid points to make the grid spacing uneven in x and y. this first plot is comparing with RPY after pushing a particle in each direction (e.g. push a particle in x and then look at the x component of velocity). this plot is done using the original procedure for the betas where we pick the smallest and use that for the kernel

and then here's the same plot but where we use a different beta for x and y depending on the grid spacing in each direction

the x-axis is h/a so I didn't care about the disagreement below 1, but it seems like using the beta that corresponds to the h in each direction helps a lot. to turn it into a proper test, I think we can set up a non-square domain and compare it to RPY and make the domain large enough that we're around a 1% error when comparing to RPY for a couple different non-square sizes

[rykerfish]

@RaulPPelaez I think this is done with two caveats. First is that we'll need to merge this PR into UAMMD and bump the version here. Second is that I kinda trashed the initialization of DPStokes while working in the beta initialization. I'm happy to rewrite it to be cleaner, just didn't get to it today and wanted to put the working code online for you to look at.

Take a look at the tests in tests_dpstokes.py. I think they're pretty thorough.

[RaulPPelaez]

You should be able to set v3.0.0 as the UAMMD version now.

[rykerfish]

Looks good to me aside from the comments I made. I'm inclined to create issues for anything outstanding that doesn't have to do with non-square domains and handle them in new PRs.