unify the conductivity unit tests (C++ and Fortran) into a single test

unit_test/test_conductivity/GNUmakefile

@@ -10,12 +10,15 @@ COMP	   = gnu
 USE_MPI    = FALSE
 USE_OMP    = FALSE
 
-USE_REACT = TRUE
-USE_EXTRA_THERMO = TRUE
+USE_REACT = FALSE
 USE_CONDUCTIVITY = TRUE
 
 EBASE = main
 
+USE_EXTRA_THERMO = TRUE
+
+USE_CXX_EOS = TRUE
+
 # define the location of the CASTRO top directory
 MICROPHYSICS_HOME  := ../..
 

unit_test/test_conductivity/Make.package

@@ -1,8 +1,12 @@
 CEXE_sources += main.cpp
 
-FEXE_headers += test_conductivity_F.H
-CEXE_headers += test_conductivity.H
+FEXE_headers += test_cond_F.H
+CEXE_headers += test_cond.H
+
+CEXE_sources += variables.cpp
+CEXE_sources += conductivity_util.cpp
+CEXE_headers += variables.H
 
-F90EXE_sources += variables.F90
-F90EXE_sources += conductivity_util.F90
 f90EXE_sources += unit_test.f90
+F90EXE_sources += variables_F.F90
+F90EXE_sources += conductivity_util_F.F90

unit_test/test_conductivity/_parameters

@@ -5,7 +5,5 @@ temp_max      real       1.d12
 
 metalicity_max  real     0.1d0
 
-test_set      character  "gr0_3d"
-
 small_temp    real        1.e4
 small_dens    real        1.e-4

unit_test/test_conductivity/conductivity_util.cpp

@@ -0,0 +1,67 @@
+#include <AMReX_PlotFileUtil.H>
+#include <AMReX_ParmParse.H>
+#include <AMReX_Print.H>
+
+#include <AMReX_Geometry.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_BCRec.H>
+
+#include <variables.H>
+#include <network.H>
+#include <eos.H>
+#include <conductivity.H>
+
+#include <cmath>
+
+using namespace amrex;
+
+void cond_test_C(const Box& bx,
+                 const Real dlogrho, const Real dlogT, const Real dmetal,
+                 const plot_t vars,
+                 Array4<Real> const sp) {
+
+  const int ih1 = network_spec_index("hydrogen-1");
+  const int ihe4 = network_spec_index("helium-4");
+
+  AMREX_PARALLEL_FOR_3D(bx, i, j, k,
+  {
+
+    // set the composition -- approximately solar
+    Real metalicity = 0.0 + static_cast<Real> (k) * dmetal;
+
+    eos_t eos_state;
+
+    for (int n = 0; n < NumSpec; n++) {
+      eos_state.xn[n] = metalicity/(NumSpec - 2);
+    }
+    eos_state.xn[ih1] = 0.75 - 0.5*metalicity;
+    eos_state.xn[ihe4] = 0.25 - 0.5*metalicity;
+
+    Real temp_zone = std::pow(10.0, std::log10(temp_min) + static_cast<Real>(j)*dlogT);
+    eos_state.T = temp_zone;
+
+    Real dens_zone = std::pow(10.0, std::log10(dens_min) + static_cast<Real>(i)*dlogrho);
+    eos_state.rho = dens_zone;
+
+    // store default state
+    sp(i, j, k, vars.irho) = eos_state.rho;
+    sp(i, j, k, vars.itemp) = eos_state.T;
+    for (int n = 0; n < NumSpec; n++) {
+      sp(i, j, k, vars.ispec+n) = eos_state.xn[n];
+    }
+
+    // call the EOS using rho, T
+    eos(eos_input_rt, eos_state);
+
+    conductivity(eos_state);
+
+    sp(i, j, k, vars.ih) = eos_state.h;
+    sp(i, j, k, vars.ie) = eos_state.e;
+    sp(i, j, k, vars.ip) = eos_state.p;
+    sp(i, j, k, vars.is) = eos_state.s;
+
+    sp(i, j, k, vars.iconductivity) = eos_state.conductivity;
+
+  });
+
+}

unit_test/test_conductivity/conductivity_util.F90 → unit_test/test_conductivity/conductivity_util_F.F90

@@ -1,5 +1,6 @@
 subroutine do_conductivity(lo, hi, &
-                           sp, s_lo, s_hi, npts) bind(C, name="do_conductivity")
+                           dlogrho, dlogT, dmetal, &
+                           sp, s_lo, s_hi) bind(C, name="do_conductivity")
 
   use variables
   use network
