From d802798da1bf2812c27bfc0787ee7814fc812a66 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 21 Aug 2024 12:54:09 -0700
Subject: [PATCH 1/4] clean up src_chemistry by moving chemistry-testing(only)
 specific code into chemistry_testing

---
 exec/chemistry_testing/Make.package           |   3 +-
 .../chemistry_testing_functions.H             |  17 ++
 .../chemistry_testing_functions.cpp           | 254 ++++++++++++++++++
 exec/chemistry_testing/main.cpp               |   5 +-
 exec/chemistry_testing/myfunc.H               |   6 -
 src_chemistry/CLE_EM.cpp                      |  40 ---
 src_chemistry/CLE_RK3.cpp                     | 108 --------
 src_chemistry/Make.package                    |   5 -
 src_chemistry/chemistry_functions.H           |  12 -
 src_chemistry/chemistry_functions.cpp         | 105 --------
 10 files changed, 276 insertions(+), 279 deletions(-)
 create mode 100644 exec/chemistry_testing/chemistry_testing_functions.H
 create mode 100644 exec/chemistry_testing/chemistry_testing_functions.cpp
 delete mode 100644 exec/chemistry_testing/myfunc.H
 delete mode 100644 src_chemistry/CLE_EM.cpp
 delete mode 100644 src_chemistry/CLE_RK3.cpp

diff --git a/exec/chemistry_testing/Make.package b/exec/chemistry_testing/Make.package
index fa194a35d..6283b5088 100644
--- a/exec/chemistry_testing/Make.package
+++ b/exec/chemistry_testing/Make.package
@@ -1,2 +1,3 @@
 CEXE_sources += main.cpp
