[SPH][Setup] Stretch mapping implementation

[DavidFang03]

The goal is to have a way to initialize a distribution of particles that follows a given density profile $\rho(r)$. We follow the stretch mapping method in Price 2018 (3.2): starting from a lattice (HCP here), each position is modified using a Newton-Raphson method, or a bisection method if necessary. In order to generate an spherical object at the end, the HCP lattice is initially cropped into a sphere and is stretched afterwards.

I've extended the HCP generator:
GeneratorLatticeHCP.hpp -> GeneratorLatticeHCP_smap_sphere.hpp
crystalLattice.hpp -> crystalLattice_smap_sphere.hpp

Usage:

setup = model.get_setup()

center = [0, 0, 0]
xmax = 1
# gen = setup.make_generator_lattice_hcp(dr, bmin, bmax)
# but instead:
gen = setup.make_generator_lattice_hcp_smap_sphere(dr, center, xmax, rhoprofile)

setup.apply_setup(gen)

The given profile has to be a function of $r$ only and needs to be defined at $r=0$ because there's always one particle at the lattice center.
Tested with

rhoprofile = lambda r:1/(r+0.1)**2
rhoprofile = np.sinc # but very slow here (because the profile is less steep?)

3 remaining issues:

• Before the first timestep, the density seems to match the given one correctly (because hpart is initialized with the expected value during the generation), but then the next computed smoothing length differs close to the origin and the corresponding density no longer matches there 1r2.pdf
• There seem to be more particles at large radii, even though the density profile decreases with radius.
• Very slow (especially with flatter profiles?)

[tdavidcl]

Very nice to see, great job !
Yona is starting to review this one so i won't interfere. I just had a few comments.
Also from what i see this is only in the spherical case right ? could we generalize it using a lambda to get the position and a lambda with the metric terms (1, 2pi r, 4 pi r^2, ...) such that we can apply it for any axis ?

[DavidFang03]

Fixes done!

Very nice to see, great job ! Yona is starting to review this one so i won"t interfere. I just had a few comments. Also from what i see this is only in the spherical case right ? could we generalize it using a lambda to get the position and a lambda with the metric terms (1, 2pi r, 4 pi r^2, ...) such that we can apply it for any axis ?

Do you mean something like

gen = setup.make_generator_lattice_hcp_smap(dr, center, xmax, rhoprofile, system="spherical", coord="r") # spherical along r
gen = setup.make_generator_lattice_hcp_smap(dr, center, xmax, rhoprofile, system="cart", coord="y") # cartesian along y

with system="cart", "spherical", cylindrical" and coord="x","y","z","r"?

And then for non-isotropic objects, maybe we can allow something like

gen = setup.make_generator_lattice_hcp_smap(dr, center, [xmax, ymax, zmax], [rhoprofile_x, rhoprofile_y], system="cart", coord=["x", "y"]) 
gen = setup.make_generator_lattice_hcp_smap(dr, center, [xmax, ymax, zmax], [rhoprofile_r, rhoprofile_z], system="cylindrical", coord=["r", "z"]) # for a disc ?

[tdavidcl]

Fixes done!

Very nice to see, great job ! Yona is starting to review this one so i won"t interfere. I just had a few comments. Also from what i see this is only in the spherical case right ? could we generalize it using a lambda to get the position and a lambda with the metric terms (1, 2pi r, 4 pi r^2, ...) such that we can apply it for any axis ?

Do you mean something like

gen = setup.make_generator_lattice_hcp_smap(dr, center, xmax, rhoprofile, system="spherical", coord="r") # spherical along r
gen = setup.make_generator_lattice_hcp_smap(dr, center, xmax, rhoprofile, system="cart", coord="y") # cartesian along y

with system="cart", "spherical", cylindrical" and coord="x","y","z","r"?

Yeah that could be the we want it to look like in the end. For now i was thinking about the stretchindiv function where one could supply rhoprofile as well as S in such a way that it becomes applicable to any direction (regardless of geometry).

Something like (with r renamed to a to mean that it is anything, x,y,z,r,theta, ...)

