Draft: copy/past boiler-plate + some changes, does not compile at the moment #422

dholladay00 · 2024-10-09T16:41:59Z

PR Summary

early beginnings of the pt space pte solver. It's mostly a copy-paste for now. Hopefully I can get things more cleaned up soon. Feel free to take a look @jhp-lanl

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

… but need to make sure where to put derivatives

jhp-lanl

Great start! Let me know when you want further review

singularity-eos/closure/mixed_cell_models.hpp

jhp-lanl · 2024-10-14T20:10:35Z

singularity-eos/closure/mixed_cell_models.hpp

+    // calculate an initil equilibrium pressure
+    Real psum = 0.0;
+    Real vsum = 0.0;
+    for (int m = 0; m < nmat; ++m) {
+      psum += vfrac[m]*press[m];
+      vsum += vfrac[m];
+    }
+    // all normalized temps set to 1 so no averaging necessary here.
+    Tequil = temp[0];
+    Pequil = psum / vsum; 


This machinery has always confused me. We pass in a Tguess that then gets saved as Tnorm (unless it isn't passed in) but it seems like we always use temp[0] as the initial starting point for the solver, yes?

We could take the same approach with press[0], but your approach does provide a reasonable guess if the pressure is evolved by the host code in some way. We will just need to note that press needs to contain sane values of some sort.

While I'm at it, I'll create an issue to rename Tguess and note that temp requires at least one sane value.

Recall that temp is in arbitrary normalized temperature units and so it will always be multiplied by Tnorm.

Don't we also normalize the pressures too?

Yep, I think they are normalized by an internal energy.

I'm fine with either duplicating the temperature approach (i.e. Pequil = press[0]) or staying with this. Whatever we do though, it needs to be noted in the documentation what values need to be sane.

The pressures are initialized from a lookup so they are potentially different whereas there is 1 temperature guess calculated and it is folded into the temperature normalization and the temperature vector is initialized to 1 for every element.

The pressures are initialized from a lookup so they are potentially different

Can you point to this code?

https://github.com/lanl/singularity-eos/blob/dholladay00/pt_solver/singularity-eos/closure/mixed_cell_models.hpp#L328

Ah okay, I'm with you now. We should definitely document this.

To summarize my understanding:

The initial temperature guess for the PTE solver is Tequil = temp[0]

The provided pressure can be nonsense since it is populated by InitBase()

InitBase() calculates the initial pressure from Tscale and the input density, NOT Tequil

You're volume fraction averaging those pressures to get Pequil, which is the initial PTE state

jhp-lanl · 2024-10-14T20:46:06Z

singularity-eos/closure/mixed_cell_models.hpp

+  void Residual() const {
+    Real tsum = 0.0;
+    Real psum = 0.0;
+    Real vsum = 0.0;
+    // is volume averaging the right thing here? 
+    for (int m = 0; m < nmat; ++m) {
+      tsum += vfrac[m]*temp[m];
+      psum += vfrac[m]*press[m];
+      vsum += vfrac[m];
+    }
+    residual[0] = Tequil - tsum/vsum;
+    residual[1] = Pequil - psum/vsum; 


I think it's worth discussing possible convergence criteria, but your suggestion may end up being the best approach.

The residual equations are formulated based off of the density and energy sums, which is why they're included in the residuals for the other solvers. Of course in the PTESolverRhoT solver, we also populate the rest of the neq - 2 entries with pressure differences, but that's because we require additional constraints in addition to the primary residual equations. Not including the energy/volume sums in the residual would mean that we're determining convergence separate from the actual residual equations themselves. This isn't necessarily bad, but we need to be intentional about it.

If we imagine a situation with a small volume of metal surrounded by gas, the pressure is highly sensitive to changes in the volume fraction of metal in the mixture. The $\left(\frac{\partial P}{\partial e}\right)_V = \rho \Gamma$ derivative for the metal is likely near 10-ish, so this implies that the energy variability will be maybe 10% of the pressure variability due to volume fraction changes. If we assume the volume fraction is maybe order 1e-05, then the mass fraction might be order 1e-02, and since the energy sums with the mass fraction, I would estimate the variability to be around 1e-03 for a doubling of the pressure. The total volume fraction variability would be near zero since the bulk modulus is large.

With your residuals, the temperature residual variability would be similar to the energy residual variability above, but the pressure residual variability would be on the order of the volume fraction, or 1e-05. So this formulation might be less sensitive to small volume fractions of stiff materials. The plus side of course is less iterations, but the downside might be a solid density inconsistent with the actual PTE pressure. I'm not sure that's very bad since small volume fractions of solid/metal are always annoying, and being robust to them is a good thing.

Maybe this is an opportunity though. Do you think you could add a switching variable to the class that switches between a traditional residual (energy and volume fraction sums) and this one you're proposing? I.e. keep what you have here but use something more like the first two entries of PTESolverRhoT when a switch is on/off?

Should some of the checks you're talking about in the alternate scheme just be placed in the checkpte function? The residual function as I understand it is really about the unknowns of the system and the checkpte function exists for checking additional constraints AFAIK.

In that case, it might make sense to add the Tequil - tsum/vsum and Pequil - psum/vsum checks to the checkpte function since they're not the actual residual equations being solved. The information is valuable, but I'm concerned that these aren't the actual residuals of the residual equations we're solving.

from my reading of the document (equation 12), these exrpessions form the residual of [T* - T, P*-P].

I'd consider equation 13 to be the residual. Equation 12 arises from the linearization of the nonlinear system, and is really just an update vector for the Newton solve.

In my experience, the convergence criterion is primarily a function of some norm of the residual itself, but it can also sometimes be related to the norm of the update when satisfying the residual equations isn't as important as finding a fixed point. In this case though, I'm concerned that volume/energy conservation is more important than all the materials being at the same temperature/pressure.

yep, I'm all sorts of confused here. You are correct. I'll fix this when I get a chance.

jhp-lanl · 2024-10-14T20:50:18Z

singularity-eos/closure/mixed_cell_models.hpp

+  bool CheckPTE() const {
+    using namespace mix_params;
+    Real mean_p = vfrac[0] * press[0];
+    Real error_p = 0.0;
+    for (int m = 1; m < nmat; ++m) {
+      mean_p += vfrac[m] * press[m];
+      error_p += residual[m + 1] * residual[m + 1];
+    }
+    error_p = std::sqrt(error_p);
+    Real error_u = std::abs(residual[1]);
+    // Check for convergence
+    bool converged_p = (error_p < pte_rel_tolerance_p * std::abs(mean_p) ||
+                        error_p < pte_abs_tolerance_p);
+    bool converged_u = (error_u < pte_rel_tolerance_e || error_u < pte_abs_tolerance_e);
+    return converged_p && converged_u;
+  }


It looks like this code needs to be updated to reflect the structure of the residual.

My suggestion above to include two residual formulations would complicate this function (perhaps requiring two implementations of the function actually).

jhp-lanl · 2024-10-14T20:51:33Z

singularity-eos/closure/mixed_cell_models.hpp

+    jacobian[0] = rtor2_dr_dT_p_sum;
+    jacobian[1] = rtor2_dr_dP_T_sum;


Make sure the sign is correct here, but let me know if there are any issues with my math too.

I'll check again, but IIRC there is a double negative that cancels.

…ate a test soon.

dholladay00 added 4 commits October 9, 2024 10:32

copy/past boiler-plate + some changes, does not compile at the moment

c3206d7

further down the road in terms of init and residual. Jacobian is next…

5ebda78

… but need to make sure where to put derivatives

fix ctor.

c5ab4d6

Some progress on jacobian calc. Not there yet, probably doesnt compile

8237aa3

jhp-lanl reviewed Oct 14, 2024

View reviewed changes

dholladay00 added 3 commits October 25, 2024 10:54

Fix residual and number of equations.

efc1033

More work on check pte, etc. Not quite there yet, should probably cre…

0e32674

…ate a test soon.

Merge branch 'main' into dholladay00/pt_solver

331e5a1

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Draft: copy/past boiler-plate + some changes, does not compile at the moment #422

Draft: copy/past boiler-plate + some changes, does not compile at the moment #422

dholladay00 commented Oct 9, 2024

jhp-lanl left a comment

jhp-lanl Oct 14, 2024

jhp-lanl Oct 14, 2024

dholladay00 Oct 16, 2024

jhp-lanl Oct 16, 2024

dholladay00 Oct 18, 2024

jhp-lanl Oct 21, 2024

dholladay00 Oct 21, 2024

jhp-lanl Oct 22, 2024

dholladay00 Oct 22, 2024

jhp-lanl Oct 22, 2024

jhp-lanl Oct 14, 2024

dholladay00 Oct 18, 2024

jhp-lanl Oct 21, 2024

dholladay00 Oct 21, 2024

jhp-lanl Oct 22, 2024

dholladay00 Oct 22, 2024

jhp-lanl Oct 14, 2024

jhp-lanl Oct 14, 2024

dholladay00 Oct 18, 2024

		jacobian[0] = rtor2_dr_dT_p_sum;
		jacobian[1] = rtor2_dr_dP_T_sum;

Draft: copy/past boiler-plate + some changes, does not compile at the moment #422

Are you sure you want to change the base?

Draft: copy/past boiler-plate + some changes, does not compile at the moment #422

Conversation

dholladay00 commented Oct 9, 2024

PR Summary

PR Checklist

jhp-lanl left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment