diff --git a/Source/Hydro/CAMR_construct_hydro_source.cpp b/Source/Hydro/CAMR_construct_hydro_source.cpp index 7be65b7..c749735 100644 --- a/Source/Hydro/CAMR_construct_hydro_source.cpp +++ b/Source/Hydro/CAMR_construct_hydro_source.cpp @@ -104,12 +104,14 @@ CAMR::construct_hydro_source (const MultiFab& S, BL_PROFILE_VAR("ctoprim()", ctop); const Real small_num = CAMRConstants::small_num; const Real dual_energy_eta = CAMR::dual_energy_eta1; + int l_allow_negative_energy = CAMR::allow_negative_energy; ParallelFor( qbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { #ifdef AMREX_USE_EB if (!flag_arr(i,j,k).isCovered()) { #endif - hydro_ctoprim(i, j, k, sarr, qarr, qauxar, *lpmap, small_num, dual_energy_eta); + hydro_ctoprim(i, j, k, sarr, qarr, qauxar, *lpmap, + small_num, dual_energy_eta, l_allow_negative_energy); #ifdef AMREX_USE_EB } else { for (int n=0; n const& q, amrex::Array4< amrex::Real> const& qa, PassMap const& pmap, - amrex::Real small_num, - amrex::Real dual_energy_eta1) + amrex::Real l_small_num, + amrex::Real l_dual_energy_eta1, + const int l_allow_negative_energy) { const amrex::Real rho = u(i, j, k, URHO); @@ -51,9 +52,11 @@ hydro_ctoprim( // of the separately updated internal energy equation. // Otherwise, we'll set e = E - K. - AMREX_ALWAYS_ASSERT(u(i,j,k,UEDEN) > 0.); + if (!l_allow_negative_energy) { + AMREX_ALWAYS_ASSERT(u(i,j,k,UEDEN) > 0.); + } - if ( (u(i,j,k,UEDEN) - kineng) / u(i,j,k,UEDEN) > dual_energy_eta1) { + if ( (u(i,j,k,UEDEN) - kineng) / u(i,j,k,UEDEN) > l_dual_energy_eta1) { q(i,j,k,QREINT) = (u(i,j,k,UEDEN) - kineng) * rhoinv; } else { q(i,j,k,QREINT) = u(i,j,k,UEINT) * rhoinv; @@ -61,7 +64,9 @@ hydro_ctoprim( const amrex::Real e = q(i,j,k,QREINT); - AMREX_ALWAYS_ASSERT(u(i,j,k,UEINT) > 0.); + if (!l_allow_negative_energy) { + AMREX_ALWAYS_ASSERT(u(i,j,k,UEINT) > 0.); + } amrex::Real T = u(i, j, k, UTEMP); amrex::Real massfrac[NUM_SPECIES]; @@ -91,7 +96,7 @@ hydro_ctoprim( qa(i, j, k, QGAMC) = gam1; cs = std::sqrt(gam1*p/rho); qa(i, j, k, QC) = cs; - qa(i, j, k, QCSML) = std::max(small_num, small_num * cs); + qa(i, j, k, QCSML) = std::max(l_small_num, l_small_num * cs); } } diff --git a/Source/Params/_cpp_parameters b/Source/Params/_cpp_parameters index 9904fab..f2e1ea3 100644 --- a/Source/Params/_cpp_parameters +++ b/Source/Params/_cpp_parameters @@ -91,7 +91,7 @@ dual_energy_eta2 Real 1.0e-4 use_pslope int 0 # Whether or not to allow internal energy to be less than zero -allow_negative_energy int 1 +allow_negative_energy int 0 # Whether or not to allow the internal energy to be less than the # internal energy corresponding to small\_temp diff --git a/Source/Params/param_includes/CAMR_defaults.H b/Source/Params/param_includes/CAMR_defaults.H index 7c6cb1f..3b450e2 100644 --- a/Source/Params/param_includes/CAMR_defaults.H +++ b/Source/Params/param_includes/CAMR_defaults.H @@ -24,7 +24,7 @@ int CAMR::dual_energy_update_E_from_e = 1; amrex::Real CAMR::dual_energy_eta1 = 1.0e0; amrex::Real CAMR::dual_energy_eta2 = 1.0e-4; int CAMR::use_pslope = 0; -int CAMR::allow_negative_energy = 1; +int CAMR::allow_negative_energy = 0; int CAMR::allow_small_energy = 1; int CAMR::transverse_reset_density = 1; int CAMR::eb_weights_type = 2;