diff --git a/src/NSCBC/channel.cpp b/src/NSCBC/channel.cpp index c607ab1fd..0e5461e89 100644 --- a/src/NSCBC/channel.cpp +++ b/src/NSCBC/channel.cpp @@ -21,6 +21,7 @@ #include "AMReX_MultiFab.H" #include "AMReX_ParallelContext.H" #include "AMReX_ParallelDescriptor.H" +#include "AMReX_ParmParse.H" #include "AMReX_REAL.H" #include "AMReX_SPACE.H" #include "AMReX_TableData.H" @@ -30,6 +31,7 @@ #include "HydroState.hpp" #include "RadhydroSimulation.hpp" #include "channel.hpp" +#include "fundamental_constants.H" #include "hydro_system.hpp" #include "physics_info.hpp" #include "physics_numVars.hpp" @@ -42,8 +44,8 @@ struct Channel { }; // dummy type to allow compile-type polymorphism via template specialization template <> struct quokka::EOS_Traits { - static constexpr double gamma = 5. / 3.; // default value - static constexpr double mean_molecular_weight = C::m_p + C::m_e; + static constexpr double gamma = 1.1; + static constexpr double mean_molecular_weight = 28.96 * C::m_u; // air static constexpr double boltzmann_constant = C::k_B; }; @@ -56,18 +58,17 @@ template <> struct Physics_Traits { static constexpr bool is_radiation_enabled = false; }; -constexpr Real Tgas0 = 300; // K -constexpr Real nH0 = 1.0; // cm^-3 -constexpr Real rho0 = nH0 * (C::m_p + C::m_e); // g cm^-3 - -// global variables needed for Dirichlet boundary condition +// global variables needed for Dirichlet boundary condition and initial conditions namespace { -AMREX_GPU_MANAGED Real u_inflow = 1.0e4; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) -AMREX_GPU_MANAGED Real v_inflow = 0; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) -AMREX_GPU_MANAGED Real w_inflow = 0; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) -AMREX_GPU_MANAGED Real P_inflow = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) -}; // namespace +Real rho0 = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +Real Tgas0 = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +Real u0 = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +AMREX_GPU_MANAGED Real u_inflow = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +AMREX_GPU_MANAGED Real v_inflow = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +AMREX_GPU_MANAGED Real w_inflow = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +AMREX_GPU_MANAGED Real P_inflow = NAN; // NOLINT(cppcoreguidelines-avoid-non-const-global-variables) +}; // namespace template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka::grid grid_elem) { @@ -77,11 +78,10 @@ template <> void RadhydroSimulation::setInitialConditionsOnGrid(quokka: amrex::GpuArray prob_hi = grid_elem.prob_hi_; const amrex::Box &indexRange = grid_elem.indexRange_; const amrex::Array4 &state_cc = grid_elem.array_; - amrex::Real const u_inflow = ::u_inflow; amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { Real const rho = rho0; - Real const xmom = 0; + Real const xmom = rho0 * u0; Real const ymom = 0; Real const zmom = 0; Real const Eint = quokka::EOS::ComputeEintFromTgas(rho0, Tgas0); @@ -146,10 +146,10 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_inflow_x1_lower(quokka::valarray< const Real c = quokka::EOS::ComputeSoundSpeed(rho, P); const Real M = std::sqrt(u * u + v * v + w * w) / c; - const Real eta_2 = -0.278; - const Real eta_3 = 1.0; - const Real eta_4 = 1.0; - const Real eta_5 = 0.278; + const Real eta_2 = 2.; //-0.278; + const Real eta_3 = 2.; //1.0; + const Real eta_4 = 2.; //1.0; + const Real eta_5 = 2.; //0.278; // see SymPy notebook for derivation quokka::valarray dQ_dx{}; @@ -181,12 +181,12 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void AMRSimulation::setCustomBounda { auto [i, j, k] = iv.dim3(); amrex::Box const &box = geom.Domain(); - const Real Lx = 1.0; const auto &domain_lo = box.loVect3d(); const auto &domain_hi = box.hiVect3d(); const int ilo = domain_lo[0]; const int ihi = domain_hi[0]; const Real dx = geom.CellSize(0); + const Real Lx = geom.prob_domain.length(0); const Real P_inflow = ::P_inflow; const Real u_inflow = ::u_inflow; @@ -285,8 +285,19 @@ auto problem_main() -> int RadhydroSimulation sim(BCs_cc); + amrex::ParmParse pp("channel"); + // initial condition parameters + pp.query("rho0", ::rho0); // initial density [g/cc] + pp.query("Tgas0", ::Tgas0); // initial temperature [K] + pp.query("u0", ::u0); // initial velocity [cm/s] + // boundary condition parameters + pp.query("u_inflow", ::u_inflow); // inflow velocity along x-axis [cm/s] + pp.query("v_inflow", ::v_inflow); // transverse inflow velocity (v_y) [cm/s] + pp.query("w_inflow", ::w_inflow); // transverse inflow velocity (v_z) [cm/s] + // compute derived parameters const Real Eint0 = quokka::EOS::ComputeEintFromTgas(rho0, Tgas0); ::P_inflow = quokka::EOS::ComputePressure(rho0, Eint0); + amrex::Print() << "Derived inflow pressure is " << ::P_inflow << " erg/cc.\n"; // Set initial conditions sim.setInitialConditions(); diff --git a/tests/NSCBC_Channel.in b/tests/NSCBC_Channel.in index d2975e6cc..ed637ad01 100644 --- a/tests/NSCBC_Channel.in +++ b/tests/NSCBC_Channel.in @@ -1,9 +1,9 @@ # ***************************************************************** # Problem size and geometry # ***************************************************************** -geometry.prob_lo = 0.0 0.0 0.0 -geometry.prob_hi = 1.0 0.25 0.25 -geometry.is_periodic = 0 1 1 +geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_hi = 100.0 12.5 12.5 +geometry.is_periodic = 0 1 1 # ***************************************************************** # VERBOSITY @@ -13,27 +13,33 @@ amr.v = 1 # verbosity in Amr # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 64 16 16 +amr.n_cell = 256 32 32 amr.max_level = 0 # number of levels = max_level + 1 -amr.blocking_factor_x = 64 -amr.blocking_factor_y = 16 -amr.blocking_factor_z = 16 # grid size must be divisible by this -amr.max_grid_size = 64 +amr.blocking_factor = 32 +amr.max_grid_size = 128 # ***************************************************************** # Quokka options # ***************************************************************** -cfl = 0.3 +cfl = 0.2 do_reflux = 1 do_subcycle = 1 -max_timesteps = 5000 -stop_time = 1.0e-4 +max_timesteps = 25000 +stop_time = 0.05 checkpoint_interval = 5000 -plotfile_interval = 20 +plotfile_interval = 100 ascent_interval = -1 hydro.rk_integrator_order = 2 hydro.reconstruction_order = 3 hydro.use_dual_energy = 0 hydro.low_level_debugging_output = 0 + +channel.rho0 = 1.1e-3 # g/cc +channel.Tgas0 = 320.8398927398333 # K +channel.u0 = 2.0e3 # cm/s + +channel.u_inflow = 4.0e3 # cm/s +channel.v_inflow = 0. # cm/s +channel.w_inflow = 0. # cm/s