Fix for internal energy scaling in PTE closure

[pdmullen]

PR Summary

The total volumetric internal energy

$$ u_\mathrm{total} = \sum_m \bar{\rho}_m e_m $$

can be a negative quantity. I believe that the internal energy scaling routines in singularity/main mixed_cell_models.hpp do not respect this. For example, here: https://github.com/lanl/singularity-eos/blob/main/singularity-eos/closure/mixed_cell_models.hpp#L313, an abs is applied to the (potentially negative) sie_total to compute uscale. In the case of a negative sie_total, the sum of the resultant, scaled u[m] can be -1. Inside residual calculations (e.g., https://github.com/lanl/singularity-eos/blob/main/singularity-eos/closure/mixed_cell_models.hpp#L594), a positive 1 is assumed for the sum over scaled material internal energies.

I believe the fix here is to just remove the abs (as here: c843e16). However, after some discussion with folks, this means that the sign of uscale can depend on the sign of sie_tot... an alternative choice is to introduce uscale_total = (sie_tot > 0.0) ? 1.0 : -1.0

PR Checklist

• Adds a test for any bugs fixed. Adds tests for new features.
• Format your changes by using the make format command after configuring with cmake.
• Document any new features, update documentation for changes made.
• Make sure the copyright notice on any files you modified is up to date.
• After creating a pull request, note it in the CHANGELOG.md file.
• LANL employees: make sure tests pass both on the github CI and on the Darwin CI

If preparing for a new release, in addition please check the following:

• Update the version in cmake.
• Move the changes in the CHANGELOG.md file under a new header for the new release, and reset the categories.
• Ensure that any when='@main' dependencies are updated to the release version in the package.py

[Yurlungur]

Thanks for the fix! This was a subtle one.

[dholladay00]

This does bring about a larger discussion of uscale in general. Generically, it is defined by sie_total at the moment, though the Fixed[P,T] solvers define it internal to the solver. My understanding is it is a quantity meant to stabilize the solver by capturing the rough order of magnitude of the energies and scaling them brings quantities close to O(1). One could think of it as a very simple preconditioner (like diagonal scaling but different, though @jdolence could provide better context). When the signs are allowed to differ, a sum may not represent such a quantity. It seems like uscale (this may be beyond the scope of this PR btw) should be defined by some norm or measure of the internal energy vector or some other scheme (shifting all energies by a constant such that all have the same sign and taking an average).

[dholladay00]

at the very least for this PR, perhaps changing the uscale sum to be a sum of absolute values of the energies should also help those solvers.

[Yurlungur]

at the very least for this PR, perhaps changing the uscale sum to be a sum of absolute values of the energies should also help those solvers.

I think this was what it was originally and it caused problems. @pdmullen can comment though.

[pdmullen]

at the very least for this PR, perhaps changing the uscale sum to be a sum of absolute values of the energies should also help those solvers.

I think this was what it was originally and it caused problems. @pdmullen can comment though.

Maybe a few thoughts here...

First, I tend to agree that there are multiple ways of handling the scaling of internal energies. As @dholladay00 argues, if we were careful, we could simply apply a shift so that the nuance of negative numbers don't bite us. @jhp-lanl on mattermost also thought about the scenario wherein the scaled internal energies summed to exactly zero, and how to best handle that situation. Lots of opportunity here.

Second, just to say it aloud... this PR is fixing a bug. The bug was that the scaled internal energies could sum to -1 when sie_tot < 0, despite mixed_cell_models.hpp assuming a sum to 1 elsewhere. Maybe I don't quite understand your comment regarding re-introducing absolute values, as my first take is that this would just reintroduce the original issue, but I suspect I am misreading your comment.

[dholladay00]

My comment wasn't necessarily related to the new introduced parameter as it seems to do the right thing. But seeing this PR with the fact that sie's can have opposite sign at the forefront of my mind, I think we should look again at uscale at some point. It is unlikely though not impossible that uscale could be 0 or near 0, which I don't think is the intent of the usage of uscale. My understanding of its usage is that it is supposed to bring sie-ish quantities close to O(1) by having uscale provide "a rough order of magnitude" of an internal energy value. Adding up values with opposing signs will not accomplish this task.

[pdmullen]

