From 819a67913ecb6b658444040117d1ed67ba859772 Mon Sep 17 00:00:00 2001
From: Zhi Chen <62574124+zhichen3@users.noreply.github.com>
Date: Tue, 28 Nov 2023 07:24:57 -0500
Subject: [PATCH] Enforce chabrier1998 or null screening with NSE_NET (#1398)

---
 Make.Microphysics_extern          | 7 +++++++
 nse_solver/make_table/GNUmakefile | 2 ++
 nse_solver/nse_solver.H           | 8 ++++----
 unit_test/test_nse/GNUmakefile    | 2 ++
 4 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/Make.Microphysics_extern b/Make.Microphysics_extern
index 787482f9f5..063e029934 100644
--- a/Make.Microphysics_extern
+++ b/Make.Microphysics_extern
@@ -145,6 +145,13 @@ ifeq ($(USE_NSE_TABLE), TRUE)
   endif
 endif
 
+# Only allow null chabrier1998 screening with NSE_NET
+ifeq ($(USE_NSE_NET),TRUE)
+  ifneq ($(or $(filter-out null chabrier1998,$(SCREEN_METHOD)),),)
+    $(error USE_NSE_NET can only be used with null or chabrier1998 screening)
+  endif
+endif
+
 clean::
 	@if [ -L helm_table.dat ]; then rm -f helm_table.dat; fi
 	@if [ -L reaclib_rate_metadata.dat ]; then rm -f reaclib_rate_metadata.dat; fi
diff --git a/nse_solver/make_table/GNUmakefile b/nse_solver/make_table/GNUmakefile
index 0b48798f98..71b95bfc2a 100644
--- a/nse_solver/make_table/GNUmakefile
+++ b/nse_solver/make_table/GNUmakefile
@@ -24,6 +24,8 @@ EOS_DIR     := helmholtz
 # This sets the network directory
 NETWORK_DIR := aprox21
 
+SCREEN_METHOD := chabrier1998
+
 CONDUCTIVITY_DIR := stellar
 
 INTEGRATOR_DIR =  VODE
diff --git a/nse_solver/nse_solver.H b/nse_solver/nse_solver.H
index a7a4ea9ff4..b357f2c16b 100644
--- a/nse_solver/nse_solver.H
+++ b/nse_solver/nse_solver.H
@@ -89,10 +89,10 @@ void apply_nse_exponent(T& nse_state) {
   Real T_in = nse_state.T_fixed > 0.0_rt ? nse_state.T_fixed : nse_state.T;
   Real exponent;
 
-  // if we use a non-null screening method.
+  // if we use chabrier1998 screening
   // Get the required terms to calculate coulomb correction term, u_c
 
-#if SCREEN_METHOD != 0
+#if SCREEN_METHOD == 4
 
   // Find n_e for original state;
   const amrex::Real n_e = nse_state.rho * nse_state.y_e / C::m_u;
@@ -112,8 +112,8 @@ void apply_nse_exponent(T& nse_state) {
 #endif
       // term for calculating u_c
 
-      // if use non-null screening, calculate the coulomb correction term
-#if SCREEN_METHOD != 0
+      // if use chabrier1998 screening, calculate the coulomb correction term
+#if SCREEN_METHOD == 4
       gamma = std::pow(zion[n], 5.0_rt/3.0_rt) * Gamma_e;
 
       // chemical potential for coulomb correction
diff --git a/unit_test/test_nse/GNUmakefile b/unit_test/test_nse/GNUmakefile
index 42590e8a9f..34e21f6078 100644
--- a/unit_test/test_nse/GNUmakefile
+++ b/unit_test/test_nse/GNUmakefile
@@ -26,6 +26,8 @@ EOS_DIR     := helmholtz
 # This sets the network directory
 NETWORK_DIR := aprox21
 
+SCREEN_METHOD := chabrier1998
+
 CONDUCTIVITY_DIR := stellar
 
 INTEGRATOR_DIR =  VODE