-CEXE_headers += myfunc.H
+CEXE_sources += chemistry_testing_functions.cpp
+CEXE_headers += chemistry_testing_functions.H
diff --git a/exec/chemistry_testing/chemistry_testing_functions.H b/exec/chemistry_testing/chemistry_testing_functions.H
new file mode 100644
index 000000000..9b8f29d23
--- /dev/null
+++ b/exec/chemistry_testing/chemistry_testing_functions.H
@@ -0,0 +1,17 @@
+#ifndef MYFUNC_H_
+#define MYFUNC_H_
+
+void main_main(const char* argv);
+
+void EMstep_chem_only(amrex::MultiFab& rho_old, amrex::MultiFab& rho_new,
+                      const amrex::Geometry geom, const amrex::Real dt);
+
+void RK3step_chem_only(amrex::MultiFab& rho_old, amrex::MultiFab& rho_new,
+                       const amrex::Geometry geom, const amrex::Real dt);
+
+void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV, MultiFab& mf_in, int startComp_in,
+                                    MultiFab& source, int startComp_out);
+
+void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV, MultiFab& mf_in, int startComp_in,
+                                    MultiFab& source, int startComp_out, MultiFab& ranchem);
+#endif
diff --git a/exec/chemistry_testing/chemistry_testing_functions.cpp b/exec/chemistry_testing/chemistry_testing_functions.cpp
new file mode 100644
index 000000000..d64d162c2
--- /dev/null
+++ b/exec/chemistry_testing/chemistry_testing_functions.cpp
@@ -0,0 +1,254 @@
+#include "common_functions.H"
+#include "chemistry_functions.H"
+#include "rng_functions.H"
+#include "chemistry_testing_functions.H"
+
+void EMstep_chem_only(MultiFab& rho_old, MultiFab& rho_new,
+                       const amrex::Geometry geom, const amrex::Real dt)
+{
+    if (reaction_type!=1) amrex::Abort("EMstep_chem_only assumes reaction_type=1");
+
+    const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
+
+    BoxArray ba = rho_old.boxArray();
+    DistributionMapping dm = rho_old.DistributionMap();
+
+    MultiFab source(ba,dm,nspecies,0);;
+    
+    MultiFab ranchem(ba,dm,nreaction,0);
+
+    // initialize white noise field
+    for (int m=0;m<nreaction;m++) {
+        MultiFabFillRandom(ranchem,m,1.,geom);
+    }
+
+    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rho_old,0,source,0,ranchem);
+
+    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        const Box& bx = mfi.tilebox();
+
+        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
+        const Array4<Real> & rho_new_fab = rho_new.array(mfi);
+        const Array4<Real> & source_fab = source.array(mfi);
+
+        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+        {
+            rho_new_fab(i,j,k,n) = rho_old_fab(i,j,k,n) + dt*source_fab(i,j,k,n);
+        });
+    }
+
+}
+
+
+void RK3step_chem_only(MultiFab& rho_old, MultiFab& rho_new,
+                       const amrex::Geometry geom, const amrex::Real dt)
+{
+    if (reaction_type!=1) amrex::Abort("RK3step_chem_only assumes reaction_type=1");
+
+    const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
+
+    BoxArray ba = rho_old.boxArray();
+    DistributionMapping dm = rho_old.DistributionMap();
+
+    MultiFab rhop(ba,dm,nspecies,0);
+    MultiFab rhop2(ba,dm,nspecies,0);
+    MultiFab rhop3(ba,dm,nspecies,0);;
+
+    MultiFab source(ba,dm,nspecies,0);;
+
+    MultiFab ranchem(ba,dm,nreaction,0);
+    MultiFab ranchem_A(ba,dm,nreaction,0);
+    MultiFab ranchem_B(ba,dm,nreaction,0);
+    
+    // weights for stochastic fluxes; swgt2 changes each stage
+    amrex::Real swgt1, swgt2;
+    swgt1 = 1.;
+
+    // initialize white noise fields
+    for (int m=0;m<nreaction;m++) {
+        MultiFabFillRandom(ranchem_A,m,1.,geom);
+        MultiFabFillRandom(ranchem_B,m,1.,geom);
+    }
+
+    // stage1
+    swgt2 = (2.*std::sqrt(2.)+std::sqrt(3.))/5.;
+
+    MultiFab::LinComb(ranchem,
+                swgt1, ranchem_A, 0,
+                swgt2, ranchem_B, 0,
+                0, nreaction, 0);
+
+    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rho_old,0,source,0,ranchem);
+
+    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        const Box& bx = mfi.tilebox();
+
+        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
+        const Array4<Real> & rhop_fab = rhop.array(mfi);
+        const Array4<Real> & source_fab = source.array(mfi);
+
+        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+        {
+            rhop_fab(i,j,k,n) = rho_old_fab(i,j,k,n) + dt*source_fab(i,j,k,n);
+        });
+    }
+
+    // stage2
+    swgt2 = (-4.*std::sqrt(2.)+3.*std::sqrt(3.))/5.;
+
+    MultiFab::LinComb(ranchem,
+                swgt1, ranchem_A, 0,
+                swgt2, ranchem_B, 0,
+                0, nreaction, 0);
+
+    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rhop,0,source,0,ranchem);
+
+    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        const Box& bx = mfi.tilebox();
+
+        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
+        const Array4<Real> & rhop_fab = rhop.array(mfi);
+        const Array4<Real> & rhop2_fab = rhop2.array(mfi);
+        const Array4<Real> & source_fab = source.array(mfi);
+
+        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+        {
+            rhop2_fab(i,j,k,n) = 0.25*( 3.*rho_old_fab(i,j,k,n) + rhop_fab(i,j,k,n) + dt*source_fab(i,j,k,n) );
+        });
+    }
+
+    // stage3
+    swgt2 = (std::sqrt(2.)-2.*std::sqrt(3.))/10.;
+
+    MultiFab::LinComb(ranchem,
+                swgt1, ranchem_A, 0,
+                swgt2, ranchem_B, 0,
+                0, nreaction, 0);
+
+    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rhop2,0,source,0,ranchem);
+
+    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+        const Box& bx = mfi.tilebox();
+
+        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
+        const Array4<Real> & rho_new_fab = rho_new.array(mfi);
+        const Array4<Real> & rhop2_fab = rhop2.array(mfi);
+        const Array4<Real> & source_fab = source.array(mfi);
+
+        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+        {
+            rho_new_fab(i,j,k,n) = (2./3.)*( 0.5*rho_old_fab(i,j,k,n) + rhop2_fab(i,j,k,n) + dt*source_fab(i,j,k,n) );
+        });
+    }
+}
+
+
+
+void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV,
+                              MultiFab& mf_in, int startComp_in,
+                              MultiFab& source, int startComp_out)
+// mf_in: input MultiFab containing mass densitities rho1, rho2, ..., rho_nspecies
+// startComp_in: position of rho1 in mf_in
+// source: output MultiFab containing source terms corresponding to rho1, rho2, ..., rho_nspecies
+// startComp_out: position of the first source term corresponding to rho1 in MultiFab source
+{
+    if (reaction_type<0 || reaction_type>2) amrex::Abort("ERROR: invalid reaction_type");
+    
+    if (reaction_type==2) amrex::Abort("ERROR: reaction_type=2 not implemented yet");
+    
+    GpuArray<amrex::Real,MAX_SPECIES> m_s;
+    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
+
+    for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.validbox();
+
+        const Array4<Real>& rho_arr = mf_in.array(mfi);
+        const Array4<Real>& source_arr = source.array(mfi);
+
+        amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, RandomEngine const& engine) noexcept
+        {
+            GpuArray<amrex::Real,MAX_SPECIES> n_dens;
+            for (int n=0; n<nspecies; n++) n_dens[n] = rho_arr(i,j,k,n+startComp_in)/m_s[n];
+            
+            GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
+            compute_reaction_rates(n_dens,avg_react_rate);
+
+            GpuArray<amrex::Real,MAX_SPECIES> sourceArr;
+            for (int n=0; n<nspecies; n++) sourceArr[n] = 0.;
+
+            for (int m=0; m<nreaction; m++)
+            {    
+                avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
+                for (int n=0; n<nspecies; n++)
+                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*avg_react_rate[m];
+            }
+            
+            if (reaction_type==1)
+            {
+                for (int m=0; m<nreaction; m++)
+                {
+                    amrex::Real W = RandomNormal(0.,1.,engine)/sqrt(dt*dV);
+                    for (int n=0; n<nspecies; n++)
+                        sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
+                }
+            }
+
+            for (int n=0; n<nspecies; n++) source_arr(i,j,k,n+startComp_out) = sourceArr[n];
+        });
+    }
+}
+
+void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV,
+                              MultiFab& mf_in, int startComp_in,
+                              MultiFab& source, int startComp_out, MultiFab& ranchem)
+// mf_in: input MultiFab containing mass densitities rho1, rho2, ..., rho_nspecies
+// startComp_in: position of rho1 in mf_in
+// source: output MultiFab containing source terms corresponding to rho1, rho2, ..., rho_nspecies
+// startComp_out: position of the first source term corresponding to rho1 in MultiFab source
+{
+    if (reaction_type!=1) amrex::Abort("ERROR: compute_chemistry_source_CLE assumes reaction_type=1");
+
+    GpuArray<amrex::Real,MAX_SPECIES> m_s;
+    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
+
+    for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.validbox();
+
+        const Array4<Real>& rho_arr = mf_in.array(mfi);
+        const Array4<Real>& source_arr = source.array(mfi);
+        const Array4<Real>& ranchem_arr = ranchem.array(mfi);
+
+        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            GpuArray<amrex::Real,MAX_SPECIES> n_dens;
+            for (int n=0; n<nspecies; n++) n_dens[n] = rho_arr(i,j,k,n+startComp_in)/m_s[n];
+
+            GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
+            compute_reaction_rates(n_dens,avg_react_rate);
+
+            GpuArray<amrex::Real,MAX_SPECIES> sourceArr;
+            for (int n=0; n<nspecies; n++) sourceArr[n] = 0.;
+
+            for (int m=0; m<nreaction; m++)
+            {
+                avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
+
+                amrex::Real W = ranchem_arr(i,j,k,m)/sqrt(dt*dV);
+
+                for (int n=0; n<nspecies; n++)
+                {
+                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*avg_react_rate[m];
+                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
+                }
+            }
+
+            for (int n=0; n<nspecies; n++) source_arr(i,j,k,n+startComp_out) = sourceArr[n];
+        });
+    }
+}
diff --git a/exec/chemistry_testing/main.cpp b/exec/chemistry_testing/main.cpp
index a17d4f476..c44c16290 100644
--- a/exec/chemistry_testing/main.cpp
+++ b/exec/chemistry_testing/main.cpp
@@ -2,10 +2,11 @@
 #include <AMReX_ParmParse.H>
 #include <AMReX_Vector.H>
 
-#include "myfunc.H"
 #include "chemistry_functions.H"
 #include "common_functions.H"
 
