From 16d8018be285ae5ef182d8feab0dec0183a2db94 Mon Sep 17 00:00:00 2001
From: Andy Nonaka <AJNonaka@lbl.gov>
Date: Mon, 19 Aug 2024 16:01:39 -0700
Subject: [PATCH] address a long-standing issue where sometimes its preferable
 to specify the avogadro number instead of Runiv.

If avogadro is supplied, then Runiv is set to avogadro * k_b

Default behavior is preserved as all codes presently set Runiv

The code will abort if you try to specify both.
---
 exec/chemistry_testing/main.cpp             | 10 +++++-----
 src_chemistry/chemistry_functions.cpp       |  6 +++---
 src_common/common_functions.cpp             |  8 ++++++++
 src_common/common_namespace.H               |  1 +
 src_compressible/compressible_functions.cpp |  4 ++--
 src_compressible/main_driver.cpp            |  2 +-
 src_compressible_stag/reservoirStag.cpp     |  4 +---
 7 files changed, 21 insertions(+), 14 deletions(-)

diff --git a/exec/chemistry_testing/main.cpp b/exec/chemistry_testing/main.cpp
index 2577074eb..a17d4f476 100644
--- a/exec/chemistry_testing/main.cpp
+++ b/exec/chemistry_testing/main.cpp
@@ -201,14 +201,14 @@ void main_main(const char* argv)
 
     for (int n=0; n<nspecies; n++)
     {
-        //amrex::Print() << ComputeSpatialMean(rho_old,n)*(Runiv/k_B)/molmass[n] << " ";
+        //amrex::Print() << ComputeSpatialMean(rho_old,n)*(avogadro)/molmass[n] << " ";
         amrex::Print() << ComputeSpatialMean(rho_old,n) << " ";
     }
 
 
     for (int n=0; n<nspecies; n++)
     {
-        //amrex::Print() << ComputeSpatialVariance(rho_old,n)*((Runiv/k_B)/molmass[n])*((Runiv/k_B)/molmass[n]) << " ";
+        //amrex::Print() << ComputeSpatialVariance(rho_old,n)*((avogadro)/molmass[n])*((avogadro)/molmass[n]) << " ";
         amrex::Print() << ComputeSpatialVariance(rho_old,n) << " ";
     }
 
