diff --git a/nse_solver/nse_compatibility/GNUmakefile b/nse_solver/nse_compatibility/GNUmakefile new file mode 100644 index 0000000000..e9ffffeb5a --- /dev/null +++ b/nse_solver/nse_compatibility/GNUmakefile @@ -0,0 +1,44 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 3 + +COMP = gnu + +USE_MPI = FALSE +USE_OMP = FALSE + +USE_REACT = TRUE + +USE_NSE_NET = TRUE +EBASE = main + +# define the location of the Microphysics top directory +MICROPHYSICS_HOME ?= ../.. + +# This sets the EOS directory +EOS_DIR := helmholtz + +# This sets the network directory +NETWORK_DIR := subch_base + +SCREEN_METHOD = chabrier1998 + +CONDUCTIVITY_DIR := stellar + +INTEGRATOR_DIR = VODE + +ifeq ($(USE_CUDA), TRUE) + INTEGRATOR_DIR := VODE +endif + +EXTERN_SEARCH += . + +Bpack := ./Make.package +Blocs := . + +include $(MICROPHYSICS_HOME)/unit_test/Make.unit_test + + diff --git a/nse_solver/nse_compatibility/Make.package b/nse_solver/nse_compatibility/Make.package new file mode 100644 index 0000000000..d87765520d --- /dev/null +++ b/nse_solver/nse_compatibility/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp +CEXE_headers += nse_network_compatibility.H diff --git a/nse_solver/nse_compatibility/README.md b/nse_solver/nse_compatibility/README.md new file mode 100644 index 0000000000..58f99d6fe4 --- /dev/null +++ b/nse_solver/nse_compatibility/README.md @@ -0,0 +1,14 @@ +# NSE Network Checking Script + +This is a script for checking whether integrated mass fractions of +a reaction network will eventually match with the mass fractions +calculated via NSE equations. Currently, this is only valid for networks +that do not involve weak rates. + +* Note that the integration doesn't achieve NSE even though we have +USE_NSE_NET=TRUE is due to the reference size of the cell. By setting +nse.nse_dx_independet=1 will fix this. + +Script will print out all burn_state cases when nse_dx_independent is +not enabled. It will only print out burn_state cases that didn't enter NSE +when nse_dx_independent=1 \ No newline at end of file diff --git a/nse_solver/nse_compatibility/_parameters b/nse_solver/nse_compatibility/_parameters new file mode 100644 index 0000000000..2b816bf783 --- /dev/null +++ b/nse_solver/nse_compatibility/_parameters @@ -0,0 +1,22 @@ +@namespace: unit_test + +run_prefix character "" + +# the final time to integrate to +tmax real 1.e3 + +# first output time -- we will output in nsteps logarithmically spaced steps between [tfirst, tmax] +tfirst real 0.0 + +# number of steps (logarithmically spaced) +nsteps integer 100 + +rho_min real 1.e7 +rho_max real 1.e9 + +nrho integer 4 + +T_min real 6.e9 +T_max real 8.e9 + +nT integer 4 diff --git a/nse_solver/nse_compatibility/inputs_nse_compatibility b/nse_solver/nse_compatibility/inputs_nse_compatibility new file mode 100644 index 0000000000..14be0cce4c --- /dev/null +++ b/nse_solver/nse_compatibility/inputs_nse_compatibility @@ -0,0 +1,26 @@ +# this inputs is used for checking the compatibility of networks with NSE_NET + +unit_test.run_prefix = "test_nse_network_compatibility_" + +integrator.burner_verbose = 0 + +# Set which jacobian to use +# 1 = analytic jacobian +# 2 = numerical jacobian + +integrator.jacobian = 1 + +integrator.renormalize_abundances = 0 + +integrator.rtol_spec = 1.0e-6 +integrator.rtol_enuc = 1.0e-6 +integrator.atol_spec = 1.0e-10 +integrator.atol_enuc = 1.0e-6 + +unit_test.tmax = 1.e3 +unit_test.tfirst = 1.e-10 +unit_test.nsteps = 100 + +integrator.integrate_energy = 0 +integrator.call_eos_in_rhs = 0 +integrator.do_species_clip = 0 \ No newline at end of file diff --git a/nse_solver/nse_compatibility/main.cpp b/nse_solver/nse_compatibility/main.cpp new file mode 100644 index 0000000000..2581910048 --- /dev/null +++ b/nse_solver/nse_compatibility/main.cpp @@ -0,0 +1,34 @@ +#include +#include +#include + +#include +#include +using namespace amrex; + +#include +#include +#include +#include +#include + +int main(int argc, char *argv[]) { + + amrex::Initialize(argc, argv); + + std::cout << "starting the single zone burn..." << std::endl; + + ParmParse ppa("amr"); + + init_unit_test(); + + // C++ EOS initialization (must be done after Fortran eos_init and init_extern_parameters) + eos_init(); + + // C++ Network, RHS, screening, rates initialization + network_init(); + + nse_network_compatibility(); + + amrex::Finalize(); +} diff --git a/nse_solver/nse_compatibility/nse_network_compatibility.H b/nse_solver/nse_compatibility/nse_network_compatibility.H new file mode 100644 index 0000000000..2d3471d95e --- /dev/null +++ b/nse_solver/nse_compatibility/nse_network_compatibility.H @@ -0,0 +1,196 @@ +#ifndef NSE_NETWORK_COMPATIBILITY_H +#define NSE_NETWORK_COMPATIBILITY_H + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace unit_test_rp; + +void burn_cell_c(const eos_t& eos_state) +{ + // Given eos_state with density, temperature and mass fractions + + burn_t burn_state; + + burn_state.rho = eos_state.rho; + burn_state.T = eos_state.T; + for (int n = 0; n < NumSpec; ++n) { + burn_state.xn[n] = eos_state.xn[n]; + } + burn_state.y_e = eos_state.y_e; + + burn_state.i = 0; + burn_state.j = 0; + burn_state.k = 0; + burn_state.T_fixed = -1.0_rt; + + // normalize -- just in case + normalize_abundances_burn(burn_state); + + // call the EOS to set initial e -- it actually doesn't matter to + // the burn but we need to keep track of e to get a valid + // temperature for the burn if we substep + + eos(eos_input_rt, burn_state); + + // we will divide the total integration time into nsteps that are + // logarithmically spaced + + if (tfirst == 0.0_rt) { + if (nsteps == 1) { + tfirst = tmax; + } else { + tfirst = tmax / nsteps; + } + } + Real dlogt = 0.0_rt; + if (nsteps == 1) { + dlogt = (std::log10(tmax) - std::log10(tfirst)); + } else { + dlogt = (std::log10(tmax) - std::log10(tfirst)) / (nsteps - 1); + } + + // save the initial state -- we'll use this to determine + // how much things changed over the entire burn + + burn_t burn_state_in = burn_state; + + Real t = 0.0; + + // store the initial internal energy -- we'll update this after + // each substep + + Real energy_initial = burn_state.e; + + // loop over steps, burn, and output the current state + + for (int n = 0; n < nsteps; n++){ + + // compute the time we wish to integrate to + + Real tend = std::pow(10.0_rt, std::log10(tfirst) + dlogt * n); + Real dt = tend - t; + + burner(burn_state, dt); + + if (! burn_state.success) { + amrex::Error("integration failed"); + } + + // state.e represents the change in energy over the burn (for + // just this sybcycle), so turn it back into a physical energy + + burn_state.e += energy_initial; + + // reset the initial energy for the next subcycle + + energy_initial = burn_state.e; + + t += dt; + } + + // Now find the mass fraction via NSE calculation + burn_state.mu_p = -3.0_rt; + burn_state.mu_n = -12.0_rt; + auto nse_state = get_actual_nse_state(burn_state); + + if (nse_dx_independent != 1 || (nse_dx_independent == 1 && !burn_state.nse)) { + std::cout << std::scientific << std::setprecision(3); + std::cout << "------------------------------------" << std::endl; + std::cout << " - added e = " << burn_state.e - burn_state_in.e << std::endl; + std::cout << " Initial vs. Final T = " << burn_state_in.T << " " + << burn_state.T << std::endl; + std::cout << " Initial vs. Final density = " << burn_state_in.rho << " " + << burn_state.rho << std::endl; + std::cout << " Initial vs. Final vs. NSE y_e = " << burn_state_in.y_e << " " + << burn_state.y_e << " " << nse_state.y_e << std::endl; + std::cout << "------------------------------------" << std::endl; + std::cout << "Element" << std::setw(14) << "Burn Xs" << std::setw(20) << "NSE Xs " + << std::setw(20) << "abs err" << std::setw(20) << "rel err" << std::endl; + for (int n = 0; n < NumSpec; ++n) { + const std::string& element = short_spec_names_cxx[n]; + Real abs_err = std::abs(burn_state.xn[n] - nse_state.xn[n]); + Real rel_err = abs_err / burn_state.xn[n]; + std::cout << " " << element << std::setw(20-element.size()) + << burn_state.xn[n] << std::setw(20) + << nse_state.xn[n] << std::setw(20) << abs_err << std::setw(20) + << rel_err << std::endl; + } + std::cout << "------------------------------------" << std::endl; + } +} + +void nse_network_compatibility() +{ + // This function compares the integrated equilibrium mass fractions + // with the results from NSE calculations given any network. + // Note that this network should not have any weak rates. + + Real dlogrho = (std::log10(rho_max) - std::log10(rho_min))/static_cast(nrho-1); + Real dlogT = (std::log10(T_max) - std::log10(T_min))/static_cast(nT-1); + + // Create massfraction 2d array to hold 3 different sets of initial mass fractions + // Corresponding to ye < 0.5, ye = 0.5, and ye > 0.5 + + amrex::Array2D massfractions; + + bool condition1_met = false; + bool condition2_met = false; + bool condition3_met = false; + + for (int n = 0; n < NumSpec; ++n) { + if (aion[n] == 2 * zion[n] && !condition1_met) { + + // Set the first encountered nuclei with Z/A = 0.5 + massfractions(1, n+1) = 0.7_rt; + massfractions(2, n+1) = 1.0_rt; + massfractions(3, n+1) = 0.7_rt; + condition1_met = true; + continue; + } + if (aion[n] > 2 * zion[n] && !condition2_met) { + + // This is for achieving ye < 0.5 + massfractions(1, n+1) = 0.3_rt; + condition2_met = true; + continue; + } + if (aion[n] < 2 * zion[n] && !condition3_met) { + + // This is for achieving ye > 0.5 + massfractions(3, n+1) = 0.3_rt; + condition3_met = true; + continue; + } + massfractions(1, n+1) = 0.0_rt; + massfractions(2, n+1) = 0.0_rt; + massfractions(3, n+1) = 0.0_rt; + } + + eos_t eos_state; + + for (int iye = 0; iye < 3; ++iye) { + for (int irho = 0; irho < nrho; ++irho) { + for (int itemp = 0; itemp < nT; ++itemp) { + + Real T = std::pow(10.0, std::log10(T_min) + itemp * dlogT); + Real rho = std::pow(10.0, std::log10(rho_min) + irho * dlogrho); + + eos_state.T = T; + eos_state.rho = rho; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = massfractions(iye+1, n+1); + } + eos(eos_input_rt, eos_state); + burn_cell_c(eos_state); + } + } + } +} +#endif