+#include "chemistry_testing_functions.H"
+
 using namespace amrex;
 
 int main (int argc, char* argv[])
@@ -267,7 +268,7 @@ void main_main(const char* argv)
 
                 const Array4<Real>& sourceArr = source.array(mfi);
 
-                amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, RandomEngine const& engine) noexcept
+                amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
                 {
                     for (int n=0; n<nspecies; n++) rhoNew(i,j,k,n) = rhoOld(i,j,k,n) + dt*sourceArr(i,j,k,n);
                 });
diff --git a/exec/chemistry_testing/myfunc.H b/exec/chemistry_testing/myfunc.H
deleted file mode 100644
index 42226b716..000000000
--- a/exec/chemistry_testing/myfunc.H
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef MYFUNC_H_
-#define MYFUNC_H_
-
-void main_main(const char* argv);
-
-#endif
diff --git a/src_chemistry/CLE_EM.cpp b/src_chemistry/CLE_EM.cpp
deleted file mode 100644
index 8f8122527..000000000
--- a/src_chemistry/CLE_EM.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#include "common_functions.H"
-#include "chemistry_functions.H"
-#include "rng_functions.H"
-
-void EMstep_chem_only(MultiFab& rho_old, MultiFab& rho_new,
-                       const amrex::Geometry geom, const amrex::Real dt)
-{
-    if (reaction_type!=1) amrex::Abort("EMstep_chem_only assumes reaction_type=1");
-
-    const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
-
-    BoxArray ba = rho_old.boxArray();
-    DistributionMapping dm = rho_old.DistributionMap();
-
-    MultiFab source(ba,dm,nspecies,0);;
-    
-    MultiFab ranchem(ba,dm,nreaction,0);
-
-    // initialize white noise field
-    for (int m=0;m<nreaction;m++) {
-        MultiFabFillRandom(ranchem,m,1.,geom);
-    }
-
-    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rho_old,0,source,0,ranchem);
-
-    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-
-        const Box& bx = mfi.tilebox();
-
-        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
-        const Array4<Real> & rho_new_fab = rho_new.array(mfi);
-        const Array4<Real> & source_fab = source.array(mfi);
-
-        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-        {
-            rho_new_fab(i,j,k,n) = rho_old_fab(i,j,k,n) + dt*source_fab(i,j,k,n);
-        });
-    }
-
-}
diff --git a/src_chemistry/CLE_RK3.cpp b/src_chemistry/CLE_RK3.cpp
deleted file mode 100644
index dc5783667..000000000
--- a/src_chemistry/CLE_RK3.cpp
+++ /dev/null
@@ -1,108 +0,0 @@
-#include "common_functions.H"
-#include "chemistry_functions.H"
-#include "rng_functions.H"
-
-void RK3step_chem_only(MultiFab& rho_old, MultiFab& rho_new,
-                       const amrex::Geometry geom, const amrex::Real dt)
-{
-    if (reaction_type!=1) amrex::Abort("RK3step_chem_only assumes reaction_type=1");
-
-    const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
-
-    BoxArray ba = rho_old.boxArray();
-    DistributionMapping dm = rho_old.DistributionMap();
-
-    MultiFab rhop(ba,dm,nspecies,0);
-    MultiFab rhop2(ba,dm,nspecies,0);
-    MultiFab rhop3(ba,dm,nspecies,0);;
-
-    MultiFab source(ba,dm,nspecies,0);;
-
-    MultiFab ranchem(ba,dm,nreaction,0);
-    MultiFab ranchem_A(ba,dm,nreaction,0);
-    MultiFab ranchem_B(ba,dm,nreaction,0);
-    
-    // weights for stochastic fluxes; swgt2 changes each stage
-    amrex::Real swgt1, swgt2;
-    swgt1 = 1.;
-
-    // initialize white noise fields
-    for (int m=0;m<nreaction;m++) {
-        MultiFabFillRandom(ranchem_A,m,1.,geom);
-        MultiFabFillRandom(ranchem_B,m,1.,geom);
-    }
-
-    // stage1
-    swgt2 = (2.*std::sqrt(2.)+std::sqrt(3.))/5.;
-
-    MultiFab::LinComb(ranchem,
-                swgt1, ranchem_A, 0,
-                swgt2, ranchem_B, 0,
-                0, nreaction, 0);
-
-    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rho_old,0,source,0,ranchem);
-
-    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-
-        const Box& bx = mfi.tilebox();
-
-        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
-        const Array4<Real> & rhop_fab = rhop.array(mfi);
-        const Array4<Real> & source_fab = source.array(mfi);
-
-        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-        {
-            rhop_fab(i,j,k,n) = rho_old_fab(i,j,k,n) + dt*source_fab(i,j,k,n);
-        });
-    }
-
-    // stage2
-    swgt2 = (-4.*std::sqrt(2.)+3.*std::sqrt(3.))/5.;
-
-    MultiFab::LinComb(ranchem,
-                swgt1, ranchem_A, 0,
-                swgt2, ranchem_B, 0,
-                0, nreaction, 0);
-
-    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rhop,0,source,0,ranchem);
-
-    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-
-        const Box& bx = mfi.tilebox();
-
-        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
-        const Array4<Real> & rhop_fab = rhop.array(mfi);
-        const Array4<Real> & rhop2_fab = rhop2.array(mfi);
-        const Array4<Real> & source_fab = source.array(mfi);
-
-        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-        {
-            rhop2_fab(i,j,k,n) = 0.25*( 3.*rho_old_fab(i,j,k,n) + rhop_fab(i,j,k,n) + dt*source_fab(i,j,k,n) );
-        });
-    }
-
-    // stage3
-    swgt2 = (std::sqrt(2.)-2.*std::sqrt(3.))/10.;
-
-    MultiFab::LinComb(ranchem,
-                swgt1, ranchem_A, 0,
-                swgt2, ranchem_B, 0,
-                0, nreaction, 0);
-
-    compute_chemistry_source_CLE_2(dt,dx[0]*dx[1]*dx[2],rhop2,0,source,0,ranchem);
-
-    for ( MFIter mfi(rho_old,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-
-        const Box& bx = mfi.tilebox();
-
-        const Array4<Real> & rho_old_fab = rho_old.array(mfi);
-        const Array4<Real> & rho_new_fab = rho_new.array(mfi);
-        const Array4<Real> & rhop2_fab = rhop2.array(mfi);
-        const Array4<Real> & source_fab = source.array(mfi);
-
-        amrex::ParallelFor(bx, nspecies, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-        {
-            rho_new_fab(i,j,k,n) = (2./3.)*( 0.5*rho_old_fab(i,j,k,n) + rhop2_fab(i,j,k,n) + dt*source_fab(i,j,k,n) );
-        });
-    }
-}
diff --git a/src_chemistry/Make.package b/src_chemistry/Make.package
index 899282e94..f6ffd66cd 100644
--- a/src_chemistry/Make.package
+++ b/src_chemistry/Make.package
@@ -1,8 +1,3 @@
 CEXE_sources   += chemistry_functions.cpp
 CEXE_headers   += chemistry_functions.H
