Skip to content

Commit

Permalink
some more cleaning
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale committed Dec 28, 2023
1 parent 046a42a commit ca00655
Show file tree
Hide file tree
Showing 2 changed files with 169 additions and 129 deletions.
18 changes: 18 additions & 0 deletions massive_star/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,3 +21,21 @@ ulimit -s 32768
```

since the arrays are put on the stack.

## Overall algorithm

The basic HSE algorithm proceeds as:

* Map the raw MESA model onto a uniform grid

* Reset the composition of any zones that are in NSE using the NSE
table

* (optional) If the first MESA r is at a larger radius than the
innermost uniform model zone, integrate inward from the first MESA
point at constant entropy to find the central density in HSE at our
model resolution. This is only necessary if the MESA model is very
low resolution.

* Integrate outward from the center taking T and composition from the
MESA model and enforcing HSE and NSE with our EOS.
280 changes: 151 additions & 129 deletions massive_star/init_1d.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@ AMREX_INLINE void init_1d() {

}

// make sure that the species (mass fractions) summ to 1
// make sure that the species (mass fractions) sum to 1

Real summ = 0.0_rt;
for (int n = 0; n < NumSpec; ++n) {
Expand Down Expand Up @@ -251,11 +251,10 @@ AMREX_INLINE void init_1d() {
Real ye;
Real xn[NumSpec];

// because the MESA model likely begins at a larger radius than
// our first HSE model zone, simple interpolation will not do a
// good job. We want to integrate in from the zone that best
// matches the first MESA model zone, assuming HSE and constant
// entropy.
// the MESA model may begin at a larger radius than our first HSE
// model zone, so simple interpolation will not do a good job. We
// want to integrate in from the zone that best matches the first
// MESA model zone, assuming HSE and constant entropy.

// find the zone in the uniformly gridded model that corresponds to the
// first zone of the original model
Expand All @@ -269,186 +268,209 @@ AMREX_INLINE void init_1d() {
}
}

// store the central density. We will iterate until the central density
// converges
std::cout << "ibegin = " << ibegin << std::endl;

Real central_density = model_mesa_hse(0, model::idens);
if (ibegin > 0) {

std::cout << "interpolated central density = " << central_density << std::endl;
// store the central density. We will iterate until the central density
// converges

bool converged_central_density = false;
Real central_density = model_mesa_hse(0, model::idens);

for (int iter_dens = 0; iter_dens < MAX_ITER; ++iter_dens) {
std::cout << "interpolated central density = " << central_density << std::endl;

// compute the enclosed mass
bool converged_central_density = false;

M_enclosed(0) = (4.0_rt/3.0_rt) * M_PI * std::pow(delx, 3) * model_mesa_hse(0, model::idens);
for (int iter_dens = 0; iter_dens < MAX_ITER; ++iter_dens) {

for (int i = 1; i <= ibegin; ++i) {
M_enclosed(i) = M_enclosed(i-1) +
(4.0_rt/3.0_rt) * M_PI * (xznr(i) - xznl(i)) *
(std::pow(xznr(i), 2) + xznl(i) * xznr(i) + std::pow(xznl(i), 2)) *
model_mesa_hse(i, model::idens);
}
// compute the enclosed mass

// now start at ibegin and integrate inward
M_enclosed(0) = (4.0_rt/3.0_rt) * M_PI * std::pow(delx, 3) * model_mesa_hse(0, model::idens);

eos_state.T = model_mesa_hse(ibegin, model::itemp);
eos_state.rho = model_mesa_hse(ibegin, model::idens);
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = model_mesa_hse(ibegin, model::ispec+n);
}
for (int i = 1; i <= ibegin; ++i) {
M_enclosed(i) = M_enclosed(i-1) +
(4.0_rt/3.0_rt) * M_PI * (xznr(i) - xznl(i)) *
(std::pow(xznr(i), 2) + xznl(i) * xznr(i) + std::pow(xznl(i), 2)) *
model_mesa_hse(i, model::idens);
}

eos_state.aux[AuxZero::iye] = model_mesa_hse(ibegin, model::iyef);
// now start at ibegin and integrate inward

set_aux(eos_state);
eos_state.T = model_mesa_hse(ibegin, model::itemp);
eos_state.rho = model_mesa_hse(ibegin, model::idens);
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = model_mesa_hse(ibegin, model::ispec+n);
}

eos(eos_input_rt, eos_state);
eos_state.aux[AuxZero::iye] = model_mesa_hse(ibegin, model::iyef);

model_mesa_hse(ibegin, model::ipres) = eos_state.p;
set_aux(eos_state);

for (int i = 0; i < problem_rp::nx; ++i) {
entropy_want(i) = eos_state.s;
}
eos(eos_input_rt, eos_state);

for (int i = ibegin-1; i >= 0; --i) {
model_mesa_hse(ibegin, model::ipres) = eos_state.p;

// as the initial guess for the temperature and density, use
// the previous zone
for (int i = 0; i < problem_rp::nx; ++i) {
entropy_want(i) = eos_state.s;
}

dens_zone = model_mesa_hse(i+1, model::idens);
temp_zone = model_mesa_hse(i+1, model::itemp);
for (int n = 0; n < NumSpec; ++n) {
xn[n] = model_mesa_hse(i, model::ispec+n);
}
for (int i = ibegin-1; i >= 0; --i) {

ye = model_mesa_hse(i, model::iyef);
// as the initial guess for the temperature and density, use
// the previous zone

// compute the gravitational acceleration on the interface between zones
// i and i+1
dens_zone = model_mesa_hse(i+1, model::idens);
temp_zone = model_mesa_hse(i+1, model::itemp);
for (int n = 0; n < NumSpec; ++n) {
xn[n] = model_mesa_hse(i, model::ispec+n);
}

Real g_zone = -C::Gconst * M_enclosed(i) / (xznr(i) * xznr(i));
ye = model_mesa_hse(i, model::iyef);

// iteration loop
// compute the gravitational acceleration on the interface between zones
// i and i+1

// start off the Newton loop by saying that the zone has not converged
Real g_zone = -C::Gconst * M_enclosed(i) / (xznr(i) * xznr(i));

bool converged_hse = false;
// iteration loop

Real p_want;
Real drho;
Real dtemp;
// start off the Newton loop by saying that the zone has not converged

for (int iter = 0; iter < MAX_ITER; ++iter) {
bool converged_hse = false;

p_want = model_mesa_hse(i+1, model::ipres) -
delx * 0.5_rt * (dens_zone + model_mesa_hse(i+1, model::idens)) * g_zone;
Real p_want;
Real drho;
Real dtemp;

// now we have two functions to zero:
// A = p_want - p(rho,T)
// B = entropy_want - s(rho,T)
// We use a two dimensional Taylor expansion and find the
// deltas for both density and temperature
for (int iter = 0; iter < MAX_ITER; ++iter) {

// (t, rho) -> (p, s)
p_want = model_mesa_hse(i+1, model::ipres) -
delx * 0.5_rt * (dens_zone + model_mesa_hse(i+1, model::idens)) * g_zone;

eos_state.T = temp_zone;
eos_state.rho = dens_zone;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = xn[n];
}
// now we have two functions to zero:
// A = p_want - p(rho,T)
// B = entropy_want - s(rho,T)
// We use a two dimensional Taylor expansion and find the
// deltas for both density and temperature

eos_state.aux[AuxZero::iye] = ye;
// (t, rho) -> (p, s)

set_aux(eos_state);
eos_state.T = temp_zone;
eos_state.rho = dens_zone;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = xn[n];
}

eos(eos_input_rt, eos_state);
eos_state.aux[AuxZero::iye] = ye;

entropy = eos_state.s;
pres_zone = eos_state.p;
set_aux(eos_state);

Real dpT = eos_state.dpdT;
Real dpd = eos_state.dpdr;
Real dsT = eos_state.dsdT;
Real dsd = eos_state.dsdr;
eos(eos_input_rt, eos_state);

Real A = p_want - pres_zone;
Real B = entropy_want(i) - entropy;
entropy = eos_state.s;
pres_zone = eos_state.p;

Real dAdT = -dpT;
Real dAdrho = -0.5_rt * delx * g_zone - dpd;
Real dBdT = -dsT;
Real dBdrho = -dsd;
Real dpT = eos_state.dpdT;
Real dpd = eos_state.dpdr;
Real dsT = eos_state.dsdT;
Real dsd = eos_state.dsdr;

dtemp = (B - (dBdrho / dAdrho) * A) /
((dBdrho / dAdrho) * dAdT - dBdT);
Real A = p_want - pres_zone;
Real B = entropy_want(i) - entropy;

drho = -(A + dAdT * dtemp) / dAdrho;
Real dAdT = -dpT;
Real dAdrho = -0.5_rt * delx * g_zone - dpd;
Real dBdT = -dsT;
Real dBdrho = -dsd;

dens_zone =
amrex::max(0.9_rt * dens_zone,
amrex::min(dens_zone + drho, 1.1_rt * dens_zone));
dtemp = (B - (dBdrho / dAdrho) * A) /
((dBdrho / dAdrho) * dAdT - dBdT);

temp_zone =
amrex::max(0.9_rt * temp_zone,
amrex::min(temp_zone + dtemp, 1.1_rt * temp_zone));
drho = -(A + dAdT * dtemp) / dAdrho;

if (std::abs(drho) < TOL * dens_zone && std::abs(dtemp) < TOL * temp_zone) {
converged_hse = true;
break;
}
dens_zone =
amrex::max(0.9_rt * dens_zone,
amrex::min(dens_zone + drho, 1.1_rt * dens_zone));

}
temp_zone =
amrex::max(0.9_rt * temp_zone,
amrex::min(temp_zone + dtemp, 1.1_rt * temp_zone));

if (! converged_hse) {
std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl;
std::cout << "integrate down" << std::endl;
std::cout << "dens_zone, temp_zone = " << dens_zone << " " << temp_zone << std::endl;
std::cout << "p_want = " << p_want << std::endl;
std::cout << "drho = " << drho << std::endl;
amrex::Error("Error: HSE non-convergence");
}
if (std::abs(drho) < TOL * dens_zone && std::abs(dtemp) < TOL * temp_zone) {
converged_hse = true;
break;
}

// call the EOS one more time for this zone and then go on to the next
// (t, rho) -> (p, s)
}

eos_state.T = temp_zone;
eos_state.rho = dens_zone;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = xn[n];
}
if (! converged_hse) {
std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl;
std::cout << "integrate down" << std::endl;
std::cout << "dens_zone, temp_zone = " << dens_zone << " " << temp_zone << std::endl;
std::cout << "p_want = " << p_want << std::endl;
std::cout << "drho = " << drho << std::endl;
amrex::Error("Error: HSE non-convergence");
}

eos_state.aux[AuxZero::iye] = ye;
// call the EOS one more time for this zone and then go on to the next
// (t, rho) -> (p, s)

set_aux(eos_state);
eos_state.T = temp_zone;
eos_state.rho = dens_zone;
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = xn[n];
}

eos(eos_input_rt, eos_state);
eos_state.aux[AuxZero::iye] = ye;

pres_zone = eos_state.p;
set_aux(eos_state);

// update the thermodynamics in this zone
model_mesa_hse(i, model::idens) = dens_zone;
model_mesa_hse(i, model::itemp) = temp_zone;
model_mesa_hse(i, model::ipres) = pres_zone;
model_mesa_hse(i, model::ientr) = eos_state.s;
eos(eos_input_rt, eos_state);

}
pres_zone = eos_state.p;

if (std::abs(model_mesa_hse(0, model::idens) - central_density) < TOL*central_density) {
converged_central_density = true;
break;
}
// update the thermodynamics in this zone
model_mesa_hse(i, model::idens) = dens_zone;
model_mesa_hse(i, model::itemp) = temp_zone;
model_mesa_hse(i, model::ipres) = pres_zone;
model_mesa_hse(i, model::ientr) = eos_state.s;

central_density = model_mesa_hse(0, model::idens);
}

}
if (std::abs(model_mesa_hse(0, model::idens) - central_density) < TOL*central_density) {
converged_central_density = true;
break;
}

if (! converged_central_density) {
amrex::Error("Error: non-convergence of central density");
}
central_density = model_mesa_hse(0, model::idens);

std::cout << "converged central density = " << model_mesa_hse(0, model::idens) << std::endl << std::endl;
}

if (! converged_central_density) {
amrex::Error("Error: non-convergence of central density");
}

std::cout << "converged central density = " << model_mesa_hse(0, model::idens) << std::endl << std::endl;

} else {

// we will integrate from zone 0, so make sure that is
// thermodynamically consistent

eos_state.T = model_mesa_hse(0, model::itemp);
eos_state.rho = model_mesa_hse(0, model::idens);
for (int n = 0; n < NumSpec; ++n) {
eos_state.xn[n] = model_mesa_hse(0, model::ispec+n);
}

eos_state.aux[AuxZero::iye] = model_mesa_hse(0, model::iyef);

set_aux(eos_state);
eos(eos_input_rt, eos_state);

model_mesa_hse(0, model::ipres) = eos_state.p;

}

// compute the full HSE model using our new central density and
// temperature, and the temperature structure as dictated by the
Expand Down

0 comments on commit ca00655

Please sign in to comment.