Skip to content
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.

Sign up for GitHub

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

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

Closed
wants to merge 7 commits into from
Closed

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

Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
c3206d7
copy/past boiler-plate + some changes, does not compile at the moment
dholladay00 Oct 9, 2024
5ebda78
further down the road in terms of init and residual. Jacobian is next…
dholladay00 Oct 9, 2024
c5ab4d6
fix ctor.
dholladay00 Oct 10, 2024
8237aa3
Some progress on jacobian calc. Not there yet, probably doesnt compile
dholladay00 Oct 10, 2024
efc1033
Fix residual and number of equations.
dholladay00 Oct 25, 2024
0e32674
More work on check pte, etc. Not quite there yet, should probably cre…
dholladay00 Oct 25, 2024
331e5a1
Merge branch 'main' into dholladay00/pt_solver
dholladay00 Oct 25, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
236 changes: 236 additions & 0 deletions singularity-eos/closure/mixed_cell_models.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -741,6 +741,242 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
Real Tequil, Ttemp;
};

// PT space solver
inline int PTESolverPTRequiredScratch(const int nmat) {
int neq = 2;
return neq * neq // jacobian
+ 4 * neq // dx, residual, and sol_scratch
+ 3 * nmat // all the nmat sized arrays
+ MAX_NUM_LAMBDAS * nmat; // the cache
}
inline size_t PTESolverPTRequiredScratchInBytes(const int nmat) {
return PTESolverRhoTRequiredScratch(nmat) * sizeof(Real);
}

template <typename EOSIndexer, typename RealIndexer, typename LambdaIndexer>
class PTESolverPT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::InitBase;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::AssignIncrement;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::nmat;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::neq;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::ResidualNorm;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::vfrac_total;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::sie_total;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::eos;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rho;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::vfrac;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::sie;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::temp;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::press;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rho_total;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::uscale;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::utotal_scale;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::MatIndex;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::TryIdealPTE;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::jacobian;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::residual;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::dx;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::u;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::rhobar;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::Cache;
using mix_impl::PTESolverBase<EOSIndexer, RealIndexer>::Tnorm;

public:
// template the ctor to get type deduction/universal references prior to c++17
template <typename EOS_t, typename Real_t, typename Lambda_t>
PORTABLE_INLINE_FUNCTION
PTESolverPT(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real sie_tot,
Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, Real_t &&press,
Lambda_t &&lambda, Real *scratch, const Real Tguess = 0.0)
: mix_impl::PTESolverBase<EOSIndexer, RealIndexer>(nmat, nmat + 1, eos, vfrac_tot,
sie_tot, rho, vfrac, sie, temp,
press, scratch, Tguess) {
vtemp = AssignIncrement(scratch, nmat);
// TODO(JCD): use whatever lambdas are passed in
/*for (int m = 0; m < nmat; m++) {
if (!variadic_utils::is_nullptr(lambda[m])) Cache[m] = lambda[m];
}*/
}

PORTABLE_INLINE_FUNCTION
Real Init() {
InitBase();
Residual();
// 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;
Comment on lines +820 to +829
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See #423

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't we also normalize the pressures too?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Can you point to this code?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

// Leave this in for now, but comment out because I'm not sure it's a good idea
// TryIdealPTE(this);
// Set the current guess for the equilibrium temperature. Note that this is already
// scaled.
return ResidualNorm();
}

PORTABLE_INLINE_FUNCTION
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;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

}

PORTABLE_INLINE_FUNCTION
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;
}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).


