From 6a9d232cf8456ac3803f3c885873deeae257587d Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 8 Oct 2025 13:27:54 -0400 Subject: [PATCH 01/20] add additional model dimension and allow 2d model state --- Exec/Make.Castro | 5 +++-- Util/model_parser/model_parser_data.H | 7 +++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/Exec/Make.Castro b/Exec/Make.Castro index 9f6990c173..927eafcf32 100644 --- a/Exec/Make.Castro +++ b/Exec/Make.Castro @@ -177,10 +177,11 @@ endif MAX_NPTS_MODEL ?= 10000 NUM_MODELS ?= 1 +DIM_MODEL ?= 1 ifeq ($(USE_MODEL_PARSER), TRUE) Bdirs += Util/model_parser - DEFINES += -DMODEL_PARSER -DNPTS_MODEL=$(MAX_NPTS_MODEL) -DNUM_MODELS=$(NUM_MODELS) + DEFINES += -DMODEL_PARSER -DNPTS_MODEL=$(MAX_NPTS_MODEL) -DNUM_MODELS=$(NUM_MODELS) -DDIM_MODEL=$(DIM_MODEL) endif ifeq ($(USE_RNG_STATE_INIT), TRUE) @@ -195,7 +196,7 @@ ifeq ($(USE_MHD), TRUE) endif ifeq ($(USE_GRAV), TRUE) - Bdirs += Source/gravity Source/scf + Bdirs += Source/gravity Source/scf DEFINES += -DGRAVITY USE_MLMG = TRUE endif diff --git a/Util/model_parser/model_parser_data.H b/Util/model_parser/model_parser_data.H index f187921316..46cef7d5c0 100644 --- a/Util/model_parser/model_parser_data.H +++ b/Util/model_parser/model_parser_data.H @@ -30,8 +30,15 @@ namespace model extern AMREX_GPU_MANAGED bool initialized; struct initial_model_t { +#if DIM_MODEL == 1 amrex::Array2D state; amrex::Array1D r; +#elif DIM_MODEL == 2 + amrex::Array3D state; + amrex::Array1D x; + amrex::Array1D y; +#endif }; // Tolerance used for getting the total star mass equal to the desired mass. From 163cd8800d4c7a581f1428aa6d9d3c98748ab134 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 13 Oct 2025 13:07:52 -0400 Subject: [PATCH 02/20] temporary update 2d model --- Exec/science/xrb_spherical/initial_model_2d.H | 557 ++++++++++++++++++ 1 file changed, 557 insertions(+) create mode 100644 Exec/science/xrb_spherical/initial_model_2d.H diff --git a/Exec/science/xrb_spherical/initial_model_2d.H b/Exec/science/xrb_spherical/initial_model_2d.H new file mode 100644 index 0000000000..56a2789d5a --- /dev/null +++ b/Exec/science/xrb_spherical/initial_model_2d.H @@ -0,0 +1,557 @@ +#ifndef INITIAL_MODEL_H +#define INITIAL_MODEL_H + +#include + +// Create a 1-d hydrostatic, atmosphere with an isothermal region +// (T_star) representing the NS, a hyperbolic tangent rise to a +// peak temperature (T_hi) representing the base of an accreted +// layer, an isoentropic profile down to a lower temperature (T_lo), +// and then isothermal. This can serve as an initial model for a +// nova or XRB. +// +// The temperature profile is: +// +// ^ +// T_hi + ^ . +// | / \ . +// | / \ . +// | / . \ . +// T_star +---------+ \ . +// | . . \ . +// | \ . +// | . . \ . +// T_lo + +----------- . +// | . . . +// +---------+---+---------------> r . +// | \ / +// | atm_delta +// |< H_star>| +// +// +// ^ +// | +// +-- dens_base +// +// dens_base is the density at a height H_star -- just below the rise +// in T up to the peak T_hi. The composition is "ash" in the lower +// isothermal region and "fuel" in the isentropic and upper +// isothermal regions. In the transition region, we apply the same +// hyperbolic tangent profile to interpolate the composition. +// +// The fuel and ash compositions are specified by the fuel?_name, +// fuel?_frac and ash?_name, ash?_frac parameters (name of the species +// and mass fraction). Where ? = 1,2,3. +// +// The model is placed into HSE by the following differencing: +// +// (1/dr) [

_i -

_{i-1} ] = (1/2) [ _i + _{i-1} ] g +// +// This will be iterated over in tandem with the EOS call, +// P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1) +// + +// this version allows for multiple initial models + +namespace fw { + constexpr amrex::Real MAX_ITER = 250; + constexpr amrex::Real TOL = 1.e-10_rt; +} + +struct model_t { + + amrex::Real xn_base[NumSpec]; + amrex::Real xn_star[NumSpec]; + amrex::Real xn_perturb[NumSpec]; + + amrex::Real dens_base; + amrex::Real T_star; + amrex::Real T_hi; + amrex::Real T_lo; + + amrex::Real H_star; + amrex::Real atm_delta; + + amrex::Real low_density_cutoff; +}; + + +// Evaluate tanh using the exponential form to workaround a PGI bug on Power9 + +AMREX_INLINE +amrex::Real evaluate_tanh(const amrex::Real z) { + + amrex::Real t; + if (std::abs(z) <= 4.0_rt) { + t = (std::exp(z) - std::exp(-z))/(std::exp(z) + std::exp(-z)); + } else if (z < -4.0_rt) { + t = -1.0_rt; + } else { + t = 1.0_rt; + } + + return t; +} + + +AMREX_INLINE +void +generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amrex::Real xmax, + const int npts_model_y, const amrex::Real ymin, const amrex::Real + const model_t& model_params, const int model_num) +{ + + + // Create a 2-d uniform grid that is identical to the mesh that we are + // mapping onto, and then we want to force it into HSE on that mesh. + + // we actually require that the number of points is the same for each + // model, so we'll just set it each time + + model::npts = npts_model_x; + model::initialized = true; + + if (npts_model_x > NPTS_MODEL) { + amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_MODEL"); + } + if (npts_model_y > NPTS_MODEL) { + amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_MODEL"); + } + + // create the grid -- cell centers + + amrex::Real dx = (xmax - xmin) / static_cast(npts_model_x); + amrex::Real dy = (ymax - ymin) / static_cast(npts_model_y); + + for (int i = 0; i < npts_model_x; i++) { + model::profile(model_num).x(i) = + xmin + (static_cast(i) + 0.5_rt) * dx; + } + + for (int j = 0; j < npts_model_y; j++) { + model::profile(model_num).y(j) = + ymin + (static_cast(j) + 0.5_rt) * dy; + } + + // find the index of the base height + // x is assumed to be the r-dir + + int index_base = -1; + for (int i = 0; i < npts_model_x; i++) { + if (model::profile(model_num).x(i) >= xmin + model_params.H_star) { + index_base = i+1; + break; + } + } + + if (index_base == -1) { + amrex::Error("ERROR: invalid base_height"); + } + + // + // put the model onto our new uniform grid + // This is the initial guess conditions. + // + + // determine the conditions at the base -- this is below the atmosphere + // this point is assumed to be at the rotating pole, θ=0 + + eos_t eos_state; + eos_state.T = model_params.T_star; + eos_state.rho = model_params.dens_base; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model_params.xn_star[n]; + } + + eos(eos_input_rt, eos_state); + + // store the conditions at the base -- we'll use the entropy later + // to constrain the isentropic layer + + amrex::Real pres_base = eos_state.p; + + // set an initial temperature profile and composition for the 2D grid + + for (int i = 0; i < npts_model_x; i++) { + + amrex::Real xc = model::profile(model_num).x(i) - + (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; + + for (int j = 0; j < npts_model_y; j++) { + + // hyperbolic tangent transition: + + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) = + model_params.xn_star[n] + + 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // force them to sum to 1 + + amrex::Real sumX = 0.0_rt; + for (int n = 0; n < NumSpec; n++) { + sumX += model::profile(model_num).state(i, j, model::ispec+n); + } + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) /= sumX; + } + + // temperature profile -- it is constant below the base + // It follows the hyperbolic tangent transition above the base + + if (i <= index_base) { + model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; + } else { + model::profile(model_num).state(i, j, model::itemp) = + model_params.T_star + + 0.5_rt * (model_params.T_hi - model_params.T_star) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // + // We will also have a hyperbolic tangent transition in the theta dir + // + + + + // the density and pressure will be determined via HSE, + // for now, set them to the base conditions + + model::profile(model_num).state(i, j, model::idens) = model_params.dens_base; + model::profile(model_num).state(i, j, model::ipres) = pres_base; + } + } + + // make the base thermodynamics consistent for this base point, + // this is assumed to be on θ=0 -- that is what we will integrate from! + + eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); + eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(index_base, 0, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + model::profile(model_num).state(index_base, 0, model::ipres) = eos_state.p; + + + // + // 2D HSE + entropy solve + // + + bool fluff = false; + bool isentropic = false; + bool flipped = false; + + // Start from the base layer at θ=0 and integrate upward + // We start out isothermal during the hyperbolic transition region + // then we 'flip' to isentropic + // the HSE state will be done putting creating an isentropic state until + // the temperature goes below T_lo -- then we will do isothermal. + // also, once the density goes below low_density_cutoff, we stop HSE + + amrex::Real entropy_base = -1.0_rt; + + for (int i = index_base+1; i < npts_model_x; i++) { + if (model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta && + !flipped) { + // This is to determine when did we get out of the hyperbolic tangent transition phase (isothermal) + // Then get the appropriate entropy for doing isentropic atmosphere + + if (i == index_base + 1) { + amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); + } + isentropic = true; + flipped = true; + + // now we need to know the entropy we are confining ourselves to + + eos_state.rho = model::profile(model_num).state(i-1, model::idens); + eos_state.T = model::profile(model_num).state(i-1, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(i-1, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + // Get entropy for isentropic atmosphere + + entropy_base = eos_state.s; + + amrex::Print() << "base density = " << eos_state.rho + << " and base Temperature " << eos_state.T << std::endl; + } + + HSE_integrate(i, 0); + } + + + // Start from the base layer at θ=0 and integrate upward and integrate to the right + + // Now we go through the intersection regions, and integrate upward. + // Since the integrated result should be the same regardless of the direction, + // we do an average from both direction + + // Start from the base layer at θ=0 and integrate downward + + // Now go through the intersection regions and integrate downward. + +} + +AMREX_INLINE +void HSE_integration(const int i, const int j, + const amrex::Real dx, const amrex::Real dy, + const int idir, + const bool isentropic, bool& fluff, + const amrex::Real entropy_base, + const int index_base, + const model_t& model_params, const int model_num) { + // + // Goal is to update density, temperature, and pressure of the current zone. + // + + // Get current state, this is the initial guess + + amrex::Real dens_zone = model::profile(model_num).state(i, j, model::idens); + amrex::Real temp_zone = model::profile(model_num).state(i, j, model::itemp); + amrex::Real xn[NumSpec]; + for (int n = 0; n < NumSpec; n++) { + xn[n] = model::profile(model_num).state(i, model::ispec+n); + } + + if (fluff) { + // If in fluff region, + // i.e. density of the zone is less than low_density_cutoff + // We stop HSE in this zone and set to a constant density and Temp. + + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + + // Call eos to get pressure: (T, rho) -> (p) + + 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); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, j, model::idens) = dens_zone; + model::profile(model_num).state(i, j, model::itemp) = temp_zone; + model::profile(model_num).state(i, j, model::ipres) = pres_zone; + return; + } + + // + // If we are not in fluff region, then we have to do Newton solve + // + + // start off the Newton loop by saying that the zone has not converged + + bool converged_hse = false; + + amrex::Real pres_zone; + amrex::Real entropy_zone; + amrex::Real p_want; + amrex::Real s_want; + amrex::Real drho; + amrex::Real dtemp = 0; + + int i0; + int j0; + + if (idir == 0) { + // Integrate upward + i0 = -1; + j0 = 0; + } else if (idir == 1) { + // Integrate downward + iX = +1; + jY = 0; + } else if (idir == 2) { + // integrate to the right + iX = 0; + jY = -1; + } else { + amrex::Error("Invalid integration direction!") + } + + for (int iter = 0; iter < fw::MAX_ITER; iter++) { + // Get the pressure from HSE eqn. + // If integrating in radial (x-) dir, then we have: + // + // P_{i,j} = P_{i∓1, j} ± 0.5 (ρ_{i∓1, j} + ρ_{i, j}) dx (g + 0.5 Ω^2 r_{i∓1/2} [1 - cos(2 θ_j)]) + // + // Note, we are using an average of the density of the two zones as + // an approximation of the interface value. + // Also, the reference state depends on the direction of integration. + // e.g. if integrating upward, reference state is below: {i-1, j} + // if integrating downward, reference state is above: {i+1, j} + // + // If integrating in theta (y-) dir, then we have: + // + // P_{i,j} = P_{i, j-1} + 0.25 (ρ_{i, j-1} + ρ_{i, j}) dy Ω^2 r_i sin(2 θ_{j-1/2}) + // + // We assumed integration in the right dir, + // hence reference state is always from left {i, j-1} + + // First just try integrate either in vertical or horizontal direction + + amrex::Real r = model::profile(model_num).x(i) + static_cast(i0) * 0.5_rt * dx; + amrex::Real twoTheta = 2.0_rt * (model::profile(model_num).y(j) + + static_cast(j0) * 0.5_rt * dy); + + if (idir == 0 || idir == 1) { + + // integrating in vertical direction + + p_want = model::profile(model_num).state(i+i0, j, model::ipres) - + static_cast(i0) * 0.5_rt * dx * + (dens_zone + model::profile(model_num).state(i+i0, j, model::idens)) * + (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + + } else if (idir == 2) { + + // integrating in to the right + + p_want = model::profile(model_num).state(i, j+j0, model::ipres) + + 0.25_rt * (dens_zone + model::profile(model_num).state(i, j+j0, model::idens)) * + Omega * Omega * r * r * dy * std::sin(twoTheta); + } + + // + // We can either have isentropic or isothermal case. + // + + if (isentropic) { + // isentropic case + // We have two unknowns: rho and T + // with two constraint eqns: entropy and pressure constraints + + 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); + + entropy_zone = eos_state.s; + pres_zone = eos_state.p; + + Real dpt = eos_state.dpdT; + Real dpd = eos_state.dpdr; + Real dst = eos_state.dsdT; + Real dsd = eos_state.dsdr; + + Real A = p_want - pres_zone; + Real B = entropy_base - entropy_zone; + + dtemp = ((dsd / (dpd - 0.5_rt * dx * gravity::const_grav)) * A - B) / + (dsd * dpt / (dpd - 0.5_rt * dx * gravity::const_grav) - dst); + + drho = (A - dpt * dtemp) / (dpd - 0.5_rt * dx * gravity::const_grav); + + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); + + // check if the density falls below our minimum cut-off -- + // if so, floor it + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } + + } else { + // isothermal case + // Since T is fixed, the only unknown is density + // with one constraint eqn: pressure constraint + + if (model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { + temp_zone = model_params.T_lo; + } + + // (t, rho) -> (p) + + 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); + + entropy_zone = eos_state.s; + pres_zone = eos_state.p; + + Real dpd = eos_state.dpdr; + + drho = (p_want - pres_zone) / (dpd - 0.5_rt * dx * gravity::const_grav); + + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + + if (std::abs(drho) < fw::TOL * dens_zone) { + converged_hse = true; + break; + } + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + } // End of isothermal case + + } // End of Newton solve + + if (!converged_hse) { + std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl; + std::cout << "integrate up" << std::endl; + std::cout << dens_zone << " " << temp_zone << std::endl; + std::cout << p_want << " " << entropy_base << " " << entropy_zone << std::endl; + std::cout << drho << " " << dtemp << std::endl; + amrex::Error("Error: HSE non-convergence"); + } + + // + // Once rho and T are determined, call eos to get pressure + // + + 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); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, model::idens) = dens_zone; + model::profile(model_num).state(i, model::itemp) = temp_zone; + model::profile(model_num).state(i, model::ipres) = pres_zone; +} +#endif From d6be463cec95ce0a2e71a690ef007730287fb2bb Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 16 Oct 2025 16:29:53 -0400 Subject: [PATCH 03/20] temporary update --- Exec/science/xrb_spherical/GNUmakefile | 5 +- Exec/science/xrb_spherical/initial_model.H | 648 +++++++++++++++++- Exec/science/xrb_spherical/initial_model_2d.H | 557 --------------- .../xrb_spherical/problem_initialize.H | 27 +- Source/rotation/Rotation.H | 21 +- Util/model_parser/model_parser_data.H | 20 +- Util/model_parser/model_parser_data.cpp | 7 +- 7 files changed, 695 insertions(+), 590 deletions(-) mode change 120000 => 100644 Exec/science/xrb_spherical/initial_model.H delete mode 100644 Exec/science/xrb_spherical/initial_model_2d.H diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index f48495a81f..24a6684db0 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -20,13 +20,14 @@ CASTRO_HOME ?= ../../.. USE_JACOBIAN_CACHING = TRUE USE_MODEL_PARSER = TRUE -NUM_MODELS := 2 +NUM_MODELS := 1 +DIM_MODEL := 2 # This sets the EOS directory in $(MICROPHYSICS_HOME)/EOS EOS_DIR := helmholtz # This sets the network directory in $(MICROPHYSICS_HOME)/networks -NETWORK_DIR := he-burn/he-burn-18a +NETWORK_DIR := he-burn/ase INTEGRATOR_DIR := RKC diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H deleted file mode 120000 index 9c923c3113..0000000000 --- a/Exec/science/xrb_spherical/initial_model.H +++ /dev/null @@ -1 +0,0 @@ -../flame_wave/initial_model.H \ No newline at end of file diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H new file mode 100644 index 0000000000..dbe45fc4da --- /dev/null +++ b/Exec/science/xrb_spherical/initial_model.H @@ -0,0 +1,647 @@ +#ifndef INITIAL_MODEL_H +#define INITIAL_MODEL_H + +#include + +// Create a 1-d hydrostatic, atmosphere with an isothermal region +// (T_star) representing the NS, a hyperbolic tangent rise to a +// peak temperature (T_hi) representing the base of an accreted +// layer, an isoentropic profile down to a lower temperature (T_lo), +// and then isothermal. This can serve as an initial model for a +// nova or XRB. +// +// The temperature profile is: +// +// ^ +// T_hi + ^ . +// | / \ . +// | / \ . +// | / . \ . +// T_star +---------+ \ . +// | . . \ . +// | \ . +// | . . \ . +// T_lo + +----------- . +// | . . . +// +---------+---+---------------> r . +// | \ / +// | atm_delta +// |< H_star>| +// +// +// ^ +// | +// +-- dens_base +// +// dens_base is the density at a height H_star -- just below the rise +// in T up to the peak T_hi. The composition is "ash" in the lower +// isothermal region and "fuel" in the isentropic and upper +// isothermal regions. In the transition region, we apply the same +// hyperbolic tangent profile to interpolate the composition. +// +// The fuel and ash compositions are specified by the fuel?_name, +// fuel?_frac and ash?_name, ash?_frac parameters (name of the species +// and mass fraction). Where ? = 1,2,3. +// +// The model is placed into HSE by the following differencing: +// +// (1/dr) [

_i -

