switch the test_react init to C++ (#448)

unit_test/test_eos/main.cpp

@@ -35,8 +35,6 @@ void main_main ()
 
     // AMREX_SPACEDIM: number of dimensions
     int n_cell, max_grid_size, do_cxx;
-    Vector<int> bc_lo(AMREX_SPACEDIM,0);
-    Vector<int> bc_hi(AMREX_SPACEDIM,0);
 
     // inputs parameters
     {
@@ -57,7 +55,7 @@ void main_main ()
 
     }
 
-    Vector<int> is_periodic(AMREX_SPACEDIM,0);
+    Vector<int> is_periodic(AMREX_SPACEDIM, 0);
     for (int idim=0; idim < AMREX_SPACEDIM; ++idim) {
       is_periodic[idim] = 1;
     }

unit_test/test_react/GNUmakefile

@@ -16,7 +16,7 @@ EBASE = main
 
 USE_CXX_EOS = TRUE
 
-USE_CXX_REACTIONS ?= FALSE
+USE_CXX_REACTIONS = TRUE
 
 ifeq ($(USE_CXX_REACTIONS),TRUE)
   DEFINES += -DCXX_REACTIONS

unit_test/test_react/Make.package

@@ -11,3 +11,4 @@ f90EXE_sources += util.f90
 CEXE_headers += variables.H
 CEXE_sources += variables.cpp
 CEXE_headers += react_zones.H
+CEXE_headers += react_util.H

unit_test/test_react/main.cpp

@@ -7,7 +7,6 @@
 #include <AMReX_iMultiFab.H>
 #include <AMReX_BCRec.H>
 
-
 using namespace amrex;
 
 #include <test_react.H>
@@ -21,7 +20,10 @@ using namespace amrex;
 #include <AMReX_buildInfo.H>
 #include <variables.H>
 #include <unit_test.H>
-
+#include <react_util.H>
+#ifdef NSE_THERMO
+#include <nse.H>
+#endif
 int main (int argc, char* argv[])
 {
     amrex::Initialize(argc, argv);
@@ -37,8 +39,6 @@ void main_main ()
 
     // AMREX_SPACEDIM: number of dimensions
     int n_cell, max_grid_size, print_every_nrhs;
-    Vector<int> bc_lo(AMREX_SPACEDIM,0);
-    Vector<int> bc_hi(AMREX_SPACEDIM,0);
 
     std::string prefix = "plt";
 
@@ -149,13 +149,55 @@ void main_main ()
     MultiFab state(ba, dm, vars.n_plot_comps, Nghost);
 
     // Initialize the state
+
+    Real dlogrho;
+    Real dlogT;
+
+    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);
+    } else {
+        dlogrho = 0.0_rt;
+        dlogT   = 0.0_rt;
+    }
+
+    init_t comp_data = setup_composition(n_cell);
+
     for ( MFIter mfi(state); mfi.isValid(); ++mfi )
     {
         const Box& bx = mfi.validbox();
 
-        init_state(AMREX_ARLIM_ANYD(bx.loVect()), AMREX_ARLIM_ANYD(bx.hiVect()),
-                   BL_TO_FORTRAN_ANYD(state[mfi]), &n_cell);
-
+        auto state_arr = state.array(mfi);
+
+        amrex::ParallelFor(bx,
+        [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
+        {
+
+            state_arr(i, j, k, vars.itemp) =
+                std::pow(10.0_rt, (std::log10(temp_min) + static_cast<Real>(j)*dlogT));
+            state_arr(i, j, k, vars.irho) =
+                std::pow(10.0_rt, (std::log10(dens_min) + static_cast<Real>(i)*dlogrho));
+
+            Real xn[NumSpec];
+            get_xn(k, comp_data, xn);
+
+            for (int n = 0; n < NumSpec; n++) {
+                state_arr(i, j, k, vars.ispec_old+n) =
+                    amrex::max(xn[n], 1.e-10_rt);
+            }
+
+            // initialize the auxillary state (in particular, for NSE)
+#ifdef NSE_THERMO
+            eos_t eos_state;
+            for (int n = 0; n < NumSpec; n++) {
+                eos_state.xn[n] = xn[n];
+            }
+            set_nse_aux_from_X(eos_state);
+            for (int n = 0; n < NumAux; n++) {
+                state_arr(i, j, k, vars.iaux_old+n) = eos_state.aux[n];
+            }
+#endif
+        });
     }
 
     // allocate a multifab for the number of RHS calls

unit_test/test_react/react_util.H

@@ -0,0 +1,143 @@
+#ifndef REACT_UTIL_H
+#define REACT_UTIL_H
+
+using namespace amrex;
+
+struct init_t {
+    int nprim;
+    int is1;
+    int is2;
+    int is3;
+
+    int n1;
+    int n2;
+    int n3;
+
+    Real Xp_min;
+    Real Xp_max;
+};
+
+AMREX_INLINE AMREX_GPU_HOST_DEVICE
+init_t setup_composition(const int nz) {
+
+    // get the primary species indices
+    init_t comp_data;
+
+    comp_data.nprim = 0;
+
+    comp_data.is1 = network_spec_index(primary_species_1);
+    if (comp_data.is1 >= 0) {
+        comp_data.nprim++;
+    }
+
+    comp_data.is2 = network_spec_index(primary_species_2);
+    if (comp_data.is2 >= 0) {
+        comp_data.nprim++;
+    }
+
+    comp_data.is3 = network_spec_index(primary_species_3);
+    if (comp_data.is3 >= 0) {
+        comp_data.nprim++;
+    }
+
+    if (comp_data.nprim == 0) {
+        amrex::Error("ERROR: no primary species set");
+    }
+
+    // figure out how many zones to allocate to the each of the primary
+    // species and the extrema for the primary species
+
+    if (comp_data.nprim == 1) {
+        comp_data.n1 = nz;
+        comp_data.n2 = 0;
+        comp_data.n3 = 0;
+
+        comp_data.Xp_min = 0.2_rt;
+        comp_data.Xp_max = 0.9_rt;
+
+    } else if (comp_data.nprim == 2) {
+        comp_data.n1 = nz/2;
+        comp_data.n2 = nz - comp_data.n1;
+        comp_data.n3 = 0;
+
+        comp_data.Xp_min = 0.2_rt;
+        comp_data.Xp_max = 0.8_rt;
+
+    } else if (comp_data.nprim == 3) {
+        comp_data.n1 = nz/3;
+        comp_data.n2 = nz/3;
+        comp_data.n3 = nz - comp_data.n1 - comp_data.n2;
+
+        comp_data.Xp_min = 0.2_rt;
+        comp_data.Xp_max = 0.7_rt;
+
+    }
+
+    return comp_data;
+}
+
+
+AMREX_INLINE AMREX_GPU_HOST_DEVICE
+void get_xn(const int k, const init_t cd, Real *xn_zone) {
+
+    for (int n = 0; n < NumSpec; n++) {
+        xn_zone[n] = 0.0_rt;
+    }
+
+    if (k < cd.n1) {
+        if (cd.nprim >= 2) xn_zone[cd.is2] = cd.Xp_min/2;
+        if (cd.nprim >= 3) xn_zone[cd.is3] = cd.Xp_min/2;
+
+        Real dX = (cd.Xp_max - cd.Xp_min) /
+            amrex::max((cd.n1 - 1), 1);
+        xn_zone[cd.is1] = cd.Xp_min + k * dX;
+
+    } else if (cd.nprim >= 2 && k < cd.n1 + cd.n2) {
+        xn_zone[cd.is1] = cd.Xp_min/2;
+        if (cd.nprim >= 3) {
+            xn_zone[cd.is3] = cd.Xp_min/2;
+        }
+        Real dX = (cd.Xp_max - cd.Xp_min) /
+            amrex::max((cd.n2 - 1), 1);
+        xn_zone[cd.is2] = cd.Xp_min + (k - cd.n1) * dX;
+
+    } else {
+        xn_zone[cd.is1] = cd.Xp_min/2;
+        xn_zone[cd.is2] = cd.Xp_min/2;
+
+        Real dX = (cd.Xp_max - cd.Xp_min) /
+            amrex::max((cd.n3 - 1), 1);
+        xn_zone[cd.is3] = cd.Xp_min + (k - (cd.n1 + cd.n2)) * dX;
+
+    }
+
+    Real excess = 0.0_rt;
+    for (int n = 0; n < NumSpec; n++) {
+        excess += xn_zone[n];
+    }
+    excess = 1.0_rt - excess;
+
+    for (int n = 0; n < NumSpec; n++) {
+        if (n == cd.is1 ||
+            (n == cd.is2 && cd.nprim >= 2) ||
+            (n == cd.is3 && cd.nprim >= 3)) {
+            continue;
+        }
+
+        xn_zone[n] = excess / (NumSpec - cd.nprim);
+    }
+
+
+    // normalize -- just in case
+
+    Real sum_X = 0.0_rt;
+    for (int n = 0; n < NumSpec; n++) {
+        sum_X += xn_zone[n];
+    }
+    for (int n = 0; n < NumSpec; n++) {
+        xn_zone[n] = xn_zone[n] / sum_X;
+    }
+
+}
+#endif
+

unit_test/test_react/react_zones.F90

@@ -15,74 +15,6 @@ module react_zones_module
 
 contains
 
-  subroutine init_state(lo, hi, &
-                        state, s_lo, s_hi, npts) bind(C, name="init_state")
-
-    integer, intent(in) :: lo(3), hi(3)
-    integer, intent(in) :: s_lo(3), s_hi(3)
-    real(rt), intent(inout) :: state(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) :: npts
-
-    real(rt) :: dlogrho, dlogT
-    real(rt), allocatable :: xn_zone(:,:)
-    real(rt) :: xn(nspec)
-
-    integer :: ii, jj, kk
-    real(rt) :: sum_X, mu_e
-
-    if (npts > 1) then
-       dlogrho = (log10(dens_max) - log10(dens_min))/(npts - 1)
-       dlogT   = (log10(temp_max) - log10(temp_min))/(npts - 1)
-    else
-       dlogrho = ZERO
-       dlogT   = ZERO
-    endif
-
-    allocate(xn_zone(nspec, 0:npts-1))   ! this assumes that lo(3) = 0
-
-    call get_xn(npts, xn_zone)
-
-    ! normalize -- just in case
-    do kk = lo(3), hi(3)
-       sum_X = sum(xn_zone(:, kk))
-       xn_zone(:, kk) = xn_zone(:, kk)/sum_X
-    enddo
-
-    do kk = lo(3), hi(3)   ! xn loop
-       do jj = lo(2), hi(2)   ! T loop
-          do ii = lo(1), hi(1)   ! rho loop
-
-             state(ii, jj, kk, p % itemp) = 10.0_rt**(log10(temp_min) + dble(jj)*dlogT)
-             state(ii, jj, kk, p % irho)  = 10.0_rt**(log10(dens_min) + dble(ii)*dlogrho)
-             state(ii, jj, kk, p % ispec_old:p % ispec_old+nspec-1) = max(xn_zone(:, kk), 1.e-10_rt)
-
-          enddo
-       enddo
-    enddo
-
-    ! initialize the auxillary state (in particular, for NSE)
-#ifdef NSE_THERMO
-    if (naux > 0) then
-       do kk = lo(3), hi(3)   ! xn loop
-          do jj = lo(2), hi(2)   ! T loop
-             do ii = lo(1), hi(1)   ! rho loop
-
-                xn(:) = state(ii,jj,kk, p % ispec_old:p % ispec_old+nspec-1)
-
-                mu_e = 1.0_rt / sum(xn(:) * zion(:) * aion_inv(:))
-                state(ii,jj,kk, p %  iaux_old+iye-1) = 1.0_rt / mu_e
-                state(ii,jj,kk, p %  iaux_old+iabar-1) = ONE / (sum(xn(:) * aion_inv(:)))
-                state(ii,jj,kk, p %  iaux_old+ibea-1) = (sum(xn(:) * bion(:) * aion_inv(:)))
-
-             end do
-          end do
-       end do
-    end if
-#endif
-
-  end subroutine init_state
-
-
   subroutine print_nrhs(lo, hi, state, s_lo, s_hi) bind(C, name="print_nrhs")
 
     implicit none

unit_test/test_react/test_react_F.H

@@ -16,10 +16,6 @@ extern "C"
                 amrex::Real* state, const int* s_lo, const int* s_hi,
                 int* n_rhs, const int* n_rhs_lo, const int* n_rhs_hi);
 
-  void init_state(const int* lo, const int* hi,
-                  amrex::Real* state, const int* s_lo, const int* s_hi,
-                  const int* npts);
-
   void print_nrhs(const int* lo, const int* hi, const int* state, const int* s_lo, const int* s_hi);
 
 #ifdef __cplusplus