Fixing Density floor in euler.F

[Andrew Emerick]

Hi all,

After a conversation with Greg, it looks there is a bug in euler.F that has been around for some time (possibly a known bug?). Fixing this is quite simple, but is invasive and will change code behavior. I'd like to hear people's thoughts about this first before I make the change and issue a pull request.

At the moment euler.F applies a hard-coded minimum density threshold of 1.0E-3 (in code units) when DualEnergyFormalism == 1 and GravityOn == 1. In addition to a hard coded density floor being scary and not obvious to the user, it applies the floor incorrectly, updating only the cell density value and ignoring the inverse density value that is used further down in the routine.  (relevant piece of code here: http://pastebin.com/aArUXuML )

My proposed changes are the following:

1) Add in user supplied parameters to turn on / off the density floor and a parameter to supply the density floor (in code units). (UseDensityFloor and DensityFloor respectively)

2) Correctly apply the density floor, updating both the density and inverse density values

Although this change will correct the code, it will have a consequence: previous versions of the code will not necessarily behave the same way, and it may be hard to actually determine which (if any) simulations this will have an effect on.

And before I do this, I wanted to see what people thought the best "default" behavior would be. A few options / questions:

1) UseDensityFloor = 0 and DensityFloor = 1.0E-3 , only turning it on if user requests

2) UseDensityFloor = 0 and DensityFloor = 1.0E-3, turning on if requested OR turning on if user doesn't specify and DualEnergyFormalism == 1 and GravityOn == 1 (similar to current situation)

3) Same as above, but turn it on it if HydroSolver is set to PPM rather than the DualEnergyFormalism check (since I believe that is the actual intended behavior anyway).

4) Is 1.0E-3 a reasonable default value for the density floor? If its a silly number, my vote is to make it something sensible rather than leaving it as 1.0E-3 for historical reasons.

Thoughts?

Andrew

---

NSF Graduate Research Fellow

Columbia University

Department of Astronomy

[Matthew Turk]

Hi Andrew,

Wow, thanks for the deep dive into this, and for finding/fixing it.

On Wed, May 18, 2016 at 11:25 AM, Andrew Emerick <aemer...@gmail.com> wrote:

Hi all,

After a conversation with Greg, it looks there is a bug in euler.F that has been around for some time (possibly a known bug?). Fixing this is quite simple, but is invasive and will change code behavior. I'd like to hear people's thoughts about this first before I make the change and issue a pull request.

At the moment euler.F applies a hard-coded minimum density threshold of 1.0E-3 (in code units) when DualEnergyFormalism == 1 and GravityOn == 1. In addition to a hard coded density floor being scary and not obvious to the user, it applies the floor incorrectly, updating only the cell density value and ignoring the inverse density value that is used further down in the routine.  (relevant piece of code here: http://pastebin.com/aArUXuML )

My proposed changes are the following:

1) Add in user supplied parameters to turn on / off the density floor and a parameter to supply the density floor (in code units). (UseDensityFloor and DensityFloor respectively)

2) Correctly apply the density floor, updating both the density and inverse density values

This sounds pretty good to me; I have a +0 on making it compile-time, but I don't have good reasons for why.

 

Although this change will correct the code, it will have a consequence: previous versions of the code will not necessarily behave the same way, and it may be hard to actually determine which (if any) simulations this will have an effect on.

And before I do this, I wanted to see what people thought the best "default" behavior would be. A few options / questions:

1) UseDensityFloor = 0 and DensityFloor = 1.0E-3 , only turning it on if user requests

2) UseDensityFloor = 0 and DensityFloor = 1.0E-3, turning on if requested OR turning on if user doesn't specify and DualEnergyFormalism == 1 and GravityOn == 1 (similar to current situation)

3) Same as above, but turn it on it if HydroSolver is set to PPM rather than the DualEnergyFormalism check (since I believe that is the actual intended behavior anyway).

I'm fine with #1, although I say that without having explored the consequences empirically.

 

4) Is 1.0E-3 a reasonable default value for the density floor? If its a silly number, my vote is to make it something sensible rather than leaving it as 1.0E-3 for historical reasons.

This is in the normalized units, so IIRC this is a reasonably small floor for cosmology simulations, which I assume is where most of the DEF=1 & GO=1 simulations were when the floor was inserted.

 

Thoughts?

Andrew

---

NSF Graduate Research Fellow

Columbia University

Department of Astronomy

--
You received this message because you are subscribed to the Google Groups "enzo-dev" group.
To unsubscribe from this group and stop receiving emails from it, send an email to enzo-dev+u...@googlegroups.com.
To post to this group, send email to enzo...@googlegroups.com.
Visit this group at https://groups.google.com/group/enzo-dev.
For more options, visit https://groups.google.com/d/optout.

[Nathan Goldbaum]

On Wed, May 18, 2016 at 11:25 AM, Andrew Emerick <aemer...@gmail.com> wrote:

I just want to say that John Forbes and I made a similar modification for our galaxy simulations. See: 

https://bitbucket.org/ngoldbaum/enzo-dev/commits/87cf1669150802ec5b8c43442e70da64dc058c59

I agree that 1.0E-3 is a silly (possibly catastrophic) default depending on the units used in the simulation.

 

Thoughts?

Andrew

---

NSF Graduate Research Fellow

Columbia University

Department of Astronomy

--

[Tom Abel]

Dear Andrew,

Yes it would be great to clean that up. 

You might want to reuse the parameters that the RK solvers in 

enzo/hydro_rk

are using for the same purpose.

See around line 1100 in ReadParameterFile.C :

    ret += sscanf(line, "SmallRho = %g", &SmallRho);

    ret += sscanf(line, "SmallP = %g", &SmallP);

    ret += sscanf(line, "SmallT = %g", &SmallT);

SmallRho is the density floor there. 

I would perhaps avoid an extra parameter UseDensityFloor and just have the

reasonable assumption that if SmallRho = 0  (or negative) nothing is done. 

Hope this is useful to you. 

Best,

Tom

[Britton Smith]

Hi Andrew, all,

My preference would be for the default behavior to be to do the correct thing, even if that means changing the default behavior.  In my experience, knowledge of new parameters doesn't disseminate very well.  Most people configure new simulations based on their own old parameter files, so I can see most simulations not getting run with the correct behavior.

Britton

[Cameron Hummels]

I agree that the defaults should be correct as opposed to "doing the thing that it has been doing forever", so I'm in favor of #1.  Great job catching this!!

Cameron

Cameron Hummels

NSF Postdoctoral Fellow

Department of Astronomy

California Institute of Technology

http://chummels.org

[David Collins]

I also agree with Tom-- there are already a number of floor parameters floating around in the code that make doing meaningful comparisons across solvers difficult.  Re-use of existing parameters would be good, and not having to turn on two switches will reduce the number of bear traps for users.

Thanks for your hard work!

d.

-- Sent from a computer.

[Andrew Emerick]

Hi everyone,

Thank you all for the advice. I've submitted the pull request with my changes, which fixes the implementation of the density floor without adding any new parameters. As Tom advised, I'm using SmallRho as the value of the density floor, turning it on if it is > 0 . I leave the default value of SmallRho the same, 1E-30 (much different than the hard coded 1E-3 value that was there before).

Andrew

---

NSF Graduate Research Fellow

Columbia University

Department of Astronomy