_{i-1} ] = (1/2) [ _i + _{i-1} ] g +// +// This will be iterated over in tandem with the EOS call, +// P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1) +// + +// this version allows for multiple initial models + +namespace fw { + constexpr amrex::Real MAX_ITER = 250; + constexpr amrex::Real TOL = 1.e-10_rt; +} + +struct model_t { + + amrex::Real xn_base[NumSpec]; + amrex::Real xn_star[NumSpec]; + amrex::Real xn_perturb[NumSpec]; + + amrex::Real dens_base; + amrex::Real T_star; + amrex::Real T_hi; + amrex::Real T_lo; + + amrex::Real H_star; + amrex::Real atm_delta; + + amrex::Real low_density_cutoff; +}; + + +// Evaluate tanh using the exponential form to workaround a PGI bug on Power9 + +AMREX_INLINE +amrex::Real evaluate_tanh(const amrex::Real z) { + + amrex::Real t; + if (std::abs(z) <= 4.0_rt) { + t = (std::exp(z) - std::exp(-z))/(std::exp(z) + std::exp(-z)); + } else if (z < -4.0_rt) { + t = -1.0_rt; + } else { + t = 1.0_rt; + } + + return t; +} + + +AMREX_INLINE +void HSE_integrate(const int i, const int j, + const amrex::Real dx, const amrex::Real dy, + const int idir, + const bool isentropic, bool& fluff, + const amrex::Real entropy_base, + const model_t& model_params, const int model_num) { + // + // Goal is to update density, temperature, and pressure of the current zone. + // Also changes the `fluff` condition of the zone if density drops below + // `low_density_cutoff`. + // + + // Get current state, this is the initial guess + + amrex::Real dens_zone = model::profile(model_num).state(i, j, model::idens); + amrex::Real temp_zone = model::profile(model_num).state(i, j, model::itemp); + amrex::Real xn[NumSpec]; + for (int n = 0; n < NumSpec; n++) { + xn[n] = model::profile(model_num).state(i, j, model::ispec+n); + } + amrex::Real pres_zone; + + if (fluff) { + // If in fluff region, + // i.e. density of the zone is less than low_density_cutoff + // We stop HSE in this zone and set to a constant density and temperature. + + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + + // Call eos to get pressure: (T, rho) -> (p) + + eos_t 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); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, j, model::idens) = dens_zone; + model::profile(model_num).state(i, j, model::itemp) = temp_zone; + model::profile(model_num).state(i, j, model::ipres) = pres_zone; + return; + } + + // + // If we are not in fluff region, then we have to do Newton solve + // + + // start off the Newton loop by saying that the zone has not converged + + bool converged_hse = false; + + amrex::Real entropy_zone; + amrex::Real p_want; + amrex::Real s_want; + amrex::Real drho; + amrex::Real dtemp = 0.0_rt; + + // Set index offset based on integration direction + + int i0; + int j0; + + if (idir == 0) { + // Integrate upward + i0 = -1; + j0 = 0; + } else if (idir == 1) { + // Integrate downward + i0 = +1; + j0 = 0; + } else if (idir == 2) { + // integrate to the right + i0 = 0; + j0 = -1; + } else { + amrex::Error("Invalid integration direction!"); + } + + // Do Newton solve + + for (int iter = 0; iter < fw::MAX_ITER; iter++) { + // Get the pressure from HSE eqn. + // If integrating in radial (x-) dir, then we have: + // + // P_{i,j} = P_{i∓1, j} ± 0.5 (ρ_{i∓1, j} + ρ_{i, j}) dx (g + 0.5 Ω^2 r_{i∓1/2} [1 - cos(2 θ_j)]) + // + // Note, we are using an average of the density of the two zones as + // an approximation of the interface value. + // Also, the reference state depends on the direction of integration. + // e.g. if integrating upward, reference state is below: {i-1, j} + // if integrating downward, reference state is above: {i+1, j} + // + // If integrating in theta (y-) dir, then we have: + // + // P_{i,j} = P_{i, j-1} + 0.25 (ρ_{i, j-1} + ρ_{i, j}) dy Ω^2 r_i sin(2 θ_{j-1/2}) + // + // We assumed integration is always in the right dir, + // hence reference state is always from left {i, j-1} + + // First just try integrate either in vertical or horizontal direction + + amrex::Real r = model::profile(model_num).x(i) + static_cast(i0) * 0.5_rt * dx; + amrex::Real twoTheta = 2.0_rt * (model::profile(model_num).y(j) + + static_cast(j0) * 0.5_rt * dy); + amrex::Real Omega = 0.0_rt; + if (castro::rotational_period > 0.0_rt && + castro::rotation_include_centrifugal) { + Omega = 2.0_rt * M_PI / castro::rotational_period; + } + + if (idir == 0 || idir == 1) { + + // integrating in vertical direction + + p_want = model::profile(model_num).state(i+i0, j, model::ipres) - + static_cast(i0) * 0.5_rt * dx * + (dens_zone + model::profile(model_num).state(i+i0, j, model::idens)) * + (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + + } else if (idir == 2) { + + // integrating in to the right + + p_want = model::profile(model_num).state(i, j+j0, model::ipres) - + static_cast(j0) * 0.25_rt * dy * + (dens_zone + model::profile(model_num).state(i, j+j0, model::idens)) * + Omega * Omega * r * r * std::sin(twoTheta); + } + + // + // We can either have isentropic or isothermal case. + // + + if (isentropic) { + // isentropic case + // We have two unknowns: rho and T + // with two constraint eqns: entropy and pressure constraints + + // Check if we have a valid base entropy when doing isentropic atmosphere + + if (entropy_base <= 0.0_rt) { + amrex::Error("isentropic atmosphere with negative entropy!"); + } + + eos_t 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); + + entropy_zone = eos_state.s; + pres_zone = eos_state.p; + + amrex::Real A = p_want - pres_zone; + amrex::Real B = entropy_base - entropy_zone; + + amrex::Real dAdT = -eos_state.dpdT; + amrex::Real dAdrho = - eos_state.dpdr; + if (idir == 0 || idir == 1) { + // integrate vertically + dAdrho += -static_cast(i0) * 0.5_rt * dx * + (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + } else if (idir == 2) { + // integrate horizontally + dAdrho += -static_cast(j0) * 0.25_rt * dy * + Omega * Omega * r * r * std::sin(twoTheta); + } + amrex::Real dBdT = -eos_state.dsdT; + amrex::Real dBdrho = -eos_state.dsdr; + + dtemp = (B - (dBdrho / dAdrho) * A) / + ((dBdrho / dAdrho) * dAdT - dBdT); + + drho = -(A + dAdT * dtemp) / dAdrho; + + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); + + // check if the density falls below our minimum cut-off -- + // if so, floor it + // Also detect whether we've entered the fluff region. + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } + + } else { + // isothermal case + // Since T is fixed, the only unknown is density + // with one constraint eqn: pressure constraint + + // (t, rho) -> (p) + + eos_t 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); + + entropy_zone = eos_state.s; + pres_zone = eos_state.p; + + amrex::Real A = p_want - pres_zone; + amrex::Real dAdrho = - eos_state.dpdr; + + if (idir == 0 || idir == 1) { + // integrate vertically + dAdrho += -static_cast(i0) * 0.5_rt * dx * + (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + } else if (idir == 2) { + // integrate horizontally + dAdrho += -static_cast(j0) * 0.25_rt * dy * + Omega * Omega * r * r * std::sin(twoTheta); + } + + drho = - A / dAdrho; + + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + if (std::abs(drho) < fw::TOL * dens_zone) { + converged_hse = true; + break; + } + + } // End of isothermal case + + } // End of Newton solve + + if (!converged_hse) { + std::cout << "Error zone " << i << ", " << j << " did not converge in 2D HSE" << std::endl; + std::cout << dens_zone << " " << temp_zone << std::endl; + std::cout << p_want << " " << entropy_base << " " << entropy_zone << std::endl; + std::cout << drho << " " << dtemp << std::endl; + amrex::Error("Error: HSE non-convergence"); + } + + // + // Once rho and T are determined, call eos to get pressure + // + + eos_t 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); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, j, model::idens) = dens_zone; + model::profile(model_num).state(i, j, model::itemp) = temp_zone; + model::profile(model_num).state(i, j, model::ipres) = pres_zone; +} + + +AMREX_INLINE +void +generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amrex::Real xmax, + const int npts_model_y, const amrex::Real ymin, const amrex::Real ymax, + const model_t& model_params, const int model_num) +{ + + + // Create a 2-d uniform grid that is identical to the mesh that we are + // mapping onto, and then we want to force it into HSE on that mesh. + + // we actually require that the number of points is the same for each + // model, so we'll just set it each time + + model::npts = npts_model_x < npts_model_y ? npts_model_x : npts_model_y; + model::initialized = true; + + if (npts_model_x > NPTS_MODEL) { + amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_MODEL"); + } + if (npts_model_y > NPTS_MODEL) { + amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_MODEL"); + } + + // create the grid -- cell centers + + amrex::Real dx = (xmax - xmin) / static_cast(npts_model_x); + amrex::Real dy = (ymax - ymin) / static_cast(npts_model_y); + + for (int i = 0; i < npts_model_x; i++) { + model::profile(model_num).x(i) = + xmin + (static_cast(i) + 0.5_rt) * dx; + } + + for (int j = 0; j < npts_model_y; j++) { + model::profile(model_num).y(j) = + ymin + (static_cast(j) + 0.5_rt) * dy; + } + + // find the index of the base height + // x is assumed to be the r-dir + + int index_base = -1; + for (int i = 0; i < npts_model_x; i++) { + if (model::profile(model_num).x(i) >= xmin + model_params.H_star) { + index_base = i+1; + break; + } + } + + if (index_base == -1) { + amrex::Error("ERROR: invalid base_height"); + } + + // + // put the model onto our new uniform grid + // This is the initial guess conditions. + // + + // determine the conditions at the base -- this is below the atmosphere + // this point is assumed to be at the rotating pole, θ=0 + + eos_t eos_state; + eos_state.T = model_params.T_star; + eos_state.rho = model_params.dens_base; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model_params.xn_star[n]; + } + + eos(eos_input_rt, eos_state); + + // store the conditions at the base -- we'll use the entropy later + // to constrain the isentropic layer + + amrex::Real pres_base = eos_state.p; + + // set an initial temperature profile and composition for the 2D grid + + for (int i = 0; i < npts_model_x; i++) { + + amrex::Real xc = model::profile(model_num).x(i) - + (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; + + for (int j = 0; j < npts_model_y; j++) { + + // vertical hyperbolic tangent transition + + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) = + model_params.xn_star[n] + + 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // force them to sum to 1 + + amrex::Real sumX = 0.0_rt; + for (int n = 0; n < NumSpec; n++) { + sumX += model::profile(model_num).state(i, j, model::ispec+n); + } + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) /= sumX; + } + + // temperature profile -- it is constant below the base + // It follows the hyperbolic tangent transition above the base + + if (i <= index_base) { + model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; + } else { + model::profile(model_num).state(i, j, model::itemp) = + model_params.T_star + + 0.5_rt * (model_params.T_hi - model_params.T_star) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // + // We will also have a hyperbolic tangent transition in the theta dir + // but for now just assume theres is no hot/cold model. + // + + + + // the density and pressure will be determined via HSE, + // for now, set them to the base conditions + + model::profile(model_num).state(i, j, model::idens) = model_params.dens_base; + model::profile(model_num).state(i, j, model::ipres) = pres_base; + } + } + + // make the base thermodynamics consistent for this base point, + // this is assumed to be on θ=0 -- that is what we will integrate from! + + eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); + eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(index_base, 0, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + model::profile(model_num).state(index_base, 0, model::ipres) = eos_state.p; + + + // + // 2D HSE + entropy solve + // + + bool fluff = false; + bool isentropic = false; + bool flipped = false; + amrex::Real entropy_base = -1.0_rt; + + // First integrate to the right from θ=0 at the base layer. + // The base layer is isothermal. + + for (int j = 0; j < npts_model_y; ++j) { + HSE_integrate(index_base, j, dx, dy, 2, + isentropic, fluff, + entropy_base, + model_params, model_num); + } + + // Then integrate vertically for all θ. + + for (int j = 0; j < npts_model_y; ++j) { + + // Integrate vertically downward from base + // Regions below the base is isothermal + + // Reset conditions for each θ + + fluff = false; + isentropic = false; + flipped = false; + entropy_base = -1.0_rt; + + for (int i = index_base-1; i >= 0; --i) { + HSE_integrate(i, j, dx, dy, 1, + isentropic, fluff, + entropy_base, + model_params, model_num); + } + + // Start from the base layer and integrate upward + // We start out isothermal during the hyperbolic transition region + // then we 'flip' to isentropic + // the HSE state will be done putting creating an isentropic state until + // the temperature goes below T_lo -- then we will do isothermal. + // also, once the density goes below low_density_cutoff, we stop HSE + + // Reset conditions for each θ + + fluff = false; + isentropic = false; + flipped = false; + entropy_base = -1.0_rt; + + for (int i = index_base+1; i < npts_model_x; i++) { + + // If above hyperbolic transition region, we can either be isentropic or isothermal, + + if (model::profile(model_num).x(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { + + if (flipped && !isentropic) { + model::profile(model_num).state(i, j, model::itemp) = model_params.T_lo; + + } else if (!flipped) { + // This is to determine when did we get out of the hyperbolic tangent transition phase (isothermal) + // Then get the appropriate entropy for doing isentropic atmosphere + if (i == index_base + 1) { + amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); + } + isentropic = true; + flipped = true; + + // now we need to know the entropy we are confining ourselves to + + eos_state.rho = model::profile(model_num).state(i-1, j, model::idens); + eos_state.T = model::profile(model_num).state(i-1, j, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(i-1, j, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + // Get entropy for isentropic atmosphere + + entropy_base = eos_state.s; + + amrex::Print() << "base density = " << eos_state.rho + << " and base Temperature " << eos_state.T << std::endl; + } + } + HSE_integrate(i, j, dx, dy, 0, + isentropic, fluff, + entropy_base, + model_params, model_num); + } + + } + + + + + + // Start from the base layer at θ=0 and integrate upward and integrate to the right + + // Now we go through the intersection regions, and integrate upward. + // Since the integrated result should be the same regardless of the direction, + // we do an average from both direction + + // Start from the base layer at θ=0 and integrate downward + + // Now go through the intersection regions and integrate downward. + +} +#endif diff --git a/Exec/science/xrb_spherical/initial_model_2d.H b/Exec/science/xrb_spherical/initial_model_2d.H deleted file mode 100644 index 56a2789d5a..0000000000 --- a/Exec/science/xrb_spherical/initial_model_2d.H +++ /dev/null @@ -1,557 +0,0 @@ -#ifndef INITIAL_MODEL_H -#define INITIAL_MODEL_H - -#include - -// Create a 1-d hydrostatic, atmosphere with an isothermal region -// (T_star) representing the NS, a hyperbolic tangent rise to a -// peak temperature (T_hi) representing the base of an accreted -// layer, an isoentropic profile down to a lower temperature (T_lo), -// and then isothermal. This can serve as an initial model for a -// nova or XRB. -// -// The temperature profile is: -// -// ^ -// T_hi + ^ . -// | / \ . -// | / \ . -// | / . \ . -// T_star +---------+ \ . -// | . . \ . -// | \ . -// | . . \ . -// T_lo + +----------- . -// | . . . -// +---------+---+---------------> r . -// | \ / -// | atm_delta -// |< H_star>| -// -// -// ^ -// | -// +-- dens_base -// -// dens_base is the density at a height H_star -- just below the rise -// in T up to the peak T_hi. The composition is "ash" in the lower -// isothermal region and "fuel" in the isentropic and upper -// isothermal regions. In the transition region, we apply the same -// hyperbolic tangent profile to interpolate the composition. -// -// The fuel and ash compositions are specified by the fuel?_name, -// fuel?_frac and ash?_name, ash?_frac parameters (name of the species -// and mass fraction). Where ? = 1,2,3. -// -// The model is placed into HSE by the following differencing: -// -// (1/dr) [

_i -