-
 CEXE_headers   += chemistry_namespace.H
-
-CEXE_sources   += CLE_RK3.cpp
-
-CEXE_sources   += CLE_EM.cpp
diff --git a/src_chemistry/chemistry_functions.H b/src_chemistry/chemistry_functions.H
index dfee3589d..57b62bc17 100644
--- a/src_chemistry/chemistry_functions.H
+++ b/src_chemistry/chemistry_functions.H
@@ -22,18 +22,6 @@ void compute_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
 AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<amrex::Real,MAX_SPECIES>& n_dens,
                                                   GpuArray<amrex::Real,MAX_REACTION>& a_r);
 
-void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV, MultiFab& mf_in, int startComp_in,
-                                    MultiFab& source, int startComp_out);
-
-void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV, MultiFab& mf_in, int startComp_in,
-                                    MultiFab& source, int startComp_out, MultiFab& ranchem);
-
-void EMstep_chem_only(MultiFab& rho_old, MultiFab& rho_new,
-                      const amrex::Geometry geom, const amrex::Real dt);
-
-void RK3step_chem_only(MultiFab& rho_old, MultiFab& rho_new,
-                       const amrex::Geometry geom, const amrex::Real dt);
-
 AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
                                                      GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt);
 
diff --git a/src_chemistry/chemistry_functions.cpp b/src_chemistry/chemistry_functions.cpp
index 8bdf0e485..ea90e8faf 100644
--- a/src_chemistry/chemistry_functions.cpp
+++ b/src_chemistry/chemistry_functions.cpp
@@ -163,111 +163,6 @@ AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<Real,MAX_SPECIES>& n_
     return;
 }
 
