[SPH][Setup] Stretch mapping implementation #1372
Draft
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
y-lapeyre
changed the title
tdavidcl
reviewed
Nov 22, 2025
Member
tdavidcl
left a comment
tdavidcl
left a comment
There was a problem hiding this comment.
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 ?
src/shammodels/sph/include/shammodels/sph/modules/setup/GeneratorLatticeHCP_smap_sphere.hpp
Outdated
Show resolved
Hide resolved
Collaborator
Author
|
Fixes done!
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 ywith 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 ? |
Member
Yeah that could be the we want it to look like in the end. For now i was thinking about the 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 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); |
Member
|
also for authorship related error in CI you can do |
y-lapeyre
requested changes
Nov 23, 2025
Collaborator
y-lapeyre
left a comment
y-lapeyre
left a comment
There was a problem hiding this comment.
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.
Member
|
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:) |
Member
|
if you update the branch the mail thing should be fixed |
Collaborator
Author
Okay, that's what I understood, it's a good lesson |
Collaborator
Author
|
I'm currently trying to migrate ModifierApplyStretchMapping into SolverGraph, I think I'll need a little feedback once I manage to compile something haha |
Member
|
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 |
Member
|
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 ? |
Collaborator
Author
|
No I haven't committed anything yet. |
Member
|
as you wish I'm ok with both options take whichever you prefer (or take less time). |
Member
|
also i've invited you on the repo to skip the workflow approval, they will start right away |
Collaborator
Author
|
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. |
…th center and radius parameters instead of box boundaries + rename
…d object is coherent with the required geometry
…ecause the final nb of particles can be difficult to estimate in advance
Collaborator
Author
|
TODO: Currently it runs on one CPU only |
Workflow reportworkflow report corresponding to commit 8cba392 Light CI is enabled. This will only run the basic tests and not the full tests. Pre-commit check reportPre-commit check: ✅ Test pipeline can run. Doxygen diff with
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
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.hppcrystalLattice.hpp->crystalLattice_smap_sphere.hppUsage:
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
3 remaining issues:
hpartis 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