C++ screening

networks/aprox13/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 24

networks/aprox19/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 33

networks/aprox21/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 35

networks/breakout/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/general_null/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/general_null/network_header.template

@@ -12,6 +12,8 @@
 constexpr int NumSpec = @@NSPEC@@;
 constexpr int NumAux = @@NAUX@@;
 
+@@PROPERTIES@@
+
 static AMREX_GPU_MANAGED amrex::Real aion[NumSpec] = {
    @@AION@@
   };

networks/general_null/write_network.py

@@ -145,7 +145,7 @@ def abort(outfile):
 
 
 def write_network(network_template, header_template,
-                  net_file,
+                  net_file, properties_file,
                   network_file, header_file):
     """read through the list of species and output the new out_file
 
@@ -164,6 +164,17 @@ def write_network(network_template, header_template,
         abort(out_file)
 
 
+    properties = {}
+    try:
+        with open(properties_file) as f:
+            for line in f:
+                key, value = line.split()
+                properties[key] = value
+    except FileNotFoundError:
+        print("no NETWORK_PROPERTIES found, skipping...")
+
+
+
     #-------------------------------------------------------------------------
     # write out the Fortran and C++ files based on the templates
     #-------------------------------------------------------------------------
@@ -265,6 +276,12 @@ def write_network(network_template, header_template,
                         for n, aux in enumerate(aux_vars):
                             fout.write("{}\"{}\",   // {} \n".format(indent, aux.name, n))
 
+                elif keyword == "PROPERTIES":
+                    if lang == "C++":
+                        for p in properties:
+                            print(p)
+                            fout.write("{}constexpr int {} = {};\n".format(indent, p, properties[p]))
+
             else:
                 fout.write(line)
 
@@ -285,14 +302,16 @@ def main():
                         help="C++ header output file name")
     parser.add_argument("-s", type=str, default="",
                         help="network file name")
+    parser.add_argument("--other_properties", type=str, default="",
+                        help="a NETWORK_PROPERTIES file with other network properties")
 
     args = parser.parse_args()
 
     if args.t == "" or args.o == "":
         sys.exit("write_probin.py: ERROR: invalid calling sequence")
 
     write_network(args.t, args.header_template,
-                  args.s,
+                  args.s, args.other_properties,
                   args.o, args.header_output)
 
 if __name__ == "__main__":

networks/ignition_chamulak/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/ignition_simple/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/iso7/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 10

networks/powerlaw/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/rprox/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/triple_alpha_plus_cago/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

networks/update_headers.py

@@ -27,6 +27,7 @@ def main():
                                 "general_null/network_header.template")
 
     net_file = os.path.join(micro_path, "networks", net, "{}.net".format(net))
+    properties_file = os.path.join(micro_path, "networks", net, "NETWORK_PROPERTIES")
 
     f90_name = os.path.join(args.odir, "network_properties.F90")
     cxx_name = os.path.join(args.odir, "network_properties.H")
@@ -37,7 +38,7 @@ def main():
         pass
 
     write_network.write_network(fortran_template, cxx_template,
-                                net_file,
+                                net_file, properties_file,
                                 f90_name, cxx_name)
 
 if __name__ == "__main__":

networks/xrb_simple/NETWORK_PROPERTIES

@@ -0,0 +1 @@
+NSCREEN 0

rates/aprox_rates.H

@@ -2201,4 +2201,4 @@ void ecapnuc(const Real etakep, const Real temp,
     }
 }
 
-#endif
+#endif

rates/aprox_rates_data.cpp

@@ -12,4 +12,5 @@ AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,4> rfd2;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,12> tfdm;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,12> tfd0;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,12> tfd1;
-AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,12> tfd2;
+AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real,12> tfd2;
+

screening/Make.package

@@ -1 +1,6 @@
 F90EXE_sources += screen.F90
+CEXE_headers += screen.H
+CEXE_headers += screen_data.H
+CEXE_sources += screening.cpp
+CEXE_sources += screen_data.cpp
+CEXE_headers += screening.H

screening/screen.F90

@@ -232,7 +232,7 @@ subroutine fill_plasma_state(state, temp, dens, y)
     ! Local variables
 
     real(rt)         :: abar, zbar, z2bar
-    real(rt)         :: ytot, rr, tempi, dtempi, deni
+    real(rt)         :: ytot, rr, tempi, dtempi
     real(rt)         :: pp, qq, dppdt, xni
 !    real(rt)         :: dppdd
 
@@ -246,7 +246,7 @@ subroutine fill_plasma_state(state, temp, dens, y)
     rr               = dens * ytot
     tempi            = ONE / temp
     dtempi           = -tempi * tempi
-    deni             = ONE / dens
+    !deni             = ONE / dens
 
     pp               = sqrt(rr*tempi*(z2bar + zbar))
     qq               = HALF/pp *(z2bar + zbar)
@@ -303,7 +303,7 @@ subroutine screen5(state,jscreen,scor,scordt,scordd)
 
 
     ! local variables
-    real(rt)         :: z1, a1, z2, a2
+    real(rt)         :: z1, z2
 
     real(rt)         :: bb,cc,dccdt, &
                         qq,dqqdt,rr,drrdt, &
@@ -325,9 +325,7 @@ subroutine screen5(state,jscreen,scor,scordt,scordd)
     ! Get the ion data based on the input index
 
     z1 = z1scr(jscreen)
-    a1 = a1scr(jscreen)
     z2 = z2scr(jscreen)
-    a2 = a2scr(jscreen)
 
     ! calculate individual screening factors
     bb       = z1 * z2

screening/screen.H

@@ -0,0 +1,99 @@
+#ifndef _screen_H_
+#define _screen_H_
+
+#include <AMReX.H>
+#include <AMReX_Array.H>
+#include <AMReX_Vector.H>
+#include <AMReX_REAL.H>
+#include <network_properties.H>
+
+#include <cmath>
+
+using namespace amrex;
+
+struct plasma_state_t {
+
+  Real qlam0z;
+  Real qlam0zdt;
+  //Real qlam0zdd;
+
+  Real taufac;
+  Real taufacdt;
+
+  Real aa;
+  Real daadt;
+  //Real daadd;
+};
+
+inline
+std::ostream& operator<< (std::ostream& o, plasma_state_t const& pstate)
+{
+  o << "qlam0z   = " << pstate.qlam0z << std::endl;
+  o << "qlam0zdt = " << pstate.qlam0zdt << std::endl;
+  o << "taufac   = " << pstate.taufac << std::endl;
+  o << "taufacdt = " << pstate.taufacdt << std::endl;
+  o << "aa       = " << pstate.aa << std::endl;
+  o << "daadt    = " << pstate.daadt << std::endl;
+
+  return o;
+}
+
+
+
+namespace scrn {
+  const Real fact = 1.25992104989487e0_rt;
+  const Real co2        = (1.0_rt/3.0_rt) * 4.248719e3_rt;
+  const Real gamefx     = 0.3e0_rt;
+  const Real gamefs     = 0.8e0_rt;
+  const Real h12_max    = 300.e0_rt;
+};
+
+
+AMREX_FORCE_INLINE
+void
+fill_plasma_state(plasma_state_t& state, const Real temp, const Real dens, Real* y) {
+
+  Real sum = 0.0_rt;
+  for (int n = 0; n < NumSpec; n++) {
+    sum += y[n];
+  }
+  Real abar = 1.0_rt / sum;
+  Real ytot = sum;
+
+  sum = 0.0_rt;
+  Real sum2 = 0.0_rt;
+  for (int n = 0; n < NumSpec; n++) {
+    sum += zion[n]*y[n];
+    sum2 += zion[n]*zion[n]*y[n];
+  }
+
+  Real zbar = sum * abar;
+  Real z2bar = sum2 * abar;
+
+  Real rr = dens * ytot;
+  Real tempi = 1.0_rt / temp;
+  Real dtempi = -tempi * tempi;
+  //Real deni = 1.0_rt / dens;
+
+  Real pp = std::sqrt(rr*tempi*(z2bar + zbar));
+  Real qq = 0.5_rt/pp *(z2bar + zbar);
+  Real dppdt = qq*rr*dtempi;
+  //Real dppdd = qq * ytot * tempi;
+
+  state.qlam0z = 1.88e8_rt * tempi * pp;
+  state.qlam0zdt = 1.88e8_rt * (dtempi*pp + tempi*dppdt);
+  //state.qlam0zdd = 1.88e8_rt * tempi * dppdd;
+
+  state.taufac = scrn::co2 * std::pow(tempi, 1.0_rt/3.0_rt);
+  state.taufacdt = -(1.0_rt/3.0_rt) * state.taufac * tempi;
+
+  qq = rr * zbar;
+  Real xni = std::pow(qq, 1.0_rt/3.0_rt);
+  //dxnidd = (1.0_rt/3.0_rt) * xni * deni;
+
+  state.aa = 2.27493e5_rt * tempi * xni;
+  state.daadt = 2.27493e5_rt * dtempi * xni;
+  //state.daadd = 2.27493e5_rt * tempi * dxnidd;
+}
+
+#endif

screening/screen_data.H

@@ -0,0 +1,47 @@
+#ifndef _screen_data_H_
+#define _screen_data_H_
+
+#include <AMReX.H>
+#include <AMReX_Array.H>
+#include <AMReX_Vector.H>
+#include <AMReX_REAL.H>
+#include <network_properties.H>
+
+class screen_factors_t {
+
+public:
+
+  amrex::Real z1 = -1;
+  amrex::Real z2 = -1;
+  amrex::Real a1 = -1;
+  amrex::Real a2 = -1;
+
+  // zs13    = (z1+z2)**(1./3.)
+  // zhat    = combination of z1 and z2 raised to the 5/3 power
+  // zhat2   = combination of z1 and z2 raised to the 5/12 power
+  // lzav    = log of effective charge
+  // aznut   = combination of a1,z1,a2,z2 raised to 1/3 power
+
+  amrex::Real zs13 = 0.0;
+  amrex::Real zs13inv = 0.0;
+  amrex::Real zhat = 0.0;
+  amrex::Real zhat2 = 0.0;
+  amrex::Real lzav = 0.0;
+  amrex::Real aznut = 0.0;
+
+  bool validate_nuclei(const amrex::Real z1_pass, const amrex::Real a1_pass,
+                       const amrex::Real z2_pass, const amrex::Real a2_pass) {
+    // a simple function for unit testing / debug runs to
+    // ensure that we are accessing the proper screening info
+
+    return (z1_pass == z1) &&
+           (z2_pass == z2) &&
+           (a1_pass == a1) &&
+           (a2_pass == a2);
+  }
+};
+
+extern AMREX_GPU_MANAGED amrex::GpuArray<screen_factors_t, NSCREEN> scn_facs;
+
+
+#endif

screening/screen_data.cpp

@@ -0,0 +1,5 @@
+#include <screen_data.H>
+
+using namespace amrex;
+
+AMREX_GPU_MANAGED amrex::GpuArray<screen_factors_t, NSCREEN> scn_facs;

screening/screening.H

@@ -0,0 +1,20 @@
+#include <screen.H>
+
+#ifndef _screening_H_
+#define _screening_H_
+
+void screening_init();
+
+void screening_finalize();
+
+void add_screening_factor(const int i,
+                          const amrex::Real z1, const amrex::Real a1,
+                          const amrex::Real z2, const amrex::Real a2);
+
+AMREX_GPU_HOST_DEVICE
+void screen5(const plasma_state_t state,
+             const int jscreen,
+             amrex::Real& scor, amrex::Real& scordt, amrex::Real& scordd);
+
+
+#endif