-void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV,
-                              MultiFab& mf_in, int startComp_in,
-                              MultiFab& source, int startComp_out)
-// mf_in: input MultiFab containing mass densitities rho1, rho2, ..., rho_nspecies
-// startComp_in: position of rho1 in mf_in
-// source: output MultiFab containing source terms corresponding to rho1, rho2, ..., rho_nspecies
-// startComp_out: position of the first source term corresponding to rho1 in MultiFab source
-{
-    if (reaction_type<0 || reaction_type>2) amrex::Abort("ERROR: invalid reaction_type");
-    
-    if (reaction_type==2) amrex::Abort("ERROR: reaction_type=2 not implemented yet");
-    
-    GpuArray<amrex::Real,MAX_SPECIES> m_s;
-    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
-
-    for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
-    {
-        const Box& bx = mfi.validbox();
-
-        const Array4<Real>& rho_arr = mf_in.array(mfi);
-        const Array4<Real>& source_arr = source.array(mfi);
-
-        amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, RandomEngine const& engine) noexcept
-        {
-            GpuArray<amrex::Real,MAX_SPECIES> n_dens;
-            for (int n=0; n<nspecies; n++) n_dens[n] = rho_arr(i,j,k,n+startComp_in)/m_s[n];
-            
-            GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
-            compute_reaction_rates(n_dens,avg_react_rate);
-
-            GpuArray<amrex::Real,MAX_SPECIES> sourceArr;
-            for (int n=0; n<nspecies; n++) sourceArr[n] = 0.;
-
-            for (int m=0; m<nreaction; m++)
-            {    
-                avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
-                for (int n=0; n<nspecies; n++)
-                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*avg_react_rate[m];
-            }
-            
-            if (reaction_type==1)
-            {
-                for (int m=0; m<nreaction; m++)
-                {
-                    amrex::Real W = RandomNormal(0.,1.,engine)/sqrt(dt*dV);
-                    for (int n=0; n<nspecies; n++)
-                        sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
-                }
-            }
-
-            for (int n=0; n<nspecies; n++) source_arr(i,j,k,n+startComp_out) = sourceArr[n];
-        });
-    }
-}
-
-void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV,
-                              MultiFab& mf_in, int startComp_in,
-                              MultiFab& source, int startComp_out, MultiFab& ranchem)
-// mf_in: input MultiFab containing mass densitities rho1, rho2, ..., rho_nspecies
-// startComp_in: position of rho1 in mf_in
-// source: output MultiFab containing source terms corresponding to rho1, rho2, ..., rho_nspecies
-// startComp_out: position of the first source term corresponding to rho1 in MultiFab source
-{
-    if (reaction_type!=1) amrex::Abort("ERROR: compute_chemistry_source_CLE assumes reaction_type=1");
-
-    GpuArray<amrex::Real,MAX_SPECIES> m_s;
-    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
-
-    for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
-    {
-        const Box& bx = mfi.validbox();
-
-        const Array4<Real>& rho_arr = mf_in.array(mfi);
-        const Array4<Real>& source_arr = source.array(mfi);
-        const Array4<Real>& ranchem_arr = ranchem.array(mfi);
-
-        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-        {
-            GpuArray<amrex::Real,MAX_SPECIES> n_dens;
-            for (int n=0; n<nspecies; n++) n_dens[n] = rho_arr(i,j,k,n+startComp_in)/m_s[n];
-
-            GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
-            compute_reaction_rates(n_dens,avg_react_rate);
-
-            GpuArray<amrex::Real,MAX_SPECIES> sourceArr;
-            for (int n=0; n<nspecies; n++) sourceArr[n] = 0.;
-
-            for (int m=0; m<nreaction; m++)
-            {
-                avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
-
-                amrex::Real W = ranchem_arr(i,j,k,m)/sqrt(dt*dV);
-
-                for (int n=0; n<nspecies; n++)
-                {
-                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*avg_react_rate[m];
-                    sourceArr[n] += m_s[n]*stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
-                }
-            }
-
-            for (int n=0; n<nspecies; n++) source_arr(i,j,k,n+startComp_out) = sourceArr[n];
-        });
-    }
-}
-
 AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
                                                      GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt)
 {

From f2433ed6c74fbdd3be9070cf437523520ee34900 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 21 Aug 2024 14:19:14 -0700
Subject: [PATCH 2/4] make MAX_REAC a compile-time option in the GNUmakefile
 instead of hard-coded in chemistry_namespace

error checking to make sure MAX_REAC and MAX_SPEC are larger than nreaction and nspecies
---
 exec/chemistry_testing/GNUmakefile     | 6 +++++-
 exec/compressible/GNUmakefile          | 4 ++++
 exec/compressible_stag/GNUmakefile     | 4 ++++
 exec/compressible_stag_mui/GNUmakefile | 3 +++
 src_chemistry/chemistry_functions.cpp  | 3 +++
 src_chemistry/chemistry_namespace.H    | 2 --
 src_common/common_functions.cpp        | 3 +++
 7 files changed, 22 insertions(+), 3 deletions(-)

diff --git a/exec/chemistry_testing/GNUmakefile b/exec/chemistry_testing/GNUmakefile
index b95ef559b..13bc1b55d 100644
--- a/exec/chemistry_testing/GNUmakefile
+++ b/exec/chemistry_testing/GNUmakefile
@@ -10,7 +10,8 @@ USE_OMP      = FALSE
 USE_CUDA     = FALSE
 COMP         = gnu
 DIM          = 3
-MAX_SPEC      = 8
+MAX_SPEC     = 8
+MAX_REAC     = 5
 
 TINY_PROFILE = FALSE
 
@@ -40,3 +41,6 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.rules
 MAXSPECIES := $(strip $(MAX_SPEC))
 DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
 
+MAXREACTION := $(strip $(MAX_REAC))
+DEFINES += -DMAX_REACTION=$(MAXREACTION)
+
diff --git a/exec/compressible/GNUmakefile b/exec/compressible/GNUmakefile
index eddabafdd..9a7e498c5 100644
--- a/exec/compressible/GNUmakefile
+++ b/exec/compressible/GNUmakefile
@@ -10,6 +10,7 @@ COMP          = gnu
 DIM           = 3
 TINY_PROFILE  = FALSE
 MAX_SPEC      = 8
+MAX_REAC     = 5
 
 USE_PARTICLES = FALSE
 
@@ -60,3 +61,6 @@ endif
 MAXSPECIES := $(strip $(MAX_SPEC))
 DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
 
+MAXREACTION := $(strip $(MAX_REAC))
+DEFINES += -DMAX_REACTION=$(MAXREACTION)
+
diff --git a/exec/compressible_stag/GNUmakefile b/exec/compressible_stag/GNUmakefile
index e754e7097..f549f704a 100644
--- a/exec/compressible_stag/GNUmakefile
+++ b/exec/compressible_stag/GNUmakefile
@@ -11,6 +11,7 @@ COMP          = gnu
 DIM           = 3
 TINY_PROFILE  = FALSE
 MAX_SPEC      = 8
+MAX_REAC     = 5
 
 USE_PARTICLES = FALSE
 DO_TURB 	  = FALSE
@@ -74,3 +75,6 @@ endif
 
 MAXSPECIES := $(strip $(MAX_SPEC))
 DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
+
+MAXREACTION := $(strip $(MAX_REAC))
+DEFINES += -DMAX_REACTION=$(MAXREACTION)
diff --git a/exec/compressible_stag_mui/GNUmakefile b/exec/compressible_stag_mui/GNUmakefile
index 2b6ca0815..d75e86fd9 100644
--- a/exec/compressible_stag_mui/GNUmakefile
+++ b/exec/compressible_stag_mui/GNUmakefile
@@ -15,6 +15,7 @@ USE_MUI       = TRUE
 USE_SLURM     = FALSE
 USE_AMREX_MPMD = TRUE
 MAX_SPEC      = 8
+MAX_REAC     = 5
 
 USE_PARTICLES = FALSE
 
@@ -127,3 +128,5 @@ endif
 MAXSPECIES := $(strip $(MAX_SPEC))
 DEFINES += -DMAX_SPECIES=$(MAXSPECIES)
 
+MAXREACTION := $(strip $(MAX_REAC))
+DEFINES += -DMAX_REACTION=$(MAXREACTION)
diff --git a/src_chemistry/chemistry_functions.cpp b/src_chemistry/chemistry_functions.cpp
index ea90e8faf..c942c7daa 100644
--- a/src_chemistry/chemistry_functions.cpp
+++ b/src_chemistry/chemistry_functions.cpp
@@ -22,6 +22,9 @@ void InitializeChemistryNamespace()
     nreaction = 0;
     // get number of reactions
     pp.query("nreaction",nreaction);
+    if (nreaction > MAX_REACTION) {
+        Abort("nreaction > MAX_REACTION; recompile with a new MAX_REAC in the GNUmakefile");
+    }
 
     // if nreaction is set to zero or not defined in the inputs file, quit the routine
     if (nreaction==0) return;
diff --git a/src_chemistry/chemistry_namespace.H b/src_chemistry/chemistry_namespace.H
index 169d17b88..df42071c2 100644
--- a/src_chemistry/chemistry_namespace.H
+++ b/src_chemistry/chemistry_namespace.H
@@ -1,6 +1,4 @@
 namespace chemistry {
-#define MAX_REACTION 5
-
     extern AMREX_GPU_MANAGED int nreaction;
 
     extern AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> rate_const; 
diff --git a/src_common/common_functions.cpp b/src_common/common_functions.cpp
index 76f45eba0..2c2db9421 100644
--- a/src_common/common_functions.cpp
+++ b/src_common/common_functions.cpp
@@ -648,6 +648,9 @@ void InitializeCommonNamespace() {
     // pp.getarr and queryarr("string",inputs,start_indx,count); can be used for arrays
 
     pp.query("nspecies",nspecies);
+    if (nspecies > MAX_SPECIES) {
+        Abort("nspecies > MAX_SPECIES; recompile with a new MAX_SPEC in the GNUmakefile");
+    }
     pp.query("nbonds",nbonds);
     
     if (pp.queryarr("prob_lo",temp)) {

From bdbae1a8beba8084a88100a3ef8a9f62edb75c5e Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 21 Aug 2024 14:21:35 -0700
Subject: [PATCH 3/4] spacing

---
 exec/compressible/GNUmakefile          | 2 +-
 exec/compressible_stag/GNUmakefile     | 2 +-
 exec/compressible_stag_mui/GNUmakefile | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/exec/compressible/GNUmakefile b/exec/compressible/GNUmakefile
index 9a7e498c5..cd7ee153f 100644
--- a/exec/compressible/GNUmakefile
+++ b/exec/compressible/GNUmakefile
@@ -10,7 +10,7 @@ COMP          = gnu
 DIM           = 3
 TINY_PROFILE  = FALSE
 MAX_SPEC      = 8
-MAX_REAC     = 5
+MAX_REAC      = 5
 
 USE_PARTICLES = FALSE
 
diff --git a/exec/compressible_stag/GNUmakefile b/exec/compressible_stag/GNUmakefile
index f549f704a..24dfa9492 100644
--- a/exec/compressible_stag/GNUmakefile
+++ b/exec/compressible_stag/GNUmakefile
@@ -11,7 +11,7 @@ COMP          = gnu
 DIM           = 3
 TINY_PROFILE  = FALSE
 MAX_SPEC      = 8
-MAX_REAC     = 5
+MAX_REAC      = 5
 
 USE_PARTICLES = FALSE
 DO_TURB 	  = FALSE
diff --git a/exec/compressible_stag_mui/GNUmakefile b/exec/compressible_stag_mui/GNUmakefile
index d75e86fd9..a04b97c26 100644
--- a/exec/compressible_stag_mui/GNUmakefile
+++ b/exec/compressible_stag_mui/GNUmakefile
@@ -15,7 +15,7 @@ USE_MUI       = TRUE
 USE_SLURM     = FALSE
 USE_AMREX_MPMD = TRUE
 MAX_SPEC      = 8
-MAX_REAC     = 5
+MAX_REAC      = 5
 
 USE_PARTICLES = FALSE
 

From fca4acdcaaf86e7b3a66ffc4e68bab9a18e66f65 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Wed, 21 Aug 2024 15:41:08 -0700
Subject: [PATCH 4/4] prep code for rewrite of new chemistry routines for the
 low Mach and reactDiff codes

---
 .../chemistry_testing_functions.H             |  14 ++
 .../chemistry_testing_functions.cpp           | 109 ++++++++++++
 src_chemistry/chemistry_functions.H           |  19 +-
 src_chemistry/chemistry_functions.cpp         | 164 ++++--------------
 src_compressible/timeStep.cpp                 |   6 +-
 src_compressible_stag/timeStepStag.cpp        |   6 +-
 6 files changed, 165 insertions(+), 153 deletions(-)

diff --git a/exec/chemistry_testing/chemistry_testing_functions.H b/exec/chemistry_testing/chemistry_testing_functions.H
index 9b8f29d23..715d01c3d 100644
--- a/exec/chemistry_testing/chemistry_testing_functions.H
+++ b/exec/chemistry_testing/chemistry_testing_functions.H
@@ -14,4 +14,18 @@ void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV, MultiFab& mf
 
 void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV, MultiFab& mf_in, int startComp_in,
                                     MultiFab& source, int startComp_out, MultiFab& ranchem);
+
+AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<amrex::Real,MAX_SPECIES>& n_dens,
+                                                  GpuArray<amrex::Real,MAX_REACTION>& a_r);
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt);
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_CLE_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
+                                                     amrex::Real dt,amrex::Real dV,RandomEngine const& engine);
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_SSA_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
+                                                     amrex::Real dt,amrex::Real dV,RandomEngine const& engine);
 #endif