@@ -15,9 +16,8 @@ subroutine do_conductivity(lo, hi, &
   integer, intent(in) :: lo(3), hi(3)
   integer, intent(in) :: s_lo(3), s_hi(3)
   real(rt), intent(inout) :: sp(s_lo(1):s_hi(1), s_lo(2):s_hi(2), s_lo(3):s_hi(3), p % n_plot_comps)
-  integer, intent(in), value :: npts
+  real(rt), intent(in), value :: dlogrho, dlogT, dmetal
 
-  real(rt) :: dlogrho, dlogT, dmetal
   real(rt) :: metalicity, temp_zone, dens_zone
   real(rt) :: xn_zone(nspec)
 
@@ -27,10 +27,6 @@ subroutine do_conductivity(lo, hi, &
   integer :: ii, jj, kk
 
   !$gpu
-
-  dlogrho = (log10(dens_max) - log10(dens_min))/(npts - 1)
-  dlogT   = (log10(temp_max) - log10(temp_min))/(npts - 1)
-  dmetal    = (metalicity_max  - ZERO)/(npts - 1)
 
   do kk = lo(3), hi(3)
 

unit_test/test_conductivity/main.cpp

@@ -9,10 +9,16 @@
 
 using namespace amrex;
 
-#include "test_conductivity.H"
-#include "test_conductivity_F.H"
+#include "test_cond.H"
+#include "test_cond_F.H"
 #include "AMReX_buildInfo.H"
 
+#include <network.H>
+#include <eos.H>
+#include <conductivity.H>
+#include <variables.H>
+
+#include <cmath>
 
 int main (int argc, char* argv[])
 {
@@ -28,7 +34,7 @@ void main_main ()
 {
 
     // AMREX_SPACEDIM: number of dimensions
-    int n_cell, max_grid_size;
+    int n_cell, max_grid_size, do_cxx;
     Vector<int> bc_lo(AMREX_SPACEDIM,0);
     Vector<int> bc_hi(AMREX_SPACEDIM,0);
 
@@ -45,6 +51,10 @@ void main_main ()
         max_grid_size = 32;
         pp.query("max_grid_size", max_grid_size);
 
+        // do_cxx = 1 for C++ EOS, 0 for Fortran EOS
+        do_cxx = 0;
+        pp.query("do_cxx", do_cxx);
+
     }
 
     Vector<int> is_periodic(AMREX_SPACEDIM,0);
@@ -98,24 +108,45 @@ void main_main ()
 
     init_unit_test(probin_file_name.dataPtr(), &probin_file_length);
 
-    // Ncomp = number of components for each array
+    init_extern_parameters();
+
+    eos_init();
+    conductivity_init();
+
     int Ncomp = -1;
-    init_variables();
-    get_ncomp(&Ncomp);
 
-    int name_len = -1;
-    get_name_len(&name_len);
+    // for C++
+    plot_t vars;
 
-    // get the variable names
+    // for F90
     Vector<std::string> varnames;
 
-    for (int i=0; i<Ncomp; i++) {
-      char* cstring[name_len+1];
-      get_var_name(cstring, &i);
-      std::string name(*cstring);
-      varnames.push_back(name);
+    if (do_cxx == 1) {
+      // C++ test
+      vars = init_variables();
+      Ncomp = vars.n_plot_comps;
+
+    } else {
+      // Fortran test
+
+      // Ncomp = number of components for each array
+      init_variables_F();
+      get_ncomp(&Ncomp);
+
+      int name_len = -1;
+      get_name_len(&name_len);
+
+      // get the variable names
+      for (int i=0; i<Ncomp; i++) {
+        char* cstring[name_len+1];
+        get_var_name(cstring, &i);
+        std::string name(*cstring);
+        varnames.push_back(name);
+      }
     }
 
+    std::cout << "Ncomp = " << Ncomp << std::endl;
+
     // time = starting time in the simulation
     Real time = 0.0;
 
@@ -125,19 +156,42 @@ void main_main ()
     // we allocate our main multifabs
     MultiFab state(ba, dm, Ncomp, Nghost);
 
+    // Initialize the state to zero; we will fill
+    // it in below in do_eos.
+    state.setVal(0.0);
+
     // What time is it now?  We'll use this to compute total run time.
     Real strt_time = ParallelDescriptor::second();
 
+
+    Real dlogrho = 0.0e0_rt;
+    Real dlogT   = 0.0e0_rt;
+    Real dmetal  = 0.0e0_rt;
+
+    if (n_cell > 1) {
+        dlogrho = (std::log10(dens_max) - std::log10(dens_min))/(n_cell - 1);
+        dlogT   = (std::log10(temp_max) - std::log10(temp_min))/(n_cell - 1);
+        dmetal  = (metalicity_max  - 0.0)/(n_cell - 1);
+    }
+
     // Initialize the state and compute the different thermodynamics
     // by inverting the EOS
     for ( MFIter mfi(state); mfi.isValid(); ++mfi )
     {
-        const Box& bx = mfi.validbox();
+      const Box& bx = mfi.validbox();
+
+      Array4<Real> const sp = state.array(mfi);
+
+      if (do_cxx == 1) {
+        cond_test_C(bx, dlogrho, dlogT, dmetal, vars, sp);
 
+      } else {
 #pragma gpu
         do_conductivity(AMREX_INT_ANYD(bx.loVect()), AMREX_INT_ANYD(bx.hiVect()),
-                        BL_TO_FORTRAN_ANYD(state[mfi]), n_cell);
+                        dlogrho, dlogT, dmetal,
+                        BL_TO_FORTRAN_ANYD(state[mfi]));
 
+      }
     }
 
     // Call the timer again and compute the maximum difference between
@@ -146,19 +200,18 @@ void main_main ()
     const int IOProc = ParallelDescriptor::IOProcessorNumber();
     ParallelDescriptor::ReduceRealMax(stop_time, IOProc);
 
-    std::string name = "test_conductivity.";
-
-    // get the name of the conductivity routine
-    int cond_len = -1;
-    get_cond_len(&cond_len);
 
-    char* cond_string[cond_len+1];
-    get_cond_name(cond_string);
-    std::string cond(*cond_string);
+    std::string name = "test_conductivity.";
+    std::string language = do_cxx == 1 ? ".cxx" : "";
 
     // Write a plotfile
-    WriteSingleLevelPlotfile(name + cond, state, varnames, geom, time, 0);
+    if (do_cxx == 1) {
+      WriteSingleLevelPlotfile(name + cond_name + language, state, vars.names, geom, time, 0);
+    } else {
+      WriteSingleLevelPlotfile(name + cond_name + language, state, varnames, geom, time, 0);
+    }
 
     // Tell the I/O Processor to write out the "run time"
     amrex::Print() << "Run time = " << stop_time << std::endl;
+
 }

unit_test/test_conductivity/test_cond.H

@@ -0,0 +1,14 @@
+#ifndef TEST_COND_H
+#define TEST_COND_H
+
+#include "extern_parameters.H"
+#include "variables.H"
+
+void main_main();
+
+void cond_test_C(const Box& bx,
+                 const Real dlogrho, const Real dlogT, const Real dmetal,
+                 const plot_t vars,
+                 Array4<Real> const sp);
+
+#endif

unit_test/test_conductivity/test_conductivity_F.H → unit_test/test_conductivity/test_cond_F.H

@@ -1,5 +1,5 @@
-#ifndef TEST_EOS_F_H_
-#define TEST_EOS_F_H_
+#ifndef TEST_COND_F_H
+#define TEST_COND_F_H
 
 #include <AMReX_BLFort.H>
 
@@ -8,7 +8,8 @@
 extern "C"
 {
 #endif
-  void init_variables();
+
+  void init_variables_F();
 
   void get_ncomp(int* ncomp);
 
@@ -23,8 +24,8 @@ extern "C"
   void init_unit_test(const int* name, const int* namlen); 
 
   void do_conductivity(const int* lo, const int* hi,
-                       amrex::Real* state, const int* s_lo, const int* s_hi,
-                       const int npts);
+                       const Real dlogrho, const Real dlogT, const Real dmetal,
+                       amrex::Real* state, const int* s_lo, const int* s_hi);
 
 #ifdef __cplusplus
 }

unit_test/test_conductivity/variables.F90 → unit_test/test_conductivity/variables_F.F90

@@ -57,7 +57,7 @@ function get_next_plot_index(this, num) result (next)
     return
   end function get_next_plot_index
 
-  subroutine init_variables() bind(C, name="init_variables")
+  subroutine init_variables_F() bind(C, name="init_variables_F")
 
     integer :: n
 
@@ -86,7 +86,7 @@ subroutine init_variables() bind(C, name="init_variables")
        p % names(p % ispec + n) = "X_" // adjustl(trim(spec_names(n+1)))
     enddo
 
-  end subroutine init_variables
+  end subroutine init_variables_F
 
   subroutine get_ncomp(ncomp_in) bind(C, name="get_ncomp")