Sorry folks... I have reverted to my original design for this fix.

In contrast to the fixed P solver implementation for example (see https://github.com/lanl/singularity-eos/blob/main/singularity-eos/closure/mixed_cell_models.hpp#L1020), when sie_tot is provided, I feel the desired scaling is set by the bulk internal energy constraint, not the initial guess.

I have gone the route of simply eliminating the abs again, which means that uscale can be either positive or negative, depending on the sign of sie_tot. If you feel too strongly against this, I think the only other reasonable strategy in the short-term is introducing a uscale_total = (sie_tot > 0.0) ? 1.0 : -1.0.

EDIT: Just for completeness, in my previous implementation (if anyone looks there), I missed a re-normalization.

[jhp-lanl]

Sorry folks... I have reverted to my original design for this fix.

In contrast to the fixed P solver implementation for example (see https://github.com/lanl/singularity-eos/blob/main/singularity-eos/closure/mixed_cell_models.hpp#L1020), when sie_tot is provided, I feel the desired scaling is set by the bulk internal energy constraint, not the initial guess.

I have gone the route of simply eliminating the abs again, which means that uscale can be either positive or negative, depending on the sign of sie_tot. If you feel too strongly against this, I think the only other reasonable strategy in the short-term is introducing a uscale_total = (sie_tot > 0.0) ? 1.0 : -1.0.

EDIT: Just for completeness, in my previous implementation (if anyone looks there), I missed a re-normalization.

I think what you had originally (in this MR at least) is the right approach and I don't think you've provided enough justification for why that approach doesn't work yet. For this MR to go in, we need strong justification for why the uscale_total approach doesn't work and why this is "the only reasonable strategy".

I am strongly advocating for the uscale_total approach because the discussion has already brought up several edge cases. Your MR shouldn't address these edge cases, but if you just remove the abs() then you're also not helping other developers address the edge cases in the future. These fixes could be setting a floor on uscale, summing the absolute values, or even just using the maximum absolute value. Regardless of what the fix is, I think it makes sense to set uscale_total the way you were setting it and not hard-code 1 as the sum convergence criterion.

[jhp-lanl]

Can you please add the case that causes the PTE solver to fail in the unit tests?
https://github.com/lanl/singularity-eos/blob/main/test/test_closure_pte.cpp

[pdmullen]

I think what you had originally (in this MR at least) is the right approach and I don't think you've provided enough justification for why that approach doesn't work yet. For this MR to go in, we need strong justification for why the uscale_total approach doesn't work and why this is "the only reasonable strategy".

I am strongly advocating for the uscale_total approach because the discussion has already brought up several edge cases. Your MR shouldn't address these edge cases, but if you just remove the abs() then you're also not helping other developers address the edge cases in the future. These fixes could be setting a floor on uscale, summing the absolute values, or even just using the maximum absolute value. Regardless of what the fix is, I think it makes sense to set uscale_total the way you were setting it and not hard-code 1 as the sum convergence criterion.

The implementation in 9c35a24 contained ill-motivated, erroneous scaling and residual estimations. This was my fault... apologies; and I could have better described the flaw. The original, actual bug fix (i.e., simply eliminating abs) is the correct patch, hence my revert. However, as you allude to, we can still consider different designs for applying the exact same patch. In particular, we can accomplish it through the introduction of a utotal_scale, which might be a more extendible design when considering different edges in the future (c.f., comments above from @dholladay00).

I suspect that earlier we weren't on the same page about what was going on in with my failed attempt in 9c35a24. If okay with all, let's forget about that implementation and reconsider afresh the most recent commit. For this reason, I am re-triggering review requests.

Can you please add the case that causes the PTE solver to fail in the unit tests?
https://github.com/lanl/singularity-eos/blob/main/test/test_closure_pte.cpp

Will try my luck at this tomorrow.

[jhp-lanl]

I suspect that earlier we weren't on the same page about what was going on in with my failed attempt in 9c35a24. If okay with all, let's forget about that implementation and reconsider afresh the most recent commit. For this reason, I am re-triggering review requests.

Yes that was my fault. I failed to appreciate what the sum over u[m] actually meant.

[jhp-lanl]

Approved when minor change is made

[jhp-lanl]

Merging. LANL gitlab tests pass