calculate con2prim enthalpy bound by logalexw · Pull Request #241 · lanl/phoebus

logalexw · 2026-03-11T17:33:41Z

PR Summary

This is to address a bug I encountered when working on the CCSNe problem; the root find in con2prim_robust.cpp isn't bracketed when initializing to the gr grid. This occurs when using a Stellar Collapse EOS (SFHo) for a progenitor at infall.

In Kastaun et al. 2021 there’s a required lower bound (h0) for relativistic enthalpy (h) s.t. 0 < h0 ≤ h. This also serves as the upper bound for the unbracketed root find, i.e. root find using regula falsi from (0, 1/h0].

h0 is currently hardcoded to 1.0, which breaks our upper bound condition if h < 1, which does happen if you’re using a StellarCollapse EOS (see fig. 1 below), espeically at Ye values that are relevant for SNe. When looking at the evolution of the root find upper bound for CCSNe progenitor, it's visible that the unbracketed root find (red dashed line) occurs when 1/h ≥ 1/h0 1 (see fig. 2 below).

Goal is to fix with something that's eos-conscious; there's existing primitives in con2prim that should be able to calculate a reasonable bound using singularity-EOS (the existing method in compute_upper_bound() isn't sufficient for CCSNe problem).

Fig. 1. Minimum enthalpy for the SFHo table at relevant electron fraction (Ye) range.

Fig. 2. Evolution of upper bound (1/h0) for a CCSNe progenitor, with location of unbracketed root find marked in red. h has been scaled so that the evolution is visible in the same range as the specific internal energy.

PR Checklist

Adds a test for any bugs fixed. Adds tests for new features.
Format your changes by calling scripts/bash/format.sh.
Explain what you did.
Make any necessary changes to the documentation.

Yurlungur · 2026-03-11T18:48:24Z

Good catch! Should have known we would need to worry about this--nubhlight also has to handle this case, as does the torus initial conditions.

c-prather · 2026-03-12T18:33:26Z

Oh, I would benefit from stealing this too I think. Though iirc I set the brackets following AthenaK's implementation, so maybe I ended up with an enthalpy bound from that.
It might be worth your time looking at KHARMA's stolen & modded version of this function, we have some downstream tweaks that really improve stability.

Yurlungur · 2026-03-25T15:51:40Z

@logalexw make sure to remove "WIP" and ping when you're ready for review. No rush though!

logalexw · 2026-03-25T17:13:44Z

I've refactored the original fix to add an parameter that stors the enthalpy lower bound and is able to be utilized in conserved to primitive reconstruction (both initialization and fixup)!

Unfortunately, there seems to be another bug (absolute value??) that causes epsilon to be strictly positive. The new updates to con2prim_robust.cpp produce the dashed line, which approximates our expected initial SIE well, but phoebus produces the "mirrored" SIE curve (in blue) on the zeroth timestep. Any ideas where this is happening? I've had no luck finding the source in src/fixup so far, which is what @AstroBarker originally suspected.

Yurlungur · 2026-03-25T17:19:53Z

Is epsilon here specific internal energy? It's like the floors.

Yurlungur · 2026-03-25T17:20:27Z

An energy floor is an optional feature in phoebus. You can try disabling or making the floor very small in the input deck.

logalexw · 2026-03-25T17:23:58Z

I thought it was the floors as well, I've set them to a reasonable negative value (for the SC EOS tables + CCSNe problem)– I was debugging with verbose output and the floors (sie, u, rho) were all at expected values (SFHo table minimum for sie, very small values for rho, etc.)

Yurlungur · 2026-03-25T17:27:25Z

Ah ok. Well if it's not floors, the next thing I would check would just to grep through the code for calls to std::abs or fabs. There will be a lot... but hopefully you can quickly narrow it down. I agree an absolute value seems like a suspect.

AstroBarker · 2026-03-25T18:34:10Z

Good work @logalexw ! When you have time could you share your workflow to reproduce this (cmake build line for example). And this flipped behavior is present in the 00000.phdf output? And you're plotting the sie output directly? And your fixup block looks like this:

<fixup>
enable_floors = true
enable_ceilings = true
enable_flux_fixup = true
enable_c2p_fixup = true
floor_type = ConstantRhoSie
# updated floors!
rho0_floor = 1.0e-10
sie0_floor = -1.0e2 # << changed back to original oom.
ceiling_type = ConstantGamSie
sie0_ceiling = 1.0e3
gam0_ceiling = 1.0e4

I am asking as I am having some trouble reproducing the behavior using this branch.

AstroBarker · 2026-03-25T18:36:09Z

And I am seeing that the docs action is failing with a sphinx-multiversion issue.. This thing is so unstable, they changed the version of something under the hood in runner. Happened to me on another project. I'll check for the fix before we merge.

Yurlungur · 2026-03-25T18:41:54Z

And I am seeing that the docs action is failing with a sphinx-multiversion issue.. This thing is so unstable, they changed the version of something under the hood in runner. Happened to me on another project. I'll check for the fix before we merge.

Yeah just need to pin sphinx. Here's the fix in singularity-eos
https://github.com/lanl/singularity-eos/blob/main/.github/workflows/docs.yml#L34

logalexw · 2026-03-25T18:48:58Z

Good work @logalexw ! When you have time could you share your workflow to reproduce this (cmake build line for example). And this flipped behavior is present in the 00000.phdf output? And you're plotting the sie output directly? And your fixup block looks like this:
<fixup>
enable_floors = true
enable_ceilings = true
enable_flux_fixup = true
enable_c2p_fixup = true
floor_type = ConstantRhoSie
# updated floors!
rho0_floor = 1.0e-10
sie0_floor = -1.0e2 # << changed back to original oom.
ceiling_type = ConstantGamSie
sie0_ceiling = 1.0e3
gam0_ceiling = 1.0e4
I am asking as I am having some trouble reproducing the behavior using this branch.

Yes, these are the correct parameters, and yes, the flip happens in 00000.phdf - I get the same result plotting both sie directly and using p.energy/p.density.

I've built this on UA HPC using cmake .. -DPHOEBUS_GEOMETRY=MonopoleSph; make -j8

logalexw · 2026-03-25T19:00:44Z

Another thing I just saw - is it intended that both drho, du and other derived quantities are always positive in fixup.cpp? Allowing for negative energy floors breaks that assumption.

Yurlungur · 2026-03-25T19:03:51Z

Another thing I just saw - is it intended that both drho, du and other derived quantities are always positive in fixup.cpp? Allowing for negative energy floors breaks that assumption.

I'm not sure I understand... shouldn't they still be always positive? If u < uflr, and |u| > |uflr|, then |du| > 0 still to add energy, right?

logalexw · 2026-03-25T23:30:18Z

previous issue with mirrored sie is resolved (thanks to @AstroBarker for the catch)!

AstroBarker

Nice work fixing this! Looks good, a couple of small suggestions for mostly cleanliness and organization.

src/fluid/con2prim_robust.hpp

AstroBarker · 2026-03-26T16:07:13Z

src/microphysics/eos_phoebus/eos_phoebus.cpp

 #include "phoebus_utils/unit_conversions.hpp"
 #include "phoebus_utils/variables.hpp"

+// inspired by SPACELOOP in geometry utils, util for global lower bound for enthalpy


I would suggest either removing this and writing the for loops out explicitly below or moving it to a better home, maybe somewhere under phoebus_utils and including the relevant header in this file. Either is fine.

I think it's fine here... It's local to this implementation file, so it can clearly be seen what it's used for. That said, the loops in geometry are written to convey intent... the bounds are always the same. If this macro is to be generalized it needs a more useful name than LOOP, like LOOP1D or something. Or PHOEBUSFOR or similar.

Personally I think this isn't really saving us enough characters that I feel the need to factor it out to utils.

I agree with Brandon; the macro originally had more functionality before I simplified the logic, but in this form the intent would probably be just as clear as a set of explicit for loops (and it's not really saving that much room).

Yurlungur

Nice work! I have a few comments about the enthalpy finding loop but my comments are non-blocking.

Yurlungur · 2026-03-26T16:20:43Z

src/microphysics/eos_phoebus/eos_phoebus.cpp

 #include "phoebus_utils/unit_conversions.hpp"
 #include "phoebus_utils/variables.hpp"

+// inspired by SPACELOOP in geometry utils, util for global lower bound for enthalpy


I think it's fine here... It's local to this implementation file, so it can clearly be seen what it's used for. That said, the loops in geometry are written to convey intent... the bounds are always the same. If this macro is to be generalized it needs a more useful name than LOOP, like LOOP1D or something. Or PHOEBUSFOR or similar.

Personally I think this isn't really saving us enough characters that I feel the need to factor it out to utils.

Yurlungur · 2026-03-26T17:46:09Z

src/microphysics/eos_phoebus/eos_phoebus.cpp

+    LOOP(y, n / 2) {
+      lambda[0] = ye;
+      LOOP(r, n) {
+        LOOP(t, n) {
+          eps = eos_sc.InternalEnergyFromDensityTemperature(rho, T, lambda) / sie_unit;
+          P = eos_sc.PressureFromDensityTemperature(rho, T, lambda) / press_unit;
+          h_min = std::min(h_min, 1 + eps + robust::ratio(P, rho));
+          T += dT;
+        }
+        rho += drho;
+      }
+      ye += dye;
+    }


This is fine for now but I think it can probably be done better... For example, one could do a Newton method to find the point of minimal enthalpy... Or by having the EOS report this itself internally. Food for thought for later.

@carlnotsagan and I were talking about this as a later wishlist item for the EOS as well. I'm definitely open to doing something more efficient in the future as well, i don't think this is the best fix per se.

Performancewise it probably doesn't matter, since it's just initial conditions. But I worry it's a bit brittle.

Yurlungur · 2026-03-26T17:47:48Z

src/microphysics/eos_phoebus/eos_phoebus.cpp

+    // spacing for each dimension (i.e. np.linspace)
+    dT = (T_max - T_min) / n;
+    drho = (rho_max - rho_min) / n;
+    dye = (ye_max - ye_min) / (n / 2); // we don't need to resolve ye as much
+    // initial conditions
+    T = T_min;
+    rho = rho_min;
+    ye = ye_min;


The table spans many orders of magnitude. I think you can use a coarser spacing if you used a log spacing in density and temperature. Or, realistically, you can simply only look at relatively low densities and temperatures, as high densities and temperatures simply won't have a small enthalpy.

logalexw · 2026-04-01T20:33:37Z

final formatting updates from the comments above are complete! if there are no further adjustments/suggestions, this should be ready to merge.

Yurlungur · 2026-04-01T20:50:58Z

Merged. Thanks for the fix @logalexw !

todos plus initial context

db15073

logalexw self-assigned this Mar 11, 2026

logalexw added bug Something isn't working enhancement New feature or request labels Mar 11, 2026

Logan Alexandria White and others added 7 commits March 13, 2026 10:10

first draft of h0 calculation

44c59d0

fixed scope; added method to Residual

f70c932

verbose basics for testing

37ee67e

works to initialize ccsne (yay)

8a535d0

new bounds condition from kastaun+21, highly rel. case

0768965

refactor, global lower bound for enthalpy

b3bcb5a

general clean up, etc

64eca66

removing verbose testing, comments

9d22331

logalexw changed the title ~~WIP: calculate con2prim enthalpy bound~~ calculate con2prim enthalpy bound Mar 25, 2026

logalexw requested a review from AstroBarker March 25, 2026 23:31

fix: pin sphinx version for docs

8ef8a95

clang formatting

bba4e30

AstroBarker approved these changes Mar 26, 2026

View reviewed changes

Yurlungur approved these changes Mar 26, 2026

View reviewed changes

minor bug fix, additional format cleanup

f3cb584

Yurlungur approved these changes Apr 1, 2026

View reviewed changes

Yurlungur merged commit 7704c35 into main Apr 1, 2026
3 checks passed

Yurlungur deleted the law/enthalpy-bound branch April 1, 2026 20:51

Conversation

logalexw commented Mar 11, 2026

PR Summary

PR Checklist

Uh oh!

Yurlungur commented Mar 11, 2026

Uh oh!

c-prather commented Mar 12, 2026

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

logalexw commented Mar 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

logalexw commented Mar 25, 2026

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

AstroBarker commented Mar 25, 2026

Uh oh!

AstroBarker commented Mar 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

logalexw commented Mar 25, 2026

Uh oh!

logalexw commented Mar 25, 2026

Uh oh!

Yurlungur commented Mar 25, 2026

Uh oh!

logalexw commented Mar 25, 2026

Uh oh!

AstroBarker left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Yurlungur left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

logalexw commented Apr 1, 2026

Uh oh!

Uh oh!

Yurlungur commented Apr 1, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

logalexw commented Mar 25, 2026 •

edited

Loading

AstroBarker commented Mar 25, 2026 •

edited

Loading

AstroBarker left a comment •

edited

Loading