diff --git a/exec/chemistry_testing/chemistry_testing_functions.cpp b/exec/chemistry_testing/chemistry_testing_functions.cpp
index d64d162c2..a27c7bd6f 100644
--- a/exec/chemistry_testing/chemistry_testing_functions.cpp
+++ b/exec/chemistry_testing/chemistry_testing_functions.cpp
@@ -252,3 +252,112 @@ void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV,
         });
     }
 }
+
+AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<Real,MAX_SPECIES>& n_dens,
+                                                  GpuArray<Real,MAX_REACTION>& a_r)
+{
+    for (int m=0; m<nreaction; m++)
+    {
+        a_r[m] = rate_const[m];
+        
+        for (int n=0; n<nspecies; n++)
+            a_r[m] *= pow(n_dens[n],stoich_coeffs_R(m,n));
+    }
+
+    return;
+}
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt)
+{
+    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
+
+    GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
+    compute_reaction_rates(n_new,avg_react_rate);
+
+    for (int m=0; m<nreaction; m++)
+    {
+        avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
+        for (int n=0; n<nspecies; n++)
+            n_new[n] += dt*stoich_coeffs_PR(m,n)*avg_react_rate[m];
+    }
+
+    return;
+}
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_CLE_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
+                                                     amrex::Real dt, amrex::Real dV,
+                                                     RandomEngine const& engine)
+{
+    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
+
+    GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
+    compute_reaction_rates(n_new,avg_react_rate);
+
+    for (int m=0; m<nreaction; m++)
+    {
+        avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
+
+        amrex::Real W = sqrt(dt/dV)*RandomNormal(0.,1.,engine);
+
+        for (int n=0; n<nspecies; n++)
+        {
+            n_new[n] += dt*stoich_coeffs_PR(m,n)*avg_react_rate[m];
+            n_new[n] += stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
+        }
+    }
+
+    return;
+}
+
+AMREX_GPU_HOST_DEVICE void advance_reaction_SSA_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
+                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
+                                                     amrex::Real dt, amrex::Real dV,
+                                                     RandomEngine const& engine)
+{
+    amrex::Real t_local = 0.;
+
+    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
+
+    while(true)
+    {
+        GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
+        compute_reaction_rates(n_new,avg_react_rate);
+
+        amrex::Real rTotal = 0.;
+        for (int m=0; m<nreaction; m++)
+        {
+            // convert reation rates to propensities
+            avg_react_rate[m] = std::max(0.,avg_react_rate[m]*dV);
+            rTotal += avg_react_rate[m];
+        }
+
+        if (rTotal==0.) break;
+
+        amrex::Real u1 = amrex::Random(engine);
+        amrex::Real tau = -log(1-u1)/rTotal;
+        t_local += tau; // update t_local
+
+        if (t_local > dt) break;
+
+        amrex::Real u2 = amrex::Random(engine);
+        u2 *= rTotal;
+
+        // find which reaction has occured
+        int which_reaction=0;
+        amrex::Real rSum = 0.;
+        for (int m=0; m<nreaction; m++)
+        {
+            rSum = rSum + avg_react_rate[m];
+            which_reaction = m;
+            if (rSum >= u2) break;
+        }
+
+        // update number densities for the reaction that has occured
+        for (int n=0; n<nspecies; n++)
+            n_new[n] += stoich_coeffs_PR(which_reaction,n)/dV;
+    }
+
+    return;
+}
diff --git a/src_chemistry/chemistry_functions.H b/src_chemistry/chemistry_functions.H
index 57b62bc17..41a5703f7 100644
--- a/src_chemistry/chemistry_functions.H
+++ b/src_chemistry/chemistry_functions.H
@@ -16,20 +16,7 @@ using namespace common;
 
 void InitializeChemistryNamespace();
 
