diff --git a/massive_star/README.md b/massive_star/README.md index 385f2cb..634c002 100644 --- a/massive_star/README.md +++ b/massive_star/README.md @@ -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. diff --git a/massive_star/init_1d.H b/massive_star/init_1d.H index 2bbcca7..0158abc 100644 --- a/massive_star/init_1d.H +++ b/massive_star/init_1d.H @@ -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) { @@ -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 @@ -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