@@ -235,7 +235,7 @@ void main_main(const char* argv)
                 {
                     GpuArray<amrex::Real,MAX_SPECIES> n_old;
                     GpuArray<amrex::Real,MAX_SPECIES> n_new;
-                    for (int n=0; n<nspecies; n++) n_old[n] = rhoOld(i,j,k,n)*(Runiv/k_B)/molmass[n];
+                    for (int n=0; n<nspecies; n++) n_old[n] = rhoOld(i,j,k,n)*(avogadro)/molmass[n];
                     switch(reaction_type){
                         case 0: // deterministic case
                             advance_reaction_det_cell(n_old,n_new,dt);
@@ -300,12 +300,12 @@ void main_main(const char* argv)
         amrex::Print()  << dt*step << " ";
         for (int n=0; n<nspecies; n++)
         {
-            //amrex::Print()  << ComputeSpatialMean(rho_new,n)*(Runiv/k_B)/molmass[n] << " ";
+            //amrex::Print()  << ComputeSpatialMean(rho_new,n)*(avogadro)/molmass[n] << " ";
             amrex::Print()  << ComputeSpatialMean(rho_new,n) << " ";
         }
         for (int n=0; n<nspecies; n++)
         {
-            //amrex::Print()  << ComputeSpatialVariance(rho_new,n)*((Runiv/k_B)/molmass[n])*((Runiv/k_B)/molmass[n]) << " ";
+            //amrex::Print()  << ComputeSpatialVariance(rho_new,n)*((avogadro)/molmass[n])*((avogadro)/molmass[n]) << " ";
             amrex::Print()  << ComputeSpatialVariance(rho_new,n) << " ";
         }
         amrex::Print() << "\n";
diff --git a/src_chemistry/chemistry_functions.cpp b/src_chemistry/chemistry_functions.cpp
index c9afc1cdf..8bdf0e485 100644
--- a/src_chemistry/chemistry_functions.cpp
+++ b/src_chemistry/chemistry_functions.cpp
@@ -105,7 +105,7 @@ void compute_chemistry_source_CLE(amrex::Real dt, amrex::Real dV,
             amrex::Real T = prim_arr(i,j,k,4);
             amrex::Real T0 = T0_chem;
             amrex::Real ctot = pres/Runiv/T;
-            amrex::Real Navo = Runiv/k_B;
+            amrex::Real Navo = avogadro;
 
             GpuArray<amrex::Real,MAX_SPECIES> ck;   // molar concentrations
             for (int n=0; n<nspecies; n++) ck[n] = prim_arr(i,j,k,6+nspecies+n)*ctot;
@@ -176,7 +176,7 @@ void compute_chemistry_source_CLE_1(amrex::Real dt, amrex::Real dV,
     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]/(Runiv/k_B);
+    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
 
     for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
     {
@@ -229,7 +229,7 @@ void compute_chemistry_source_CLE_2(amrex::Real dt, amrex::Real dV,
     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]/(Runiv/k_B);
+    for (int n=0; n<nspecies; n++) m_s[n] = molmass[n]/(avogadro);
 
     for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
     {
diff --git a/src_common/common_functions.cpp b/src_common/common_functions.cpp
index d51fb1107..76f45eba0 100644
--- a/src_common/common_functions.cpp
+++ b/src_common/common_functions.cpp
@@ -68,6 +68,7 @@ AMREX_GPU_MANAGED amrex::Real common::variance_coef_ener;
 AMREX_GPU_MANAGED amrex::Real common::k_B;
 AMREX_GPU_MANAGED amrex::Real common::h_bar;
 AMREX_GPU_MANAGED amrex::Real common::Runiv;
+AMREX_GPU_MANAGED amrex::Real common::avogadro;
 AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> common::T_init;
 AMREX_GPU_MANAGED int      common::algorithm_type;
 int                        common::barodiffusion_type;
@@ -764,6 +765,13 @@ void InitializeCommonNamespace() {
     pp.query("k_B",k_B);
     pp.query("h_bar",h_bar);
     pp.query("Runiv",Runiv);
+    avogadro = Runiv / k_B;
+    if (pp.query("avogadro",avogadro) ) {
+        Runiv = k_B * avogadro;
+    }
+    if (pp.query("Runiv",Runiv) && pp.query("avogadro",avogadro)) {
+        Abort("Cannot specify both Runiv and avogadro");
+    }
     if (pp.queryarr("T_init",temp)) {
         for (int i=0; i<nspecies; ++i) {
             T_init[i] = temp[i];
diff --git a/src_common/common_namespace.H b/src_common/common_namespace.H
index 12b865c78..5f8c68487 100644
--- a/src_common/common_namespace.H
+++ b/src_common/common_namespace.H
@@ -86,6 +86,7 @@ namespace common {
     extern AMREX_GPU_MANAGED amrex::Real k_B;
     extern AMREX_GPU_MANAGED amrex::Real h_bar;
     extern AMREX_GPU_MANAGED amrex::Real Runiv;
+    extern AMREX_GPU_MANAGED amrex::Real avogadro;
     extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, MAX_SPECIES> T_init;
 
     // Algorithm control / selection
diff --git a/src_compressible/compressible_functions.cpp b/src_compressible/compressible_functions.cpp
index 8ab3611a7..96bb6f502 100644
--- a/src_compressible/compressible_functions.cpp
+++ b/src_compressible/compressible_functions.cpp
@@ -229,13 +229,13 @@ void InitConsVar(MultiFab& cons,
                 Real pscale = 884.147e3;
 
                 // rho0 = m * p / (k_B * T)
-                Real rhoscale = molmass[0] / (Runiv / k_B) * pscale / (k_B * T_init[0]);
+                Real rhoscale = molmass[0] / avogadro * pscale / (k_B * T_init[0]);
 
                 // compute pressure (needed to compute density)
                 Real pres = pscale+rhoscale*velscale*velscale*cos(2.*pi*x/Lf)*cos(4.*pi*y/Lf)*(cos(4.*pi*z/Lf)+2.);
 
                 // density
-                cu(i,j,k,0) = (molmass[0] / (Runiv / k_B)) * pres / (k_B * T_init[0]);
+                cu(i,j,k,0) = (molmass[0] / avogadro) * pres / (k_B * T_init[0]);
                
                 // momentum
                 cu(i,j,k,1) =  velscale*cu(i,j,k,0)*sin(2.*pi*x/Lf)*cos(2.*pi*y/Lf)*cos(2.*pi*z/Lf);
diff --git a/src_compressible/main_driver.cpp b/src_compressible/main_driver.cpp
index f682dc08b..504ae956a 100644
--- a/src_compressible/main_driver.cpp
+++ b/src_compressible/main_driver.cpp
@@ -337,7 +337,7 @@ void main_driver(const char* argv)
     
         p0 = 884.147e3;
         dProb = (AMREX_SPACEDIM==2) ? 1./(n_cells[0]*n_cells[1]) : 1./(n_cells[0]*n_cells[1]*n_cells[2]);
-        rho0 = molmass[0] / (Runiv / k_B) * p0 / (k_B * T_init[0]);
+        rho0 = molmass[0] / avogadro * p0 / (k_B * T_init[0]);
         nu0 = 0.185;
     }
         
diff --git a/src_compressible_stag/reservoirStag.cpp b/src_compressible_stag/reservoirStag.cpp
index b59bf5627..54a0807c4 100644
--- a/src_compressible_stag/reservoirStag.cpp
+++ b/src_compressible_stag/reservoirStag.cpp
@@ -21,11 +21,9 @@ ComputeFluxMomReservoir(const MultiFab& cons0_in, const MultiFab& prim0_in,
     const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();
     Box dom(geom.Domain());
 
-    Real N_A = Runiv/k_B; // Avagadro's number
-    
     GpuArray<Real,MAX_SPECIES> mass;
     for (int l=0;l<nspecies;++l) {
-        mass[l] = molmass[l]/(N_A);
+        mass[l] = molmass[l]/avogadro;
     }
 
     for (int d=0; d<AMREX_SPACEDIM; ++d) {