static Tscal stretchindiv(
            Tscal a,
            std::function<Tscal(Tscal)> rhoprofile, //rho(a)
            std::function<Tscal(Tscal)> S,                 //surface element (a)
            Tscal integral_profile,
            Tscal rmin,
            Tscal rmax,
            Tvec center) {
....
            auto rhodS = [&rhoprofile, &S](Tscal a) -> Tscal {
                return S(a)*rhoprofile(a);
            };

If i'm right it is sufficient to generalize it in such a way that can supply anything (even theta/phi streching).

Also maybe we can add a std::function<Tscal(Tvec)> a_from_pos, std::function<Tvec(Tscal, Tvec)> a_to_pos that perform the conversion supplied by the user. Then in the generator you could do

        Tscal a = a_from_pos(r_a);
        ....
        Tscal new_a = stretchindiv(a, rhoprofile, S, integral_profile, rmin, rmax, center);
        ...
        Tvec new_r_a = a_to_pos(a, r_a);

[tdavidcl]

also for authorship related error in CI you can do python3 ../buildbot/update_authors.py from your build directory

[y-lapeyre]

Nicely done :)
Could you generalize this setup by turning it into a modifier that takes whatever lattice generator in input and applies the stretch transformation? You can name it something like ModifierApplyStretchMapping for instance. Inspiration can be found in ModifierOffset.hpp / cpp for instance.

[tdavidcl]

yeah you can expect the fail with kernel_call since a std::function can not be mapped in the GPU kernel (it is a CPU function pointer which can not exist in the GPU world :sad:)

[tdavidcl]

if you update the branch the mail thing should be fixed

[DavidFang03]

yeah you can expect the fail with kernel_call since a std::function can not be mapped in the GPU kernel (it is a CPU function pointer which can not exist in the GPU world :sad:)

Okay, that's what I understood, it's a good lesson

[DavidFang03]

I'm currently trying to migrate ModifierApplyStretchMapping into SolverGraph, I think I'll need a little feedback once I manage to compile something haha

[tdavidcl]

honestly as you prefer for this one, the solvergraph is really made for the solver itself so it is quite quirky when it comes to the setups. I haven't really started integrating it in that part of the code so don't worry if it is not used it there

[tdavidcl]

i don't see any solvergraph related thing in the current state of the PR have you committed it or did you meant the setup graph ?

[DavidFang03]

No I haven't committed anything yet.
What I'm trying to do is write a node (src/shammodels/sph/src/modules/StretchMap.cpp) which takes positions as edges and I would write the stretch mapping process inside. It would also determine the particles that should be killed if they're outside the integration boundaries (typically it would be equivalent to crop a sphere before stretch mapping, in the spherical case), similarly to GetParticlesOutsideSphere
Then in the Modifier, the sequence would be stretchmap,killparticles,set_hpart or sth like that. But maybe it's not that useful, what do you think? If I want to kill these particles I'll still need to write sth with solvergraph though (or not?).

[tdavidcl]

as you wish I'm ok with both options take whichever you prefer (or take less time).
Also could you add a example in the doc on how to use the strechmapping, we can do a call if you need any help on that side.

[tdavidcl]

also i've invited you on the repo to skip the workflow approval, they will start right away

[DavidFang03]

Okay, I'll leave the SolverGraph aside for the moment I think. I'll do it if I have some time at the end of my project. For the doc I can do a call today at 5:30 pm, or Monday if that's too late.
Thanks for the invitation :)

[github-actions]

Workflow report

workflow report corresponding to commit 494dcad
Commiter email is [email protected]

Pre-commit check report

Some failures were detected in base source checks checks.
Check the On PR / Linting / Base source checks (pull_request) job in the tests for more detailled output

❌ Authorship update required

The following files had their author headers updated by the author update script.

Please run the script again (python3 buildbot/update_authors.py) and commit these changes.

Note: The list below is only partial. Only the first 10 files are shown.

• /src/shammath/include/shammath/integrator.hpp

Suggested changes

Detailed changes :

diff --git a/src/shammath/include/shammath/integrator.hpp b/src/shammath/include/shammath/integrator.hpp
index 9f4adebf..dc7a6bcb 100644
--- a/src/shammath/include/shammath/integrator.hpp
+++ b/src/shammath/include/shammath/integrator.hpp
@@ -11,6 +11,7 @@
 
 /**
  * @file integrator.hpp
+ * @author David Fang ([email protected])
  * @author Timothée David--Cléris ([email protected])
  * @brief
  *