PORTABLE_INLINE_FUNCTION
void Jacobian() const {
using namespace mix_params;
// sum_m d e_m / dT )_P
Real dedT_P_sum = 0.0;
// sum_m d e_m / dP )_T
Real dedP_T_sum = 0.0;
// - sum_m rho_tot / rho_m^2 * d rho_m / dT )_P
Real rtor2_dr_dT_p_sum = 0.0;
// - sum_m rho_tot / rho_m^2 d rho_m / dP )_T
Real rtor2_dr_dP_T_sum = 0.0;

for (int m = 0; m < nmat; m++) {
//////////////////////////////
// perturb volume fractions
//////////////////////////////
Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps;
const Real vf_pert = vfrac[m] + dv;
const Real rho_pert = robust::ratio(rhobar[m], vf_pert);

Real p_pert{};
Real e_pert =
eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]);
p_pert =
robust::ratio(this->GetPressureFromPreferred(eos[m], rho_pert, Tnorm * Tequil,
e_pert, Cache[m], false),
uscale);
Real dpdv = robust::ratio((p_pert - press[m]), dv);
Real dedv = robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dv);
//////////////////////////////
// perturb temperature
//////////////////////////////
Real dT = Tequil * derivative_eps;
e_pert = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * (Tequil + dT),
Cache[m]);
p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m],
Tnorm * (Tequil + dT), e_pert,
Cache[m], false),
uscale);
Real dpdT = robust::ratio((p_pert - press[m]), dT);
rtor2_dr_dT_p_sum += rho_tot / (rho[m]*rho[m]) * dpdT / dpdv;
dedT_sum += robust::ratio(rhobar[m] * robust::ratio(e_pert, uscale) - u[m], dT);
}

// Fill in the Jacobian
for (int i = 0; i < neq * neq; ++i)
jacobian[i] = 0.0;
jacobian[0] = rtor2_dr_dT_p_sum;
jacobian[1] = rtor2_dr_dP_T_sum;
Comment on lines +908 to +909
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

jacobian[2] = dedT_P_sum;
jacobian[3] = dedP_T_sum;
}

PORTABLE_INLINE_FUNCTION
Real ScaleDx() const {
using namespace mix_params;
// Each check reduces the scale further if necessary
Real scale = 1.0;
// control how big of a step toward vfrac = 0 is allowed
for (int m = 0; m < nmat; ++m) {
if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) {
scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]);
}
}
const Real Tnew = Tequil + scale * dx[nmat];
// control how big of a step toward rho = rho(Pmin) is allowed
for (int m = 0; m < nmat; m++) {
const Real rho_min =
std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew));
const Real alpha_max = robust::ratio(rhobar[m], rho_min);
if (alpha_max < vfrac[m]) {
// Despite our best efforts, we're already in the unstable regime (i.e.
// dPdV_T > 0) so we would actually want to *increase* the step instead
// of decreasing it. As a result, this code doesn't work as intended and
// should be skipped.
continue;
}
if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) {
scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]);
}
}
// control how big of a step toward T = 0 is allowed
if (scale * dx[nmat] < -0.95 * Tequil) {
scale = robust::ratio(-0.95 * Tequil, dx[nmat]);
}
// Now apply the overall scaling
for (int i = 0; i < neq; ++i)
dx[i] *= scale;
return scale;
}

// Update the solution and return new residual. Possibly called repeatedly with
// different scale factors as part of a line search
PORTABLE_INLINE_FUNCTION
Real TestUpdate(const Real scale, bool const cache_state = false) {
if (cache_state) {
// Store the current state in temp variables for first iteration of line
// search
Ttemp = Tequil;
for (int m = 0; m < nmat; ++m)
vtemp[m] = vfrac[m];
}
Tequil = Ttemp + scale * dx[nmat];
for (int m = 0; m < nmat; ++m) {
vfrac[m] = vtemp[m] + scale * dx[m];
rho[m] = robust::ratio(rhobar[m], vfrac[m]);
u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature(
rho[m], Tnorm * Tequil, Cache[m]);
sie[m] = robust::ratio(u[m], rhobar[m]);
u[m] = robust::ratio(u[m], uscale);
temp[m] = Tequil;
press[m] =
robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * Tequil,
sie[m], Cache[m], false),
uscale);
}
Residual();
return ResidualNorm();
}

private:
Real *vtemp;
Real Tequil, Ttemp, Pequil;
};

// fixed temperature solver
inline int PTESolverFixedTRequiredScratch(const int nmat) {
int neq = nmat;
Expand Down
Loading