-void compute_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
-                                  MultiFab& prim, MultiFab& source, MultiFab& ranchem);
-
-AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<amrex::Real,MAX_SPECIES>& n_dens,
-                                                  GpuArray<amrex::Real,MAX_REACTION>& a_r);
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt);
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_CLE_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
-                                                     amrex::Real dt,amrex::Real dV,RandomEngine const& engine);
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_SSA_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
-                                                     amrex::Real dt,amrex::Real dV,RandomEngine const& engine);
+// used in compressible code only
+void compute_compressible_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
+                                               MultiFab& prim, MultiFab& source, MultiFab& ranchem);
 #endif
diff --git a/src_chemistry/chemistry_functions.cpp b/src_chemistry/chemistry_functions.cpp
index c942c7daa..1bcd46524 100644
--- a/src_chemistry/chemistry_functions.cpp
+++ b/src_chemistry/chemistry_functions.cpp
@@ -4,16 +4,24 @@
 AMREX_GPU_MANAGED int chemistry::nreaction;
 
 AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> chemistry::rate_const;
-AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> chemistry::alpha_param;
-AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> chemistry::beta_param;
-AMREX_GPU_MANAGED amrex::Real chemistry::T0_chem;
 
-AMREX_GPU_MANAGED Array2D<int,0, MAX_REACTION,0, MAX_SPECIES> chemistry::stoich_coeffs_R; 
+// from the fortran code, stoich_coeffs_R = stoichiometric_factors(spec,1,reac)
+// from the fortran code, stoich_coeffs_P = stoichiometric_factors(spec,2,reac)
+// stoich_coeffs_PR = stoich_coeffs_P - stoich_coeffs_R
+AMREX_GPU_MANAGED Array2D<int,0, MAX_REACTION,0, MAX_SPECIES> chemistry::stoich_coeffs_R;
 AMREX_GPU_MANAGED Array2D<int,0, MAX_REACTION,0, MAX_SPECIES> chemistry::stoich_coeffs_P;
 AMREX_GPU_MANAGED Array2D<int,0, MAX_REACTION,0, MAX_SPECIES> chemistry::stoich_coeffs_PR;
 
+// from the fortran code this was use_Poisson_rng (0=CLE; 1=tau leaping; -1=deterministic; 2=SSA)
+// here it's being used as reaction_type (0=deterministic; 1=CLE; 2=SSA; 3=tau leap)
 AMREX_GPU_MANAGED int chemistry::reaction_type;
 
+// specific to compressible codes
+AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> chemistry::alpha_param;
+AMREX_GPU_MANAGED GpuArray<amrex::Real, MAX_REACTION> chemistry::beta_param;
+AMREX_GPU_MANAGED amrex::Real chemistry::T0_chem;
+
+
 void InitializeChemistryNamespace()
 {
     // extract inputs parameters
@@ -34,20 +42,6 @@ void InitializeChemistryNamespace()
     pp.getarr("rate_const",k_tmp,0,nreaction);
     for (int m=0; m<nreaction; m++) rate_const[m] = k_tmp[m];
 
-    // get alpha parameter
-    std::vector<amrex::Real> alpha_tmp(MAX_REACTION);
-    pp.queryarr("alpha_param",alpha_tmp,0,nreaction);
-    for (int m=0; m<nreaction; m++) alpha_param[m] = alpha_tmp[m];
-    
-    // get beta parameter
-    std::vector<amrex::Real> beta_tmp(MAX_REACTION);
-    pp.queryarr("beta_param",beta_tmp,0,nreaction);
-    for (int m=0; m<nreaction; m++) beta_param[m] = beta_tmp[m];
-    
-    T0_chem = 0.;
-    // get temperature T0 for rate constants
-    pp.query("T0_chem",T0_chem);
-
     // get stoich coeffs for reactants
     for (int m=0; m<nreaction; m++)
     {
@@ -81,16 +75,33 @@ void InitializeChemistryNamespace()
         for (int n=0; n<nspecies; n++)
             stoich_coeffs_PR(m,n) = stoich_coeffs_P(m,n)-stoich_coeffs_R(m,n);
 
-    // get reaction type: Deterministic, CLE or SSA
+    // get reaction type (0=deterministic; 1=CLE; 2=SSA; 3=tau leap)
     pp.get("reaction_type",reaction_type);
+
+    // get alpha parameter for compressible code
+    std::vector<amrex::Real> alpha_tmp(MAX_REACTION);
+    pp.queryarr("alpha_param",alpha_tmp,0,nreaction);
+    for (int m=0; m<nreaction; m++) alpha_param[m] = alpha_tmp[m];
+    
+    // get beta parameter for compressible code
+    std::vector<amrex::Real> beta_tmp(MAX_REACTION);
+    pp.queryarr("beta_param",beta_tmp,0,nreaction);
+    for (int m=0; m<nreaction; m++) beta_param[m] = beta_tmp[m];
     
+    T0_chem = 0.;
+    // get temperature T0 for rate constants for compressible code
+    pp.query("T0_chem",T0_chem);
+
     return;
 }
 
-void compute_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
-                                  MultiFab& prim, MultiFab& source, MultiFab& ranchem)
+// used in compressible code only
+void compute_compressible_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
+                                               MultiFab& prim, MultiFab& source, MultiFab& ranchem)
 {
-    if (reaction_type!=0 && reaction_type!=1) amrex::Abort("ERROR: compute_chemistry_source_CLE only works for reaction_type = 0 or 1");
+    if (reaction_type!=0 && reaction_type!=1) {
+        amrex::Abort("ERROR: compute_compressible_chemistry_source_CLE only works for reaction_type = 0 or 1");
+    }
 
     if (T0_chem<=0.) amrex::Abort("ERROR: T0_chem>0 expected");
 
@@ -151,112 +162,3 @@ void compute_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
         });
     }
 }