_{i-1} ] = (1/2) [ _i + _{i-1} ] g -// -// This will be iterated over in tandem with the EOS call, -// P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1) -// - -// this version allows for multiple initial models - -namespace fw { - constexpr amrex::Real MAX_ITER = 250; - constexpr amrex::Real TOL = 1.e-10_rt; -} - -struct model_t { - - amrex::Real xn_base[NumSpec]; - amrex::Real xn_star[NumSpec]; - amrex::Real xn_perturb[NumSpec]; - - amrex::Real dens_base; - amrex::Real T_star; - amrex::Real T_hi; - amrex::Real T_lo; - - amrex::Real H_star; - amrex::Real atm_delta; - - amrex::Real low_density_cutoff; -}; - - -// Evaluate tanh using the exponential form to workaround a PGI bug on Power9 - -AMREX_INLINE -amrex::Real evaluate_tanh(const amrex::Real z) { - - amrex::Real t; - if (std::abs(z) <= 4.0_rt) { - t = (std::exp(z) - std::exp(-z))/(std::exp(z) + std::exp(-z)); - } else if (z < -4.0_rt) { - t = -1.0_rt; - } else { - t = 1.0_rt; - } - - return t; -} - - -AMREX_INLINE -void -generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amrex::Real xmax, - const int npts_model_y, const amrex::Real ymin, const amrex::Real - const model_t& model_params, const int model_num) -{ - - - // Create a 2-d uniform grid that is identical to the mesh that we are - // mapping onto, and then we want to force it into HSE on that mesh. - - // we actually require that the number of points is the same for each - // model, so we'll just set it each time - - model::npts = npts_model_x; - model::initialized = true; - - if (npts_model_x > NPTS_MODEL) { - amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_MODEL"); - } - if (npts_model_y > NPTS_MODEL) { - amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_MODEL"); - } - - // create the grid -- cell centers - - amrex::Real dx = (xmax - xmin) / static_cast(npts_model_x); - amrex::Real dy = (ymax - ymin) / static_cast(npts_model_y); - - for (int i = 0; i < npts_model_x; i++) { - model::profile(model_num).x(i) = - xmin + (static_cast(i) + 0.5_rt) * dx; - } - - for (int j = 0; j < npts_model_y; j++) { - model::profile(model_num).y(j) = - ymin + (static_cast(j) + 0.5_rt) * dy; - } - - // find the index of the base height - // x is assumed to be the r-dir - - int index_base = -1; - for (int i = 0; i < npts_model_x; i++) { - if (model::profile(model_num).x(i) >= xmin + model_params.H_star) { - index_base = i+1; - break; - } - } - - if (index_base == -1) { - amrex::Error("ERROR: invalid base_height"); - } - - // - // put the model onto our new uniform grid - // This is the initial guess conditions. - // - - // determine the conditions at the base -- this is below the atmosphere - // this point is assumed to be at the rotating pole, θ=0 - - eos_t eos_state; - eos_state.T = model_params.T_star; - eos_state.rho = model_params.dens_base; - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model_params.xn_star[n]; - } - - eos(eos_input_rt, eos_state); - - // store the conditions at the base -- we'll use the entropy later - // to constrain the isentropic layer - - amrex::Real pres_base = eos_state.p; - - // set an initial temperature profile and composition for the 2D grid - - for (int i = 0; i < npts_model_x; i++) { - - amrex::Real xc = model::profile(model_num).x(i) - - (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; - - for (int j = 0; j < npts_model_y; j++) { - - // hyperbolic tangent transition: - - for (int n = 0; n < NumSpec; n++) { - model::profile(model_num).state(i, j, model::ispec+n) = - model_params.xn_star[n] + - 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * - (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); - } - - // force them to sum to 1 - - amrex::Real sumX = 0.0_rt; - for (int n = 0; n < NumSpec; n++) { - sumX += model::profile(model_num).state(i, j, model::ispec+n); - } - for (int n = 0; n < NumSpec; n++) { - model::profile(model_num).state(i, j, model::ispec+n) /= sumX; - } - - // temperature profile -- it is constant below the base - // It follows the hyperbolic tangent transition above the base - - if (i <= index_base) { - model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; - } else { - model::profile(model_num).state(i, j, model::itemp) = - model_params.T_star + - 0.5_rt * (model_params.T_hi - model_params.T_star) * - (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); - } - - // - // We will also have a hyperbolic tangent transition in the theta dir - // - - - - // the density and pressure will be determined via HSE, - // for now, set them to the base conditions - - model::profile(model_num).state(i, j, model::idens) = model_params.dens_base; - model::profile(model_num).state(i, j, model::ipres) = pres_base; - } - } - - // make the base thermodynamics consistent for this base point, - // this is assumed to be on θ=0 -- that is what we will integrate from! - - eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); - eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model::profile(model_num).state(index_base, 0, model::ispec+n); - } - - eos(eos_input_rt, eos_state); - - model::profile(model_num).state(index_base, 0, model::ipres) = eos_state.p; - - - // - // 2D HSE + entropy solve - // - - bool fluff = false; - bool isentropic = false; - bool flipped = false; - - // Start from the base layer at θ=0 and integrate upward - // We start out isothermal during the hyperbolic transition region - // then we 'flip' to isentropic - // the HSE state will be done putting creating an isentropic state until - // the temperature goes below T_lo -- then we will do isothermal. - // also, once the density goes below low_density_cutoff, we stop HSE - - amrex::Real entropy_base = -1.0_rt; - - for (int i = index_base+1; i < npts_model_x; i++) { - if (model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta && - !flipped) { - // This is to determine when did we get out of the hyperbolic tangent transition phase (isothermal) - // Then get the appropriate entropy for doing isentropic atmosphere - - if (i == index_base + 1) { - amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); - } - isentropic = true; - flipped = true; - - // now we need to know the entropy we are confining ourselves to - - eos_state.rho = model::profile(model_num).state(i-1, model::idens); - eos_state.T = model::profile(model_num).state(i-1, model::itemp); - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model::profile(model_num).state(i-1, model::ispec+n); - } - - eos(eos_input_rt, eos_state); - - // Get entropy for isentropic atmosphere - - entropy_base = eos_state.s; - - amrex::Print() << "base density = " << eos_state.rho - << " and base Temperature " << eos_state.T << std::endl; - } - - HSE_integrate(i, 0); - } - - - // Start from the base layer at θ=0 and integrate upward and integrate to the right - - // Now we go through the intersection regions, and integrate upward. - // Since the integrated result should be the same regardless of the direction, - // we do an average from both direction - - // Start from the base layer at θ=0 and integrate downward - - // Now go through the intersection regions and integrate downward. - -} - -AMREX_INLINE -void HSE_integration(const int i, const int j, - const amrex::Real dx, const amrex::Real dy, - const int idir, - const bool isentropic, bool& fluff, - const amrex::Real entropy_base, - const int index_base, - const model_t& model_params, const int model_num) { - // - // Goal is to update density, temperature, and pressure of the current zone. - // - - // Get current state, this is the initial guess - - amrex::Real dens_zone = model::profile(model_num).state(i, j, model::idens); - amrex::Real temp_zone = model::profile(model_num).state(i, j, model::itemp); - amrex::Real xn[NumSpec]; - for (int n = 0; n < NumSpec; n++) { - xn[n] = model::profile(model_num).state(i, model::ispec+n); - } - - if (fluff) { - // If in fluff region, - // i.e. density of the zone is less than low_density_cutoff - // We stop HSE in this zone and set to a constant density and Temp. - - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; - - // Call eos to get pressure: (T, rho) -> (p) - - 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); - - pres_zone = eos_state.p; - - // update the thermodynamics in this zone - - model::profile(model_num).state(i, j, model::idens) = dens_zone; - model::profile(model_num).state(i, j, model::itemp) = temp_zone; - model::profile(model_num).state(i, j, model::ipres) = pres_zone; - return; - } - - // - // If we are not in fluff region, then we have to do Newton solve - // - - // start off the Newton loop by saying that the zone has not converged - - bool converged_hse = false; - - amrex::Real pres_zone; - amrex::Real entropy_zone; - amrex::Real p_want; - amrex::Real s_want; - amrex::Real drho; - amrex::Real dtemp = 0; - - int i0; - int j0; - - if (idir == 0) { - // Integrate upward - i0 = -1; - j0 = 0; - } else if (idir == 1) { - // Integrate downward - iX = +1; - jY = 0; - } else if (idir == 2) { - // integrate to the right - iX = 0; - jY = -1; - } else { - amrex::Error("Invalid integration direction!") - } - - for (int iter = 0; iter < fw::MAX_ITER; iter++) { - // Get the pressure from HSE eqn. - // If integrating in radial (x-) dir, then we have: - // - // P_{i,j} = P_{i∓1, j} ± 0.5 (ρ_{i∓1, j} + ρ_{i, j}) dx (g + 0.5 Ω^2 r_{i∓1/2} [1 - cos(2 θ_j)]) - // - // Note, we are using an average of the density of the two zones as - // an approximation of the interface value. - // Also, the reference state depends on the direction of integration. - // e.g. if integrating upward, reference state is below: {i-1, j} - // if integrating downward, reference state is above: {i+1, j} - // - // If integrating in theta (y-) dir, then we have: - // - // P_{i,j} = P_{i, j-1} + 0.25 (ρ_{i, j-1} + ρ_{i, j}) dy Ω^2 r_i sin(2 θ_{j-1/2}) - // - // We assumed integration in the right dir, - // hence reference state is always from left {i, j-1} - - // First just try integrate either in vertical or horizontal direction - - amrex::Real r = model::profile(model_num).x(i) + static_cast(i0) * 0.5_rt * dx; - amrex::Real twoTheta = 2.0_rt * (model::profile(model_num).y(j) + - static_cast(j0) * 0.5_rt * dy); - - if (idir == 0 || idir == 1) { - - // integrating in vertical direction - - p_want = model::profile(model_num).state(i+i0, j, model::ipres) - - static_cast(i0) * 0.5_rt * dx * - (dens_zone + model::profile(model_num).state(i+i0, j, model::idens)) * - (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); - - } else if (idir == 2) { - - // integrating in to the right - - p_want = model::profile(model_num).state(i, j+j0, model::ipres) + - 0.25_rt * (dens_zone + model::profile(model_num).state(i, j+j0, model::idens)) * - Omega * Omega * r * r * dy * std::sin(twoTheta); - } - - // - // We can either have isentropic or isothermal case. - // - - if (isentropic) { - // isentropic case - // We have two unknowns: rho and T - // with two constraint eqns: entropy and pressure constraints - - 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); - - entropy_zone = eos_state.s; - pres_zone = eos_state.p; - - Real dpt = eos_state.dpdT; - Real dpd = eos_state.dpdr; - Real dst = eos_state.dsdT; - Real dsd = eos_state.dsdr; - - Real A = p_want - pres_zone; - Real B = entropy_base - entropy_zone; - - dtemp = ((dsd / (dpd - 0.5_rt * dx * gravity::const_grav)) * A - B) / - (dsd * dpt / (dpd - 0.5_rt * dx * gravity::const_grav) - dst); - - drho = (A - dpt * dtemp) / (dpd - 0.5_rt * dx * gravity::const_grav); - - dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); - temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); - - // check if the density falls below our minimum cut-off -- - // if so, floor it - - if (dens_zone < model_params.low_density_cutoff) { - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; - converged_hse = true; - fluff = true; - break; - } - - // if (A < TOL .and. B < ETOL) then - if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { - converged_hse = true; - break; - } - - } else { - // isothermal case - // Since T is fixed, the only unknown is density - // with one constraint eqn: pressure constraint - - if (model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { - temp_zone = model_params.T_lo; - } - - // (t, rho) -> (p) - - 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); - - entropy_zone = eos_state.s; - pres_zone = eos_state.p; - - Real dpd = eos_state.dpdr; - - drho = (p_want - pres_zone) / (dpd - 0.5_rt * dx * gravity::const_grav); - - dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); - - if (std::abs(drho) < fw::TOL * dens_zone) { - converged_hse = true; - break; - } - - if (dens_zone < model_params.low_density_cutoff) { - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; - converged_hse = true; - fluff = true; - break; - } - - } // End of isothermal case - - } // End of Newton solve - - if (!converged_hse) { - std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl; - std::cout << "integrate up" << std::endl; - std::cout << dens_zone << " " << temp_zone << std::endl; - std::cout << p_want << " " << entropy_base << " " << entropy_zone << std::endl; - std::cout << drho << " " << dtemp << std::endl; - amrex::Error("Error: HSE non-convergence"); - } - - // - // Once rho and T are determined, call eos to get pressure - // - - 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); - - pres_zone = eos_state.p; - - // update the thermodynamics in this zone - - model::profile(model_num).state(i, model::idens) = dens_zone; - model::profile(model_num).state(i, model::itemp) = temp_zone; - model::profile(model_num).state(i, model::ipres) = pres_zone; -} -#endif diff --git a/Exec/science/xrb_spherical/problem_initialize.H b/Exec/science/xrb_spherical/problem_initialize.H index cfb11741ff..e6dfea4f61 100644 --- a/Exec/science/xrb_spherical/problem_initialize.H +++ b/Exec/science/xrb_spherical/problem_initialize.H @@ -151,10 +151,12 @@ void problem_initialize () auto dx = fine_geom.CellSizeArray(); auto dx_model = dx[0]; + auto dy_model = dx[1]; int nx_model = static_cast((probhi[0] - problo[0]) / dx_model); - + int ny_model = static_cast((probhi[1] - problo[1]) / + dy_model); int ng = 4; // now generate the initial models @@ -172,32 +174,35 @@ void problem_initialize () generate_initial_model(nx_model + ng, problo[0] - ng * dx_model, probhi[0], + ny_model + ng, + problo[1] - ng * dy_model, + probhi[1], model_params, 0); // now create a perturbed model -- we want the same base conditions // a hotter temperature - model_params.T_hi = model_params.T_hi + problem::dtemp; + // model_params.T_hi = model_params.T_hi + problem::dtemp; - generate_initial_model(nx_model + ng, - problo[0] - ng * dx_model, - probhi[0], - model_params, 1); + // generate_initial_model(nx_model + ng, + // problo[0] - ng * dx_model, + // probhi[0], + // model_params, 1); // set center for (int d = 0; d < AMREX_SPACEDIM; d++) { - // problem::center[d] = 0.5_rt * (problo[d] + probhi[d]); problem::center[d] = 0.0_rt; } - // set the ambient state for the upper boundary condition + // set the ambient state for the upper boundary condition using upper most r zone at θ=0 + + ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, ng, model::idens); - ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, model::idens); - ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, model::itemp); + ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, ng, model::itemp); for (int n = 0; n < NumSpec; n++) { ambient::ambient_state[UFS+n] = - ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, model::ispec+n); + ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, ng, model::ispec+n); } ambient::ambient_state[UMX] = 0.0_rt; diff --git a/Source/rotation/Rotation.H b/Source/rotation/Rotation.H index 1585b71ed8..d7da5ffafe 100644 --- a/Source/rotation/Rotation.H +++ b/Source/rotation/Rotation.H @@ -93,15 +93,12 @@ rotational_acceleration(const GpuArray& r, const GpuArray& v, bool c2 = (castro::rotation_include_coriolis == 1 && coriolis) ? true : false; - auto r_vec = r; - - GpuArray omega_cross_v; - cross_product(omega, v, omega_cross_v); - if (c1) { // For Spherical coordinate, the physical position vector is just r rhat, // but the input r will be in the format of (r,theta,0). So change it to (r,0,0) + auto r_vec = r; + if (coord == 2) { r_vec[1] = 0.0_rt; r_vec[2] = 0.0_rt; @@ -118,6 +115,8 @@ rotational_acceleration(const GpuArray& r, const GpuArray& v, } if (c2) { + GpuArray omega_cross_v; + cross_product(omega, v, omega_cross_v); for (int idir = 0; idir < 3; idir++) { Sr[idir] -= 2.0_rt * omega_cross_v[idir]; } @@ -134,15 +133,15 @@ rotational_potential(const GpuArray& r, const GpuArray& omega, // Real phi = 0.0_rt; - auto r_vec = r; - - if (coord == 2) { - r_vec[1] = 0.0_rt; - r_vec[2] = 0.0_rt; - } if (castro::rotation_include_centrifugal == 1) { + auto r_vec = r; + if (coord == 2) { + r_vec[1] = 0.0_rt; + r_vec[2] = 0.0_rt; + } + GpuArray omega_cross_r; cross_product(omega, r, omega_cross_r); diff --git a/Util/model_parser/model_parser_data.H b/Util/model_parser/model_parser_data.H index 46cef7d5c0..dedb98876e 100644 --- a/Util/model_parser/model_parser_data.H +++ b/Util/model_parser/model_parser_data.H @@ -26,18 +26,24 @@ namespace model constexpr int iaux = -1; #endif - extern AMREX_GPU_MANAGED int npts; +#if DIM_MODEL >= 1 + extern AMREX_GPU_MANAGED int npts_x; +#endif +#if DIM_MODEL >= 2 + extern AMREX_GPU_MANAGED int npts_y; +#endif + extern AMREX_GPU_MANAGED bool initialized; struct initial_model_t { #if DIM_MODEL == 1 - amrex::Array2D state; - amrex::Array1D r; + amrex::Array2D state; + amrex::Array1D x; #elif DIM_MODEL == 2 - amrex::Array3D state; - amrex::Array1D x; - amrex::Array1D y; + amrex::Array3D state; + amrex::Array1D x; + amrex::Array1D y; #endif }; diff --git a/Util/model_parser/model_parser_data.cpp b/Util/model_parser/model_parser_data.cpp index cfbc475cbe..b83158c19a 100644 --- a/Util/model_parser/model_parser_data.cpp +++ b/Util/model_parser/model_parser_data.cpp @@ -4,7 +4,12 @@ namespace model { - AMREX_GPU_MANAGED int npts; +#if DIM_MODEL >= 1 + AMREX_GPU_MANAGED int npts_x; +#endif +#if DIM_MODEL >= 2 + AMREX_GPU_MANAGED int npts_y; +#endif AMREX_GPU_MANAGED bool initialized; AMREX_GPU_MANAGED amrex::Array1D profile; From 7b1679881fa1570b7a91b50ebca59e8b6f182a18 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 18:19:09 -0400 Subject: [PATCH 04/20] update initial model to work with 2d --- Exec/Make.Castro | 28 +- Util/model_parser/model_parser.H | 349 ++++++++++++++++++++---- Util/model_parser/model_parser_data.H | 12 + Util/model_parser/model_parser_data.cpp | 6 +- 4 files changed, 336 insertions(+), 59 deletions(-) diff --git a/Exec/Make.Castro b/Exec/Make.Castro index 9f6990c173..777b8fbac2 100644 --- a/Exec/Make.Castro +++ b/Exec/Make.Castro @@ -175,14 +175,36 @@ ifeq ($(USE_ALL_CASTRO), TRUE) Source/problems Source/sources endif -MAX_NPTS_MODEL ?= 10000 NUM_MODELS ?= 1 +DIM_MODEL ?= 1 + +# Only 1 or 2 dimensions are supported for now. +ifeq ($(filter $(DIM_MODEL),1 2),) + $(error Invalid DIM_MODEL=$(DIM_MODEL). Must be 1 or 2.) +endif + +ifeq ($(DIM_MODEL), 1) + MAX_NPTS_MODEL ?= 10000 +else ifeq ($(DIM_MODEL), 2) + MAX_NPTS_X_MODEL ?= 10000 + MAX_NPTS_Y_MODEL ?= 10000 +endif ifeq ($(USE_MODEL_PARSER), TRUE) Bdirs += Util/model_parser - DEFINES += -DMODEL_PARSER -DNPTS_MODEL=$(MAX_NPTS_MODEL) -DNUM_MODELS=$(NUM_MODELS) + DEFINES += -DMODEL_PARSER \ + -DNUM_MODELS=$(NUM_MODELS) \ + -DDIM_MODEL=$(DIM_MODEL) + + ifeq ($(DIM_MODEL), 1) + DEFINES += -DNPTS_MODEL=$(MAX_NPTS_MODEL) + else ifeq ($(DIM_MODEL), 2) + DEFINES += -DNPTS_X_MODEL=$(MAX_NPTS_X_MODEL) \ + -DNPTS_Y_MODEL=$(MAX_NPTS_Y_MODEL) + endif endif + ifeq ($(USE_RNG_STATE_INIT), TRUE) DEFINES += -DRNG_STATE_INIT endif @@ -195,7 +217,7 @@ ifeq ($(USE_MHD), TRUE) endif ifeq ($(USE_GRAV), TRUE) - Bdirs += Source/gravity Source/scf + Bdirs += Source/gravity Source/scf DEFINES += -DGRAVITY USE_MLMG = TRUE endif diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index b96e4fd86f..098889068f 100644 --- a/Util/model_parser/model_parser.H +++ b/Util/model_parser/model_parser.H @@ -37,6 +37,7 @@ using namespace amrex; /// Note that the first data column is the coordinate, but that is not /// listed as a variable. /// +/// /// The new format has a single line header, labeling the columns, /// including the coordinate. This format plays nicely with NumPy /// genfromtxt to read in the header @@ -52,10 +53,19 @@ using namespace amrex; /// density, temperature, pressure and composition. /// /// composition is assumed to be in terms of mass fractions +/// +/// For 2D initial models the second column data is assumed to be +/// the other coordinate. It is assumed that the first coordinate +/// repeats and the cycles through the second coordinate. i.e. +/// +/// 195312.5000 140523.5000 5437711139. 8805500.952 .4695704813E+28 0.3 0.7 0 +/// 195312.5000 531148.5000 5410152416. 8816689.836 0.4663923963E+28 0.3 0.7 0 +/// ... 921773.5000 5410153244. 8826843.234 0.4621245962E+28 0.3 0.7 0 +/// 585937.5000 140523.5000 5423405349. 8813452.123 .4682340123E+28 0.3 0.7 0 +/// // remove whitespace -- from stackoverflow - namespace model_string { inline std::string& ltrim(std::string& s) @@ -79,7 +89,7 @@ namespace model_string } } - +#if DIM_MODEL == 1 /// /// return the index into the model coordinate, loc, such that /// model::profile(model_index).r(loc) < r < model::profile(model_index).r(loc+1) @@ -89,20 +99,20 @@ namespace model_string /// AMREX_INLINE AMREX_GPU_HOST_DEVICE int -locate(const Real r, const int model_index) { +locate(const amrex::Real r, const int model_index) { int loc; if (r <= model::profile(model_index).r(0)) { loc = 0; - } else if (r > model::profile(model_index).r(model::npts-2)) { - loc = model::npts-2; + } else if (r > model::profile(model_index).r(model::npts - 2)) { + loc = model::npts - 2; } else { int ilo = 0; - int ihi = model::npts-2; + int ihi = model::npts - 2; while (ilo+1 != ihi) { int imid = (ilo + ihi) / 2; @@ -120,17 +130,88 @@ locate(const Real r, const int model_index) { return loc; } +#elif DIM_MODEL == 2 +AMREX_INLINE AMREX_GPU_HOST_DEVICE +int +locate_x(const amrex::Real x, const int model_index) { + int loc; + + if (x <= model::profile(model_index).x(0)) { + loc = 0; + + } else if (x > model::profile(model_index).x(model::npts_x - 2)) { + loc = model::npts_x - 2; + + } else { + + int ilo = 0; + int ihi = model::npts_x - 2; + + while (ilo+1 != ihi) { + int imid = (ilo + ihi) / 2; + + if (x <= model::profile(model_index).x(imid)) { + ihi = imid; + } else { + ilo = imid; + } + } + + loc = ilo; + } + + return loc; +} + +AMREX_INLINE AMREX_GPU_HOST_DEVICE +int +locate_y(const amrex::Real y, const int model_index) { + + int loc; + + if (y <= model::profile(model_index).y(0)) { + loc = 0; + + } else if (y > model::profile(model_index).y(model::npts_y - 2)) { + loc = model::npts_y - 2; + } else { + + int ilo = 0; + int ihi = model::npts_y - 2; + + while (ilo+1 != ihi) { + int imid = (ilo + ihi) / 2; + + if (x <= model::profile(model_index).y(imid)) { + ihi = imid; + } else { + ilo = imid; + } + } + + loc = ilo; + } + + return loc; +} +#endif + + +#if DIM_MODEL == 1 +// +// Linear 1D interpolation +// AMREX_INLINE AMREX_GPU_HOST_DEVICE -Real -interpolate(const Real r, const int var_index, const int model_index=0) { +amrex::Real +interpolate(const amrex::Real r, const int var_index, const int model_index=0) { // this gives us an index such that profile.r(id) < r < profile.r(id+1) int id = locate(r, model_index); - Real slope; - Real interp; + amrex::Real slope; + amrex::Real interp; slope = (model::profile(model_index).state(id+1, var_index) - model::profile(model_index).state(id, var_index)) / @@ -141,10 +222,10 @@ interpolate(const Real r, const int var_index, const int model_index=0) { // safety check to make sure interp lies within the bounding points. We don't // do this at the lower boundary, which usually corresponds to the center of the star. if (r >= model::profile(model_index).r(0)) { - Real minvar = std::min(model::profile(model_index).state(id+1, var_index), - model::profile(model_index).state(id, var_index)); - Real maxvar = std::max(model::profile(model_index).state(id+1, var_index), - model::profile(model_index).state(id, var_index)); + amrex::Real minvar = std::min(model::profile(model_index).state(id+1, var_index), + model::profile(model_index).state(id, var_index)); + amrex::Real maxvar = std::max(model::profile(model_index).state(id+1, var_index), + model::profile(model_index).state(id, var_index)); interp = amrex::Clamp(interp, minvar, maxvar); } @@ -152,13 +233,83 @@ interpolate(const Real r, const int var_index, const int model_index=0) { } +#elif DIM_MODEL == 2 +// +// Linear 2D interpolation +// +AMREX_INLINE AMREX_GPU_HOST_DEVICE +amrex::Real +interpolate(const amrex::Real x, const amrex::Real y, + const int var_index, const int model_index=0) { + + // this gives us an index such that profile.x(id) < x < profile.x(id+1) + + int idx = locate_x(x, model_index); + int idy = locate_y(y, model_index); + + amrex::Real dx = model::profile(model_index).x(idx+1) - + model::profile(model_index).x(idx); + + amrex::Real dy = model::profile(model_index).y(idy+1) - + model::profile(model_index).y(idy); + + amrex::Real D = model::profile(model_index).state(idx, idy, var_index); + amrex::Real C = (model::profile(model_index).state(idx, idy+1, var_index) - + model::profile(model_index).state(idx, idy, var_index)) / dy; + amrex::Real B = (model::profile(model_index).state(idx+1, idy, var_index) - + model::profile(model_index).state(idx, idy, var_index)) / dx; + amrex::Real A = (model::profile(model_index).state(idx+1, idy+1, var_index) - + B * dx - C * dy - D) / (dx * dy); + + amrex::Real interp = A * (x - model::profile(model_index).x(idx)) * (y - model::profile(model_index).y(idy)) + + B * (x - model::profile(model_index).x(idx)) + C * (y - model::profile(model_index).y(idy)) + D; + + // safety check to make sure interp lies within the bounding points. We don't + // do this at the lower boundary, which usually corresponds to the center of the star. + if (x >= model::profile(model_index).x(0) && y >= model::profile(model_index).y(0)) { + amrex::Real minvar = amrex::min(model::profile(model_index).state(idx+1, idy+1, var_index) + model::profile(model_index).state(idx+1, idy , var_index), + model::profile(model_index).state(idx , idy+1, var_index), + model::profile(model_index).state(idx , idy , var_index)); + + amrex::Real maxvar = amrex::max(model::profile(model_index).state(idx+1, idy+1, var_index) + model::profile(model_index).state(idx+1, idy , var_index), + model::profile(model_index).state(idx , idy+1, var_index), + model::profile(model_index).state(idx , idy , var_index)); + + interp = amrex::Clamp(interp, minvar, maxvar); + + } else if (x >= model::profile(model_index).x(0)) { + amrex::Real minvar = amrex::min(model::profile(model_index).state(idx+1, idy, var_index), + model::profile(model_index).state(idx , idy, var_index)); + + amrex::Real maxvar = amrex::max(model::profile(model_index).state(idx+1, idy, var_index), + model::profile(model_index).state(idx , idy, var_index)); + + interp = amrex::Clamp(interp, minvar, maxvar); + + } else if (y >= model::profile(model_index).y(0)) { + amrex::Real minvar = amrex::min(model::profile(model_index).state(idx, idy+1, var_index), + model::profile(model_index).state(idx, idy , var_index)); + + amrex::Real maxvar = amrex::max(model::profile(model_index).state(idx, idy+1, var_index), + model::profile(model_index).state(idx, idy , var_index)); + + interp = amrex::Clamp(interp, minvar, maxvar); + + } + + return interp; + +} +#endif /// /// Subsample the interpolation to get an averaged profile. For this we need to know the /// 3D coordinate (relative to the model center) and cell size. /// AMREX_GPU_HOST_DEVICE AMREX_INLINE -Real interpolate_3d (const Real* loc, const Real* dx, int var_index, int nsub = 1, int model_index = 0) +amrex::Real interpolate_3d (const amrex::Real* loc, const amrex::Real* dx, int var_index, int nsub = 1, int model_index = 0) { // We perform a sub-grid-scale interpolation, where // nsub determines the number of intervals we split the zone into. @@ -168,18 +319,18 @@ Real interpolate_3d (const Real* loc, const Real* dx, int var_index, int nsub = // k = 1 --> z = loc(3) (halfway between left and right edge) // k = 2 --> z = loc(3) + dx(3) / 3 (1/6 of way from right edge of zone) - Real interp = 0.0_rt; + amrex::Real interp = 0.0_rt; for (int k = 0; k < nsub; ++k) { - Real z = loc[2] + (static_cast(k) + 0.5_rt * (1 - nsub)) * dx[2] / nsub; + amrex::Real z = loc[2] + (static_cast(k) + 0.5_rt * (1 - nsub)) * dx[2] / nsub; for (int j = 0; j < nsub; ++j) { - Real y = loc[1] + (static_cast(j) + 0.5_rt * (1 - nsub)) * dx[1] / nsub; + amrex::Real y = loc[1] + (static_cast(j) + 0.5_rt * (1 - nsub)) * dx[1] / nsub; for (int i = 0; i < nsub; ++i) { - Real x = loc[0] + (static_cast(i) + 0.5_rt * (1 - nsub)) * dx[0] / nsub; + amrex::Real x = loc[0] + (static_cast(i) + 0.5_rt * (1 - nsub)) * dx[0] / nsub; - Real dist = std::sqrt(x * x + y * y + z * z); + amrex::Real dist = std::sqrt(x * x + y * y + z * z); interp += interpolate(dist, var_index, model_index); } @@ -205,9 +356,9 @@ Real interpolate_3d (const Real* loc, const Real* dx, int var_index, int nsub = // we switch from the initial model to ambient material. AMREX_INLINE -void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, - const Real core_comp[NumSpec], Real temperature, Real dx, - Real envelope_mass, const Real envelope_comp[NumSpec], int model_index = 0) +void establish_hse (amrex::Real& mass_want, amrex::Real& central_density_want, amrex::Real& radius, + const amrex::Real core_comp[NumSpec], amrex::Real temperature, amrex::Real dx, + amrex::Real envelope_mass, const amrex::Real envelope_comp[NumSpec], int model_index = 0) { // Note that if central_density > 0, then this initial model generator will use it in calculating // the model. If mass is also provided in this case, we assume it is an estimate used for the purpose of @@ -231,9 +382,9 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, const int max_hse_iter = 250; int max_mass_iter; - Real rho_c, rho_c_old, drho_c; - Real mass, mass_old; - Real p_want, drho; + amrex::Real rho_c, rho_c_old, drho_c; + amrex::Real mass, mass_old; + amrex::Real p_want, drho; if (central_density_want > 0.0_rt) { @@ -307,10 +458,10 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, // Keep track of the mass enclosed below the current zone. - Real rl = 0.0_rt; - Real rr = rl + dx; + amrex::Real rl = 0.0_rt; + amrex::Real rr = rl + dx; - Real M_enclosed = (4.0_rt / 3.0_rt) * M_PI * (std::pow(rr, 3) - std::pow(rl, 3)) * model.state(0, model::idens); + amrex::Real M_enclosed = (4.0_rt / 3.0_rt) * M_PI * (std::pow(rr, 3) - std::pow(rl, 3)) * model.state(0, model::idens); mass = M_enclosed; //------------------------------------------------------------------------- @@ -341,7 +492,7 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, set_aux_comp_from_X(eos_state); #endif - Real g = -C::Gconst * M_enclosed / (std::pow(rl, 2)); + amrex::Real g = -C::Gconst * M_enclosed / (std::pow(rl, 2)); //---------------------------------------------------------------------- @@ -366,7 +517,7 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, // We difference HSE about the interface between the current // zone and the one just inside. - Real rho_avg = 0.5_rt * (model.state(i, model::idens) + model.state(i-1, model::idens)); + amrex::Real rho_avg = 0.5_rt * (model.state(i, model::idens) + model.state(i-1, model::idens)); p_want = model.state(i-1, model::ipres) + dx * rho_avg * g; eos(eos_input_rt, eos_state); @@ -412,7 +563,7 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, // Discretize the mass enclosed as (4 pi / 3) * rho * dr * (rl**2 + rl * rr + rr**2). - Real dM = (4.0_rt / 3.0_rt) * M_PI * model.state(i, model::idens) * dx * + amrex::Real dM = (4.0_rt / 3.0_rt) * M_PI * model.state(i, model::idens) * dx * (rr * rr + rl * rr + rl * rl); M_enclosed += dM; @@ -470,6 +621,7 @@ void establish_hse (Real& mass_want, Real& central_density_want, Real& radius, AMREX_INLINE void read_model_file(std::string& model_file, const int model_index=0) { + // This is assumed to work with 1D initial model only for now. bool found_model, found_dens, found_temp, found_pres, found_velr; bool found_spec[NumSpec]; @@ -537,6 +689,11 @@ read_model_file(std::string& model_file, const int model_index=0) { std::string word; iss >> word; +#if DIM_MODEL == 2 + // For 2D model, second column name is also coordinate -- skip + iss >> word; +#endif + while (iss >> word) { varnames_stored.push_back(word); } @@ -552,26 +709,62 @@ read_model_file(std::string& model_file, const int model_index=0) { // start reading in the data - amrex::Vector vars_stored; + amrex::Vector vars_stored; vars_stored.resize(nvars_model_file); int i{0}; +#if DIM_MODEL == 2 + int j{0}; +#endif + + // To track whether we read the first line. + // relevant to keep track of indices for DIM_MODEL==2 + bool firstLineRead = false; + while (getline(initial_model_file, line)) { + std::istringstream iss(line); + if (i > NPTS_MODEL) { amrex::Error("Error: model has more than NPTS_MODEL points, Increase MAX_NPTS_MODEL"); } + iss >> model::profile(model_index).r(i); - std::istringstream iss(line); +#if DIM_MODEL == 2 + // Reset index if the current index corresponds to position at index=0 + if (firstLineRead) { + if (iss == model::profile(model_index).x(0)) { + i = 0; + } + } else { + if (i > NPTS_X_MODEL) { + amrex::Error("Error: model has more than NPTS_X_MODEL points, Increase MAX_NPTS_X_MODEL"); + } + iss >> model::profile(model_index).x(i); + } - iss >> model::profile(model_index).r(i); + if (firstLineRead) { + if (model::profile(model_index).y(j) == model::profile(model_index).y(0)) { + j = 0; + } + } else { + if (j > NPTS_Y_MODEL) { + amrex::Error("Error: model has more than NPTS_Y_MODEL points, Increase MAX_NPTS_Y_MODEL"); + } + iss >> model::profile(model_index).y(j); + } +#endif - for (int j = 0; j < nvars_model_file; j++) { - iss >> vars_stored[j]; + for (int n = 0; n < nvars_model_file; n++) { + iss >> vars_stored[n]; } - for (int j = 0; j < model::nvars; j++) { - model::profile(model_index).state(i,j) = 0.0_rt; + for (int n = 0; n < model::nvars; n++) { +#if DIM_MODEL == 1 + model::profile(model_index).state(i, n) = 0.0_rt; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i, j, n) = 0.0_rt; +#endif } // make sure that each of the variables we care about is found @@ -588,36 +781,56 @@ read_model_file(std::string& model_file, const int model_index=0) { } #endif - for (int j = 0; j < nvars_model_file; j++) { + for (int n = 0; n < nvars_model_file; n++) { // keep track of whether the current variable from the model // file is one that we care about found_model = false; - if (varnames_stored[j] == "density") { - model::profile(model_index).state(i,model::idens) = vars_stored[j]; + if (varnames_stored[n] == "density") { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::idens) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::idens) = vars_stored[n]; +#endif found_model = true; found_dens = true; - } else if (varnames_stored[j] == "temperature") { - model::profile(model_index).state(i,model::itemp) = vars_stored[j]; + } else if (varnames_stored[n] == "temperature") { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::itemp) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::itemp) = vars_stored[n]; +#endif found_model = true; found_temp = true; - } else if (varnames_stored[j] == "pressure") { - model::profile(model_index).state(i,model::ipres) = vars_stored[j]; + } else if (varnames_stored[n] == "pressure") { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::ipres) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::ipres) = vars_stored[n]; +#endif found_model = true; found_pres = true; - } else if (varnames_stored[j] == "velocity") { - model::profile(model_index).state(i,model::ivelr) = vars_stored[j]; + } else if (varnames_stored[n] == "velocity") { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::ivelr) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::ivelr) = vars_stored[n]; +#endif found_model = true; found_velr = true; } else { for (int comp = 0; comp < NumSpec; comp++) { - if (varnames_stored[j] == spec_names_cxx[comp]) { - model::profile(model_index).state(i,model::ispec+comp) = vars_stored[j]; + if (varnames_stored[n] == spec_names_cxx[comp]) { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::ispec+comp) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::ispec+comp) = vars_stored[n]; +#endif found_model = true; found_spec[comp] = true; break; @@ -626,8 +839,12 @@ read_model_file(std::string& model_file, const int model_index=0) { #if NAUX_NET > 0 if (!found_model) { for (int comp = 0; comp < NumAux; comp++) { - if (varnames_stored[j] == aux_names_cxx[comp]) { - model::profile(model_index).state(i,model::iaux+comp) = vars_stored[j]; + if (varnames_stored[n] == aux_names_cxx[comp]) { +#if DIM_MODEL == 1 + model::profile(model_index).state(i,model::iaux+comp) = vars_stored[n]; +#elif DIM_MODEL == 2 + model::profile(model_index).state(i,j,model::iaux+comp) = vars_stored[n]; +#endif found_model = true; found_aux[comp] = true; break; @@ -639,9 +856,10 @@ read_model_file(std::string& model_file, const int model_index=0) { // yell if we didn't find the current variable - if (!found_model && i == 0) { + if (!found_model && !firstLineRead + ) { amrex::Print() << Font::Bold << FGColor::Yellow - << "[WARNING] variable not found: " << varnames_stored[j] + << "[WARNING] variable not found: " << varnames_stored[n] << ResetDisplay << std::endl; } @@ -650,7 +868,7 @@ read_model_file(std::string& model_file, const int model_index=0) { // were all the variables we care about provided? - if (i == 0) { + if (!firstLineRead) { if (!found_dens) { amrex::Print() << Font::Bold << FGColor::Yellow << "[WARNING] density not provided in inputs file" @@ -694,14 +912,35 @@ read_model_file(std::string& model_file, const int model_index=0) { #endif } +#if DIM_MODEL == 1 i++; +#elif DIM_MODEL == 2 + // update index if current position is different from last position + if (model::profile(model_index).x(i) != model::profile(model_index).x(i-1)) { + i++; + } + if (model::profile(model_index).y(j) != model::profile(model_index).y(j-1)) { + j++; + } +#endif + firstLineRead = true; } // end of loop over lines in the model file +#if DIM_MODEL == 1 model::npts = i; +#elif DIM_MODEL == 2 + model::npts_x = i; + model::npts_y = j; +#endif amrex::Print() << Font::Bold << FGColor::Green +#if DIM_MODEL == 1 << model::npts << " points found in the initial model" +#elif DIM_MODEL == 2 + << model::npts_x << " x-points found in the initial model" + << model::npts_y << " y-points found in the initial model" +#endif << ResetDisplay << std::endl; initial_model_file.close(); diff --git a/Util/model_parser/model_parser_data.H b/Util/model_parser/model_parser_data.H index f187921316..25098bc714 100644 --- a/Util/model_parser/model_parser_data.H +++ b/Util/model_parser/model_parser_data.H @@ -26,12 +26,24 @@ namespace model constexpr int iaux = -1; #endif +#if DIM_MODEL == 1 extern AMREX_GPU_MANAGED int npts; +#elif DIM_MODEL == 2 + extern AMREX_GPU_MANAGED int npts_x; + extern AMREX_GPU_MANAGED int npts_y; +#endif extern AMREX_GPU_MANAGED bool initialized; struct initial_model_t { +#if DIM_MODEL == 1 amrex::Array2D state; amrex::Array1D r; +#elif DIM_MODEL == 2 + amrex::Array3D state; + amrex::Array1D x; + amrex::Array1D y; +#endif }; // Tolerance used for getting the total star mass equal to the desired mass. diff --git a/Util/model_parser/model_parser_data.cpp b/Util/model_parser/model_parser_data.cpp index cfbc475cbe..1860056363 100644 --- a/Util/model_parser/model_parser_data.cpp +++ b/Util/model_parser/model_parser_data.cpp @@ -3,8 +3,12 @@ namespace model { - +#if DIM_MODEL == 1 AMREX_GPU_MANAGED int npts; +#elif DIM_MODEL == 2 + AMREX_GPU_MANAGED int npts_x; + AMREX_GPU_MANAGED int npts_y; +#endif AMREX_GPU_MANAGED bool initialized; AMREX_GPU_MANAGED amrex::Array1D profile; From a85f7d7b1099d83401402dca335afa105fc049f2 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 18:21:35 -0400 Subject: [PATCH 05/20] delete extra comment line --- Util/model_parser/model_parser.H | 1 - 1 file changed, 1 deletion(-) diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index 098889068f..b4b53907cf 100644 --- a/Util/model_parser/model_parser.H +++ b/Util/model_parser/model_parser.H @@ -37,7 +37,6 @@ using namespace amrex; /// Note that the first data column is the coordinate, but that is not /// listed as a variable. /// -/// /// The new format has a single line header, labeling the columns, /// including the coordinate. This format plays nicely with NumPy /// genfromtxt to read in the header From 068ddb7718f0b111803fe5d4c03d9bb471518c69 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 21:01:59 -0400 Subject: [PATCH 06/20] add good define --- .github/workflows/good_defines.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/good_defines.txt b/.github/workflows/good_defines.txt index cdfaba0e22..99e88357f4 100644 --- a/.github/workflows/good_defines.txt +++ b/.github/workflows/good_defines.txt @@ -14,6 +14,7 @@ BL_LANG_FORT BL_LAZY BL_USE_SETBUF MODEL_PARSER +DIM_MODEL DIFFUSION DO_PROBLEM_POST_INIT DO_PROBLEM_POST_RESTART From eb2b86deb4caea78d0c0e11a5b4c54b8af2ad2ac Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 21:13:21 -0400 Subject: [PATCH 07/20] more guards --- Util/model_parser/model_parser.H | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index b4b53907cf..1cb2abd1e1 100644 --- a/Util/model_parser/model_parser.H +++ b/Util/model_parser/model_parser.H @@ -303,6 +303,7 @@ interpolate(const amrex::Real x, const amrex::Real y, } #endif +#if DIM_MODEL == 1 /// /// Subsample the interpolation to get an averaged profile. For this we need to know the /// 3D coordinate (relative to the model center) and cell size. @@ -616,6 +617,8 @@ void establish_hse (amrex::Real& mass_want, amrex::Real& central_density_want, a model::initialized = true; model::npts = NPTS_MODEL; } +#endif + AMREX_INLINE void From 716d39a0c3dd1bdf512c668ee3314497af9e913b Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 21:31:35 -0400 Subject: [PATCH 08/20] use npts_x --- Exec/science/xrb_spherical/problem_initialize.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/problem_initialize.H b/Exec/science/xrb_spherical/problem_initialize.H index e6dfea4f61..e480b950cb 100644 --- a/Exec/science/xrb_spherical/problem_initialize.H +++ b/Exec/science/xrb_spherical/problem_initialize.H @@ -197,12 +197,12 @@ void problem_initialize () // set the ambient state for the upper boundary condition using upper most r zone at θ=0 - ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, ng, model::idens); + ambient::ambient_state[URHO] = model::profile(0).state(model::npts_x-1, ng, model::idens); - ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, ng, model::itemp); + ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts_x-1, ng, model::itemp); for (int n = 0; n < NumSpec; n++) { ambient::ambient_state[UFS+n] = - ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, ng, model::ispec+n); + ambient::ambient_state[URHO] * model::profile(0).state(model::npts_x-1, ng, model::ispec+n); } ambient::ambient_state[UMX] = 0.0_rt; From be6fb284125671f3c207b0c1be7b86db2a7f9232 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 21:31:53 -0400 Subject: [PATCH 09/20] more fixes --- Util/model_parser/model_parser.H | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/Util/model_parser/model_parser.H b/Util/model_parser/model_parser.H index 1cb2abd1e1..dcf9faddf4 100644 --- a/Util/model_parser/model_parser.H +++ b/Util/model_parser/model_parser.H @@ -182,7 +182,7 @@ locate_y(const amrex::Real y, const int model_index) { while (ilo+1 != ihi) { int imid = (ilo + ihi) / 2; - if (x <= model::profile(model_index).y(imid)) { + if (y <= model::profile(model_index).y(imid)) { ihi = imid; } else { ilo = imid; @@ -266,12 +266,12 @@ interpolate(const amrex::Real x, const amrex::Real y, // safety check to make sure interp lies within the bounding points. We don't // do this at the lower boundary, which usually corresponds to the center of the star. if (x >= model::profile(model_index).x(0) && y >= model::profile(model_index).y(0)) { - amrex::Real minvar = amrex::min(model::profile(model_index).state(idx+1, idy+1, var_index) + amrex::Real minvar = amrex::min(model::profile(model_index).state(idx+1, idy+1, var_index), model::profile(model_index).state(idx+1, idy , var_index), model::profile(model_index).state(idx , idy+1, var_index), model::profile(model_index).state(idx , idy , var_index)); - amrex::Real maxvar = amrex::max(model::profile(model_index).state(idx+1, idy+1, var_index) + amrex::Real maxvar = amrex::max(model::profile(model_index).state(idx+1, idy+1, var_index), model::profile(model_index).state(idx+1, idy , var_index), model::profile(model_index).state(idx , idy+1, var_index), model::profile(model_index).state(idx , idy , var_index)); @@ -727,33 +727,38 @@ read_model_file(std::string& model_file, const int model_index=0) { std::istringstream iss(line); - if (i > NPTS_MODEL) { +#if DIM_MODEL == 1 + if (i >= NPTS_MODEL) { amrex::Error("Error: model has more than NPTS_MODEL points, Increase MAX_NPTS_MODEL"); } iss >> model::profile(model_index).r(i); +#elif DIM_MODEL == 2 + amrex::Real temp; + iss >> temp; -#if DIM_MODEL == 2 // Reset index if the current index corresponds to position at index=0 if (firstLineRead) { - if (iss == model::profile(model_index).x(0)) { + if (temp == model::profile(model_index).x(0)) { i = 0; } } else { - if (i > NPTS_X_MODEL) { + if (i >= NPTS_X_MODEL) { amrex::Error("Error: model has more than NPTS_X_MODEL points, Increase MAX_NPTS_X_MODEL"); } - iss >> model::profile(model_index).x(i); + model::profile(model_index).x(i) = temp; } + iss >> temp; + if (firstLineRead) { - if (model::profile(model_index).y(j) == model::profile(model_index).y(0)) { + if (temp == model::profile(model_index).y(0)) { j = 0; } } else { - if (j > NPTS_Y_MODEL) { + if (j >= NPTS_Y_MODEL) { amrex::Error("Error: model has more than NPTS_Y_MODEL points, Increase MAX_NPTS_Y_MODEL"); } - iss >> model::profile(model_index).y(j); + model::profile(model_index).y(j) = temp; } #endif From be552c9010116ef273aee73cea97b7a0c363936a Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 22:14:32 -0400 Subject: [PATCH 10/20] update initial model --- Exec/science/xrb_spherical/GNUmakefile | 2 + Exec/science/xrb_spherical/initial_model.H | 7 +-- .../problem_initialize_state_data.H | 49 +++++++++++-------- 3 files changed, 34 insertions(+), 24 deletions(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index 24a6684db0..6b61ad975d 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -22,6 +22,8 @@ USE_JACOBIAN_CACHING = TRUE USE_MODEL_PARSER = TRUE NUM_MODELS := 1 DIM_MODEL := 2 +MAX_NPTS_X_MODEL = 9900 +MAX_NPTS_Y_MODEL = 9000 # This sets the EOS directory in $(MICROPHYSICS_HOME)/EOS EOS_DIR := helmholtz diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index dbe45fc4da..fbfe7b9b47 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -400,13 +400,14 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // we actually require that the number of points is the same for each // model, so we'll just set it each time - model::npts = npts_model_x < npts_model_y ? npts_model_x : npts_model_y; + model::npts_x = npts_model_x; + model::npts_y = npts_model_y; model::initialized = true; - if (npts_model_x > NPTS_MODEL) { + if (npts_model_x > NPTS_X_MODEL) { amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_MODEL"); } - if (npts_model_y > NPTS_MODEL) { + if (npts_model_y > NPTS_Y_MODEL) { amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_MODEL"); } diff --git a/Exec/science/xrb_spherical/problem_initialize_state_data.H b/Exec/science/xrb_spherical/problem_initialize_state_data.H index 0431916a8e..495ae4e291 100644 --- a/Exec/science/xrb_spherical/problem_initialize_state_data.H +++ b/Exec/science/xrb_spherical/problem_initialize_state_data.H @@ -10,42 +10,49 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_initialize_state_data (int i, int j, int k, - Array4 const& state, + amrex::Array4 const& state, const GeometryData& geomdata) { - const Real* dx = geomdata.CellSize(); - const Real* problo = geomdata.ProbLo(); + const amrex::Real* dx = geomdata.CellSize(); + const amrex::Real* problo = geomdata.ProbLo(); - Real r = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; - Real theta = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; + amrex::Real r = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + amrex::Real theta = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; // blending factor - Real f; + // amrex::Real f; - if (theta < problem::theta_half_max) { - f = 1.0_rt; + // if (theta < problem::theta_half_max) { + // f = 1.0_rt; - } else if (theta > problem::theta_half_max + problem::theta_half_width) { - f = 0.0_rt; + // } else if (theta > problem::theta_half_max + problem::theta_half_width) { + // f = 0.0_rt; - } else { - f = -(theta - problem::theta_half_max) / problem::theta_half_width + 1.0_rt; - } + // } else { + // f = -(theta - problem::theta_half_max) / problem::theta_half_width + 1.0_rt; + // } + + // state(i,j,k,URHO) = f * interpolate(r, model::idens, 1) + + // (1.0_rt - f) * interpolate(r, model::idens, 0); - state(i,j,k,URHO) = f * interpolate(r, model::idens, 1) + - (1.0_rt - f) * interpolate(r, model::idens, 0); + // state(i,j,k,UTEMP) = f * interpolate(r, model::itemp, 1) + + // (1.0_rt - f) * interpolate(r, model::itemp, 0); - state(i,j,k,UTEMP) = f * interpolate(r, model::itemp, 1) + - (1.0_rt - f) * interpolate(r, model::itemp, 0); + // amrex::Real temppres = f * interpolate(r, model::ipres, 1) + + // (1.0_rt - f) * interpolate(r, model::ipres, 0); - Real temppres = f * interpolate(r, model::ipres, 1) + - (1.0_rt - f) * interpolate(r, model::ipres, 0); + // for (int n = 0; n < NumSpec; n++) { + // state(i,j,k,UFS+n) = f * interpolate(r, model::ispec+n, 1) + + // (1.0_rt - f) * interpolate(r, model::ispec+n, 0); + // } + state(i,j,k,URHO) = interpolate(r, theta, model::idens, 0); + state(i,j,k,UTEMP) = interpolate(r, theta, model::itemp, 0); + amrex::Real temppres = interpolate(r, theta, model::ipres, 0); for (int n = 0; n < NumSpec; n++) { - state(i,j,k,UFS+n) = f * interpolate(r, model::ispec+n, 1) + - (1.0_rt - f) * interpolate(r, model::ispec+n, 0); + state(i,j,k,UFS+n) = interpolate(r, theta, model::ispec+n, 0); } eos_t eos_state; From f7aa88966871e1b3374ffa3629376866731302a4 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 22:17:07 -0400 Subject: [PATCH 11/20] update interpolate in castro_react --- Source/reactions/Castro_react.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index a2b35d160b..34c0c98cb0 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -281,6 +281,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra rr[2] = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; #endif +#if DIM_MODEL == 1 Real dist; if (domain_is_plane_parallel) { @@ -290,7 +291,8 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra } burn_state.T_fixed = interpolate(dist, model::itemp); - +#elif DIM_MODEL == 2 + burn_state.T_fix = interpolate(rr[0], rr[1], model::itemp); } #endif @@ -632,6 +634,7 @@ Castro::react_state(Real time, Real dt) rr[2] = problo[2] + dx[2] * (static_cast(k) + 0.5_rt) - problem::center[2]; #endif +#if DIM_MODEL == 1 Real dist; if (domain_is_plane_parallel) { @@ -641,7 +644,8 @@ Castro::react_state(Real time, Real dt) } burn_state.T_fixed = interpolate(dist, model::itemp); - +#elif DIM_MODEL == 2 + burn_state.T_fix = interpolate(rr[0], rr[1], model::itemp); } #endif From 5326438c8f89a7ae8334e1297f6fbc40475f60ba Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 22:18:32 -0400 Subject: [PATCH 12/20] fix --- Source/reactions/Castro_react.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 34c0c98cb0..7ae236292a 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -292,7 +292,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra burn_state.T_fixed = interpolate(dist, model::itemp); #elif DIM_MODEL == 2 - burn_state.T_fix = interpolate(rr[0], rr[1], model::itemp); + burn_state.T_fixed = interpolate(rr[0], rr[1], model::itemp); } #endif @@ -645,7 +645,7 @@ Castro::react_state(Real time, Real dt) burn_state.T_fixed = interpolate(dist, model::itemp); #elif DIM_MODEL == 2 - burn_state.T_fix = interpolate(rr[0], rr[1], model::itemp); + burn_state.T_fixed = interpolate(rr[0], rr[1], model::itemp); } #endif From 92474ccd9ba71c7768c83b90a99473c5264654f7 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 20 Oct 2025 22:19:39 -0400 Subject: [PATCH 13/20] mroe fix --- Source/reactions/Castro_react.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index 7ae236292a..de12e5a7fd 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -293,6 +293,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra burn_state.T_fixed = interpolate(dist, model::itemp); #elif DIM_MODEL == 2 burn_state.T_fixed = interpolate(rr[0], rr[1], model::itemp); +#endif } #endif @@ -646,6 +647,7 @@ Castro::react_state(Real time, Real dt) burn_state.T_fixed = interpolate(dist, model::itemp); #elif DIM_MODEL == 2 burn_state.T_fixed = interpolate(rr[0], rr[1], model::itemp); +#endif } #endif From 3f6280c32153dc5d5fff3202ab03bfe42446324d Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 22 Oct 2025 14:04:04 -0400 Subject: [PATCH 14/20] fix initial model --- Exec/science/xrb_spherical/initial_model.H | 45 +++++++++++++--------- 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index fbfe7b9b47..15d5b6a402 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -98,7 +98,7 @@ AMREX_INLINE void HSE_integrate(const int i, const int j, const amrex::Real dx, const amrex::Real dy, const int idir, - const bool isentropic, bool& fluff, + bool& isentropic, bool& fluff, const amrex::Real entropy_base, const model_t& model_params, const int model_num) { // @@ -284,11 +284,24 @@ void HSE_integrate(const int i, const int j, dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } + + // Check if we entered isothermal region when integrating upward + if (temp_zone < model_params.T_lo && idir == 0) { + temp_zone = model_params.T_lo; + isentropic = false; + } + // check if the density falls below our minimum cut-off -- // if so, floor it // Also detect whether we've entered the fluff region. + // Note that, its only possible when integrating upward. - if (dens_zone < model_params.low_density_cutoff) { + if (dens_zone < model_params.low_density_cutoff && idir == 0) { dens_zone = model_params.low_density_cutoff; temp_zone = model_params.T_lo; converged_hse = true; @@ -296,12 +309,6 @@ void HSE_integrate(const int i, const int j, break; } - // if (A < TOL .and. B < ETOL) then - if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { - converged_hse = true; - break; - } - } else { // isothermal case // Since T is fixed, the only unknown is density @@ -338,16 +345,17 @@ void HSE_integrate(const int i, const int j, dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); - if (dens_zone < model_params.low_density_cutoff) { - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; + if (std::abs(drho) < fw::TOL * dens_zone) { converged_hse = true; - fluff = true; break; } - if (std::abs(drho) < fw::TOL * dens_zone) { + // fluff region only possible when integrate upward + if (dens_zone < model_params.low_density_cutoff && idir == 0) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; converged_hse = true; + fluff = true; break; } @@ -397,18 +405,15 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // Create a 2-d uniform grid that is identical to the mesh that we are // mapping onto, and then we want to force it into HSE on that mesh. - // we actually require that the number of points is the same for each - // model, so we'll just set it each time - model::npts_x = npts_model_x; model::npts_y = npts_model_y; model::initialized = true; if (npts_model_x > NPTS_X_MODEL) { - amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_MODEL"); + amrex::Error("Error: model has more than NPTS_MODEL points in x-dir, Increase MAX_NPTS_X_MODEL"); } if (npts_model_y > NPTS_Y_MODEL) { - amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_MODEL"); + amrex::Error("Error: model has more than NPTS_MODEL points in y-dir, Increase MAX_NPTS_Y_MODEL"); } // create the grid -- cell centers @@ -520,6 +525,8 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // make the base thermodynamics consistent for this base point, // this is assumed to be on θ=0 -- that is what we will integrate from! + // here 0-index is actually not θ=0, due to extra ghost cell, + // we might need to determine the actual index. eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); @@ -544,7 +551,7 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // First integrate to the right from θ=0 at the base layer. // The base layer is isothermal. - for (int j = 0; j < npts_model_y; ++j) { + for (int j = 1; j < npts_model_y; ++j) { HSE_integrate(index_base, j, dx, dy, 2, isentropic, fluff, entropy_base, From 03f3213ece5ac71f88738dc7037305b33595967e Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 27 Oct 2025 22:28:00 -0400 Subject: [PATCH 15/20] separate out initial condition setting into separate function --- Exec/science/xrb_spherical/initial_model.H | 186 +++++++++++---------- 1 file changed, 99 insertions(+), 87 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index 15d5b6a402..2d6786506d 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -94,6 +94,100 @@ amrex::Real evaluate_tanh(const amrex::Real z) { } +AMREX_INLINE +void set_initial_condition(const int index_base, + const int npts_model_x, const amrex::Real xmin, + const int npts_model_y, const amrex::Real ymin, + const model_t& model_params, const int model_num) { + + // Set the initial condition of the mesh + + eos_t eos_state; + eos_state.T = model_params.T_star; + eos_state.rho = model_params.dens_base; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model_params.xn_star[n]; + } + + eos(eos_input_rt, eos_state); + + // store the conditions at the base -- we'll use the entropy later + // to constrain the isentropic layer + + amrex::Real pres_base = eos_state.p; + + // set an initial temperature profile and composition for the 2D grid + + for (int i = 0; i < npts_model_x; i++) { + + amrex::Real xc = model::profile(model_num).x(i) - + (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; + + for (int j = 0; j < npts_model_y; j++) { + + // vertical hyperbolic tangent transition + + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) = + model_params.xn_star[n] + + 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // force them to sum to 1 + + amrex::Real sumX = 0.0_rt; + for (int n = 0; n < NumSpec; n++) { + sumX += model::profile(model_num).state(i, j, model::ispec+n); + } + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, j, model::ispec+n) /= sumX; + } + + // temperature profile -- it is constant below the base + // It follows the hyperbolic tangent transition above the base + + if (i <= index_base) { + model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; + } else { + model::profile(model_num).state(i, j, model::itemp) = + model_params.T_star + + 0.5_rt * (model_params.T_hi - model_params.T_star) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // + // We will also have a hyperbolic tangent transition in the theta dir + // but for now just assume theres is no hot/cold model. + // + + + + // the density and pressure will be determined via HSE, + // for now, set them to the base conditions as initial guess + + model::profile(model_num).state(i, j, model::idens) = model_params.dens_base; + model::profile(model_num).state(i, j, model::ipres) = pres_base; + } + } + + // make the base thermodynamics consistent for this base point, + // this is assumed to be on θ=0 -- that is what we will integrate from! + // here 0-index is actually not θ=0, due to extra ghost cell, + // we might need to determine the actual index. + + eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); + eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(index_base, 0, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + model::profile(model_num).state(index_base, 0, model::ipres) = eos_state.p; +} + + AMREX_INLINE void HSE_integrate(const int i, const int j, const amrex::Real dx, const amrex::Real dy, @@ -451,93 +545,10 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // This is the initial guess conditions. // - // determine the conditions at the base -- this is below the atmosphere - // this point is assumed to be at the rotating pole, θ=0 - - eos_t eos_state; - eos_state.T = model_params.T_star; - eos_state.rho = model_params.dens_base; - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model_params.xn_star[n]; - } - - eos(eos_input_rt, eos_state); - - // store the conditions at the base -- we'll use the entropy later - // to constrain the isentropic layer - - amrex::Real pres_base = eos_state.p; - - // set an initial temperature profile and composition for the 2D grid - - for (int i = 0; i < npts_model_x; i++) { - - amrex::Real xc = model::profile(model_num).x(i) - - (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; - - for (int j = 0; j < npts_model_y; j++) { - - // vertical hyperbolic tangent transition - - for (int n = 0; n < NumSpec; n++) { - model::profile(model_num).state(i, j, model::ispec+n) = - model_params.xn_star[n] + - 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * - (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); - } - - // force them to sum to 1 - - amrex::Real sumX = 0.0_rt; - for (int n = 0; n < NumSpec; n++) { - sumX += model::profile(model_num).state(i, j, model::ispec+n); - } - for (int n = 0; n < NumSpec; n++) { - model::profile(model_num).state(i, j, model::ispec+n) /= sumX; - } - - // temperature profile -- it is constant below the base - // It follows the hyperbolic tangent transition above the base - - if (i <= index_base) { - model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; - } else { - model::profile(model_num).state(i, j, model::itemp) = - model_params.T_star + - 0.5_rt * (model_params.T_hi - model_params.T_star) * - (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); - } - - // - // We will also have a hyperbolic tangent transition in the theta dir - // but for now just assume theres is no hot/cold model. - // - - - - // the density and pressure will be determined via HSE, - // for now, set them to the base conditions - - model::profile(model_num).state(i, j, model::idens) = model_params.dens_base; - model::profile(model_num).state(i, j, model::ipres) = pres_base; - } - } - - // make the base thermodynamics consistent for this base point, - // this is assumed to be on θ=0 -- that is what we will integrate from! - // here 0-index is actually not θ=0, due to extra ghost cell, - // we might need to determine the actual index. - - eos_state.rho = model::profile(model_num).state(index_base, 0, model::idens); - eos_state.T = model::profile(model_num).state(index_base, 0, model::itemp); - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model::profile(model_num).state(index_base, 0, model::ispec+n); - } - - eos(eos_input_rt, eos_state); - - model::profile(model_num).state(index_base, 0, model::ipres) = eos_state.p; - + set_initial_condition(index_base, + npts_model_x, xmin, + npts_model_y, ymin, + model_params, model_num); // // 2D HSE + entropy solve @@ -613,6 +624,7 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // now we need to know the entropy we are confining ourselves to + eos_t eos_state; eos_state.rho = model::profile(model_num).state(i-1, j, model::idens); eos_state.T = model::profile(model_num).state(i-1, j, model::itemp); for (int n = 0; n < NumSpec; n++) { From 1c967601908bd34e26f7af292ca0ab0d750cfd5b Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 27 Oct 2025 23:35:58 -0400 Subject: [PATCH 16/20] change to smaller max npts x/y model --- Exec/science/xrb_spherical/GNUmakefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index 6b61ad975d..b46bbb9a8f 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -22,8 +22,8 @@ USE_JACOBIAN_CACHING = TRUE USE_MODEL_PARSER = TRUE NUM_MODELS := 1 DIM_MODEL := 2 -MAX_NPTS_X_MODEL = 9900 -MAX_NPTS_Y_MODEL = 9000 +MAX_NPTS_X_MODEL = 3100 +MAX_NPTS_Y_MODEL = 1000 # This sets the EOS directory in $(MICROPHYSICS_HOME)/EOS EOS_DIR := helmholtz From bc3c2c0b68332e2053edb9579f389807b72fedec Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 29 Oct 2025 13:06:22 -0400 Subject: [PATCH 17/20] split out hse pres calcualtion and try averaging for integration --- Exec/science/xrb_spherical/initial_model.H | 230 ++++++++++++++------- 1 file changed, 152 insertions(+), 78 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index 2d6786506d..b525bc601e 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -188,6 +188,76 @@ void set_initial_condition(const int index_base, } +AMREX_INLINE +amrex::Real HSE_pressure(const int i, const int j, const int idir, + const amrex::Real dx, const amrex::Real dy, + const amrex::Real dens_zone, const amrex::Real Omega2, + const int model_num) { + // Evaluate the HSE pressure given density of the current zone + // using information of the surrounding zones. + + amrex::Real hse_p; + + if (idir == 0 || idir == 1) { + + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); + + if (idir == 0) { + // Integrate upward + + amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; + amrex::Real dens_im = 0.5_rt * + (dens_zone + model::profile(model_num).state(i-1, j, model::idens)); + amrex::Real dp_im1 = -dens_im * dx * (gravity::const_grav + + 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); + + hse_p = model::profile(model_num).state(i-1, j, model::ipres) - dp_im1; + + } else { + // Integrate downward + + amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; + amrex::Real dens_ip = 0.5_rt * + (dens_zone + model::profile(model_num).state(i+1, j, model::idens)); + amrex::Real dp_ip1 = dens_ip * dx * (gravity::const_grav + + 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); + + hse_p = model::profile(model_num).state(i+1, j, model::ipres) - dp_ip1; + } + + // Do an average between zone from bottom and left if not along the pole + if (j != 0) { + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + amrex::Real dens_jm = 0.5_rt * + (dens_zone + model::profile(model_num).state(i, j-1, model::idens)); + amrex::Real dp_jm1 = -0.5_rt * dy * dens_jm * Omega2 * rc * rc * sinTwoThetam; + + hse_p = 0.5_rt * (hse_p + model::profile(model_num).state(i, j-1, model::ipres) - dp_jm1); + } + + } else if (idir == 2) { + // Integrate to the right + + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + + amrex::Real dens_jm = 0.5_rt * + (dens_zone + model::profile(model_num).state(i, j-1, model::idens)); + amrex::Real dp_jm1 = -0.5_rt * dy * dens_jm * Omega2 * rc * rc * sinTwoThetam; + + hse_p = model::profile(model_num).state(i, j-1, model::ipres) - dp_jm1; + + } else { + amrex::Error("Invalid integrate direction"); + } + + return hse_p; +} + + AMREX_INLINE void HSE_integrate(const int i, const int j, const amrex::Real dx, const amrex::Real dy, @@ -254,25 +324,10 @@ void HSE_integrate(const int i, const int j, amrex::Real drho; amrex::Real dtemp = 0.0_rt; - // Set index offset based on integration direction - - int i0; - int j0; - - if (idir == 0) { - // Integrate upward - i0 = -1; - j0 = 0; - } else if (idir == 1) { - // Integrate downward - i0 = +1; - j0 = 0; - } else if (idir == 2) { - // integrate to the right - i0 = 0; - j0 = -1; - } else { - amrex::Error("Invalid integration direction!"); + amrex::Real Omega2 = 0.0_rt; + if (castro::rotational_period > 0.0_rt && + castro::rotation_include_centrifugal) { + Omega2 = std::pow(2.0_rt * M_PI / castro::rotational_period, 2); } // Do Newton solve @@ -296,35 +351,9 @@ void HSE_integrate(const int i, const int j, // We assumed integration is always in the right dir, // hence reference state is always from left {i, j-1} - // First just try integrate either in vertical or horizontal direction - - amrex::Real r = model::profile(model_num).x(i) + static_cast(i0) * 0.5_rt * dx; - amrex::Real twoTheta = 2.0_rt * (model::profile(model_num).y(j) + - static_cast(j0) * 0.5_rt * dy); - amrex::Real Omega = 0.0_rt; - if (castro::rotational_period > 0.0_rt && - castro::rotation_include_centrifugal) { - Omega = 2.0_rt * M_PI / castro::rotational_period; - } - - if (idir == 0 || idir == 1) { - - // integrating in vertical direction + // Get HSE pressure based on integration direction - p_want = model::profile(model_num).state(i+i0, j, model::ipres) - - static_cast(i0) * 0.5_rt * dx * - (dens_zone + model::profile(model_num).state(i+i0, j, model::idens)) * - (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); - - } else if (idir == 2) { - - // integrating in to the right - - p_want = model::profile(model_num).state(i, j+j0, model::ipres) - - static_cast(j0) * 0.25_rt * dy * - (dens_zone + model::profile(model_num).state(i, j+j0, model::idens)) * - Omega * Omega * r * r * std::sin(twoTheta); - } + amrex::Real p_want = HSE_pressure(i, j, idir, dx, dy, dens_zone, Omega2, model_num); // // We can either have isentropic or isothermal case. @@ -357,16 +386,46 @@ void HSE_integrate(const int i, const int j, amrex::Real B = entropy_base - entropy_zone; amrex::Real dAdT = -eos_state.dpdT; - amrex::Real dAdrho = - eos_state.dpdr; + amrex::Real dAdrho = 0.0_rt; + + // Add d HSE_p / drho contribution if (idir == 0 || idir == 1) { - // integrate vertically - dAdrho += -static_cast(i0) * 0.5_rt * dx * - (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); + + if (idir == 0) { + // Integrate upward + + amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; + dAdrho += 0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); + + } else { + // Integrate downward + + amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; + dAdrho += -0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); + } + + // Do an average between zone from bottom and left if not along the pole + if (j != 0) { + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + + dAdrho += 0.5_rt * (dAdrho + 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam); + } + } else if (idir == 2) { - // integrate horizontally - dAdrho += -static_cast(j0) * 0.25_rt * dy * - Omega * Omega * r * r * std::sin(twoTheta); + // Integrate to the right + + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + + dAdrho += 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam; } + + dAdrho += - eos_state.dpdr; amrex::Real dBdT = -eos_state.dsdT; amrex::Real dBdrho = -eos_state.dsdr; @@ -423,18 +482,47 @@ void HSE_integrate(const int i, const int j, pres_zone = eos_state.p; amrex::Real A = p_want - pres_zone; - amrex::Real dAdrho = - eos_state.dpdr; + amrex::Real dAdrho = 0.0_rt; + // Add d HSE_p / drho contribution if (idir == 0 || idir == 1) { - // integrate vertically - dAdrho += -static_cast(i0) * 0.5_rt * dx * - (gravity::const_grav + 0.5_rt * Omega * Omega * r * (1.0_rt - std::cos(twoTheta))); + + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); + + if (idir == 0) { + // Integrate upward + + amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; + dAdrho += 0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); + + } else { + // Integrate downward + + amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; + dAdrho += -0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); + } + + // Do an average between zone from bottom and left if not along the pole + if (j != 0) { + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + + dAdrho += 0.5_rt * (dAdrho + 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam); + } + } else if (idir == 2) { - // integrate horizontally - dAdrho += -static_cast(j0) * 0.25_rt * dy * - Omega * Omega * r * r * std::sin(twoTheta); + // Integrate to the right + + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); + amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + + dAdrho += 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam; } + dAdrho += - eos_state.dpdr; + drho = - A / dAdrho; dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); @@ -554,14 +642,14 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // 2D HSE + entropy solve // + // First integrate to the right from θ=0 at the base layer. + // The base layer is isothermal. + bool fluff = false; bool isentropic = false; bool flipped = false; amrex::Real entropy_base = -1.0_rt; - // First integrate to the right from θ=0 at the base layer. - // The base layer is isothermal. - for (int j = 1; j < npts_model_y; ++j) { HSE_integrate(index_base, j, dx, dy, 2, isentropic, fluff, @@ -649,19 +737,5 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr } - - - - - // Start from the base layer at θ=0 and integrate upward and integrate to the right - - // Now we go through the intersection regions, and integrate upward. - // Since the integrated result should be the same regardless of the direction, - // we do an average from both direction - - // Start from the base layer at θ=0 and integrate downward - - // Now go through the intersection regions and integrate downward. - } #endif From 7e20c281c52a846098ea9d4a0e67b8513bf7203a Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 3 Nov 2025 13:54:41 -0500 Subject: [PATCH 18/20] try to use relaxation after 1d integration --- Exec/science/xrb_spherical/initial_model.H | 979 ++++++++++++------ .../xrb_spherical/problem_initialize.H | 10 +- 2 files changed, 686 insertions(+), 303 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index b525bc601e..9bfc926b86 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -95,7 +95,7 @@ amrex::Real evaluate_tanh(const amrex::Real z) { AMREX_INLINE -void set_initial_condition(const int index_base, +void set_initial_condition(const int index_base, const int index_isentropic, const int npts_model_x, const amrex::Real xmin, const int npts_model_y, const amrex::Real ymin, const model_t& model_params, const int model_num) { @@ -149,11 +149,20 @@ void set_initial_condition(const int index_base, if (i <= index_base) { model::profile(model_num).state(i, j, model::itemp) = model_params.T_star; - } else { + } else if (i >= index_base && i <= index_isentropic) { model::profile(model_num).state(i, j, model::itemp) = model_params.T_star + 0.5_rt * (model_params.T_hi - model_params.T_star) * (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } else { + model::profile(model_num).state(i, j, model::itemp) = + model_params.T_hi - std::pow(xc, 3); + } + + // Floor temperature by T_lo + + if (model::profile(model_num).state(i, j, model::itemp) < model_params.T_lo) { + model::profile(model_num).state(i, j, model::itemp) = model_params.T_lo; } // @@ -188,89 +197,33 @@ void set_initial_condition(const int index_base, } -AMREX_INLINE -amrex::Real HSE_pressure(const int i, const int j, const int idir, - const amrex::Real dx, const amrex::Real dy, - const amrex::Real dens_zone, const amrex::Real Omega2, - const int model_num) { - // Evaluate the HSE pressure given density of the current zone - // using information of the surrounding zones. - - amrex::Real hse_p; - - if (idir == 0 || idir == 1) { - - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); - - if (idir == 0) { - // Integrate upward - - amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; - amrex::Real dens_im = 0.5_rt * - (dens_zone + model::profile(model_num).state(i-1, j, model::idens)); - amrex::Real dp_im1 = -dens_im * dx * (gravity::const_grav + - 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); - - hse_p = model::profile(model_num).state(i-1, j, model::ipres) - dp_im1; - - } else { - // Integrate downward - - amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; - amrex::Real dens_ip = 0.5_rt * - (dens_zone + model::profile(model_num).state(i+1, j, model::idens)); - amrex::Real dp_ip1 = dens_ip * dx * (gravity::const_grav + - 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); - - hse_p = model::profile(model_num).state(i+1, j, model::ipres) - dp_ip1; - } - - // Do an average between zone from bottom and left if not along the pole - if (j != 0) { - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); - amrex::Real dens_jm = 0.5_rt * - (dens_zone + model::profile(model_num).state(i, j-1, model::idens)); - amrex::Real dp_jm1 = -0.5_rt * dy * dens_jm * Omega2 * rc * rc * sinTwoThetam; - - hse_p = 0.5_rt * (hse_p + model::profile(model_num).state(i, j-1, model::ipres) - dp_jm1); - } - - } else if (idir == 2) { - // Integrate to the right - - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); - - amrex::Real dens_jm = 0.5_rt * - (dens_zone + model::profile(model_num).state(i, j-1, model::idens)); - amrex::Real dp_jm1 = -0.5_rt * dy * dens_jm * Omega2 * rc * rc * sinTwoThetam; - - hse_p = model::profile(model_num).state(i, j-1, model::ipres) - dp_jm1; - - } else { - amrex::Error("Invalid integrate direction"); - } - - return hse_p; -} - - AMREX_INLINE void HSE_integrate(const int i, const int j, const amrex::Real dx, const amrex::Real dy, - const int idir, + const int idir, const int index_isentropic, bool& isentropic, bool& fluff, - const amrex::Real entropy_base, + amrex::Real& entropy_base, const model_t& model_params, const int model_num) { // - // Goal is to update density, temperature, and pressure of the current zone. + // Goal is to update density, temperature, and pressure of the current zone based on integration // Also changes the `fluff` condition of the zone if density drops below // `low_density_cutoff`. // + // Comment out the following to get pure isothermal profile + + // if (i == index_isentropic + 1) { + // // First time reaching isentropic atmosphere from siothermal transition region + // // set isentropic to true + // isentropic = true; + // } + + // if (i > index_isentropic && !isentropic) { + // // If higher than the atmosphere region and its isothermal, + // // then we have reached the low temperature region, floor the temperature + // model::profile(model_num).state(i, j, model::itemp) = model_params.T_lo; + // } + // Get current state, this is the initial guess amrex::Real dens_zone = model::profile(model_num).state(i, j, model::idens); @@ -320,16 +273,33 @@ void HSE_integrate(const int i, const int j, amrex::Real entropy_zone; amrex::Real p_want; - amrex::Real s_want; amrex::Real drho; - amrex::Real dtemp = 0.0_rt; + amrex::Real dtemp; amrex::Real Omega2 = 0.0_rt; - if (castro::rotational_period > 0.0_rt && - castro::rotation_include_centrifugal) { + if (castro::rotational_period > 0.0_rt && castro::rotation_include_centrifugal) { Omega2 = std::pow(2.0_rt * M_PI / castro::rotational_period, 2); } + // Set index offset based on integration direction + int i0; + int j0; + if (idir == 0) { + // Integrate upward + i0 = -1; + j0 = 0; + } else if (idir == 1) { + // Integrate downward + i0 = +1; + j0 = 0; + } else if (idir == 2) { + // integrate to the right + i0 = 0; + j0 = -1; + } else { + amrex::Error("Invalid integration direction!"); + } + // Do Newton solve for (int iter = 0; iter < fw::MAX_ITER; iter++) { @@ -351,204 +321,638 @@ void HSE_integrate(const int i, const int j, // We assumed integration is always in the right dir, // hence reference state is always from left {i, j-1} - // Get HSE pressure based on integration direction + eos_t 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); - amrex::Real p_want = HSE_pressure(i, j, idir, dx, dy, dens_zone, Omega2, model_num); + entropy_zone = eos_state.s; + pres_zone = eos_state.p; - // - // We can either have isentropic or isothermal case. - // + // Get HSE pressure and the density derivative based on integration direction - if (isentropic) { - // isentropic case - // We have two unknowns: rho and T - // with two constraint eqns: entropy and pressure constraints + amrex::Real r = model::profile(model_num).x(i) + static_cast(i0) * 0.5_rt * dx; + amrex::Real twoTheta = 2.0_rt * model::profile(model_num).y(j) + static_cast(j0) * dy; + + amrex::Real dHSEp_drho; + + if (idir == 0 || idir == 1) { + // integrating in vertical direction - // Check if we have a valid base entropy when doing isentropic atmosphere + dHSEp_drho = -static_cast(i0) * 0.5_rt * dx * + (gravity::const_grav + 0.5_rt * Omega2 * r * (1.0_rt - std::cos(twoTheta))); + p_want = model::profile(model_num).state(i+i0, j, model::ipres) + + (dens_zone + model::profile(model_num).state(i+i0, j, model::idens)) * dHSEp_drho; + + } else if (idir == 2) { + // integrating in to the right + + dHSEp_drho = -static_cast(j0) * 0.25_rt * dy * + Omega2 * r * r * std::sin(twoTheta); + + p_want = model::profile(model_num).state(i, j+j0, model::ipres) + + (dens_zone + model::profile(model_num).state(i, j+j0, model::idens)) * dHSEp_drho; + } + + amrex::Real A = p_want - pres_zone; + amrex::Real dAdT = -eos_state.dpdT; + amrex::Real dAdrho = dHSEp_drho - eos_state.dpdr; + + // Find dtemp if dealing with isentropic atmosphere + + dtemp = 0.0_rt; + if (isentropic) { if (entropy_base <= 0.0_rt) { amrex::Error("isentropic atmosphere with negative entropy!"); } - eos_t 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]; - } + amrex::Real B = entropy_base - entropy_zone; + amrex::Real dBdT = -eos_state.dsdT; + amrex::Real dBdrho = -eos_state.dsdr; - eos(eos_input_rt, eos_state); + dtemp = (B - (dBdrho / dAdrho) * A) / + ((dBdrho / dAdrho) * dAdT - dBdT); + } - entropy_zone = eos_state.s; - pres_zone = eos_state.p; + drho = -(A + dAdT * dtemp) / dAdrho; - amrex::Real A = p_want - pres_zone; - amrex::Real B = entropy_base - entropy_zone; + // Update zone - amrex::Real dAdT = -eos_state.dpdT; - amrex::Real dAdrho = 0.0_rt; + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); - // Add d HSE_p / drho contribution - if (idir == 0 || idir == 1) { + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); + // Check if we entered isothermal region when integrating upward + if (isentropic && temp_zone < model_params.T_lo && idir == 0) { + temp_zone = model_params.T_lo; + isentropic = false; + } - if (idir == 0) { - // Integrate upward + // check if the density falls below our minimum cut-off -- + // if so, floor it + // Also detect whether we've entered the fluff region. + // Note that, its only possible when integrating upward. - amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; - dAdrho += 0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); + if (dens_zone < model_params.low_density_cutoff && idir == 0) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } - } else { - // Integrate downward + } // End of Newton solve - amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; - dAdrho += -0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); - } + if (!converged_hse) { + std::cout << "Error zone " << i << ", " << j << " did not converge in 2D HSE" << std::endl; + std::cout << dens_zone << " " << temp_zone << std::endl; + std::cout << p_want << " " << entropy_base << " " << entropy_zone << std::endl; + std::cout << drho << " " << dtemp << std::endl; + amrex::Error("Error: HSE non-convergence"); + } - // Do an average between zone from bottom and left if not along the pole - if (j != 0) { - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + // + // Once rho and T are determined, call eos to get pressure + // - dAdrho += 0.5_rt * (dAdrho + 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam); - } + eos_t 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]; + } - } else if (idir == 2) { - // Integrate to the right + eos(eos_input_rt, eos_state); - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + pres_zone = eos_state.p; - dAdrho += 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam; - } + // update the thermodynamics in this zone - dAdrho += - eos_state.dpdr; - amrex::Real dBdT = -eos_state.dsdT; - amrex::Real dBdrho = -eos_state.dsdr; + model::profile(model_num).state(i, j, model::idens) = dens_zone; + model::profile(model_num).state(i, j, model::itemp) = temp_zone; + model::profile(model_num).state(i, j, model::ipres) = pres_zone; - dtemp = (B - (dBdrho / dAdrho) * A) / - ((dBdrho / dAdrho) * dAdT - dBdT); + // If reaching the zone prior to the isentropic atmosphere, + // save the entropy, which will be used for isentropic integration + if (i == index_isentropic) { + entropy_base = eos_state.s; + } +} - drho = -(A + dAdT * dtemp) / dAdrho; - dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); - temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); +AMREX_INLINE +void +set_boundary_condition(const int npts_model_x, const int npts_model_y, + const amrex::Real dx, const amrex::Real dy, + const int index_base, const int index_isentropic, + bool fluff_mask[NPTS_X_MODEL][NPTS_Y_MODEL], + const model_t& model_params, const int model_num) { - // if (A < TOL .and. B < ETOL) then - if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { - converged_hse = true; - break; - } + // This function sets the boundary condition along θ_min, θ_max, + // r_min and r=H_star + 3 delta_atm. + // This defines the boundary condition for the bottom isothermal region. + + /* + Integrate the grid, and use the net result as boundary condition + or the initial guess for the smoothing procedure. + */ + + // + // integrate along r_min from θ_min (bottom boundary) + // + + bool fluff = false; + bool isentropic = false; + amrex::Real entropy_base = -1.0_rt; + + for (int j = 1; j < npts_model_y; j++) { + HSE_integrate(index_base, j, dx, dy, 2, + index_isentropic, + isentropic, fluff, entropy_base, + model_params, model_num); + fluff_mask[index_base][j] = fluff; + } + + // Integrate horizontally + + for (int j = 0; j < npts_model_y; ++j) { - // Check if we entered isothermal region when integrating upward - if (temp_zone < model_params.T_lo && idir == 0) { - temp_zone = model_params.T_lo; - isentropic = false; + fluff = false; + isentropic = false; + entropy_base = -1.0_rt; + + // First integrate downward from base + + for (int i = index_base-1; i >= 0; --i) { + HSE_integrate(i, j, dx, dy, 1, + index_isentropic, + isentropic, fluff, entropy_base, + model_params, model_num); + fluff_mask[i][j] = fluff; + } + + // Then integrate upward from base + + fluff = false; + isentropic = false; + entropy_base = -1.0_rt; + + for (int i = index_base+1; i < npts_model_x; i++) { + HSE_integrate(i, j, dx, dy, 0, + index_isentropic, + isentropic, fluff, entropy_base, + model_params, model_num); + fluff_mask[i][j] = fluff; + } + } + + // + // integrate along θ_min from r_base (left boundary) + // + // bool fluff = false; + // bool isentropic = false; + // amrex::Real entropy_base = -1.0_rt; + + // // First integrate downward from base + + // for (int i = index_base-1; i >= 0; --i) { + // HSE_integrate(i, 0, dx, dy, 1, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + // // Then integrate upward from base + + // fluff = false; + // isentropic = false; + // entropy_base = -1.0_rt; + + // for (int i = index_base+1; i < npts_model_x; i++) { + // HSE_integrate(i, 0, dx, dy, 0, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + // // + // // integrate along r_min from θ_min (bottom boundary) + // // + + // fluff = false; + // isentropic = false; + // entropy_base = -1.0_rt; + + // for (int j = 1; j < npts_model_y; j++) { + // HSE_integrate(0, j, dx, dy, 2, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + // // integrate along θ_max from r_min (right boundary) + + // fluff = false; + // isentropic = false; + // entropy_base = -1.0_rt; + + // // First integrate downward from base + + // for (int i = index_base-1; i >= 0; --i) { + // HSE_integrate(i, npts_model_y-1, dx, dy, 1, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + // // Then integrate upward from base + + // fluff = false; + // isentropic = false; + // entropy_base = -1.0_rt; + + // for (int i = index_base+1; i < npts_model_x; i++) { + // HSE_integrate(i, npts_model_y-1, dx, dy, 0, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + // // for (int i = 1; i < npts_model_x; i++) { + // // HSE_integrate(i, npts_model_y-1, dx, dy, 0, + // // index_isentropic, + // // isentropic, fluff, entropy_base, + // // model_params, model_num); + // // } + + // // Set upper boundary to fluff condition + // fluff= true; + // for (int j = 0; j < npts_model_y; ++j) { + // fluff_mask[npts_model_x-1][j] = fluff; + + // HSE_integrate(npts_model_x-1, j, dx, dy, 2, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + + // // integrate along r = H_star + 3 delta_atm (isothermal top boundary) + + // fluff = false; + // isentropic = false; + // entropy_base = -1.0_rt; + + // for (int j = 1; j < npts_model_y; j++) { + // HSE_integrate(index_isentropic+1, j, dx, dy, 2, + // index_isentropic, + // isentropic, fluff, entropy_base, + // model_params, model_num); + // } + + +} + +AMREX_INLINE +amrex::Real compute_relative_error(const int npts_model_x, const int npts_model_y, + const amrex::Real dx, const amrex::Real dy, + const bool fluff_mask[NPTS_X_MODEL][NPTS_Y_MODEL], + const model_t& model_params, const int model_num) { + + // Computes the relative error for the smoothing process + + amrex::Real Omega2 = 0.0_rt; + if (castro::rotational_period > 0.0_rt && castro::rotation_include_centrifugal) { + Omega2 = std::pow(2.0_rt * M_PI / castro::rotational_period, 2); + } + + amrex::Real residual_norm = 0.0_rt;; + amrex::Real f_norm = 0.0_rt; + + for (int i = 1; i < npts_model_x-1; ++i) { + for (int j = 1; j < npts_model_y-1; ++j) { + + // Ignore current zone if: + // 1) the surrounding zone is fluff then we're at the top boundary interface + // 2) the current zone is fluff + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j-1] || fluff_mask[i][j+1]) { + continue; } - // check if the density falls below our minimum cut-off -- - // if so, floor it - // Also detect whether we've entered the fluff region. - // Note that, its only possible when integrating upward. - - if (dens_zone < model_params.low_density_cutoff && idir == 0) { - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; - converged_hse = true; - fluff = true; - break; + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real rci = 1.0_rt / rc; + amrex::Real rc2i = 1.0_rt / (rc * rc); + amrex::Real rp = rc + 0.5_rt * dx; + amrex::Real rm = rc - 0.5_rt * dx; + amrex::Real rp2 = rp * rp; + amrex::Real rm2 = rm * rm; + + amrex::Real sinThetac = std::sin(model::profile(model_num).y(j)); + amrex::Real sinThetap = std::sin(model::profile(model_num).y(j) + 0.5_rt * dy); + amrex::Real sinThetam = std::sin(model::profile(model_num).y(j) - 0.5_rt * dy); + + amrex::Real cosThetac = std::cos(model::profile(model_num).y(j)); + amrex::Real cosThetap = std::cos(model::profile(model_num).y(j) + 0.5_rt * dy); + amrex::Real cosThetam = std::cos(model::profile(model_num).y(j) - 0.5_rt * dy); + + amrex::Real dxi = 1.0_rt / dx; + amrex::Real dx2i = 1.0_rt / (dx * dx); + amrex::Real dySinThetaci = 1.0_rt / (dy * sinThetac); + amrex::Real dy2sinThetaci = 1.0_rt / (dy * dy * sinThetac); + + // amrex::Real sinTwoThetac = std::sin(2.0_rt * model::profile(model_num).y(j)); + // amrex::Real sinTwoThetap = std::sin(2.0_rt * model::profile(model_num).y(j) + dy); + // amrex::Real sinTwoThetam = std::sin(2.0_rt * model::profile(model_num).y(j) - dy); + // amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(2.0_rt * model::profile(model_num).y(j)); + + amrex::Real psi_rc = Omega2 * rc * sinThetac * sinThetac + gravity::const_grav; + amrex::Real psi_rp = Omega2 * rp * sinThetac * sinThetac + gravity::const_grav; + amrex::Real psi_rm = Omega2 * rm * sinThetac * sinThetac + gravity::const_grav; + + amrex::Real psi_tc = Omega2 * rc * sinThetac * cosThetac; + amrex::Real psi_tp = Omega2 * rc * sinThetap * cosThetap; + amrex::Real psi_tm = Omega2 * rc * sinThetam * cosThetam; + + // amrex::Real psi_rc = 0.5_rt * Omega2 * rc * oneMinusCosTwoThetac + gravity::const_grav; + // amrex::Real psi_rp = 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac + gravity::const_grav; + // amrex::Real psi_rm = 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac + gravity::const_grav; + + // amrex::Real psi_tc = 0.5_rt * Omega2 * rc * sinTwoThetac; + // amrex::Real psi_tp = 0.5_rt * Omega2 * rc * sinTwoThetap; + // amrex::Real psi_tm = 0.5_rt * Omega2 * rc * sinTwoThetam; + + // amrex::Real psi_drho = psi_rc * (1.0_rt + (psi_tc * psi_tc) / (psi_rc * psi_rc)) * + // (model::profile(model_num).state(i+1, j, model::idens) - + // model::profile(model_num).state(i-1, j, model::idens)) / (2.0_rt * dx); + + // amrex::Real psi_drho = 0.5_rt * ((model::profile(model_num).state(i+1, j , model::idens) - + // model::profile(model_num).state(i-1, j , model::idens)) * + // psi_rc / dx + + // (model::profile(model_num).state(i , j+1, model::idens) - + // model::profile(model_num).state(i , j-1, model::idens)) * + // psi_tc / (rc * dy)); + + // amrex::Real dHSEp_drho = 2.0_rt * (Omega2 + gravity::const_grav / rc); + + // amrex::Real f = model::profile(model_num).state(i, j, model::idens) * dHSEp_drho + psi_drho; + + amrex::Real rhoip = 0.5_rt * (model::profile(model_num).state(i+1, j , model::idens) + + model::profile(model_num).state(i , j , model::idens)); + amrex::Real rhoim = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + + model::profile(model_num).state(i-1, j , model::idens)); + amrex::Real rhojp = 0.5_rt * (model::profile(model_num).state(i , j+1, model::idens) + + model::profile(model_num).state(i , j , model::idens)); + amrex::Real rhojm = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + + model::profile(model_num).state(i , j-1, model::idens)); + + amrex::Real f = (rp2 * rhoip * psi_rp - rm2 * rhoim * psi_rm) * rc2i * dxi + + (sinThetap * rhojp * psi_tp - sinThetam * rhojm * psi_tm) * rci * dySinThetaci; + + amrex::Real d2p = rc2i * + (dx2i * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + + rm2 * model::profile(model_num).state(i-1, j , model::ipres) - + (rp2 + rm2) * model::profile(model_num).state(i , j , model::ipres)) + + dy2sinThetaci * (sinThetap * model::profile(model_num).state(i , j+1, model::ipres) + + sinThetam * model::profile(model_num).state(i , j-1, model::ipres) - + (sinThetap + sinThetam) * model::profile(model_num).state(i , j , model::ipres))); + + amrex::Real residual = f - d2p; + + // Do L1 norm + if (std::abs(residual) > std::abs(f)) { + std::cout << "f in " << i << " " << j << " is " << f << std::endl; + std::cout << "d2p in " << i << " " << j << " is " << d2p << std::endl; } + f_norm += std::abs(f); + residual_norm += std::abs(residual); + } + } - } else { - // isothermal case - // Since T is fixed, the only unknown is density - // with one constraint eqn: pressure constraint + // Return relative error + return residual_norm / f_norm; +} - // (t, rho) -> (p) - eos_t 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]; - } +AMREX_INLINE +void HSE_smooth(const int i, const int j, + const amrex::Real dx, const amrex::Real dy, + bool fluff_mask[NPTS_X_MODEL][NPTS_Y_MODEL], + const model_t& model_params, const int model_num) { + // Smoothing procedure for HSE + + // Get current state, this is the initial guess + + amrex::Real dens_zone = model::profile(model_num).state(i, j, model::idens); + amrex::Real temp_zone = model::profile(model_num).state(i, j, model::itemp); + amrex::Real xn[NumSpec]; + for (int n = 0; n < NumSpec; n++) { + xn[n] = model::profile(model_num).state(i, j, model::ispec+n); + } + amrex::Real pres_zone; - eos(eos_input_rt, eos_state); + // Ignore current zone if: + // 1) the surrounding zone is fluff then we're at the top boundary interface + // 2) the current zone is fluff + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j-1] || fluff_mask[i][j+1]) { + return; + } - entropy_zone = eos_state.s; - pres_zone = eos_state.p; + // if (fluff_mask[i][j]) { + // If in fluff region, + // i.e. density of the zone is less than low_density_cutoff + // We stop HSE in this zone and set to a constant density and temperature. - amrex::Real A = p_want - pres_zone; - amrex::Real dAdrho = 0.0_rt; + // dens_zone = model_params.low_density_cutoff; + // temp_zone = model_params.T_lo; - // Add d HSE_p / drho contribution - if (idir == 0 || idir == 1) { + // // Call eos to get pressure: (T, rho) -> (p) - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(twoThetac); + // eos_t 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]; + // } - if (idir == 0) { - // Integrate upward + // eos(eos_input_rt, eos_state); - amrex::Real rm = model::profile(model_num).x(i) - 0.5_rt * dx; - dAdrho += 0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac); + // pres_zone = eos_state.p; - } else { - // Integrate downward + // // update the thermodynamics in this zone - amrex::Real rp = model::profile(model_num).x(i) + 0.5_rt * dx; - dAdrho += -0.5_rt * dx * (gravity::const_grav + 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac); - } + // model::profile(model_num).state(i, j, model::idens) = dens_zone; + // model::profile(model_num).state(i, j, model::itemp) = temp_zone; + // model::profile(model_num).state(i, j, model::ipres) = pres_zone; + // return; + // } - // Do an average between zone from bottom and left if not along the pole - if (j != 0) { - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + // + // If we are not in fluff region, then we have to do Newton solve + // - dAdrho += 0.5_rt * (dAdrho + 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam); - } + // start off the Newton loop by saying that the zone has not converged - } else if (idir == 2) { - // Integrate to the right + bool converged_hse = false; - amrex::Real rc = model::profile(model_num).x(i); - amrex::Real twoThetac = 2.0_rt * model::profile(model_num).y(j); - amrex::Real sinTwoThetam = std::sin(twoThetac - dy); + amrex::Real entropy_zone; + amrex::Real p_want; + amrex::Real drho; + amrex::Real dtemp; - dAdrho += 0.25_rt * dy * Omega2 * rc * rc * sinTwoThetam; - } + // precompute some terms for evaluating HSE pressure. + amrex::Real Omega2 = 0.0_rt; + if (castro::rotational_period > 0.0_rt && castro::rotation_include_centrifugal) { + Omega2 = std::pow(2.0_rt * M_PI / castro::rotational_period, 2); + } - dAdrho += - eos_state.dpdr; - drho = - A / dAdrho; + amrex::Real rc = model::profile(model_num).x(i); + amrex::Real rc2 = rc * rc; + amrex::Real rci = 1.0_rt / rc; + amrex::Real rc2i = 1.0_rt / (rc * rc); + amrex::Real rp = rc + 0.5_rt * dx; + amrex::Real rm = rc - 0.5_rt * dx; + amrex::Real rp2 = rp * rp; + amrex::Real rm2 = rm * rm; - dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + amrex::Real sinThetac = std::sin(model::profile(model_num).y(j)); + amrex::Real sinThetap = std::sin(model::profile(model_num).y(j) + 0.5_rt * dy); + amrex::Real sinThetam = std::sin(model::profile(model_num).y(j) - 0.5_rt * dy); - if (std::abs(drho) < fw::TOL * dens_zone) { - converged_hse = true; - break; - } + amrex::Real cosThetac = std::cos(model::profile(model_num).y(j)); + amrex::Real cosThetap = std::cos(model::profile(model_num).y(j) + 0.5_rt * dy); + amrex::Real cosThetam = std::cos(model::profile(model_num).y(j) - 0.5_rt * dy); - // fluff region only possible when integrate upward - if (dens_zone < model_params.low_density_cutoff && idir == 0) { - dens_zone = model_params.low_density_cutoff; - temp_zone = model_params.T_lo; - converged_hse = true; - fluff = true; - break; - } + amrex::Real dx2 = dx * dx; + amrex::Real dxi = 1.0_rt / dx; + amrex::Real dx2i = 1.0_rt / (dx * dx); + amrex::Real dy2sinThetac = dy * dy * sinThetac; + amrex::Real dySinThetaci = 1.0_rt / (dy * sinThetac); + amrex::Real dy2sinThetaci = 1.0_rt / (dy * dy * sinThetac); + + // amrex::Real sinTwoThetac = std::sin(2.0_rt * model::profile(model_num).y(j)); + // amrex::Real sinTwoThetap = std::sin(2.0_rt * model::profile(model_num).y(j) + dy); + // amrex::Real sinTwoThetam = std::sin(2.0_rt * model::profile(model_num).y(j) - dy); + // amrex::Real oneMinusCosTwoThetac = 1.0_rt - std::cos(2.0_rt * model::profile(model_num).y(j)); + + amrex::Real psi_rc = Omega2 * rc * sinThetac * sinThetac + gravity::const_grav; + amrex::Real psi_rp = Omega2 * rp * sinThetac * sinThetac + gravity::const_grav; + amrex::Real psi_rm = Omega2 * rm * sinThetac * sinThetac + gravity::const_grav; + + amrex::Real psi_tc = Omega2 * rc * sinThetac * cosThetac; + amrex::Real psi_tp = Omega2 * rc * sinThetap * cosThetap; + amrex::Real psi_tm = Omega2 * rc * sinThetam * cosThetam; + + // amrex::Real psi_rc = 0.5_rt * Omega2 * rc * oneMinusCosTwoThetac + gravity::const_grav; + // amrex::Real psi_rp = 0.5_rt * Omega2 * rp * oneMinusCosTwoThetac + gravity::const_grav; + // amrex::Real psi_rm = 0.5_rt * Omega2 * rm * oneMinusCosTwoThetac + gravity::const_grav; + + // amrex::Real psi_tc = 0.5_rt * Omega2 * rc * sinTwoThetac; + // amrex::Real psi_tp = 0.5_rt * Omega2 * rc * sinTwoThetap; + // amrex::Real psi_tm = 0.5_rt * Omega2 * rc * sinTwoThetam; + + // amrex::Real fac = psi_rc * (1.0_rt + (psi_tc * psi_tc) / (psi_rc * psi_rc)) * + // (model::profile(model_num).state(i+1, j, model::idens) - + // model::profile(model_num).state(i-1, j, model::idens)) / (2.0_rt * dx); + + // amrex::Real dHSEp_drho = 2.0_rt * (Omega2 + gravity::const_grav / rc); + + // Do Newton solve + + for (int iter = 0; iter < fw::MAX_ITER; iter++) { + + eos_t 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); + + entropy_zone = eos_state.s; + pres_zone = eos_state.p; + + // Get HSE pressure and the density derivative + + // amrex::Real f = dens_zone * dHSEp_drho + fac; - } // End of isothermal case + amrex::Real rhoip = 0.5_rt * (model::profile(model_num).state(i+1, j , model::idens) + + model::profile(model_num).state(i , j , model::idens)); + amrex::Real rhoim = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + + model::profile(model_num).state(i-1, j , model::idens)); + amrex::Real rhojp = 0.5_rt * (model::profile(model_num).state(i , j+1, model::idens) + + model::profile(model_num).state(i , j , model::idens)); + amrex::Real rhojm = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + + model::profile(model_num).state(i , j-1, model::idens)); + + amrex::Real f = (rp2 * rhoip * psi_rp - rm2 * rhoim * psi_rm) * rc2i * dxi + + (sinThetap * rhojp * psi_tp - sinThetam * rhojm * psi_tm) * rci * dySinThetaci; + + amrex::Real dHSEp_drho = 0.5_rt * ((rp2 * psi_rp - rm2 * psi_rm) * rc2i * dxi + + (sinThetap * psi_tp - sinThetam * psi_tm) * rci * dySinThetaci); + + p_want = ( + dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + + rm2 * model::profile(model_num).state(i-1, j , model::ipres)) + + dx2 * (sinThetap * model::profile(model_num).state(i , j+1, model::ipres) + + sinThetam * model::profile(model_num).state(i , j-1, model::ipres)) - + rc2 * dx2 * dy2sinThetac * f + ) / ( dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap + sinThetam) ); + + amrex::Real A = p_want - pres_zone; + amrex::Real dAdT = -eos_state.dpdT; + amrex::Real dAdrho = dHSEp_drho - eos_state.dpdr; + + // Find dtemp if dealing with isentropic atmosphere + + dtemp = 0.0_rt; + drho = -(A + dAdT * dtemp) / dAdrho; + + // Update zone + + dens_zone = amrex::Clamp(dens_zone + drho , 0.9_rt * dens_zone, 1.1_rt * dens_zone); + temp_zone = amrex::Clamp(temp_zone + dtemp, 0.9_rt * temp_zone, 1.1_rt * temp_zone); + + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } + + // check if the density falls below our minimum cut-off -- + // if so, floor it + // Also detect whether we've entered the fluff region. + // Note that, its only possible when integrating upward. + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff_mask[i][j] = true; + break; + } } // End of Newton solve if (!converged_hse) { std::cout << "Error zone " << i << ", " << j << " did not converge in 2D HSE" << std::endl; std::cout << dens_zone << " " << temp_zone << std::endl; - std::cout << p_want << " " << entropy_base << " " << entropy_zone << std::endl; + std::cout << p_want << " " << " " << entropy_zone << std::endl; std::cout << drho << " " << dtemp << std::endl; amrex::Error("Error: HSE non-convergence"); } @@ -614,12 +1018,11 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr } // find the index of the base height - // x is assumed to be the r-dir int index_base = -1; for (int i = 0; i < npts_model_x; i++) { if (model::profile(model_num).x(i) >= xmin + model_params.H_star) { - index_base = i+1; + index_base = i; break; } } @@ -628,114 +1031,94 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr amrex::Error("ERROR: invalid base_height"); } + // Find the index of isentropic atmosphere, + // this is to find the entropy of the last zone of the + // isothermal layer prior to the isentropic atmosphere + + int index_isentropic = -1; + for (int i = 0; i < npts_model_x; i++) { + if (model::profile(model_num).x(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { + index_isentropic = i - 1; + if (i == index_base + 1) { + amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); + } + break; + } + } + // // put the model onto our new uniform grid // This is the initial guess conditions. // - set_initial_condition(index_base, + set_initial_condition(index_base, index_isentropic, npts_model_x, xmin, npts_model_y, ymin, model_params, model_num); // - // 2D HSE + entropy solve + // set boundary conditions // - // First integrate to the right from θ=0 at the base layer. - // The base layer is isothermal. + bool fluff_mask[NPTS_X_MODEL][NPTS_Y_MODEL] = {false}; - bool fluff = false; - bool isentropic = false; - bool flipped = false; - amrex::Real entropy_base = -1.0_rt; + set_boundary_condition(npts_model_x, npts_model_y, dx, dy, + index_base, index_isentropic, fluff_mask, + model_params, model_num); - for (int j = 1; j < npts_model_y; ++j) { - HSE_integrate(index_base, j, dx, dy, 2, - isentropic, fluff, - entropy_base, - model_params, model_num); - } + amrex::Real relative_error = compute_relative_error(npts_model_x, npts_model_y, + dx, dy, fluff_mask, + model_params, model_num); + std::cout << "initial relative error is " << relative_error << std::endl; - // Then integrate vertically for all θ. + // + // Multigrid Relaxation + Newton solve + // - for (int j = 0; j < npts_model_y; ++j) { + int smooth_count = 0; + int converging_counter = 0; - // Integrate vertically downward from base - // Regions below the base is isothermal + while (converging_counter < 20) { - // Reset conditions for each θ + // Do smoothing over the interior grid - fluff = false; - isentropic = false; - flipped = false; - entropy_base = -1.0_rt; + for (int i = 1; i < npts_model_x-1; ++i) { + for (int j = 1; j < npts_model_y-1; ++j) { - for (int i = index_base-1; i >= 0; --i) { - HSE_integrate(i, j, dx, dy, 1, - isentropic, fluff, - entropy_base, - model_params, model_num); + // Ignore current zone if: + // 1) the surrounding zone is fluff then we're at the top boundary interface + // 2) the current zone is fluff + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j-1] || fluff_mask[i][j+1]) { + continue; + } + HSE_smooth(i, j, dx, dy, + fluff_mask, + model_params, model_num); + } } - // Start from the base layer and integrate upward - // We start out isothermal during the hyperbolic transition region - // then we 'flip' to isentropic - // the HSE state will be done putting creating an isentropic state until - // the temperature goes below T_lo -- then we will do isothermal. - // also, once the density goes below low_density_cutoff, we stop HSE + // Compute relative error - // Reset conditions for each θ - - fluff = false; - isentropic = false; - flipped = false; - entropy_base = -1.0_rt; + if (smooth_count % 10 == 0) { + amrex::Real new_relative_error = compute_relative_error(npts_model_x, npts_model_y, + dx, dy, fluff_mask, + model_params, model_num); + std::cout << "relative error is " << new_relative_error << std::endl; - for (int i = index_base+1; i < npts_model_x; i++) { - - // If above hyperbolic transition region, we can either be isentropic or isothermal, - - if (model::profile(model_num).x(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { - - if (flipped && !isentropic) { - model::profile(model_num).state(i, j, model::itemp) = model_params.T_lo; - - } else if (!flipped) { - // This is to determine when did we get out of the hyperbolic tangent transition phase (isothermal) - // Then get the appropriate entropy for doing isentropic atmosphere - if (i == index_base + 1) { - amrex::Error("atm_delta is too small for the grid spacing (at least one transition zone is required)"); - } - isentropic = true; - flipped = true; - - // now we need to know the entropy we are confining ourselves to - - eos_t eos_state; - eos_state.rho = model::profile(model_num).state(i-1, j, model::idens); - eos_state.T = model::profile(model_num).state(i-1, j, model::itemp); - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = model::profile(model_num).state(i-1, j, model::ispec+n); - } - - eos(eos_input_rt, eos_state); - - // Get entropy for isentropic atmosphere - - entropy_base = eos_state.s; - - amrex::Print() << "base density = " << eos_state.rho - << " and base Temperature " << eos_state.T << std::endl; - } + if (relative_error < new_relative_error) { + ++converging_counter; } - HSE_integrate(i, j, dx, dy, 0, - isentropic, fluff, - entropy_base, - model_params, model_num); + relative_error = new_relative_error; } - } + ++smooth_count; + if (smooth_count > 1e5) { + break; + } + } + std::cout << "Number of relaxation iterations took are " << smooth_count << std::endl; } #endif diff --git a/Exec/science/xrb_spherical/problem_initialize.H b/Exec/science/xrb_spherical/problem_initialize.H index e480b950cb..167cb546b1 100644 --- a/Exec/science/xrb_spherical/problem_initialize.H +++ b/Exec/science/xrb_spherical/problem_initialize.H @@ -174,8 +174,8 @@ void problem_initialize () generate_initial_model(nx_model + ng, problo[0] - ng * dx_model, probhi[0], - ny_model + ng, - problo[1] - ng * dy_model, + ny_model, + problo[1], probhi[1], model_params, 0); @@ -197,12 +197,12 @@ void problem_initialize () // set the ambient state for the upper boundary condition using upper most r zone at θ=0 - ambient::ambient_state[URHO] = model::profile(0).state(model::npts_x-1, ng, model::idens); + ambient::ambient_state[URHO] = model::profile(0).state(model::npts_x-1, 0, model::idens); - ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts_x-1, ng, model::itemp); + ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts_x-1, 0, model::itemp); for (int n = 0; n < NumSpec; n++) { ambient::ambient_state[UFS+n] = - ambient::ambient_state[URHO] * model::profile(0).state(model::npts_x-1, ng, model::ispec+n); + ambient::ambient_state[URHO] * model::profile(0).state(model::npts_x-1, 0, model::ispec+n); } ambient::ambient_state[UMX] = 0.0_rt; From cfd1756c9caa2b7d8bb36b6c3e44bd7e5899e0ba Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 3 Nov 2025 14:39:48 -0500 Subject: [PATCH 19/20] fix dp drho --- Exec/science/xrb_spherical/initial_model.H | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index 9bfc926b86..467dec517e 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -467,9 +467,8 @@ set_boundary_condition(const int npts_model_x, const int npts_model_y, or the initial guess for the smoothing procedure. */ - // - // integrate along r_min from θ_min (bottom boundary) - // + + // Integrate horizontally along the base bool fluff = false; bool isentropic = false; @@ -483,7 +482,7 @@ set_boundary_condition(const int npts_model_x, const int npts_model_y, fluff_mask[index_base][j] = fluff; } - // Integrate horizontally + // Integrate vertically for (int j = 0; j < npts_model_y; ++j) { @@ -516,6 +515,7 @@ set_boundary_condition(const int npts_model_x, const int npts_model_y, } } + // // integrate along θ_min from r_base (left boundary) // @@ -731,10 +731,10 @@ amrex::Real compute_relative_error(const int npts_model_x, const int npts_model_ amrex::Real residual = f - d2p; // Do L1 norm - if (std::abs(residual) > std::abs(f)) { - std::cout << "f in " << i << " " << j << " is " << f << std::endl; - std::cout << "d2p in " << i << " " << j << " is " << d2p << std::endl; - } + // if (std::abs(residual) > std::abs(f)) { + // std::cout << "f in " << i << " " << j << " is " << f << std::endl; + // std::cout << "d2p in " << i << " " << j << " is " << d2p << std::endl; + // } f_norm += std::abs(f); residual_norm += std::abs(residual); } @@ -903,8 +903,10 @@ void HSE_smooth(const int i, const int j, amrex::Real f = (rp2 * rhoip * psi_rp - rm2 * rhoim * psi_rm) * rc2i * dxi + (sinThetap * rhojp * psi_tp - sinThetam * rhojm * psi_tm) * rci * dySinThetaci; - amrex::Real dHSEp_drho = 0.5_rt * ((rp2 * psi_rp - rm2 * psi_rm) * rc2i * dxi + - (sinThetap * psi_tp - sinThetam * psi_tm) * rci * dySinThetaci); + amrex::Real df_drho = 0.5_rt * ((rp2 * psi_rp - rm2 * psi_rm) * rc2i * dxi + + (sinThetap * psi_tp - sinThetam * psi_tm) * rci * dySinThetaci); + amrex::Real dHSEp_drho = - rc2 * dx2 * dy2sinThetac * df_drho / + (dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap + sinThetam)); p_want = ( dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + From ef3a8bb29dd81645a7b07d0398a96e873586ca90 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 3 Nov 2025 18:13:41 -0500 Subject: [PATCH 20/20] allow left and right boundary to vary by assuming outflow condition --- Exec/science/xrb_spherical/initial_model.H | 116 ++++++++++++++++++--- 1 file changed, 99 insertions(+), 17 deletions(-) diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H index 467dec517e..c150ffb04c 100644 --- a/Exec/science/xrb_spherical/initial_model.H +++ b/Exec/science/xrb_spherical/initial_model.H @@ -717,6 +717,12 @@ amrex::Real compute_relative_error(const int npts_model_x, const int npts_model_ amrex::Real rhojm = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + model::profile(model_num).state(i , j-1, model::idens)); + if (j == 0) { + rhojm = model::profile(model_num).state(i , j , model::idens); + } else if (j == model::npts_y - 1) { + rhojp = model::profile(model_num).state(i , j , model::idens); + } + amrex::Real f = (rp2 * rhoip * psi_rp - rm2 * rhoim * psi_rm) * rc2i * dxi + (sinThetap * rhojp * psi_tp - sinThetam * rhojm * psi_tm) * rci * dySinThetaci; @@ -728,6 +734,28 @@ amrex::Real compute_relative_error(const int npts_model_x, const int npts_model_ sinThetam * model::profile(model_num).state(i , j-1, model::ipres) - (sinThetap + sinThetam) * model::profile(model_num).state(i , j , model::ipres))); + if (j == 0) { + d2p = rc2i * + ( + dx2i * (rp2 * (model::profile(model_num).state(i+1, j , model::ipres) - + model::profile(model_num).state(i , j , model::ipres)) - + rm2 * (model::profile(model_num).state(i , j , model::ipres) - + model::profile(model_num).state(i-1, j , model::ipres))) + + dy2sinThetaci * sinThetap * (model::profile(model_num).state(i , j+1, model::ipres) - + model::profile(model_num).state(i , j , model::ipres)) + ); + } else if (j == model::npts_y - 1){ + d2p = rc2i * + ( + dx2i * (rp2 * (model::profile(model_num).state(i+1, j , model::ipres) - + model::profile(model_num).state(i , j , model::ipres)) - + rm2 * (model::profile(model_num).state(i , j , model::ipres) - + model::profile(model_num).state(i-1, j , model::ipres))) - + dy2sinThetaci * sinThetam * (model::profile(model_num).state(i , j , model::ipres) - + model::profile(model_num).state(i , j-1, model::ipres)) + ); + } + amrex::Real residual = f - d2p; // Do L1 norm @@ -900,21 +928,59 @@ void HSE_smooth(const int i, const int j, amrex::Real rhojm = 0.5_rt * (model::profile(model_num).state(i , j , model::idens) + model::profile(model_num).state(i , j-1, model::idens)); + if (j == 0) { + rhojm = model::profile(model_num).state(i , j , model::idens); + } else if (j == model::npts_y - 1) { + rhojp = model::profile(model_num).state(i , j , model::idens); + } + amrex::Real f = (rp2 * rhoip * psi_rp - rm2 * rhoim * psi_rm) * rc2i * dxi + (sinThetap * rhojp * psi_tp - sinThetam * rhojm * psi_tm) * rci * dySinThetaci; amrex::Real df_drho = 0.5_rt * ((rp2 * psi_rp - rm2 * psi_rm) * rc2i * dxi + (sinThetap * psi_tp - sinThetam * psi_tm) * rci * dySinThetaci); + if (j == 0){ + df_drho = 0.5_rt * ((rp2 * psi_rp - rm2 * psi_rm) * rc2i * dxi + + (sinThetap * psi_tp - 2.0_rt * sinThetam * psi_tm) * rci * dySinThetaci); + } else if (j == model::npts_y - 1) { + df_drho = 0.5_rt * ((rp2 * psi_rp - 2.0_rt * rm2 * psi_rm) * rc2i * dxi + + (2.0_rt * sinThetap * psi_tp - sinThetam * psi_tm) * rci * dySinThetaci); + } + + if (j == 0) { + p_want = ( + dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + + rm2 * model::profile(model_num).state(i-1, j , model::ipres)) + + dx2 * (sinThetap * model::profile(model_num).state(i , j+1, model::ipres)) - + rc2 * dx2 * dy2sinThetac * f + ) / ( dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap) ); + } else if (j == model::npts_y - 1) { + p_want = ( + dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + + rm2 * model::profile(model_num).state(i-1, j , model::ipres)) + + dx2 * (sinThetam * model::profile(model_num).state(i , j-1, model::ipres)) - + rc2 * dx2 * dy2sinThetac * f + ) / ( dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetam) ); + + } else { + p_want = ( + dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + + rm2 * model::profile(model_num).state(i-1, j , model::ipres)) + + dx2 * (sinThetap * model::profile(model_num).state(i , j+1, model::ipres) + + sinThetam * model::profile(model_num).state(i , j-1, model::ipres)) - + rc2 * dx2 * dy2sinThetac * f + ) / ( dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap + sinThetam) ); + } + amrex::Real dHSEp_drho = - rc2 * dx2 * dy2sinThetac * df_drho / (dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap + sinThetam)); - - p_want = ( - dy2sinThetac * (rp2 * model::profile(model_num).state(i+1, j , model::ipres) + - rm2 * model::profile(model_num).state(i-1, j , model::ipres)) + - dx2 * (sinThetap * model::profile(model_num).state(i , j+1, model::ipres) + - sinThetam * model::profile(model_num).state(i , j-1, model::ipres)) - - rc2 * dx2 * dy2sinThetac * f - ) / ( dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap + sinThetam) ); + if (j == 0) { + dHSEp_drho = - rc2 * dx2 * dy2sinThetac * df_drho / + (dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetap)); + } else if (j == model::npts_y - 1) { + dHSEp_drho = - rc2 * dx2 * dy2sinThetac * df_drho / + (dy2sinThetac * (rp2 + rm2) + dx2 * (sinThetam)); + } amrex::Real A = p_want - pres_zone; amrex::Real dAdT = -eos_state.dpdT; @@ -1085,16 +1151,32 @@ generate_initial_model(const int npts_model_x, const amrex::Real xmin, const amr // Do smoothing over the interior grid for (int i = 1; i < npts_model_x-1; ++i) { - for (int j = 1; j < npts_model_y-1; ++j) { + for (int j = 0; j < npts_model_y; ++j) { + + // Ignore current zone if: + // 1) the surrounding zone is fluff then we're at the top boundary interface + // 2) the current zone is fluff + + if (j == 0) { + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j+1]) { + continue; + } + } else if (j == npts_model_y - 1) { + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j-1]) { + continue; + } + } else { + if (fluff_mask[i][j] || + fluff_mask[i+1][j] || fluff_mask[i-1][j] || + fluff_mask[i][j-1] || fluff_mask[i][j+1]) { + continue; + } + } - // Ignore current zone if: - // 1) the surrounding zone is fluff then we're at the top boundary interface - // 2) the current zone is fluff - if (fluff_mask[i][j] || - fluff_mask[i+1][j] || fluff_mask[i-1][j] || - fluff_mask[i][j-1] || fluff_mask[i][j+1]) { - continue; - } HSE_smooth(i, j, dx, dy, fluff_mask, model_params, model_num);