diff --git a/unit_test/test_eos/main.cpp b/unit_test/test_eos/main.cpp index 0276483165..e40b7bb0d3 100644 --- a/unit_test/test_eos/main.cpp +++ b/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 bc_lo(AMREX_SPACEDIM,0); - Vector bc_hi(AMREX_SPACEDIM,0); // inputs parameters { @@ -57,7 +55,7 @@ void main_main () } - Vector is_periodic(AMREX_SPACEDIM,0); + Vector is_periodic(AMREX_SPACEDIM, 0); for (int idim=0; idim < AMREX_SPACEDIM; ++idim) { is_periodic[idim] = 1; } diff --git a/unit_test/test_react/GNUmakefile b/unit_test/test_react/GNUmakefile index 502dbba9f7..64646e8af9 100644 --- a/unit_test/test_react/GNUmakefile +++ b/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 diff --git a/unit_test/test_react/Make.package b/unit_test/test_react/Make.package index 199be48f68..80cf07e2b0 100644 --- a/unit_test/test_react/Make.package +++ b/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 diff --git a/unit_test/test_react/main.cpp b/unit_test/test_react/main.cpp index e2e799f9dd..6fb8d5e82d 100644 --- a/unit_test/test_react/main.cpp +++ b/unit_test/test_react/main.cpp @@ -7,7 +7,6 @@ #include #include - using namespace amrex; #include @@ -21,7 +20,10 @@ using namespace amrex; #include #include #include - +#include +#ifdef NSE_THERMO +#include +#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 bc_lo(AMREX_SPACEDIM,0); - Vector 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(j)*dlogT)); + state_arr(i, j, k, vars.irho) = + std::pow(10.0_rt, (std::log10(dens_min) + static_cast(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 diff --git a/unit_test/test_react/react_util.H b/unit_test/test_react/react_util.H new file mode 100644 index 0000000000..2025d6d779 --- /dev/null +++ b/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 + diff --git a/unit_test/test_react/react_zones.F90 b/unit_test/test_react/react_zones.F90 index df0e20c85b..f03b37343c 100644 --- a/unit_test/test_react/react_zones.F90 +++ b/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 diff --git a/unit_test/test_react/test_react_F.H b/unit_test/test_react/test_react_F.H index a6aa7b7981..9e78952560 100644 --- a/unit_test/test_react/test_react_F.H +++ b/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