-
-AMREX_GPU_HOST_DEVICE void compute_reaction_rates(GpuArray<Real,MAX_SPECIES>& n_dens,
-                                                  GpuArray<Real,MAX_REACTION>& a_r)
-{
-    for (int m=0; m<nreaction; m++)
-    {
-        a_r[m] = rate_const[m];
-        
-        for (int n=0; n<nspecies; n++)
-            a_r[m] *= pow(n_dens[n],stoich_coeffs_R(m,n));
-    }
-
-    return;
-}
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_det_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,amrex::Real dt)
-{
-    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
-
-    GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
-    compute_reaction_rates(n_new,avg_react_rate);
-
-    for (int m=0; m<nreaction; m++)
-    {
-        avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
-        for (int n=0; n<nspecies; n++)
-            n_new[n] += dt*stoich_coeffs_PR(m,n)*avg_react_rate[m];
-    }
-
-    return;
-}
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_CLE_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
-                                                     amrex::Real dt, amrex::Real dV,
-                                                     RandomEngine const& engine)
-{
-    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
-
-    GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
-    compute_reaction_rates(n_new,avg_react_rate);
-
-    for (int m=0; m<nreaction; m++)
-    {
-        avg_react_rate[m] = std::max(0.,avg_react_rate[m]);
-
-        amrex::Real W = sqrt(dt/dV)*RandomNormal(0.,1.,engine);
-
-        for (int n=0; n<nspecies; n++)
-        {
-            n_new[n] += dt*stoich_coeffs_PR(m,n)*avg_react_rate[m];
-            n_new[n] += stoich_coeffs_PR(m,n)*sqrt(avg_react_rate[m])*W;
-        }
-    }
-
-    return;
-}
-
-AMREX_GPU_HOST_DEVICE void advance_reaction_SSA_cell(GpuArray<amrex::Real,MAX_SPECIES>& n_old,
-                                                     GpuArray<amrex::Real,MAX_SPECIES>& n_new,
-                                                     amrex::Real dt, amrex::Real dV,
-                                                     RandomEngine const& engine)
-{
-    amrex::Real t_local = 0.;
-
-    for (int n=0; n<nspecies; n++) n_new[n] = n_old[n];
-
-    while(true)
-    {
-        GpuArray<amrex::Real,MAX_REACTION> avg_react_rate;
-        compute_reaction_rates(n_new,avg_react_rate);
-
-        amrex::Real rTotal = 0.;
-        for (int m=0; m<nreaction; m++)
-        {
-            // convert reation rates to propensities
-            avg_react_rate[m] = std::max(0.,avg_react_rate[m]*dV);
-            rTotal += avg_react_rate[m];
-        }
-
-        if (rTotal==0.) break;
-
-        amrex::Real u1 = amrex::Random(engine);
-        amrex::Real tau = -log(1-u1)/rTotal;
-        t_local += tau; // update t_local
-
-        if (t_local > dt) break;
-
-        amrex::Real u2 = amrex::Random(engine);
-        u2 *= rTotal;
-
-        // find which reaction has occured
-        int which_reaction=0;
-        amrex::Real rSum = 0.;
-        for (int m=0; m<nreaction; m++)
-        {
-            rSum = rSum + avg_react_rate[m];
-            which_reaction = m;
-            if (rSum >= u2) break;
-        }
-
-        // update number densities for the reaction that has occured
-        for (int n=0; n<nspecies; n++)
-            n_new[n] += stoich_coeffs_PR(which_reaction,n)/dV;
-    }
-
-    return;
-}
diff --git a/src_compressible/timeStep.cpp b/src_compressible/timeStep.cpp
index 39f1b8ec9..1aa25fa7f 100644
--- a/src_compressible/timeStep.cpp
+++ b/src_compressible/timeStep.cpp
@@ -139,7 +139,7 @@ void RK3step(MultiFab& cu, MultiFab& cup, MultiFab& cup2, MultiFab& /*cup3*/,
                     stoch_weights[1], ranchem_B, 0,
                     0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
+        compute_compressible_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
     }
 
     for ( MFIter mfi(cu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
@@ -247,7 +247,7 @@ void RK3step(MultiFab& cu, MultiFab& cup, MultiFab& cup2, MultiFab& /*cup3*/,
                     stoch_weights[1], ranchem_B, 0,
                     0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
+        compute_compressible_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
     }
 
     for ( MFIter mfi(cu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
@@ -351,7 +351,7 @@ void RK3step(MultiFab& cu, MultiFab& cup, MultiFab& cup2, MultiFab& /*cup3*/,
                     stoch_weights[1], ranchem_B, 0,
                     0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
+        compute_compressible_chemistry_source_CLE(dt,dx[0]*dx[1]*dx[2],prim,source,ranchem);
     }
 
     for ( MFIter mfi(cu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
diff --git a/src_compressible_stag/timeStepStag.cpp b/src_compressible_stag/timeStepStag.cpp
index caeeefabd..335cf4acb 100644
--- a/src_compressible_stag/timeStepStag.cpp
+++ b/src_compressible_stag/timeStepStag.cpp
@@ -357,7 +357,7 @@ void RK3stepStag(MultiFab& cu,
             stoch_weights[1], ranchem_B, 0,
             0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
+        compute_compressible_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
     }
 
 #if defined(TURB)
@@ -684,7 +684,7 @@ void RK3stepStag(MultiFab& cu,
             stoch_weights[1], ranchem_B, 0,
             0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
+        compute_compressible_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
     }
     
 #if defined(TURB)
@@ -1018,7 +1018,7 @@ void RK3stepStag(MultiFab& cu,
             stoch_weights[1], ranchem_B, 0,
             0, nreaction, 0);
 
-        compute_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
+        compute_compressible_chemistry_source_CLE(dt, dx[0]*dx[1]*dx[2], prim, source, ranchem);
     }
     
 #if defined(TURB)