add a script to test network compatibility with NSE_NET (#1416)

This script is based on burn_cell. This script does constant temperature evolution and compares it to the NSE results. It prints out tables at the end that show the absolute and relative errors of mass fractions.

There are three cases for ye: <0.5, =0.5, >0.5
Both temperature and density are set at runtime.

Since we didn't set a reference cell size, we should let nse.nse_dx_independent=1, if we want to see whether the burn_state would enter NSE.

If nse.nse_dx_independent=1, then the script only prints out cases when burn_state failed to enter NSE.

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
+
+

nse_solver/nse_compatibility/Make.package

@@ -0,0 +1,2 @@
+CEXE_sources += main.cpp
+CEXE_headers += nse_network_compatibility.H

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

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

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

nse_solver/nse_compatibility/main.cpp

@@ -0,0 +1,34 @@
+#include <iostream>
+#include <cstring>
+#include <vector>
+
+#include <AMReX_ParmParse.H>
+#include <AMReX_MultiFab.H>
+using namespace amrex;
+
+#include <extern_parameters.H>
+#include <eos.H>
+#include <network.H>
+#include <nse_network_compatibility.H>
+#include <unit_test.H>
+
+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();
+}

nse_solver/nse_compatibility/nse_network_compatibility.H

@@ -0,0 +1,196 @@
+#ifndef NSE_NETWORK_COMPATIBILITY_H
+#define NSE_NETWORK_COMPATIBILITY_H
+
+#include <extern_parameters.H>
+#include <eos.H>
+#include <network.H>
+#include <burner.H>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <react_util.H>
+
+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<Real>(nrho-1);
+    Real dlogT = (std::log10(T_max) - std::log10(T_min))/static_cast<Real>(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<Real, 1, 3, 1, NumSpec+1, Order::C> 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