diff --git a/Source/BC.H b/Source/BC.H index 09a7b3a4..8344747c 100644 --- a/Source/BC.H +++ b/Source/BC.H @@ -183,9 +183,9 @@ struct BCFill amrex::GpuArray flux_of_heat_flux = { rxx_eq, ryy_eq, rzz_eq, rxy_eq, rxz_eq, ryz_eq}; - set_extended_grad_expansion_generic( + data(iv, q) = set_extended_grad_expansion_generic( two_rho_e, heat_flux, flux_of_heat_flux, m_mesh_speed, m_weights[q], - ev, m_theta0, m_zero_vec, 1.0, data(iv, q)); + ev, m_theta0, m_zero_vec, 1.0); } AMREX_GPU_DEVICE @@ -236,8 +236,8 @@ struct BCFill const amrex::Real rho_bc = (rho_out + rho_tan) / (1.0 - ndir * vel_bc[idir]); - set_equilibrium_value( - rho_bc, vel_bc, RT, m_mesh_speed, m_weights[q], ev, data(iv, q)); + data(iv, q) = set_equilibrium_value( + rho_bc, vel_bc, RT, m_mesh_speed, m_weights[q], ev); rho_bc_out = rho_bc; } @@ -287,8 +287,8 @@ struct BCFill amrex::RealVect vel_bc(AMREX_D_DECL(0.0, 0.0, 0.0)); vel_bc[idir] = ndir * (1.0 - (rho_out + rho_tan) / rho_bc); - set_equilibrium_value( - rho_bc, vel_bc, RT, m_mesh_speed, m_weights[q], ev, data(iv, q)); + data(iv, q) = set_equilibrium_value( + rho_bc, vel_bc, RT, m_mesh_speed, m_weights[q], ev); rho_bc_out = rho_out + rho_tan; } @@ -340,9 +340,9 @@ struct BCFill amrex::GpuArray flux_of_heat_flux = { rxx_eq, ryy_eq, rzz_eq, rxy_eq, rxz_eq, ryz_eq}; - set_extended_grad_expansion_generic( + data(iv, q) = set_extended_grad_expansion_generic( two_rho_e, heat_flux, flux_of_heat_flux, m_mesh_speed, m_weights[q], - ev, m_theta0, m_zero_vec, 1.0, data(iv, q)); + ev, m_theta0, m_zero_vec, 1.0); } AMREX_GPU_DEVICE diff --git a/Source/IC.H b/Source/IC.H index 67b22e5b..f0dab1e2 100644 --- a/Source/IC.H +++ b/Source/IC.H @@ -39,7 +39,7 @@ struct Constant vel = velocity; } - AMREX_GPU_DEVICE void thermal_initialise( + AMREX_GPU_DEVICE void thermal_initialize( const amrex::IntVect& /*iv*/, amrex::GeometryData const& /*geom*/, amrex::Real& rho, @@ -117,7 +117,7 @@ struct TaylorGreen } AMREX_GPU_DEVICE - void thermal_initialise( + void thermal_initialize( const amrex::IntVect& iv, amrex::GeometryData const& geom, amrex::Real& rho, @@ -201,13 +201,12 @@ struct ViscosityTest vel[0] = velocity[0] + a0 * std::sin(2.0 * M_PI * y / wave_length); vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif } - AMREX_GPU_DEVICE - void thermal_initialise( + AMREX_GPU_DEVICE void thermal_initialize( const amrex::IntVect& iv, amrex::GeometryData const& geom, amrex::Real& rho, @@ -228,9 +227,9 @@ struct ViscosityTest vel[0] = velocity[0] + a0 * std::sin(2.0 * M_PI * y / wave_length); vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif temperature = initial_temperature; R = R_u / m_bar; @@ -285,13 +284,13 @@ struct ThermalDiffusivityTest vel[0] = velocity[0]; vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif } AMREX_GPU_DEVICE - void thermal_initialise( + void thermal_initialize( const amrex::IntVect& iv, amrex::GeometryData const& geom, amrex::Real& rho, @@ -314,9 +313,9 @@ struct ThermalDiffusivityTest vel[0] = velocity[0]; vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif gamma = adiabatic_exponent; @@ -370,17 +369,16 @@ struct SodTest vel[0] = velocity[0]; vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } - +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif rho = density + 0.5 * (1.0 + std::tanh((x - x_discontinuity) * 3.0)) * (density_ratio * density - density); } AMREX_GPU_DEVICE - void thermal_initialise( + void thermal_initialize( const amrex::IntVect& iv, amrex::GeometryData const& geom, amrex::Real& rho, @@ -399,9 +397,9 @@ struct SodTest vel[0] = velocity[0]; vel[1] = velocity[1]; - if (AMREX_SPACEDIM == 3) { - vel[2] = velocity[2]; - } +#if AMREX_SPACEDIM == 3 + vel[2] = velocity[2]; +#endif gamma = adiabatic_exponent; @@ -488,13 +486,13 @@ public: ic.R_u / constants::AIR_MOLAR_MASS; amrex::Real gamma = 1.667; - ic.thermal_initialise( + ic.thermal_initialize( iv, geom, rho, vel, temperature, specific_gas_constant, gamma); - set_equilibrium_value( + f_arrs[nbx](i, j, k, q) = set_equilibrium_value( rho, vel, specific_gas_constant * temperature, - l_mesh_speed, wt, ev, f_arrs[nbx](i, j, k, q)); + l_mesh_speed, wt, ev); amrex::Real cv = specific_gas_constant / (gamma - 1.0); @@ -514,9 +512,10 @@ public: amrex::GpuArray flux_of_heat_flux = { rxx_eq, ryy_eq, rzz_eq, rxy_eq, rxz_eq, ryz_eq}; - set_extended_grad_expansion_generic( - two_rho_e, heat_flux, flux_of_heat_flux, l_mesh_speed, - wt, ev, theta0, zero_vec, 1.0, g_arrs[nbx](i, j, k, q)); + g_arrs[nbx](i, j, k, q) = + set_extended_grad_expansion_generic( + two_rho_e, heat_flux, flux_of_heat_flux, + l_mesh_speed, wt, ev, theta0, zero_vec, 1.0); }); } else { amrex::Abort("Invalid model_type"); diff --git a/Source/LBM.cpp b/Source/LBM.cpp index 04285ae9..f6d7d116 100644 --- a/Source/LBM.cpp +++ b/Source/LBM.cpp @@ -669,34 +669,32 @@ void LBM::macrodata_to_equilibrium(const int lev) const auto& ev = evs[q]; - amrex::Real temperature = + const amrex::Real temperature = md_arr(iv, constants::TEMPERATURE_IDX); - amrex::Real omega = + const amrex::Real omega = 1.0 / (nu / (specific_gas_constant * temperature * dt) + 0.5); - amrex::Real omega_one = + const amrex::Real omega_one = 1.0 / (alpha / (specific_gas_constant * temperature * dt) + 0.5); - amrex::Real omega_one_by_omega = omega_one / omega; - amrex::Real omega_corr = (2.0 - omega) / (2.0 * omega * rho); + const amrex::Real omega_one_by_omega = omega_one / omega; + const amrex::Real omega_corr = + (2.0 - omega) / (2.0 * omega * rho); - amrex::Real pxx_ext = + const amrex::Real pxx_ext = vel[0] * vel[0] + specific_gas_constant * temperature + dt * (omega_corr)*d_arr(iv, constants::D_Q_CORR_X_IDX); - amrex::Real pyy_ext = + const amrex::Real pyy_ext = vel[1] * vel[1] + specific_gas_constant * temperature + dt * (omega_corr)*d_arr(iv, constants::D_Q_CORR_Y_IDX); - amrex::Real pzz_ext(0.0); - if (AMREX_SPACEDIM == 3) { - pzz_ext = - vel[2] * vel[2] + specific_gas_constant * temperature + - dt * (omega_corr)*d_arr(iv, constants::D_Q_CORR_Z_IDX); - } + const amrex::Real pzz_ext = AMREX_D_PICK( + 0.0, 0.0, + vel[2] * vel[2] + specific_gas_constant * temperature + + dt * (omega_corr)*d_arr(iv, constants::D_Q_CORR_Z_IDX)); - set_extended_equilibrium_value( - rho, vel, pxx_ext, pyy_ext, pzz_ext, l_mesh_speed, wt, ev, - eq_arr(iv, q)); + eq_arr(iv, q) = set_extended_equilibrium_value( + rho, vel, pxx_ext, pyy_ext, pzz_ext, l_mesh_speed, wt, ev); amrex::Real AMREX_D_DECL(qx_eq = 0.0, qy_eq = 0.0, qz_eq = 0.0); amrex::Real rxx_eq(0.0), ryy_eq(0.0), rzz_eq(0.0), rxy_eq(0.0), @@ -755,9 +753,9 @@ void LBM::macrodata_to_equilibrium(const int lev) amrex::GpuArray flux_of_heat_flux = { rxx_eq, ryy_eq, rzz_eq, rxy_eq, rxz_eq, ryz_eq}; - set_extended_grad_expansion_generic( + eq_arr_g(iv, q) = set_extended_grad_expansion_generic( two_rho_e, heat_flux_mrt, flux_of_heat_flux, l_mesh_speed, - wt, ev, theta0, zero_vec, 1.0, eq_arr_g(iv, q)); + wt, ev, theta0, zero_vec, 1.0); } }); amrex::Gpu::synchronize(); @@ -849,16 +847,12 @@ void LBM::f_to_macrodata(const int lev) pxx += ev[0] * ev[0] * f_arr(iv, q); pyy += ev[1] * ev[1] * f_arr(iv, q); - if (AMREX_SPACEDIM == 3) { - pzz += ev[2] * ev[2] * f_arr(iv, q); - } pxy += ev[0] * ev[1] * f_arr(iv, q); - if (AMREX_SPACEDIM == 3) { - pxz += ev[0] * ev[2] * f_arr(iv, q); - } - if (AMREX_SPACEDIM == 3) { - pyz += ev[1] * ev[2] * f_arr(iv, q); - } +#if AMREX_SPACEDIM == 3 + pzz += ev[2] * ev[2] * f_arr(iv, q); + pxz += ev[0] * ev[2] * f_arr(iv, q); + pyz += ev[1] * ev[2] * f_arr(iv, q); +#endif two_rho_e += g_arr(iv, q); @@ -982,22 +976,15 @@ void LBM::compute_q_corrections(const int lev) const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); if (if_arr(iv, 0) == 1) { - const amrex::Real d_qxxx = gradient( + d_arr(iv, constants::D_Q_CORR_X_IDX) = gradient( 0, constants::Q_CORR_X_IDX, iv, idx, dbox, if_arr, md_arr); - const amrex::Real d_qyyy = gradient( + d_arr(iv, constants::D_Q_CORR_Y_IDX) = gradient( 1, constants::Q_CORR_Y_IDX, iv, idx, dbox, if_arr, md_arr); - amrex::Real d_qzzz = 0.0; - - if (AMREX_SPACEDIM == 3) { - d_qzzz = gradient( - 2, constants::Q_CORR_Z_IDX, iv, idx, dbox, if_arr, - md_arr); - } - - d_arr(iv, constants::D_Q_CORR_X_IDX) = d_qxxx; - d_arr(iv, constants::D_Q_CORR_Y_IDX) = d_qyyy; - d_arr(iv, constants::D_Q_CORR_Z_IDX) = d_qzzz; +#if AMREX_SPACEDIM == 3 + d_arr(iv, constants::D_Q_CORR_Z_IDX) = gradient( + 2, constants::Q_CORR_Z_IDX, iv, idx, dbox, if_arr, md_arr); +#endif } }); amrex::Gpu::synchronize(); @@ -1302,8 +1289,8 @@ void LBM::fill_f_inside_eb(const int lev) [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int q) noexcept { if (is_fluid_arrs[nbx](i, j, k, 0) == 0) { - set_population_zero(f_arrs[nbx](i, j, k, q)); - set_population_zero(g_arrs[nbx](i, j, k, q)); + f_arrs[nbx](i, j, k, q) = 0.0; + g_arrs[nbx](i, j, k, q) = 0.0; } }); diff --git a/Source/Utilities.H b/Source/Utilities.H index 8c56ac50..f3e7514e 100644 --- a/Source/Utilities.H +++ b/Source/Utilities.H @@ -8,14 +8,13 @@ #include "Constants.H" namespace lbm { -AMREX_GPU_DEVICE AMREX_INLINE void set_equilibrium_value( +AMREX_GPU_DEVICE AMREX_INLINE amrex::Real set_equilibrium_value( const amrex::Real& rho, const amrex::RealVect& vel, const amrex::Real& r_temperature, const amrex::Real& /*mesh_speed*/, const amrex::Real& /*wt*/, - const amrex::IntVect& ev, - amrex::Real& eq) + const amrex::IntVect& ev) { // Warning! Achtung! Strictly for D3Q27 only // This entropic equilibrium is very stable. Relaxes all 27 moments. @@ -26,50 +25,44 @@ AMREX_GPU_DEVICE AMREX_INLINE void set_equilibrium_value( const amrex::Real pxx = vel[0] * vel[0] + r_temperature; const amrex::Real pyy = vel[1] * vel[1] + r_temperature; - amrex::Real pzz = 0.0; - if (AMREX_SPACEDIM == 3) { - pzz = vel[2] * vel[2] + r_temperature; - } const amrex::Real phix = ev[0] * 0.5 * vel[0] + std::abs(ev[0]) * (1.50 * pxx - 1.0) - pxx + 1.0; const amrex::Real phiy = ev[1] * 0.5 * vel[1] + std::abs(ev[1]) * (1.50 * pyy - 1.0) - pyy + 1.0; - amrex::Real phiz = 1.0; - if (AMREX_SPACEDIM == 3) { - phiz = ev[2] * 0.5 * vel[2] + std::abs(ev[2]) * (1.50 * pzz - 1.0) - - pzz + 1.0; - } +#if AMREX_SPACEDIM == 3 + const amrex::Real pzz = vel[2] * vel[2] + r_temperature; + const amrex::Real phiz = + ev[2] * 0.5 * vel[2] + std::abs(ev[2]) * (1.50 * pzz - 1.0) - pzz + 1.0; +#endif - eq = rho * phix * phiy * phiz; + return rho * AMREX_D_TERM(phix, *phiy, *phiz); } -AMREX_GPU_DEVICE AMREX_INLINE void set_extended_equilibrium_value( +AMREX_GPU_DEVICE AMREX_INLINE amrex::Real set_extended_equilibrium_value( const amrex::Real& rho, const amrex::RealVect& vel, const amrex::Real& pxx, const amrex::Real& pyy, - const amrex::Real& pzz, + const amrex::Real& AMREX_D_PICK(, /*pzz*/, pzz), const amrex::Real& /*mesh_speed*/, const amrex::Real& /*wt*/, - const amrex::IntVect& ev, - amrex::Real& eq) + const amrex::IntVect& ev) { // Warning! Achtung! Strictly for D3Q27 only const amrex::Real phix = ev[0] * 0.5 * vel[0] + std::abs(ev[0]) * (1.50 * pxx - 1.0) - pxx + 1.0; const amrex::Real phiy = ev[1] * 0.5 * vel[1] + std::abs(ev[1]) * (1.50 * pyy - 1.0) - pyy + 1.0; - amrex::Real phiz = 1.0; - if (AMREX_SPACEDIM == 3) { - phiz = ev[2] * 0.5 * vel[2] + std::abs(ev[2]) * (1.50 * pzz - 1.0) - - pzz + 1.0; - } +#if AMREX_SPACEDIM == 3 + const amrex::Real phiz = + ev[2] * 0.5 * vel[2] + std::abs(ev[2]) * (1.50 * pzz - 1.0) - pzz + 1.0; +#endif - eq = rho * phix * phiy * phiz; + return rho * AMREX_D_TERM(phix, *phiy, *phiz); } -AMREX_GPU_DEVICE AMREX_INLINE void set_extended_grad_expansion_generic( +AMREX_GPU_DEVICE AMREX_INLINE amrex::Real set_extended_grad_expansion_generic( const amrex::Real& rho, const amrex::RealVect& momentum, const amrex::GpuArray& flux_of_momentum, @@ -84,8 +77,7 @@ AMREX_GPU_DEVICE AMREX_INLINE void set_extended_grad_expansion_generic( const amrex::IntVect& ev, const amrex::Real& theta0, const amrex::RealVect& frame_velocity, - const amrex::Real& s, - amrex::Real& f) + const amrex::Real& s) { // Note: This is lattice independent, as long as the correct lattice // temperature theta0 is passed. theta0 = \sum_i w_i cx_i * cx_i @@ -152,7 +144,7 @@ AMREX_GPU_DEVICE AMREX_INLINE void set_extended_grad_expansion_generic( uy * ((jz)-rho * uz) - uz * ((jy)-rho * uy)) * one_bystheta0 * one_bystheta0;); - f = rho + AMREX_D_TERM(a1x * ev[0], +a1y * ev[1], +a1z * ev[2]); + amrex::Real f = rho + AMREX_D_TERM(a1x * ev[0], +a1y * ev[1], +a1z * ev[2]); f += 0.5 * @@ -166,11 +158,8 @@ AMREX_GPU_DEVICE AMREX_INLINE void set_extended_grad_expansion_generic( 2.0 * (ev[1] * ev[2] - 0) * a2yz)); f *= wt; -} -AMREX_GPU_DEVICE AMREX_INLINE void set_population_zero(amrex::Real& eq) -{ - eq = 0.0; + return f; } AMREX_GPU_DEVICE AMREX_INLINE amrex::Real get_energy( @@ -262,10 +251,10 @@ AMREX_GPU_DEVICE AMREX_INLINE void get_equilibrium_moments( amrex::RealVect& qEq, amrex::Real& RxxEq, amrex::Real& RyyEq, - amrex::Real& RzzEq, + amrex::Real& AMREX_D_PICK(, /*RzzEq*/, RzzEq), amrex::Real& RxyEq, - amrex::Real& RxzEq, - amrex::Real& RyzEq) + amrex::Real& AMREX_D_PICK(, /*RxzEq*/, RxzEq), + amrex::Real& AMREX_D_PICK(, /*RyzEq*/, RyzEq)) { amrex::Real energy = total_energy / (2.0 * rho); @@ -275,22 +264,16 @@ AMREX_GPU_DEVICE AMREX_INLINE void get_equilibrium_moments( qEq[0] = 2.0 * rho * vel[0] * h; qEq[1] = 2.0 * rho * vel[1] * h; - if (AMREX_SPACEDIM == 3) { - qEq[2] = 2.0 * rho * vel[2] * h; - } RxxEq = 2.0 * rho * vel[0] * vel[0] * (h + (p / rho)) + 2.0 * p * h; RyyEq = 2.0 * rho * vel[1] * vel[1] * (h + (p / rho)) + 2.0 * p * h; - if (AMREX_SPACEDIM == 3) { - RzzEq = 2.0 * rho * vel[2] * vel[2] * (h + (p / rho)) + 2.0 * p * h; - } RxyEq = 2.0 * rho * vel[0] * vel[1] * (h + (p / rho)) + 0; - if (AMREX_SPACEDIM == 3) { - RxzEq = 2.0 * rho * vel[0] * vel[2] * (h + (p / rho)) + 0; - } - if (AMREX_SPACEDIM == 3) { - RyzEq = 2.0 * rho * vel[1] * vel[2] * (h + (p / rho)) + 0; - } +#if AMREX_SPACEDIM == 3 + qEq[2] = 2.0 * rho * vel[2] * h; + RzzEq = 2.0 * rho * vel[2] * vel[2] * (h + (p / rho)) + 2.0 * p * h; + RxzEq = 2.0 * rho * vel[0] * vel[2] * (h + (p / rho)) + 0; + RyzEq = 2.0 * rho * vel[1] * vel[2] * (h + (p / rho)) + 0; +#endif } AMREX_GPU_DEVICE AMREX_INLINE amrex::Real gradient(