-
Notifications
You must be signed in to change notification settings - Fork 35
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Particle masses bugfix #89
Particle masses bugfix #89
Conversation
…articleDescr classes
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
someone who deals with particles algo needs to look at the comment by Stefan in file enzo_EnzoMethodPmDeposit.cpp . I have also left a comment there flagging the issue.
Thanks for pointing that out. That comment was due to a misunderstanding on my part. I'll delete it now. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
src/Cello/data_ParticleDescr.cpp
Outdated
@@ -253,6 +253,24 @@ int ParticleDescr::attribute_index (int it, std::string attribute_name) const | |||
|
|||
//---------------------------------------------------------------------- | |||
|
|||
bool ParticleDescr::has_attribute(int it, std::string attribute_name) const | |||
{ | |||
ASSERT1("ParticleDescr::is_attribute", |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Mismatch, but I'd just go with is_attribute
to be consistent with #49
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I named has_attribute as such to maintain consistency with the has_constant method. I have now decided to name the methods 'is_attribute' and 'is_constant' instead.
@stefanarridge it looks like this PR has some conflicts, as well as some comments. Would mind fixing the conflicts, addressing the comments, and then letting the reviewers know it's ready for review? Thank you! |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All looks good! I only had a minor comment that isn't related to the density -> mass changes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This looks great! I've got a couple of minor suggestions. I'm a little concerned that one of your changes might have broken gas density deposition in cosmological simulations (there is a comment highlighting that in enzo_EnzoMethodPmDeposit
). But I could definitely be wrong.
Apologies for not providing these suggestions sooner. I wasn't originally intending to fully review this because I didn't really understand the Particle machinery (but James's presentation at the workshop made things a lot clearer).
if (rank >= 1) hx *= cosmo_a; | ||
if (rank >= 2) hy *= cosmo_a; | ||
if (rank >= 3) hz *= cosmo_a; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It seems like that removing this is going to affect the way that the gas density gets deposited (which happens at the end of this function).
If this is true, then it's a little disturbing that none of our automated tests seem to detect this.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Okay, I'll look into this. I was under the impression that that both gas density and particle density were comoving quantities. I suppose if a physical length scale is required for depositing the gas density, then I can define two separate variables i.e hx_com and hy_phys, or something similar.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was just trying to point out that this change was going to modify the deposititon of the gas density.
You could definitely be correct (I know basically nothing about the units of cosmological simulations). In that case the old implementation must have been wrong.
Could somebody familiar with the density units used for cosmological simulations in Enzo quickly clarify this for us? @jwise77 @bwoshea @gregbryan
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I know what you mean. It would be good to understand which quantities need to be passed to dep_grid_cic (i.e., physical or comoving). I am also concerned that the value of dtf which is passed is not correct. Currently it is set to alpha_, which is a dimensionless parameter which is 0.5 by default, rather than a quantity with dimensions of time.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In Enzo, the particle "mass" is actually a comoving density, adjusted for the level that it is on, so that when it gets deposited to the GravitatingMassField (which is the comoving density field used to calculate the gravitational acceleration), there is no multiplication by a volume. Here, I'm not quite sure what should be done -- the mass of a particle is a constant, so if you wanted to compute a comoving density, then you would want to divide by the comoving volume of the cell. Normally dx is comoving, so a factor of the scale factor is not required, but I have to admit I have not really looked at what is going on in this routine, so best to take what I say with a grain of salt.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@gregbryan , the purpose of this PR is to redefine the "mass" attribute, so that it is always treated as mass, rather than (sometimes) as density. As far as I know, in a cosmological simulation, length scales and densities are all comoving quantities, which is why I deleted the lines which multiply hx,hy,hz by "cosmo_a", so that, when the mass is multiplied by 1/(hx * hy * hz), the output is a comoving density.
The issue which @mabruzzo pointed out, is that this will affect the deposition of gas density, which is done by calling the 'dep_grid_cic' function:
`
FORTRAN_NAME(dep_grid_cic)(de,de_gas,temp,
vx, vy, vz,
&dtf, rfield, &rank,
&hxf,&hyf,&hzf,
&mx,&my,&mz,
&gxi,&gyi,&gzi,
&nxi,&nyi,&nzi,
&i0,&i0,&i0,
&nx,&ny,&nz,
&i1,&i1,&i1);`
Should hxf, hyf, hzf be comoving or physical length scales?
I'm also surprised that dtf is set equal to 'alpha_'. Should it not be the same as the 'dt' used to drift the particles, i.e.:
const double dt = alpha_ * block->dt() / cosmo_a;
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I took a look at dep_grid_cic
, and I think that it's effectively just advancing the density over dt
(i.e. the analogous operation to drifting the particles). I'm also pretty sure that the output density de_gas
has the same units as the input de
- As an aside, I'm pretty sure that a lot of these arguments that get passed have to do with complications that arise from the patch-based AMR used in the original Enzo.
If I'm remembering correctly, the fields used to store gas properties store values in comoving units for cosmological simulations, so hxf
, hyf
, and hzf
should probably be comoving.
Furthermore, since the "density_total"
field is the sum of the "density"
field (after drifting) and the "density_particle"
field, I suspect that "density_total"
, "density_particle"
and "density_particle_accumulate"
should all be comoving densities (in cosmological simulations).
I could be wrong, but this makes the most sense to me. The easiest way to check this would probably be to run some tests and compare results against Enzo.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Without doing a deep-dive back into this, it's unclear to me if you've made any changes on this front. Lmk if you have made changes, and I'll take another look
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm also surprised that dtf is set equal to 'alpha_'. Should it not be the same as the 'dt' used to drift the particles, i.e.:
const double dt = alpha_ * block->dt() / cosmo_a;
I agree that it should be the same dt as the particle drift, but in Enzo, dt is actually zero for the PPM solver but not the RK2 solver (diff from 2011).
I looked at a very old version of Enzo from 2005, and I found this comment that using a non-zero dt creates an instability. I think we should change it to zero to be consistent with Enzo.
float dt = (DepositTime - Time)/a;
dt = 0; /* the previous line is correct but creates an instability */
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, I looked this over and I think removing the multiplication by cosmo_a is the right thing to do. As others have pointed out, hx,hy,hz should be comoving, so this is not necessary. Regarding alpha_, this is for the particles, and I think that is only set to zero for the baryons, so OK to keep as is?
…prints the density rather than inv_vol
…thodPmDeposit and added comment to improve readability
…with constant mass, rather than separately for each batch
…field is now the same as that used to advance the particles
…n branch. This is likely not correct but this was done for purposes of consistency and to ensure that this branch only changes the treatment of particle masses.
…ented successful compilation.
Merged with main
… simplified the code a little bit
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
All looks good, but I will comment on the dtf = _alpha
factor in a reply to John R's comment.
Sorry for the late comments here. Firstly, I don't foresee any issues with the yt frontend. The existence or absence of the "is_gravitating" attribute makes for an easy check for whether the data was made before and after this change. I apologize if this has already been discussed and if so I'm happy to leave it. Supporting both density and mass particles seems unnecessary. Is there any benefit to defining particle masses as densities beyond skipping a multiplication/division? It seems like we are making things more complicated for ourselves by continuing to support both options. |
There is in face little to no computational benefit in using densities rather than masses: in the "density" case, the densities have to be multiplied by 2^(rank * level) (rank is the number of dimensions, level is the refinement level), while in the "mass" case, the masses have to be multiplied by the inverse cell volume. In principle, I agree with you, it would be better not to have the "density" option, but I am allowing for this option in order to make it easier for users to convert their old input parameter files or IC files. For example, if you look at the changes made to |
I think this is an opportunity for us to rip off a relatively painless bandaid. Since the density option still requires a renormalization to the local grid cell size, this doesn't even work the way it does in enzo-dev, so this could be a potential point of confusion for analysis or if more code is ported from that codebase. My personal opinion is that leaving the density option in leaves the door open for bugs, but I would like to hear other opinions. @bwoshea, @jwise77, @gregbryan, @fearmayo, @peeples, @mabruzzo ? |
I don't much (any) experience using particles, but FWIW, I'm all for ripping of this band-aid. It doesn't seem like there's much downside |
if (rank >= 1) hx *= cosmo_a; | ||
if (rank >= 2) hy *= cosmo_a; | ||
if (rank >= 3) hz *= cosmo_a; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
OK, I looked this over and I think removing the multiplication by cosmo_a is the right thing to do. As others have pointed out, hx,hy,hz should be comoving, so this is not necessary. Regarding alpha_, this is for the particles, and I think that is only set to zero for the baryons, so OK to keep as is?
|
||
enzo_float hxf = hx * cosmo_a; | ||
enzo_float hyf = hy * cosmo_a; | ||
enzo_float hzf = hz * cosmo_a; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think this cosmo_a factors should not be here, for the same reason discussed earlier.
|
||
enzo_float hxf = hx * cosmo_a; | ||
enzo_float hyf = hy * cosmo_a; | ||
enzo_float hzf = hz * cosmo_a; | ||
enzo_float dtf = alpha_; |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Definitely agree this should be alpha_ * dt.
I also think it may be simpler just to only allow mass. I always regretted the Enzo decision to make it density (the level of confusion added was not worth the small gains). I do know that you have already invested quite a bit of effort, so would completely understand punting on that for now. I think the bigger issue is that we need a test problem to verify the choices made. For example, there is a discussion about using hx or hx * cosmo_a. I think we have arrived at the right thing, but there is only one answer that is right and it would be great to be able to show it. For cosmology, one could test the growth rate of the Zeldovich pancake or the power spectrum for a cosmo run. For non-cosmology, maybe just the growth timescale for the collapse problem? |
It won't be much work to get rid of the "density" option. However, users will have to update parameter files / IC files which set a value for particle "mass", although it should be as simple as just multiplying by the root level cell volume.
I have been considering making a separate PR which will address the "gas density deposition" issue, as well as a few other issues with the gravity method which I think need to be addressed. The issue I have is that I am not completely certain what the
The way that I think about it is that it is essentially depositing mass in the same way as is done for the particles, only this time, the "particle positions" are the cell centers, the "particle velocities" are the values of the velocity field, and the "particle masses" are the cell masses (density times cell volume). If this is indeed how its done, it would be relatively straightforward to write this down in C++ code in `src/Enzo/enzo_EnzoMethodPmDeposit.cpp", which would simplify things a lot. Obviously, how to test whether its "correct" is a separate question which I'm not sure how to answer. |
34d8923
to
dc5bbd1
Compare
…"Particle:mass_is_mass" parameter, which will nearly always have the value "dummy", but there is a small chance it will have a different value. The purpose of this parameter is to be a flag for the Enzo-E frontend of yt, indicating that particle masses should be treated as quantities with dimensions of mass, rather than quantities with dimensions of density. It is only important that this parameter exists, it does not matter what value it has.
dc5bbd1
to
af6e42f
Compare
…e a "density" attribute / constant.
…ave a "mass" constant rather than a "density" constant, and the value of the "mass" constant is the value of the previous "density" constant multiplied by the root level cell volume.
Update: I have now made the following changes:
Again, I suggest this doesn't get merged in until the yt Enzo-E frontend is updated, which @brittonsmith is working on. I welcome comments and questions in the meantime. |
Do we need the "mass_is_mass" attribute? I guess we use it for now and then silently end support for it in due course? |
This is what I'm using on the yt side to detect whether the data was made before of after this change, which determines whether the density->mass conversion happens. The yt PR with the relevant changes can be found here. |
Previously the drift time for the gas density was set to the value of `alpha` (taken from `Method:pm_deposit:alpha`). As @stefanarridge pointed out in PR enzo-project#89 this didn't make a lot of sense. I also ran into issues with introducing a Stable Jeans Wave regression test in PR enzo-project#186 that required me to set `Method:pm_deposit:alpha` to 0. For comparison, the drift time for the gravitating mass particles is set to `alpha * dt / cosmo_a`, in which: - `alpha` again comes from `Method:pm_deposit:alpha` - `dt` is the time-step during the current cycle - `cosmo_a` is the scale factor computed at `time + alpha * dt` (where `time` is the time at the start of the current cycle) To investigate this issue I checked the original code from enzo-dev in `Grid_DepositBaryons.C`. I concluded that enzo-dev sets the gas density's drift-timestep should get set in exactly the same way as the Particle's drift timestep. However, when using the PPM and Zeus Hydro-solvers, the drift timestep is set to 0 (this is consistent with a comment @johnwise made during his review of PR enzo-project#189). Since our 2 primary hydro-solvers in Enzo-E (Ppm and VL+CT) currently handle these source terms with the same temporal order as these solvers, we now force the gas-deposition drift timestep to be 0 in `EnzoMethodPmDeposit`. We will need to revisit this exact behavior when we eventually modify the VL+CT solver to handle self-gravity in a higher temporal-order manner.
…posit Previously the drift time for the gas density was set to the value of `alpha` (taken from `Method:pm_deposit:alpha`). As @stefanarridge pointed out in PR enzo-project#89 this didn't make a lot of sense. I also ran into issues with introducing a Stable Jeans Wave regression test in PR enzo-project#186 that required me to set `Method:pm_deposit:alpha` to 0. For comparison, the drift time for the gravitating mass particles is set to `alpha * dt / cosmo_a`, in which: - `alpha` again comes from `Method:pm_deposit:alpha` - `dt` is the time-step during the current cycle - `cosmo_a` is the scale factor computed at `time + alpha * dt` (where `time` is the time at the start of the current cycle) To investigate this issue I checked the original code from enzo-dev in `Grid_DepositBaryons.C`. I concluded that enzo-dev sets the gas density's drift-timestep should get set in exactly the same way as the Particle's drift timestep. However, when using the PPM and Zeus Hydro-solvers, the drift timestep is set to 0 (this is consistent with a comment @johnwise made during his review of PR enzo-project#189). Since our 2 primary hydro-solvers in Enzo-E (Ppm and VL+CT) currently handle these source terms with the same temporal order as these solvers, we now force the gas-deposition drift timestep to be 0 in `EnzoMethodPmDeposit`. We will need to revisit this exact behavior when we eventually modify the VL+CT solver to handle self-gravity in a higher temporal-order manner.
This changes
EnzoMethodPmDeposit
so that that all particles in the 'has_mass' group are taken into account when calculating the mass field, and the particles' mass attribute/constant is always interpreted as a quantity with dimensions of mass, rather than sometimes density. This fixes issue #83.This may cause some issues with some of the Collapse tests, since those were written assuming that 'mass' is really density.
Update (1 March 2022):
I have made the following changes:
1 - The
"has_mass"
group of particle types is now called"is_gravitating"
.2 - Particle types can have either an attribute or a constant called either
"mass"
or"density"
. Thepm_deposit
method handles the four separate cases."mass"
is treated as a quantity with dimensions of mass, and"density"
is treated as a quantity with dimensions of density.I have updated some of the parameter files to reflect these changes, mostly by replacing
"mass"
with"density"
and"has_mass"
with"is_gravitating"
in appropriate places. Users will have to update their parameter files likewise. I have not changed theIsolatedGalaxy
parameter files, since it would appear (based oninput/IsolatedGalaxy/README
) that the"mass"
attribute was already being treated as a mass quantity. I have also left alone the parameter files ininput/MergeStars
since the functionality of these tests do not depend on whether"mass"
is a mass or density quantity.Regarding the issue concerning the depositing of the gas density field raised by @mabruzzo , I have undone the changes which I previously made to this part of the
pm_deposit
method. However, since I believe that the way this is done currently is incorrect, and were not fixed by the changes I made in any case, I leave this for a future PR.Update (8 March 2022):
This PR is now ready for review. I am re-requesting a review from @jwise77 since there have been substantial changes since he previously approved it.