From f8094cb7679bbad89e931f98c3985a401e468213 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 14:11:03 +0200 Subject: [PATCH 01/56] compute pressure from microphysics eos --- src/hydro_system.hpp | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 33c787291..54ee9b9d3 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -28,6 +28,16 @@ #include "radiation_system.hpp" #include "valarray.hpp" +// Microphysics headers +#include "burn_type.H" +#include "eos.H" +#include "extern_parameters.H" + +#ifdef PRIMORDIAL_CHEM +#include "actual_eos_data.H" +#endif + + // this struct is specialized by the user application code // template struct HydroSystem_Traits { @@ -300,7 +310,31 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; - const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); + +#ifdef PRIMORDIAL_CHEM + auto massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + if (massScalars) { + + eos_t chemstate; // cannot use burn_t here; it does not contain pressure + chemstate.rho = rho; + chemstate.e = thermal_energy / rho; + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + eos(eos_input_re, chemstate); + const auto P = chemstate.p; + } +#else + chem_eos_t estate; + estate.rho = rho; + estate.e = thermal_energy / rho; + estate.mu = quokka::EOS_Traits::mean_molecular_weight / C::m_u; + eos(eos_input_re, estate); + const auto P = estate.p; +#endif + + // const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); return P; } From 289e1605fb58610813e118d384d9cd0b0177315a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 26 Jul 2023 12:12:11 +0000 Subject: [PATCH 02/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro_system.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 54ee9b9d3..5817d8afd 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -37,7 +37,6 @@ #include "actual_eos_data.H" #endif - // this struct is specialized by the user application code // template struct HydroSystem_Traits { @@ -317,7 +316,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure eos_t chemstate; // cannot use burn_t here; it does not contain pressure chemstate.rho = rho; - chemstate.e = thermal_energy / rho; + chemstate.e = thermal_energy / rho; const auto &massArray = *massScalars; for (int nn = 0; nn < nmscalars_; ++nn) { chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) From 343429c1a9ac79defa437b441c800df408db50e5 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 20:19:27 +0200 Subject: [PATCH 03/56] call EOS to compute pressure --- src/hydro_system.hpp | 26 +++----------------------- 1 file changed, 3 insertions(+), 23 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 5817d8afd..6ac38edcb 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -309,31 +309,11 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; + // const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); -#ifdef PRIMORDIAL_CHEM - auto massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); - if (massScalars) { - - eos_t chemstate; // cannot use burn_t here; it does not contain pressure - chemstate.rho = rho; - chemstate.e = thermal_energy / rho; - const auto &massArray = *massScalars; - for (int nn = 0; nn < nmscalars_; ++nn) { - chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) - } - eos(eos_input_re, chemstate); - const auto P = chemstate.p; - } -#else - chem_eos_t estate; - estate.rho = rho; - estate.e = thermal_energy / rho; - estate.mu = quokka::EOS_Traits::mean_molecular_weight / C::m_u; - eos(eos_input_re, estate); - const auto P = estate.p; -#endif + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); - // const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); return P; } From 642b3afc4aa86c161d0d18d3394a9484839c0ba7 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 20:20:08 +0200 Subject: [PATCH 04/56] add function to compute pressure --- src/EOS.hpp | 40 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/EOS.hpp b/src/EOS.hpp index 7dbf8f75d..f101fafc0 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -48,6 +48,8 @@ template class EOS [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputeEintTempDerivative(amrex::Real rho, amrex::Real Tgas, const std::optional> massScalars = {}) -> amrex::Real; + [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto + ComputePressure(amrex::Real rho, amrex::Real Eint, const std::optional> massScalars = {}) -> amrex::Real; private: static constexpr amrex::Real gamma_ = EOS_Traits::gamma; @@ -176,6 +178,44 @@ EOS::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re return dEint_dT; } +template +AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(amrex::Real rho, amrex::Real Eint, + const std::optional> massScalars) + -> amrex::Real +{ + // return pressure for an ideal gas + amrex::Real P = NAN; +#ifdef PRIMORDIAL_CHEM + eos_t chemstate; // cannot use burn_t because it does not have pressure + chemstate.rho = rho; + chemstate.e = Eint / rho; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } + + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } + + eos(eos_input_re, chemstate); + P = chemstate.p; +#else + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.e = Eint / rho; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_re, estate); + P = estate.p; + } +#endif + return P; +} + } // namespace quokka #endif // EOS_HPP_ From 3db466e3f9b3063d926a8ceda7cf10f52b1a0801 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 26 Jul 2023 18:21:48 +0000 Subject: [PATCH 05/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/EOS.hpp | 58 ++++++++++++++++++++++---------------------- src/hydro_system.hpp | 4 +-- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index f101fafc0..2febb2a72 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -180,40 +180,40 @@ EOS::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re template AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(amrex::Real rho, amrex::Real Eint, - const std::optional> massScalars) + const std::optional> massScalars) -> amrex::Real { - // return pressure for an ideal gas - amrex::Real P = NAN; + // return pressure for an ideal gas + amrex::Real P = NAN; #ifdef PRIMORDIAL_CHEM - eos_t chemstate; // cannot use burn_t because it does not have pressure - chemstate.rho = rho; - chemstate.e = Eint / rho; - // initialize array of number densities - for (int ii = 0; ii < NumSpec; ++ii) { - chemstate.xn[ii] = -1.0; - } - - if (massScalars) { - const auto &massArray = *massScalars; - for (int nn = 0; nn < nmscalars_; ++nn) { - chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) - } - } - - eos(eos_input_re, chemstate); - P = chemstate.p; + eos_t chemstate; // cannot use burn_t because it does not have pressure + chemstate.rho = rho; + chemstate.e = Eint / rho; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } + + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } + + eos(eos_input_re, chemstate); + P = chemstate.p; #else - if constexpr (gamma_ != 1.0) { - chem_eos_t estate; - estate.rho = rho; - estate.e = Eint / rho; - estate.mu = mean_molecular_weight_ / C::m_u; - eos(eos_input_re, estate); - P = estate.p; - } + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.e = Eint / rho; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_re, estate); + P = estate.p; + } #endif - return P; + return P; } } // namespace quokka diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 6ac38edcb..f73241ae9 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -311,8 +311,8 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto thermal_energy = E - kinetic_energy; // const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); - amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); return P; } From aa5a3f60461df09d06849ca3bdb1cf3919ad516f Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 20:26:17 +0200 Subject: [PATCH 06/56] return pressure times density --- src/EOS.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 2febb2a72..640cbf896 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -202,7 +202,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am } eos(eos_input_re, chemstate); - P = chemstate.p; + P = chemstate.p * rho; #else if constexpr (gamma_ != 1.0) { chem_eos_t estate; @@ -210,7 +210,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am estate.e = Eint / rho; estate.mu = mean_molecular_weight_ / C::m_u; eos(eos_input_re, estate); - P = estate.p; + P = estate.p * rho; } #endif return P; From 9e28e431af02cfab15bdb4707020a9a9c3ba2c79 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 20:32:49 +0200 Subject: [PATCH 07/56] updated microp --- .gitmodules | 2 +- extern/Microphysics | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.gitmodules b/.gitmodules index f61e8946a..3cee8997f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -10,4 +10,4 @@ [submodule "extern/Microphysics"] path = extern/Microphysics url = https://github.com/psharda/Microphysics - branch = development + branch = compute-pressure diff --git a/extern/Microphysics b/extern/Microphysics index ea4204146..8a627847f 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit ea4204146eb2269605764e432a022cb1ae090398 +Subproject commit 8a627847f891ab8c0eb8ce859c56c95a9e41bde8 From 75baed9b39d918c279be82cf8cc0f7fcd7c48640 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 21:02:19 +0200 Subject: [PATCH 08/56] change mean mol weight from NAN to arbitrary +ve value to get finite pressure from microp gamma_law --- src/HydroHighMach/test_hydro_highmach.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HydroHighMach/test_hydro_highmach.cpp b/src/HydroHighMach/test_hydro_highmach.cpp index 82075656a..9b8d74fb3 100644 --- a/src/HydroHighMach/test_hydro_highmach.cpp +++ b/src/HydroHighMach/test_hydro_highmach.cpp @@ -30,7 +30,7 @@ struct HighMachProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; From 580f3e51fcba8b061953641a80b79327f2229fcd Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 21:30:53 +0200 Subject: [PATCH 09/56] use gamma_law to compute pressure --- src/hydro_system.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index f73241ae9..245b600b2 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -159,7 +159,7 @@ template void HydroSystem::ConservedToPrimitive( const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto Eint_cons = E - kinetic_energy; - const amrex::Real Pgas = Eint_cons * (HydroSystem::gamma_ - 1.0); + const amrex::Real Pgas = ComputePressure(cons[bx], i, j, k); const amrex::Real eint_cons = Eint_cons / rho; const amrex::Real eint_aux = Eint_aux / rho; @@ -213,7 +213,7 @@ template auto HydroSystem::maxSignalSpeedLocal(a } else { const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); const auto Eint = Etot - kinetic_energy; - const auto P = Eint * (HydroSystem::gamma_ - 1.0); + const auto P = ComputePressure(cons[bx], i, j, k); // const auto P = Eint * (HydroSystem::gamma_ - 1.0); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } return {cs + abs_vel}; From c4a36f0380ae712e1420b5c375fa8fba5d89da67 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 26 Jul 2023 19:31:07 +0000 Subject: [PATCH 10/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro_system.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 245b600b2..2d8b44a12 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -213,7 +213,8 @@ template auto HydroSystem::maxSignalSpeedLocal(a } else { const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); const auto Eint = Etot - kinetic_energy; - const auto P = ComputePressure(cons[bx], i, j, k); // const auto P = Eint * (HydroSystem::gamma_ - 1.0); + const auto P = + ComputePressure(cons[bx], i, j, k); // const auto P = Eint * (HydroSystem::gamma_ - 1.0); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } return {cs + abs_vel}; From 380fec701191486fcc6cf2725d106144fbf106cb Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 22:04:57 +0200 Subject: [PATCH 11/56] define a finite mean molecular weight --- src/HydroShocktube/test_hydro_shocktube.cpp | 2 +- src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp | 2 +- src/HydroVacuum/test_hydro_vacuum.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/HydroShocktube/test_hydro_shocktube.cpp b/src/HydroShocktube/test_hydro_shocktube.cpp index a1e6036f2..fbb5417ea 100644 --- a/src/HydroShocktube/test_hydro_shocktube.cpp +++ b/src/HydroShocktube/test_hydro_shocktube.cpp @@ -28,7 +28,7 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp b/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp index fa7f921f1..a21f5a502 100644 --- a/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp +++ b/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp @@ -33,7 +33,7 @@ template <> struct SimulationData { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/HydroVacuum/test_hydro_vacuum.cpp b/src/HydroVacuum/test_hydro_vacuum.cpp index 42b8fcf99..bb02e7fcb 100644 --- a/src/HydroVacuum/test_hydro_vacuum.cpp +++ b/src/HydroVacuum/test_hydro_vacuum.cpp @@ -27,7 +27,7 @@ struct ShocktubeProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; From fb67ff572b95bd182b5167aad8bf65ee1e96908f Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 22:09:19 +0200 Subject: [PATCH 12/56] more EOS calls --- src/hydro_system.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 2d8b44a12..2069f8942 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -213,8 +213,7 @@ template auto HydroSystem::maxSignalSpeedLocal(a } else { const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); const auto Eint = Etot - kinetic_energy; - const auto P = - ComputePressure(cons[bx], i, j, k); // const auto P = Eint * (HydroSystem::gamma_ - 1.0); + const auto P = ComputePressure(cons[bx], i, j, k); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } return {cs + abs_vel}; @@ -247,7 +246,7 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4::gamma_ - 1.0); + const auto P = ComputePressure(cons, i, j, k); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } AMREX_ASSERT(cs > 0.); @@ -275,7 +274,7 @@ template auto HydroSystem::CheckStatesValid(amre const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; - const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); + const auto P = ComputePressure(cons[bx], i, j, k); bool negativeDensity = (rho <= 0.); bool negativePressure = (P <= 0.); @@ -310,7 +309,6 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; - // const auto P = thermal_energy * (HydroSystem::gamma_ - 1.0); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); From e94edecaa145cadb747e767c1d51ebe741084125 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 22:11:51 +0200 Subject: [PATCH 13/56] dont need microp headers --- src/hydro_system.hpp | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 2069f8942..70ed13291 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -28,15 +28,6 @@ #include "radiation_system.hpp" #include "valarray.hpp" -// Microphysics headers -#include "burn_type.H" -#include "eos.H" -#include "extern_parameters.H" - -#ifdef PRIMORDIAL_CHEM -#include "actual_eos_data.H" -#endif - // this struct is specialized by the user application code // template struct HydroSystem_Traits { From a70fe90ed6f26a1f92494306112fedab21aa130b Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 22:19:37 +0200 Subject: [PATCH 14/56] removed unused var --- src/hydro_system.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 70ed13291..86da237d3 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -203,7 +203,6 @@ template auto HydroSystem::maxSignalSpeedLocal(a cs = cs_iso_; } else { const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); - const auto Eint = Etot - kinetic_energy; const auto P = ComputePressure(cons[bx], i, j, k); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } From 6cea067923240dc92268ffe144d8529ae5961eb7 Mon Sep 17 00:00:00 2001 From: psharda Date: Wed, 26 Jul 2023 22:27:38 +0200 Subject: [PATCH 15/56] remove unused vars --- src/hydro_system.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 86da237d3..eeeb102f3 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -202,7 +202,6 @@ template auto HydroSystem::maxSignalSpeedLocal(a if constexpr (is_eos_isothermal()) { cs = cs_iso_; } else { - const auto Etot = cons[bx](i, j, k, HydroSystem::energy_index); const auto P = ComputePressure(cons[bx], i, j, k); cs = std::sqrt(HydroSystem::gamma_ * P / rho); } @@ -235,7 +234,6 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4::gamma_ * P / rho); } From c9ffd1c75384ffcb9a76f2fe1c829013a7fbc28d Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 00:44:31 +0200 Subject: [PATCH 16/56] remove unused var; assert P --- src/hydro_system.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index eeeb102f3..4d6aa1770 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -231,10 +231,8 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4::gamma_ * P / rho); } AMREX_ASSERT(cs > 0.); From 95db06916c669ee310ebaacb55a8c154ae17fb76 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 12:27:07 +0200 Subject: [PATCH 17/56] fixed warnings --- src/RadMarshak/test_radiation_marshak.cpp | 2 +- src/RadMarshakCGS/test_radiation_marshak_cgs.cpp | 2 +- src/RadMatterCoupling/test_radiation_matter_coupling.cpp | 2 +- .../test_radiation_matter_coupling_rsla.cpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/RadMarshak/test_radiation_marshak.cpp b/src/RadMarshak/test_radiation_marshak.cpp index c4f4ec104..b7986e508 100644 --- a/src/RadMarshak/test_radiation_marshak.cpp +++ b/src/RadMarshak/test_radiation_marshak.cpp @@ -31,7 +31,7 @@ constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson; constexpr double T_initial = 1.0e-2; template <> struct quokka::EOS_Traits { - static constexpr double mean_molecular_mass = 1.0; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = 1.0; static constexpr double gamma = 5. / 3.; }; diff --git a/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp b/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp index 39cdbf0ad..15e208280 100644 --- a/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp +++ b/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp @@ -33,7 +33,7 @@ constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson; constexpr double T_initial = 1.0e4; // K template <> struct quokka::EOS_Traits { - static constexpr double mean_molecular_mass = C::m_u; + static constexpr double mean_molecular_weight = C::m_u; static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; diff --git a/src/RadMatterCoupling/test_radiation_matter_coupling.cpp b/src/RadMatterCoupling/test_radiation_matter_coupling.cpp index 3e554b767..1c2582e50 100644 --- a/src/RadMatterCoupling/test_radiation_matter_coupling.cpp +++ b/src/RadMatterCoupling/test_radiation_matter_coupling.cpp @@ -31,7 +31,7 @@ template <> struct SimulationData { }; template <> struct quokka::EOS_Traits { - static constexpr double mean_molecular_mass = C::m_u; + static constexpr double mean_molecular_weight = C::m_u; static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; diff --git a/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp b/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp index 3c7f2e97c..0d774bfc9 100644 --- a/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp +++ b/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp @@ -33,7 +33,7 @@ template <> struct SimulationData { }; template <> struct quokka::EOS_Traits { - static constexpr double mean_molecular_mass = C::m_u; + static constexpr double mean_molecular_weight = C::m_u; static constexpr double boltzmann_constant = C::k_B; static constexpr double gamma = 5. / 3.; }; From dac7af8af1fe4300eb20cd0d775fcba3bb231293 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 16:02:04 +0200 Subject: [PATCH 18/56] use EOS compute pressure directly --- src/hydro_system.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 4d6aa1770..45c431772 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -824,8 +824,11 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu // (pressure_index is actually eint) const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); - P_L = rho_L * eint_L * (gamma_ - 1.0); - P_R = rho_R * eint_R * (gamma_ - 1.0); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); + // P_L = rho_L * eint_L * (gamma_ - 1.0); + // P_R = rho_R * eint_R * (gamma_ - 1.0); + P_L = quokka::EOS::ComputePressure(rho_L, eint_L, massScalars); + P_R = quokka::EOS::ComputePressure(rho_R, eint_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); From 483d3b7ec21d51ef43b0fbee792893105c997692 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 17:49:02 +0200 Subject: [PATCH 19/56] updated microp --- .gitmodules | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitmodules b/.gitmodules index 3cee8997f..0f61cc08f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -10,4 +10,4 @@ [submodule "extern/Microphysics"] path = extern/Microphysics url = https://github.com/psharda/Microphysics - branch = compute-pressure + branch = compute-cs From 7db199bdda31a14493269ce84a4b55f90c3fcee3 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 17:54:10 +0200 Subject: [PATCH 20/56] updated microp --- extern/Microphysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/Microphysics b/extern/Microphysics index 8a627847f..ed47c5c13 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit 8a627847f891ab8c0eb8ce859c56c95a9e41bde8 +Subproject commit ed47c5c13980d0440d5bc2d9a22a9a91440d5ad9 From 906b4d95780f41cee37b5961e32bb3abba80357d Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:00:30 +0200 Subject: [PATCH 21/56] added ComputeSoundSpeed --- src/EOS.hpp | 52 +++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 640cbf896..5a0464488 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -50,6 +50,8 @@ template class EOS -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputePressure(amrex::Real rho, amrex::Real Eint, const std::optional> massScalars = {}) -> amrex::Real; + [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto + ComputeSoundSpeed(amrex::Real rho, amrex::Real Pres, const std::optional> massScalars = {}) -> amrex::Real; private: static constexpr amrex::Real gamma_ = EOS_Traits::gamma; @@ -66,7 +68,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeTgasFromEin amrex::Real Tgas = NAN; #ifdef PRIMORDIAL_CHEM - burn_t chemstate; + eos_t chemstate; chemstate.rho = rho; chemstate.e = Eint / rho; // initialize array of number densities @@ -106,7 +108,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromTga amrex::Real Eint = NAN; #ifdef PRIMORDIAL_CHEM - burn_t chemstate; + eos_t chemstate; chemstate.rho = rho; // Define and initialize Tgas here amrex::Real const Tgas_value = Tgas; @@ -147,7 +149,7 @@ EOS::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re amrex::Real dEint_dT = NAN; #ifdef PRIMORDIAL_CHEM - burn_t chemstate; + eos_t chemstate; chemstate.rho = rho; // we don't need Tgas to find chemstate.dedT, but we still need to initialize chemstate.T because we are using the 'rt' EOS mode chemstate.T = NAN; @@ -202,7 +204,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am } eos(eos_input_re, chemstate); - P = chemstate.p * rho; + P = chemstate.p; #else if constexpr (gamma_ != 1.0) { chem_eos_t estate; @@ -210,12 +212,52 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am estate.e = Eint / rho; estate.mu = mean_molecular_weight_ / C::m_u; eos(eos_input_re, estate); - P = estate.p * rho; + P = estate.p; } #endif return P; } + +template +AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed(amrex::Real rho, amrex::Real Pres, + const std::optional> massScalars) + -> amrex::Real +{ + // return sound speed for an ideal gas + amrex::Real cs = NAN; + +#ifdef PRIMORDIAL_CHEM + eos_t chemstate; + chemstate.rho = rho; + chemstate.p = Pres; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } + + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } + + eos(eos_input_rp, chemstate); + cs = chemstate.cs; +#else + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.p = Pres; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_rp, estate); + cs = estate.cs; + } +#endif + return cs; +} + } // namespace quokka #endif // EOS_HPP_ From 235c18de78a81f9528628e529f279a84230b328c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 16:02:53 +0000 Subject: [PATCH 22/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/EOS.hpp | 55 ++++++++++++++++++++++++++--------------------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 5a0464488..004740979 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -218,44 +218,43 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am return P; } - template AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed(amrex::Real rho, amrex::Real Pres, - const std::optional> massScalars) + const std::optional> massScalars) -> amrex::Real { - // return sound speed for an ideal gas - amrex::Real cs = NAN; + // return sound speed for an ideal gas + amrex::Real cs = NAN; #ifdef PRIMORDIAL_CHEM - eos_t chemstate; - chemstate.rho = rho; - chemstate.p = Pres; - // initialize array of number densities - for (int ii = 0; ii < NumSpec; ++ii) { - chemstate.xn[ii] = -1.0; - } + eos_t chemstate; + chemstate.rho = rho; + chemstate.p = Pres; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } - if (massScalars) { - const auto &massArray = *massScalars; - for (int nn = 0; nn < nmscalars_; ++nn) { - chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) - } - } + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } - eos(eos_input_rp, chemstate); - cs = chemstate.cs; + eos(eos_input_rp, chemstate); + cs = chemstate.cs; #else - if constexpr (gamma_ != 1.0) { - chem_eos_t estate; - estate.rho = rho; - estate.p = Pres; - estate.mu = mean_molecular_weight_ / C::m_u; - eos(eos_input_rp, estate); - cs = estate.cs; - } + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.p = Pres; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_rp, estate); + cs = estate.cs; + } #endif - return cs; + return cs; } } // namespace quokka From abf95f9c6251c330c0f9ea6657a11f06364520ef Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:06:09 +0200 Subject: [PATCH 23/56] updated microp --- extern/Microphysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/Microphysics b/extern/Microphysics index ed47c5c13..45122a3f3 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit ed47c5c13980d0440d5bc2d9a22a9a91440d5ad9 +Subproject commit 45122a3f346da0f4ddd2bc9f5b5941bf8e14bcf1 From 1112f9f4cb3592aecb455027520c8e11008da8a6 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:12:10 +0200 Subject: [PATCH 24/56] spell check --- src/EOS.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 004740979..57ce9a0f5 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -51,7 +51,7 @@ template class EOS [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputePressure(amrex::Real rho, amrex::Real Eint, const std::optional> massScalars = {}) -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto - ComputeSoundSpeed(amrex::Real rho, amrex::Real Pres, const std::optional> massScalars = {}) -> amrex::Real; + ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure, const std::optional> massScalars = {}) -> amrex::Real; private: static constexpr amrex::Real gamma_ = EOS_Traits::gamma; @@ -188,7 +188,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am // return pressure for an ideal gas amrex::Real P = NAN; #ifdef PRIMORDIAL_CHEM - eos_t chemstate; // cannot use burn_t because it does not have pressure + eos_t chemstate; chemstate.rho = rho; chemstate.e = Eint / rho; // initialize array of number densities @@ -219,7 +219,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am } template -AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed(amrex::Real rho, amrex::Real Pres, +AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed(amrex::Real rho, amrex::Real Pressure, const std::optional> massScalars) -> amrex::Real { @@ -229,7 +229,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed( #ifdef PRIMORDIAL_CHEM eos_t chemstate; chemstate.rho = rho; - chemstate.p = Pres; + chemstate.p = Pressure; // initialize array of number densities for (int ii = 0; ii < NumSpec; ++ii) { chemstate.xn[ii] = -1.0; @@ -248,7 +248,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeSoundSpeed( if constexpr (gamma_ != 1.0) { chem_eos_t estate; estate.rho = rho; - estate.p = Pres; + estate.p = Pressure; estate.mu = mean_molecular_weight_ / C::m_u; eos(eos_input_rp, estate); cs = estate.cs; From ebbfe818173155522223332d2cc89f5c51fdbdff Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:13:55 +0200 Subject: [PATCH 25/56] pass quokka gamma to microphysics eos --- src/RadhydroSimulation.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 0ded9fc67..513712d8c 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -129,6 +129,10 @@ template class RadhydroSimulation : public AMRSimulation::gamma); + // initialize Microphysics params init_extern_parameters(); // initialize Microphysics EOS @@ -142,6 +146,10 @@ template class RadhydroSimulation : public AMRSimulation::gamma); + // initialize Microphysics params init_extern_parameters(); // initialize Microphysics EOS From 27ee2139a0abca7fd9067b924ec8340cc38e1ace Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 16:14:09 +0000 Subject: [PATCH 26/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/RadhydroSimulation.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp index 513712d8c..2a2abf7fb 100644 --- a/src/RadhydroSimulation.hpp +++ b/src/RadhydroSimulation.hpp @@ -130,8 +130,8 @@ template class RadhydroSimulation : public AMRSimulation::gamma); + amrex::ParmParse eos("eos"); + eos.add("eos_gamma", quokka::EOS_Traits::gamma); // initialize Microphysics params init_extern_parameters(); @@ -146,9 +146,9 @@ template class RadhydroSimulation : public AMRSimulation::gamma); + // set gamma + amrex::ParmParse eos("eos"); + eos.add("eos_gamma", quokka::EOS_Traits::gamma); // initialize Microphysics params init_extern_parameters(); From a7648ca98ca852e77fbb2819c254c480188167a4 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:41:54 +0200 Subject: [PATCH 27/56] add function ComputeEintFromPres --- src/EOS.hpp | 47 +++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 2 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 57ce9a0f5..9cfe1f679 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -46,6 +46,8 @@ template class EOS [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputeEintFromTgas(amrex::Real rho, amrex::Real Tgas, const std::optional> massScalars = {}) -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto + ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, const std::optional> massScalars = {}) -> amrex::Real; + [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputeEintTempDerivative(amrex::Real rho, amrex::Real Tgas, const std::optional> massScalars = {}) -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto @@ -64,7 +66,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeTgasFromEin const std::optional> massScalars) -> amrex::Real { - // return temperature for an ideal gas + // return temperature for an ideal gas given density and internal energy amrex::Real Tgas = NAN; #ifdef PRIMORDIAL_CHEM @@ -104,7 +106,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromTga const std::optional> massScalars) -> amrex::Real { - // return internal energy density for a gamma-law ideal gas + // return internal energy density given density and temperature amrex::Real Eint = NAN; #ifdef PRIMORDIAL_CHEM @@ -140,6 +142,47 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromTga return Eint; } + +template +AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, + const std::optional> massScalars) + -> amrex::Real +{ + // return internal energy density given density and pressure + amrex::Real Eint = NAN; + +#ifdef PRIMORDIAL_CHEM + eos_t chemstate; + chemstate.rho = rho; + chemstate.p = Pressure; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } + + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } + + eos(eos_input_rp, chemstate); + Eint = chemstate.e * chemstate.rho; +#else + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.p = Pressure; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_rp, estate); + Eint = estate.e * rho * boltzmann_constant_ / C::k_B; + } +#endif + return Eint; +} + + template AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real Tgas, From e62a6ed8c5cae366aff2c75e72f0be62cf0832be Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 16:43:43 +0000 Subject: [PATCH 28/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/EOS.hpp | 63 ++++++++++++++++++++++++++--------------------------- 1 file changed, 31 insertions(+), 32 deletions(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 9cfe1f679..5dd288163 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -46,7 +46,8 @@ template class EOS [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputeEintFromTgas(amrex::Real rho, amrex::Real Tgas, const std::optional> massScalars = {}) -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto - ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, const std::optional> massScalars = {}) -> amrex::Real; + ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, const std::optional> massScalars = {}) + -> amrex::Real; [[nodiscard]] AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE static auto ComputeEintTempDerivative(amrex::Real rho, amrex::Real Tgas, const std::optional> massScalars = {}) -> amrex::Real; @@ -142,47 +143,45 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromTga return Eint; } - template AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromPres(amrex::Real rho, amrex::Real Pressure, - const std::optional> massScalars) + const std::optional> massScalars) -> amrex::Real { - // return internal energy density given density and pressure - amrex::Real Eint = NAN; + // return internal energy density given density and pressure + amrex::Real Eint = NAN; #ifdef PRIMORDIAL_CHEM - eos_t chemstate; - chemstate.rho = rho; - chemstate.p = Pressure; - // initialize array of number densities - for (int ii = 0; ii < NumSpec; ++ii) { - chemstate.xn[ii] = -1.0; - } - - if (massScalars) { - const auto &massArray = *massScalars; - for (int nn = 0; nn < nmscalars_; ++nn) { - chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) - } - } - - eos(eos_input_rp, chemstate); - Eint = chemstate.e * chemstate.rho; + eos_t chemstate; + chemstate.rho = rho; + chemstate.p = Pressure; + // initialize array of number densities + for (int ii = 0; ii < NumSpec; ++ii) { + chemstate.xn[ii] = -1.0; + } + + if (massScalars) { + const auto &massArray = *massScalars; + for (int nn = 0; nn < nmscalars_; ++nn) { + chemstate.xn[nn] = massArray[nn] / spmasses[nn]; // massScalars are partial densities (massFractions * rho) + } + } + + eos(eos_input_rp, chemstate); + Eint = chemstate.e * chemstate.rho; #else - if constexpr (gamma_ != 1.0) { - chem_eos_t estate; - estate.rho = rho; - estate.p = Pressure; - estate.mu = mean_molecular_weight_ / C::m_u; - eos(eos_input_rp, estate); - Eint = estate.e * rho * boltzmann_constant_ / C::k_B; - } + if constexpr (gamma_ != 1.0) { + chem_eos_t estate; + estate.rho = rho; + estate.p = Pressure; + estate.mu = mean_molecular_weight_ / C::m_u; + eos(eos_input_rp, estate); + Eint = estate.e * rho * boltzmann_constant_ / C::k_B; + } #endif - return Eint; + return Eint; } - template AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real Tgas, From 1c64f8e0039b579fdc592237a0cb85d081c60002 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:51:12 +0200 Subject: [PATCH 29/56] use EOS to find pres and cs --- src/hydro_system.hpp | 68 ++++++++++++++++++++++++++++++++------------ 1 file changed, 50 insertions(+), 18 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 45c431772..5fc88a200 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -74,6 +74,8 @@ template class HydroSystem : public HyperbolicSystem const &cons, int i, int j, int k) -> amrex::Real; + AMREX_GPU_DEVICE static auto ComputeSoundSpeed(amrex::Array4 const &cons, int i, int j, int k) -> amrex::Real; + AMREX_GPU_DEVICE static auto ComputeVelocityX1(amrex::Array4 const &cons, int i, int j, int k) -> amrex::Real; AMREX_GPU_DEVICE static auto ComputeVelocityX2(amrex::Array4 const &cons, int i, int j, int k) -> amrex::Real; @@ -203,7 +205,7 @@ template auto HydroSystem::maxSignalSpeedLocal(a cs = cs_iso_; } else { const auto P = ComputePressure(cons[bx], i, j, k); - cs = std::sqrt(HydroSystem::gamma_ * P / rho); + cs = ComputeSoundSpeed(cons[bx], i, j, k); } return {cs + abs_vel}; }); @@ -233,7 +235,7 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4::gamma_ * P / rho); + cs = ComputeSoundSpeed(cons, i, j, k); } AMREX_ASSERT(cs > 0.); @@ -302,6 +304,28 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure return P; } +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputeSoundSpeed(amrex::Array4 const &cons, int i, int j, int k) + -> amrex::Real +{ + const auto rho = cons(i, j, k, density_index); + const auto px = cons(i, j, k, x1Momentum_index); + const auto py = cons(i, j, k, x2Momentum_index); + const auto pz = cons(i, j, k, x3Momentum_index); + const auto E = cons(i, j, k, energy_index); // *total* gas energy per unit volume + const auto vx = px / rho; + const auto vy = py / rho; + const auto vz = pz / rho; + const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); + const auto thermal_energy = E - kinetic_energy; + + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + amrex::Real cs = quokka::EOS::ComputeSoundSpeed(rho, P, massScalars); + + return cs; +} + template AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputeVelocityX1(amrex::Array4 const &cons, int i, int j, int k) -> amrex::Real @@ -468,11 +492,16 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const if constexpr (reconstruct_eint) { // compute (rho e) (gamma - 1) - Pplus2 *= primVar(i + 2, j, k, primDensity_index) * (gamma_ - 1.0); - Pplus1 *= primVar(i + 1, j, k, primDensity_index) * (gamma_ - 1.0); - P *= primVar(i, j, k, primDensity_index) * (gamma_ - 1.0); - Pminus1 *= primVar(i - 1, j, k, primDensity_index) * (gamma_ - 1.0); - Pminus2 *= primVar(i - 2, j, k, primDensity_index) * (gamma_ - 1.0); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 2, j, k); + Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), primVar(i + 2, j, k, primDensity_index)*Pplus2, massScalars)); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); + Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), primVar(i + 1, j, k, primDensity_index)*Pplus1, massScalars)); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); + P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index)*P, massScalars)); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); + Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), primVar(i - 1, j, k, primDensity_index)*Pminus1, massScalars)); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); + Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), primVar(i - 2, j, k, primDensity_index)*Pminus2, massScalars)); } if constexpr (is_eos_isothermal()) { @@ -497,7 +526,7 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const const double chi_min = std::max(0., std::min(1., (beta_max - beta) / (beta_max - beta_min))); // Z is a measure of shock strength (Eq. 76 of Miller & Colella 2002) - const double K_S = gamma_ * P; // equal to \rho c_s^2 + const double K_S = BENBEN use computesoundspeed function here in place of gamma_ * P; // equal to \rho c_s^2 const double Z = std::abs(Pplus1 - Pminus1) / K_S; // check for converging flow along the normal direction DIR (Eq. 77) @@ -628,12 +657,14 @@ void HydroSystem::EnforceLimits(amrex::Real const densityFloor, amrex if (!HydroSystem::is_eos_isothermal()) { // recompute gas energy (to prevent P < 0) amrex::Real const Eint_star = Etot - 0.5 * rho_new * vsq; - amrex::Real const P_star = Eint_star * (HydroSystem::gamma_ - 1.); + + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); + amrex::Real const P_star = quokka::EOS::ComputePressure(rho_new, Eint_star, massScalars); amrex::Real P_new = P_star; if (P_star < pressureFloor) { P_new = pressureFloor; #pragma nv_diag_suppress divide_by_zero - amrex::Real const Etot_new = P_new / (HydroSystem::gamma_ - 1.) + 0.5 * rho_new * vsq; + amrex::Real const Etot_new = quokka::EOS::ComputeEintFromPres(rho_new, P_new, massScalars) + 0.5 * rho_new * vsq; state[bx](i, j, k, HydroSystem::energy_index) = Etot_new; } } @@ -825,10 +856,9 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - // P_L = rho_L * eint_L * (gamma_ - 1.0); - // P_R = rho_R * eint_R * (gamma_ - 1.0); - P_L = quokka::EOS::ComputePressure(rho_L, eint_L, massScalars); - P_R = quokka::EOS::ComputePressure(rho_R, eint_R, massScalars); + P_L = quokka::EOS::ComputePressure(rho_L, eint_L*rho_L, massScalars); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + P_R = quokka::EOS::ComputePressure(rho_R, eint_R*rho_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); @@ -843,11 +873,13 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu Eint_R = x1RightState(i, j, k, primEint_index); } - cs_L = std::sqrt(gamma_ * P_L / rho_L); - cs_R = std::sqrt(gamma_ * P_R / rho_R); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); + cs_L = quokka::EOS::ComputeSoundSpeed(rho_L, P_L, massScalars); + E_L = quokka::EOS::ComputeEintFromPres(rho_L, P_L, massScalars) + ke_L; - E_L = P_L / (gamma_ - 1.0) + ke_L; - E_R = P_R / (gamma_ - 1.0) + ke_R; + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + cs_R = quokka::EOS::ComputeSoundSpeed(rho_R, P_R, massScalars); + E_R = quokka::EOS::ComputeEintFromPres(rho_R, P_R, massScalars) + ke_R; } AMREX_ASSERT(cs_L > 0.0); From 6d31572ad388de35549291bc15defdac4404822f Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 16:51:26 +0000 Subject: [PATCH 30/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro_system.hpp | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 5fc88a200..4b3aff7db 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -308,22 +308,22 @@ template AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputeSoundSpeed(amrex::Array4 const &cons, int i, int j, int k) -> amrex::Real { - const auto rho = cons(i, j, k, density_index); - const auto px = cons(i, j, k, x1Momentum_index); - const auto py = cons(i, j, k, x2Momentum_index); - const auto pz = cons(i, j, k, x3Momentum_index); - const auto E = cons(i, j, k, energy_index); // *total* gas energy per unit volume - const auto vx = px / rho; - const auto vy = py / rho; - const auto vz = pz / rho; - const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); - const auto thermal_energy = E - kinetic_energy; - - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); - amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + const auto rho = cons(i, j, k, density_index); + const auto px = cons(i, j, k, x1Momentum_index); + const auto py = cons(i, j, k, x2Momentum_index); + const auto pz = cons(i, j, k, x3Momentum_index); + const auto E = cons(i, j, k, energy_index); // *total* gas energy per unit volume + const auto vx = px / rho; + const auto vy = py / rho; + const auto vz = pz / rho; + const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); + const auto thermal_energy = E - kinetic_energy; + + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); amrex::Real cs = quokka::EOS::ComputeSoundSpeed(rho, P, massScalars); - return cs; + return cs; } template @@ -658,8 +658,8 @@ void HydroSystem::EnforceLimits(amrex::Real const densityFloor, amrex // recompute gas energy (to prevent P < 0) amrex::Real const Eint_star = Etot - 0.5 * rho_new * vsq; - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); - amrex::Real const P_star = quokka::EOS::ComputePressure(rho_new, Eint_star, massScalars); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); + amrex::Real const P_star = quokka::EOS::ComputePressure(rho_new, Eint_star, massScalars); amrex::Real P_new = P_star; if (P_star < pressureFloor) { P_new = pressureFloor; @@ -856,9 +856,9 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - P_L = quokka::EOS::ComputePressure(rho_L, eint_L*rho_L, massScalars); + P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); - P_R = quokka::EOS::ComputePressure(rho_R, eint_R*rho_R, massScalars); + P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); From cdd18da2ab1c77714f6cfa4436c5ba45092663d7 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 18:58:14 +0200 Subject: [PATCH 31/56] replace K_S with sound spped calc --- src/hydro_system.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 4b3aff7db..6bb5fd486 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -526,7 +526,8 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const const double chi_min = std::max(0., std::min(1., (beta_max - beta) / (beta_max - beta_min))); // Z is a measure of shock strength (Eq. 76 of Miller & Colella 2002) - const double K_S = BENBEN use computesoundspeed function here in place of gamma_ * P; // equal to \rho c_s^2 + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); + const double K_S = std::pow(quokka::EOS::ComputeSoundSpeed(primVar(i, j, k, primDensity_index), P, massScalars), 2) * primVar(i, j, k, primDensity_index); const double Z = std::abs(Pplus1 - Pminus1) / K_S; // check for converging flow along the normal direction DIR (Eq. 77) From 54723a3be21cef60d201b7990912d30f25779725 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 16:58:27 +0000 Subject: [PATCH 32/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro_system.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 6bb5fd486..6140aa874 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -527,7 +527,8 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const // Z is a measure of shock strength (Eq. 76 of Miller & Colella 2002) amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); - const double K_S = std::pow(quokka::EOS::ComputeSoundSpeed(primVar(i, j, k, primDensity_index), P, massScalars), 2) * primVar(i, j, k, primDensity_index); + const double K_S = std::pow(quokka::EOS::ComputeSoundSpeed(primVar(i, j, k, primDensity_index), P, massScalars), 2) * + primVar(i, j, k, primDensity_index); const double Z = std::abs(Pplus1 - Pminus1) / K_S; // check for converging flow along the normal direction DIR (Eq. 77) From c7ca8a05dc84c94abe9f9d67f9b041db89eb3b84 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 19:06:30 +0200 Subject: [PATCH 33/56] avoid re-declaaration --- src/hydro_system.hpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 6140aa874..64adf55a1 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -493,15 +493,15 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const if constexpr (reconstruct_eint) { // compute (rho e) (gamma - 1) amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 2, j, k); - Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), primVar(i + 2, j, k, primDensity_index)*Pplus2, massScalars)); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); - Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), primVar(i + 1, j, k, primDensity_index)*Pplus1, massScalars)); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); - P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index)*P, massScalars)); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); - Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), primVar(i - 1, j, k, primDensity_index)*Pminus1, massScalars)); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); - Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), primVar(i - 2, j, k, primDensity_index)*Pminus2, massScalars)); + Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), primVar(i + 2, j, k, primDensity_index)*Pplus2, massScalars); + massScalars = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); + Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), primVar(i + 1, j, k, primDensity_index)*Pplus1, massScalars); + massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); + P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index)*P, massScalars); + massScalars = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); + Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), primVar(i - 1, j, k, primDensity_index)*Pminus1, massScalars); + massScalars = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); + Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), primVar(i - 2, j, k, primDensity_index)*Pminus2, massScalars); } if constexpr (is_eos_isothermal()) { @@ -859,7 +859,7 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu const double eint_R = x1RightState(i, j, k, pressure_index); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy @@ -879,7 +879,7 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu cs_L = quokka::EOS::ComputeSoundSpeed(rho_L, P_L, massScalars); E_L = quokka::EOS::ComputeEintFromPres(rho_L, P_L, massScalars) + ke_L; - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); cs_R = quokka::EOS::ComputeSoundSpeed(rho_R, P_R, massScalars); E_R = quokka::EOS::ComputeEintFromPres(rho_R, P_R, massScalars) + ke_R; } From 8c56cfcb9f2c1a8cbeafb85977c46abfff02a2b2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:06:43 +0000 Subject: [PATCH 34/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/hydro_system.hpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 64adf55a1..bfcd0d1b1 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -493,15 +493,19 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const if constexpr (reconstruct_eint) { // compute (rho e) (gamma - 1) amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 2, j, k); - Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), primVar(i + 2, j, k, primDensity_index)*Pplus2, massScalars); + Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), + primVar(i + 2, j, k, primDensity_index) * Pplus2, massScalars); massScalars = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); - Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), primVar(i + 1, j, k, primDensity_index)*Pplus1, massScalars); + Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), + primVar(i + 1, j, k, primDensity_index) * Pplus1, massScalars); massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); - P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index)*P, massScalars); + P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index) * P, massScalars); massScalars = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); - Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), primVar(i - 1, j, k, primDensity_index)*Pminus1, massScalars); + Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), + primVar(i - 1, j, k, primDensity_index) * Pminus1, massScalars); massScalars = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); - Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), primVar(i - 2, j, k, primDensity_index)*Pminus2, massScalars); + Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), + primVar(i - 2, j, k, primDensity_index) * Pminus2, massScalars); } if constexpr (is_eos_isothermal()) { From 6ec3096661cdc960dfdd09d5842f2d1eef26bed5 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 19:25:09 +0200 Subject: [PATCH 35/56] fix warning --- src/hydro_system.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index bfcd0d1b1..db1cce707 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -204,7 +204,6 @@ template auto HydroSystem::maxSignalSpeedLocal(a if constexpr (is_eos_isothermal()) { cs = cs_iso_; } else { - const auto P = ComputePressure(cons[bx], i, j, k); cs = ComputeSoundSpeed(cons[bx], i, j, k); } return {cs + abs_vel}; From a5d2714de9edbb4359ba3afb963a5dd2f8c83b8a Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 19:43:07 +0200 Subject: [PATCH 36/56] use microp EOS --- src/HydroContact/test_hydro_contact.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index 5576b7af3..464df5774 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -78,8 +78,8 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid( state_cc(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0.; state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0.; - state_cc(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx); - state_cc(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + state_cc(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + state_cc(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); }); } @@ -117,8 +117,8 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFa stateExact(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx); - stateExact(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + stateExact(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + stateExact(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); }); } From 8146b4eb2aa879fc50ca189db999f4e7c6a991e5 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 19:43:52 +0200 Subject: [PATCH 37/56] remove ASSERT --- src/hydro_system.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index db1cce707..4f9cac85b 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -233,7 +233,6 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4 0.); From 1fcab921855a737efab9245d05a2a05e75ab1d88 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 27 Jul 2023 17:44:16 +0000 Subject: [PATCH 38/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/HydroContact/test_hydro_contact.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index 464df5774..2d26adbbd 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -117,8 +117,9 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFa stateExact(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); - stateExact(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); + stateExact(i, j, k, HydroSystem::energy_index) = + quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + stateExact(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); }); } From 9dfb12ccb3da93f8615347dbbd0f80fbb5d1e998 Mon Sep 17 00:00:00 2001 From: psharda Date: Thu, 27 Jul 2023 19:54:20 +0200 Subject: [PATCH 39/56] remove unused var --- src/hydro_system.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 4f9cac85b..1cb82aca9 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -232,7 +232,6 @@ void HydroSystem::ComputeMaxSignalSpeed(amrex::Array4 0.); From 49c46c426e506b23b33c3f03a1cdabbf0eb33c9f Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 10:42:48 +0200 Subject: [PATCH 40/56] fix bug in EintFromPres --- src/EOS.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 5dd288163..4f16634f9 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -176,7 +176,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputeEintFromPre estate.p = Pressure; estate.mu = mean_molecular_weight_ / C::m_u; eos(eos_input_rp, estate); - Eint = estate.e * rho * boltzmann_constant_ / C::k_B; + Eint = estate.e * rho; } #endif return Eint; From f3b6ea0a76e5f44166dd22e038363606bb9f199f Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 10:52:37 +0200 Subject: [PATCH 41/56] declaring P_L as const double --- src/hydro_system.hpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 1cb82aca9..76628ff12 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -294,10 +294,8 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); - return P; } @@ -660,7 +658,6 @@ void HydroSystem::EnforceLimits(amrex::Real const densityFloor, amrex if (!HydroSystem::is_eos_isothermal()) { // recompute gas energy (to prevent P < 0) amrex::Real const Eint_star = Etot - 0.5 * rho_new * vsq; - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state[bx], i, j, k); amrex::Real const P_star = quokka::EOS::ComputePressure(rho_new, Eint_star, massScalars); amrex::Real P_new = P_star; @@ -859,9 +856,9 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); + const double P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); - P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); + const double P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); From 0bfb82eacb504ae511ae634bc4cd04e7d68aea65 Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 10:54:22 +0200 Subject: [PATCH 42/56] remove unused var --- src/HydroContact/test_hydro_contact.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index 2d26adbbd..ca63bb1d6 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -112,7 +112,6 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFa stateExact(i, j, k, n) = 0.; } - const auto gamma = quokka::EOS_Traits::gamma; stateExact(i, j, k, HydroSystem::density_index) = rho; stateExact(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; From a50ea5a6494bb284291c808dd910662f4af9327e Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 14:05:55 +0200 Subject: [PATCH 43/56] remove declaration of P_L --- src/hydro_system.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 76628ff12..d5e4b11db 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -856,9 +856,9 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - const double P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); + P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); - const double P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); + P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); From 13c623e697b697ac8fefca8f41f0083ad82b6e0a Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 15:54:59 +0200 Subject: [PATCH 44/56] updated microp: ensure exact cancellation --- extern/Microphysics | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extern/Microphysics b/extern/Microphysics index 45122a3f3..d6d91d423 160000 --- a/extern/Microphysics +++ b/extern/Microphysics @@ -1 +1 @@ -Subproject commit 45122a3f346da0f4ddd2bc9f5b5941bf8e14bcf1 +Subproject commit d6d91d423e88dc7a8ee9ede1b4101e1a4d312c88 From e07c5e4795cc16f0cd26fe680e2b7edaaaa0332f Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 16:06:19 +0200 Subject: [PATCH 45/56] remove unused var --- src/HydroContact/test_hydro_contact.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp index ca63bb1d6..a8146d602 100644 --- a/src/HydroContact/test_hydro_contact.cpp +++ b/src/HydroContact/test_hydro_contact.cpp @@ -70,7 +70,6 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid( AMREX_ASSERT(!std::isnan(rho)); AMREX_ASSERT(!std::isnan(P)); - const auto gamma = quokka::EOS_Traits::gamma; for (int n = 0; n < ncomp_cc; ++n) { state_cc(i, j, k, n) = 0.; } From dc4cb277cf424641b87350697b3c5d3a2bba1c03 Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 16:59:59 +0200 Subject: [PATCH 46/56] use finite mean_molecular_weight --- src/HydroBlast2D/test_hydro2d_blast.cpp | 2 +- src/HydroBlast3D/test_hydro3d_blast.cpp | 2 +- src/HydroKelvinHelmholz/test_hydro2d_kh.cpp | 2 +- src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp | 2 +- src/RayleighTaylor2D/test_hydro2d_rt.cpp | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/HydroBlast2D/test_hydro2d_blast.cpp b/src/HydroBlast2D/test_hydro2d_blast.cpp index e8a3988e4..094bce5e5 100644 --- a/src/HydroBlast2D/test_hydro2d_blast.cpp +++ b/src/HydroBlast2D/test_hydro2d_blast.cpp @@ -27,7 +27,7 @@ struct BlastProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/HydroBlast3D/test_hydro3d_blast.cpp b/src/HydroBlast3D/test_hydro3d_blast.cpp index 9a5e698a6..b698e8d5f 100644 --- a/src/HydroBlast3D/test_hydro3d_blast.cpp +++ b/src/HydroBlast3D/test_hydro3d_blast.cpp @@ -34,7 +34,7 @@ bool test_passes = false; // if one of the energy checks fails, set to false template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp index f561fe6f4..ac402cf59 100644 --- a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp +++ b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp @@ -24,7 +24,7 @@ struct KelvinHelmholzProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp b/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp index 05b1e5ae6..d2ba59bb3 100644 --- a/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp +++ b/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp @@ -23,7 +23,7 @@ struct RichtmeyerMeshkovProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; diff --git a/src/RayleighTaylor2D/test_hydro2d_rt.cpp b/src/RayleighTaylor2D/test_hydro2d_rt.cpp index 6255e9b8f..3800768a2 100644 --- a/src/RayleighTaylor2D/test_hydro2d_rt.cpp +++ b/src/RayleighTaylor2D/test_hydro2d_rt.cpp @@ -24,7 +24,7 @@ struct RTProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 1.4; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; From 73567968a2915e6122a286be9c4c8647877fc79a Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 17:02:31 +0200 Subject: [PATCH 47/56] use finite mol wt and use microphysics for EOS --- src/HydroQuirk/test_quirk.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/HydroQuirk/test_quirk.cpp b/src/HydroQuirk/test_quirk.cpp index 7ba3ac4b0..b5e532ba6 100644 --- a/src/HydroQuirk/test_quirk.cpp +++ b/src/HydroQuirk/test_quirk.cpp @@ -43,7 +43,7 @@ struct QuirkProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; @@ -125,8 +125,8 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(qu state_cc(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; state_cc(i, j, k, HydroSystem::x2Momentum_index) = rho * vy; state_cc(i, j, k, HydroSystem::x3Momentum_index) = rho * vz; - state_cc(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.) + 0.5 * rho * v_sq; - state_cc(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + state_cc(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * v_sq; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); }); } From 04c26d713d466beaf9ab3aa06bcd236170ada31a Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 21:11:44 +0200 Subject: [PATCH 48/56] use finite mean mol wt --- src/ShockCloud/cloud.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ShockCloud/cloud.cpp b/src/ShockCloud/cloud.cpp index a60435dd9..9fd6e54c7 100644 --- a/src/ShockCloud/cloud.cpp +++ b/src/ShockCloud/cloud.cpp @@ -39,7 +39,7 @@ constexpr double seconds_in_year = 3.154e7; template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; // default value - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; From 9dcf4912b2071180fef924be412dab321f9cb1d3 Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 21:39:25 +0200 Subject: [PATCH 49/56] use microphysics eos --- src/PassiveScalar/test_scalars.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/PassiveScalar/test_scalars.cpp b/src/PassiveScalar/test_scalars.cpp index 7f52c9859..6c8f94cce 100644 --- a/src/PassiveScalar/test_scalars.cpp +++ b/src/PassiveScalar/test_scalars.cpp @@ -78,8 +78,8 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(q state_cc(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0.; state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0.; - state_cc(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx); - state_cc(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + state_cc(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + state_cc(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); state_cc(i, j, k, HydroSystem::scalar0_index) = scalar; }); } @@ -120,8 +120,8 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFab stateExact(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.) + 0.5 * rho * (vx * vx); - stateExact(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + stateExact(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + stateExact(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); stateExact(i, j, k, HydroSystem::scalar0_index) = scalar; }); } @@ -166,7 +166,7 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFab const auto s = val_exact.at(HydroSystem::scalar0_index)[i]; const auto vx = xmom / rho; const auto Eint = E - 0.5 * rho * (vx * vx); - const auto P = (quokka::EOS_Traits::gamma - 1.) * Eint; + const auto P = quokka::EOS::ComputePressure(rho, Eint); d_exact.push_back(rho); vx_exact.push_back(vx); P_exact.push_back(P); @@ -180,7 +180,7 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFab const auto fs = values.at(HydroSystem::scalar0_index)[i]; const auto fvx = fxmom / frho; const auto fEint = fE - 0.5 * frho * (fvx * fvx); - const auto fP = (quokka::EOS_Traits::gamma - 1.) * fEint; + const auto fP = quokka::EOS::ComputePressure(frho, fEint); d_final.push_back(frho); vx_final.push_back(fvx); P_final.push_back(fP); From bed26e3d5d162ad964317dcaca81fd32666374c8 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 28 Jul 2023 19:39:51 +0000 Subject: [PATCH 50/56] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/PassiveScalar/test_scalars.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/PassiveScalar/test_scalars.cpp b/src/PassiveScalar/test_scalars.cpp index 6c8f94cce..a7a553b4a 100644 --- a/src/PassiveScalar/test_scalars.cpp +++ b/src/PassiveScalar/test_scalars.cpp @@ -120,7 +120,8 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFab stateExact(i, j, k, HydroSystem::x1Momentum_index) = rho * vx; stateExact(i, j, k, HydroSystem::x2Momentum_index) = 0.; stateExact(i, j, k, HydroSystem::x3Momentum_index) = 0.; - stateExact(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); + stateExact(i, j, k, HydroSystem::energy_index) = + quokka::EOS::ComputeEintFromPres(rho, P) + 0.5 * rho * (vx * vx); stateExact(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); stateExact(i, j, k, HydroSystem::scalar0_index) = scalar; }); From 16f441a91ab8a4dfa066dc0de1ac39eeeaa06b53 Mon Sep 17 00:00:00 2001 From: psharda Date: Fri, 28 Jul 2023 22:02:54 +0200 Subject: [PATCH 51/56] remove unused var --- src/PassiveScalar/test_scalars.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/PassiveScalar/test_scalars.cpp b/src/PassiveScalar/test_scalars.cpp index a7a553b4a..e768e1ef0 100644 --- a/src/PassiveScalar/test_scalars.cpp +++ b/src/PassiveScalar/test_scalars.cpp @@ -70,7 +70,6 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(q scalar = 0.0; } - const auto gamma = quokka::EOS_Traits::gamma; for (int n = 0; n < state_cc.nComp(); ++n) { state_cc(i, j, k, n) = 0.; } @@ -112,7 +111,6 @@ void RadhydroSimulation::computeReferenceSolution(amrex::MultiFab scalar = 0.0; } - const auto gamma = quokka::EOS_Traits::gamma; for (int n = 0; n < ncomp; ++n) { stateExact(i, j, k, n) = 0.; } From b7c3da0297dca2a124562b088d2d51a5c309ca5a Mon Sep 17 00:00:00 2001 From: psharda Date: Sat, 29 Jul 2023 14:56:58 +0200 Subject: [PATCH 52/56] if rho = 0, pres is 0 --- src/EOS.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/EOS.hpp b/src/EOS.hpp index 4f16634f9..e931217f2 100644 --- a/src/EOS.hpp +++ b/src/EOS.hpp @@ -251,7 +251,12 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS::ComputePressure(am if constexpr (gamma_ != 1.0) { chem_eos_t estate; estate.rho = rho; - estate.e = Eint / rho; + // if rho is 0, pass 0 to state.e + if (rho == 0.0) { + estate.e = 0; + } else { + estate.e = Eint / rho; + } estate.mu = mean_molecular_weight_ / C::m_u; eos(eos_input_re, estate); P = estate.p; From 6dd9326df563d6a92de6a010e382bde72e26c1a9 Mon Sep 17 00:00:00 2001 From: psharda Date: Sat, 29 Jul 2023 23:19:17 +0200 Subject: [PATCH 53/56] return finite pressure for isothermal gamma --- src/hydro_system.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index d5e4b11db..75d1ab976 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -295,7 +295,11 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); - amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + if constexpr (is_eos_isothermal()) { + amrex::Real P = rho * cs_iso_ * cs_iso_; + } else { + amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + } return P; } From 50631c5187a1ae324ba376c195ba9e6bf306fb69 Mon Sep 17 00:00:00 2001 From: psharda Date: Sat, 29 Jul 2023 23:34:33 +0200 Subject: [PATCH 54/56] @psharda return finite pressure for isothermal gamma --- src/hydro_system.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index 75d1ab976..c8831e19b 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -294,11 +294,13 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto HydroSystem::ComputePressure const auto vz = pz / rho; const auto kinetic_energy = 0.5 * rho * (vx * vx + vy * vy + vz * vz); const auto thermal_energy = E - kinetic_energy; - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + amrex::Real P = NAN; + if constexpr (is_eos_isothermal()) { - amrex::Real P = rho * cs_iso_ * cs_iso_; + P = rho * cs_iso_ * cs_iso_; } else { - amrex::Real P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(cons, i, j, k); + P = quokka::EOS::ComputePressure(rho, thermal_energy, massScalars); } return P; } From deac5f33931c8ec04e82ee5b112127cab7c2ff6d Mon Sep 17 00:00:00 2001 From: psharda Date: Sun, 30 Jul 2023 09:46:23 +0200 Subject: [PATCH 55/56] use finite mol wt --- src/SphericalCollapse/spherical_collapse.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/SphericalCollapse/spherical_collapse.cpp b/src/SphericalCollapse/spherical_collapse.cpp index 253f289b5..a26c25c8d 100644 --- a/src/SphericalCollapse/spherical_collapse.cpp +++ b/src/SphericalCollapse/spherical_collapse.cpp @@ -28,7 +28,7 @@ struct CollapseProblem { template <> struct quokka::EOS_Traits { static constexpr double gamma = 5. / 3.; - static constexpr double mean_molecular_weight = NAN; + static constexpr double mean_molecular_weight = 1.0; static constexpr double boltzmann_constant = C::k_B; }; @@ -76,13 +76,12 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid AMREX_ASSERT(!std::isnan(rho)); AMREX_ASSERT(!std::isnan(P)); - const amrex::Real gamma = quokka::EOS_Traits::gamma; state_cc(i, j, k, HydroSystem::density_index) = rho; state_cc(i, j, k, HydroSystem::x1Momentum_index) = 0; state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0; state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0; - state_cc(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.); - state_cc(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + state_cc(i, j, k, HydroSystem::energy_index) = quokka::EOS::ComputeEintFromPres(rho, P); + state_cc(i, j, k, HydroSystem::internalEnergy_index) = quokka::EOS::ComputeEintFromPres(rho, P); }); } From 77c35b1ed3a39c44798c445e6b363813ce1a0dea Mon Sep 17 00:00:00 2001 From: psharda Date: Sun, 30 Jul 2023 21:58:17 +0200 Subject: [PATCH 56/56] minor style changes --- src/hydro_system.hpp | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/hydro_system.hpp b/src/hydro_system.hpp index c8831e19b..4a56b78c9 100644 --- a/src/hydro_system.hpp +++ b/src/hydro_system.hpp @@ -493,20 +493,20 @@ void HydroSystem::ComputeFlatteningCoefficients(amrex::MultiFab const if constexpr (reconstruct_eint) { // compute (rho e) (gamma - 1) - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i + 2, j, k); + amrex::GpuArray massScalars_plus2 = RadSystem::ComputeMassScalars(primVar, i + 2, j, k); Pplus2 = quokka::EOS::ComputePressure(primVar(i + 2, j, k, primDensity_index), - primVar(i + 2, j, k, primDensity_index) * Pplus2, massScalars); - massScalars = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); + primVar(i + 2, j, k, primDensity_index) * Pplus2, massScalars_plus2); + amrex::GpuArray massScalars_plus1 = RadSystem::ComputeMassScalars(primVar, i + 1, j, k); Pplus1 = quokka::EOS::ComputePressure(primVar(i + 1, j, k, primDensity_index), - primVar(i + 1, j, k, primDensity_index) * Pplus1, massScalars); - massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); + primVar(i + 1, j, k, primDensity_index) * Pplus1, massScalars_plus1); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(primVar, i, j, k); P = quokka::EOS::ComputePressure(primVar(i, j, k, primDensity_index), primVar(i, j, k, primDensity_index) * P, massScalars); - massScalars = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); + amrex::GpuArray massScalars_minus1 = RadSystem::ComputeMassScalars(primVar, i - 1, j, k); Pminus1 = quokka::EOS::ComputePressure(primVar(i - 1, j, k, primDensity_index), - primVar(i - 1, j, k, primDensity_index) * Pminus1, massScalars); - massScalars = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); + primVar(i - 1, j, k, primDensity_index) * Pminus1, massScalars_minus1); + amrex::GpuArray massScalars_minus2 = RadSystem::ComputeMassScalars(primVar, i - 2, j, k); Pminus2 = quokka::EOS::ComputePressure(primVar(i - 2, j, k, primDensity_index), - primVar(i - 2, j, k, primDensity_index) * Pminus2, massScalars); + primVar(i - 2, j, k, primDensity_index) * Pminus2, massScalars_minus2); } if constexpr (is_eos_isothermal()) { @@ -861,10 +861,10 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu // (pressure_index is actually eint) const double eint_L = x1LeftState(i, j, k, pressure_index); const double eint_R = x1RightState(i, j, k, pressure_index); - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars); - massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); - P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars); + amrex::GpuArray massScalars_L = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); + P_L = quokka::EOS::ComputePressure(rho_L, eint_L * rho_L, massScalars_L); + amrex::GpuArray massScalars_R = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + P_R = quokka::EOS::ComputePressure(rho_R, eint_R * rho_R, massScalars_R); // auxiliary Eint is actually (auxiliary) specific internal energy Eint_L = rho_L * x1LeftState(i, j, k, primEint_index); @@ -879,13 +879,13 @@ void HydroSystem::ComputeFluxes(amrex::MultiFab &x1Flux_mf, amrex::Mu Eint_R = x1RightState(i, j, k, primEint_index); } - amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); - cs_L = quokka::EOS::ComputeSoundSpeed(rho_L, P_L, massScalars); - E_L = quokka::EOS::ComputeEintFromPres(rho_L, P_L, massScalars) + ke_L; + amrex::GpuArray massScalars_L = RadSystem::ComputeMassScalars(x1LeftState, i, j, k); + cs_L = quokka::EOS::ComputeSoundSpeed(rho_L, P_L, massScalars_L); + E_L = quokka::EOS::ComputeEintFromPres(rho_L, P_L, massScalars_L) + ke_L; - massScalars = RadSystem::ComputeMassScalars(x1RightState, i, j, k); - cs_R = quokka::EOS::ComputeSoundSpeed(rho_R, P_R, massScalars); - E_R = quokka::EOS::ComputeEintFromPres(rho_R, P_R, massScalars) + ke_R; + amrex::GpuArray massScalars_R = RadSystem::ComputeMassScalars(x1RightState, i, j, k); + cs_R = quokka::EOS::ComputeSoundSpeed(rho_R, P_R, massScalars_R); + E_R = quokka::EOS::ComputeEintFromPres(rho_R, P_R, massScalars_R) + ke_R; } AMREX_ASSERT(cs_L > 0.0);