From 5e7d611d233514c590e4aed039ee01c416f59bc5 Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Wed, 21 Jun 2023 09:39:57 -0600 Subject: [PATCH 01/16] Fix box type in abl boundary plane input for MAC vel (#860) --- amr-wind/wind_energy/ABLBoundaryPlane.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.cpp b/amr-wind/wind_energy/ABLBoundaryPlane.cpp index af48493482..b8beb70398 100644 --- a/amr-wind/wind_energy/ABLBoundaryPlane.cpp +++ b/amr-wind/wind_energy/ABLBoundaryPlane.cpp @@ -894,6 +894,11 @@ void ABLBoundaryPlane::populate_data( amrex::Abort("No inflow data at this level."); } + if (ori.isHigh()) { + amrex::Warning( + "We typically don't inflow boundary planes on the high side."); + } + const size_t nc = mfab.nComp(); #ifdef AMREX_USE_OMP @@ -902,7 +907,10 @@ void ABLBoundaryPlane::populate_data( for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const auto& sbx = mfi.growntilebox(1); + auto sbx = mfi.growntilebox(1); + if (!sbx.cellCentered()) { + sbx.enclosedCells(); + } const auto& src = m_in_data.interpolate_data(ori, lev); const auto& bx = sbx & src.box(); if (bx.isEmpty()) { From b41395c2f906b294eb3d73119e6f8f44212ba0cf Mon Sep 17 00:00:00 2001 From: "Georgios (Yorgos) Deskos" Date: Wed, 21 Jun 2023 13:51:32 -0600 Subject: [PATCH 02/16] init tke profile (#864) --- amr-wind/wind_energy/ABLFieldInit.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/amr-wind/wind_energy/ABLFieldInit.cpp b/amr-wind/wind_energy/ABLFieldInit.cpp index b74b981278..cd696d9d61 100644 --- a/amr-wind/wind_energy/ABLFieldInit.cpp +++ b/amr-wind/wind_energy/ABLFieldInit.cpp @@ -203,7 +203,7 @@ void ABLFieldInit::perturb_temperature( void ABLFieldInit::init_tke( const amrex::Geometry& geom, amrex::MultiFab& tke_fab) const { - tke_fab.setVal(m_tke_init, 1); + tke_fab.setVal(m_tke_init); if (!m_tke_init_profile) { return; @@ -219,18 +219,19 @@ void ABLFieldInit::init_tke( #endif for (amrex::MFIter mfi(tke_fab, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const auto& bx = mfi.growntilebox(1); + const auto& bx = mfi.tilebox(); const auto& tke = tke_fab.array(mfi); // Profile definition from Beare et al. (2006) amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + if (z < tke_cutoff_height) { tke(i, j, k) = tke_init_factor * - std::pow(1. - z / tke_cutoff_height, 3.0); + std::pow(1. - z / tke_cutoff_height, 3); } else { - tke(i, j, k) = 0.; + tke(i, j, k) = 1.e-20; } }); } From 5ae35331c43715818366a0cf2c89424095774b81 Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Thu, 22 Jun 2023 13:38:29 -0600 Subject: [PATCH 03/16] Fix for fill abl inflow (#862) --- amr-wind/wind_energy/ABLFillInflow.H | 2 +- amr-wind/wind_energy/ABLFillInflow.cpp | 12 +++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/amr-wind/wind_energy/ABLFillInflow.H b/amr-wind/wind_energy/ABLFillInflow.H index 8edfe3951f..6172e8a9cd 100644 --- a/amr-wind/wind_energy/ABLFillInflow.H +++ b/amr-wind/wind_energy/ABLFillInflow.H @@ -15,7 +15,7 @@ namespace amr_wind { * * \sa ABLBoundaryPlane */ -class ABLFillInflow : public FieldFillPatchOps +class ABLFillInflow : public FieldFillPatchOps { public: ABLFillInflow( diff --git a/amr-wind/wind_energy/ABLFillInflow.cpp b/amr-wind/wind_energy/ABLFillInflow.cpp index 588367fe74..bb5c85e4a8 100644 --- a/amr-wind/wind_energy/ABLFillInflow.cpp +++ b/amr-wind/wind_energy/ABLFillInflow.cpp @@ -7,7 +7,7 @@ ABLFillInflow::ABLFillInflow( const amrex::AmrCore& mesh, const SimTime& time, const ABLBoundaryPlane& bndry_plane) - : FieldFillPatchOps( + : FieldFillPatchOps( field, mesh, time, FieldInterpolator::CellConsLinear) , m_bndry_plane(bndry_plane) {} @@ -21,7 +21,8 @@ void ABLFillInflow::fillpatch( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillpatch(lev, time, mfab, nghost, fstate); + FieldFillPatchOps::fillpatch( + lev, time, mfab, nghost, fstate); m_bndry_plane.populate_data(lev, time, m_field, mfab); } @@ -33,7 +34,7 @@ void ABLFillInflow::fillpatch_from_coarse( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillpatch_from_coarse( + FieldFillPatchOps::fillpatch_from_coarse( lev, time, mfab, nghost, fstate); m_bndry_plane.populate_data(lev, time, m_field, mfab); @@ -46,7 +47,8 @@ void ABLFillInflow::fillphysbc( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillphysbc(lev, time, mfab, nghost, fstate); + FieldFillPatchOps::fillphysbc( + lev, time, mfab, nghost, fstate); m_bndry_plane.populate_data(lev, time, m_field, mfab); } @@ -87,7 +89,7 @@ void ABLFillInflow::fillpatch_sibling_fields( } } - FieldFillPatchOps::fillpatch_sibling_fields( + FieldFillPatchOps::fillpatch_sibling_fields( lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, fstate, itype); for (int i = 0; i < static_cast(mfabs.size()); i++) { From ffd72051dc1a9bf5eaa618435da2f31fcba621cc Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Thu, 22 Jun 2023 14:53:05 -0600 Subject: [PATCH 04/16] Fix MPL inflow (#865) --- amr-wind/wind_energy/ABLFillMPL.H | 2 +- amr-wind/wind_energy/ABLFillMPL.cpp | 12 +++++++----- amr-wind/wind_energy/ABLModulatedPowerLaw.cpp | 10 ++++++++-- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/amr-wind/wind_energy/ABLFillMPL.H b/amr-wind/wind_energy/ABLFillMPL.H index 9da6783823..6c9eea8aa2 100644 --- a/amr-wind/wind_energy/ABLFillMPL.H +++ b/amr-wind/wind_energy/ABLFillMPL.H @@ -15,7 +15,7 @@ namespace amr_wind { * * \sa ABLBoundaryPlane */ -class ABLFillMPL : public FieldFillPatchOps +class ABLFillMPL : public FieldFillPatchOps { public: ABLFillMPL( diff --git a/amr-wind/wind_energy/ABLFillMPL.cpp b/amr-wind/wind_energy/ABLFillMPL.cpp index fe6fd18a55..661c9b4e65 100644 --- a/amr-wind/wind_energy/ABLFillMPL.cpp +++ b/amr-wind/wind_energy/ABLFillMPL.cpp @@ -7,7 +7,7 @@ ABLFillMPL::ABLFillMPL( const amrex::AmrCore& mesh, const SimTime& time, const ABLModulatedPowerLaw& abl_mpl) - : FieldFillPatchOps( + : FieldFillPatchOps( field, mesh, time, FieldInterpolator::CellConsLinear) , m_abl_mpl(abl_mpl) {} @@ -21,7 +21,8 @@ void ABLFillMPL::fillpatch( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillpatch(lev, time, mfab, nghost, fstate); + FieldFillPatchOps::fillpatch( + lev, time, mfab, nghost, fstate); if (m_field.base_name() == "velocity") { m_abl_mpl.set_velocity(lev, time, m_field, mfab); @@ -37,7 +38,7 @@ void ABLFillMPL::fillpatch_from_coarse( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillpatch_from_coarse( + FieldFillPatchOps::fillpatch_from_coarse( lev, time, mfab, nghost, fstate); if (m_field.base_name() == "velocity") { @@ -54,7 +55,8 @@ void ABLFillMPL::fillphysbc( const amrex::IntVect& nghost, const FieldState fstate) { - FieldFillPatchOps::fillphysbc(lev, time, mfab, nghost, fstate); + FieldFillPatchOps::fillphysbc( + lev, time, mfab, nghost, fstate); if (m_field.base_name() == "velocity") { m_abl_mpl.set_velocity(lev, time, m_field, mfab); @@ -100,7 +102,7 @@ void ABLFillMPL::fillpatch_sibling_fields( } } - FieldFillPatchOps::fillpatch_sibling_fields( + FieldFillPatchOps::fillpatch_sibling_fields( lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, fstate, itype); for (int i = 0; i < static_cast(mfabs.size()); i++) { diff --git a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp index a83c893831..a5b0a21dba 100644 --- a/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp +++ b/amr-wind/wind_energy/ABLModulatedPowerLaw.cpp @@ -173,7 +173,10 @@ void ABLModulatedPowerLaw::set_velocity( : amrex::adjCellHi(domain, idir, nghost); for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) { - const auto& gbx = amrex::grow(mfi.validbox(), nghost); + auto gbx = amrex::grow(mfi.validbox(), nghost); + if (!gbx.cellCentered()) { + gbx.enclosedCells(); + } const auto& bx = gbx & dbx; if (!bx.ok()) { continue; @@ -247,7 +250,10 @@ void ABLModulatedPowerLaw::set_temperature( : amrex::adjCellHi(domain, idir, nghost); for (amrex::MFIter mfi(mfab); mfi.isValid(); ++mfi) { - const auto& gbx = amrex::grow(mfi.validbox(), nghost); + auto gbx = amrex::grow(mfi.validbox(), nghost); + if (!gbx.cellCentered()) { + gbx.enclosedCells(); + } const auto& bx = gbx & dbx; if (!bx.ok()) { continue; From 4b71037218723e0c63d54c140423ef503ac3c912 Mon Sep 17 00:00:00 2001 From: prakash <120606615+moprak-nrel@users.noreply.github.com> Date: Mon, 26 Jun 2023 13:26:58 -0600 Subject: [PATCH 05/16] add implementaiton for a log-law based wall model (#847) * add implementaiton for a log-law based wall model * make clang-format happy * make LogLaw header only and add gpu fix * linter fix * Add unit test and remove print for gpu * remove .h from cmakelist and fix gpu * Add a Schumann type model too * change default constructor for loglaw, and fix bug in utau_mean initialization * add a base amd model that can work without ABL physics to test the log-law wall-model for channel flow * Remove default constructors * Add the newly added wall-models to the doc, along with documentation for the base AMD model and Smagorisnky model * allow the specification of a reference index for the log-law instead of just the first cell * Change to a while loop, and abort if u_tau hasn't converged * make clang happy * formatting change * fix breaking changes from pr #842 * gpu fix for for the notherm model, needed to use CellSizeArray, and a GPUArray for grid spacing --------- Co-authored-by: Marc T. Henry de Frahan --- .../boundary_conditions/wall_models/LogLaw.H | 59 ++++++++ .../wall_models/ShearStressSimple.H | 40 +++++ .../wall_models/WallFunction.H | 26 +++- .../wall_models/WallFunction.cpp | 137 +++++++++++++++++- amr-wind/turbulence/LES/AMDNoTherm.H | 84 +++++++++++ amr-wind/turbulence/LES/AMDNoTherm.cpp | 80 ++++++++++ amr-wind/turbulence/LES/CMakeLists.txt | 1 + docs/sphinx/theory/theory.rst | 135 ++++++++++++++++- unit_tests/CMakeLists.txt | 1 + unit_tests/boundary_conditions/CMakeLists.txt | 4 + .../boundary_conditions/test_log_law.cpp | 47 ++++++ unit_tests/turbulence/test_turbulence_LES.cpp | 98 +++++++++++++ 12 files changed, 698 insertions(+), 14 deletions(-) create mode 100644 amr-wind/boundary_conditions/wall_models/LogLaw.H create mode 100644 amr-wind/boundary_conditions/wall_models/ShearStressSimple.H create mode 100644 amr-wind/turbulence/LES/AMDNoTherm.H create mode 100644 amr-wind/turbulence/LES/AMDNoTherm.cpp create mode 100644 unit_tests/boundary_conditions/CMakeLists.txt create mode 100644 unit_tests/boundary_conditions/test_log_law.cpp diff --git a/amr-wind/boundary_conditions/wall_models/LogLaw.H b/amr-wind/boundary_conditions/wall_models/LogLaw.H new file mode 100644 index 0000000000..dc02629da3 --- /dev/null +++ b/amr-wind/boundary_conditions/wall_models/LogLaw.H @@ -0,0 +1,59 @@ +#ifndef LogLaw_H +#define LogLaw_H + +#include "AMReX_AmrCore.H" +#include +namespace amr_wind { +struct LogLaw +{ + /* + * A simple wall model that sets the wall-shear stress + * based on computing u_tau given the horizontal velocity + * magnitude at a zref. This is akin to an explicit non-linear + * Robin boundary condition at the wall. + */ + + // Log law constants from Lee & Moser 2015 + // https://doi.org/10.1017/jfm.2015.268. + amrex::Real B{4.27}; + amrex::Real kappa{0.384}; + int max_iters = 25; // Max iterations for u_tau Newton-Raphson solve + // Reference height for log law + amrex::Real zref; + int ref_index{0}; + amrex::Real nu; // molecular viscosity + // u_tau state variable, gets updated in update_utau depending on + // the type of wall model used + amrex::Real utau_mean{1.0}; + amrex::Real wspd_mean; // mean horizontal velocity magnitude + + void update_utau_mean() { utau_mean = get_utau(wspd_mean); } + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real + get_utau(amrex::Real wspd) const + { + amrex::Real utau_iter = -1; + amrex::Real wspd_pred; + amrex::Real wspd_deriv; + amrex::Real zplus; + amrex::Real utau = utau_mean; + int iter = 0; + while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters) { + utau_iter = utau; + zplus = zref * utau / nu; + // Get wspd for a given utau from log-law + wspd_pred = utau * (std::log(zplus) / kappa + B); + wspd_deriv = (1 + std::log(zplus)) / kappa + B; // d(wspd)/d(utau) + utau = + utau - (wspd_pred - wspd) / wspd_deriv; // Newton-Raphson update + ++iter; + } + if (iter == max_iters) { + amrex::Abort(); + } + return utau; + } +}; +} /* namespace amr_wind */ + +#endif /* LogLaw_H */ diff --git a/amr-wind/boundary_conditions/wall_models/ShearStressSimple.H b/amr-wind/boundary_conditions/wall_models/ShearStressSimple.H new file mode 100644 index 0000000000..679efdc26e --- /dev/null +++ b/amr-wind/boundary_conditions/wall_models/ShearStressSimple.H @@ -0,0 +1,40 @@ +#ifndef SHEARSTRESSSIMPLE_H +#define SHEARSTRESSSIMPLE_H + +#include "amr-wind/boundary_conditions/wall_models/LogLaw.H" +#include "amr-wind/wind_energy/ShearStress.H" + +namespace amr_wind { + +struct SimpleShearSchumann +{ + explicit SimpleShearSchumann(const amr_wind::LogLaw& ll) + : utau2(ll.utau_mean * ll.utau_mean), wspd_mean(ll.wspd_mean), m_ll(ll) + {} + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real + get_shear(amrex::Real u, amrex::Real /* wspd */) const + { + return u / wspd_mean * utau2; + }; + + amrex::Real utau2; + amrex::Real wspd_mean; + const amr_wind::LogLaw m_ll; +}; +struct SimpleShearLogLaw +{ + explicit SimpleShearLogLaw(const amr_wind::LogLaw& ll) : m_ll(ll) {} + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real + get_shear(amrex::Real u, amrex::Real wspd) const + { + amrex::Real utau = m_ll.get_utau(wspd); + return utau * utau * u / wspd; + }; + + const amr_wind::LogLaw m_ll; +}; +} // namespace amr_wind + +#endif /* SHEARSTRESSSIMPLE_H */ diff --git a/amr-wind/boundary_conditions/wall_models/WallFunction.H b/amr-wind/boundary_conditions/wall_models/WallFunction.H index 5247c285a2..fff1dd1082 100644 --- a/amr-wind/boundary_conditions/wall_models/WallFunction.H +++ b/amr-wind/boundary_conditions/wall_models/WallFunction.H @@ -3,6 +3,8 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/core/FieldBCOps.H" +#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/boundary_conditions/wall_models/LogLaw.H" namespace amr_wind { @@ -16,9 +18,14 @@ namespace amr_wind { class WallFunction { public: - explicit WallFunction(const CFDSim& sim); + explicit WallFunction(CFDSim& sim); - amrex::Real utau() const { return m_utau; } + amrex::Real utau() const { return m_log_law.utau_mean; } + LogLaw log_law() const { return m_log_law; } + + //! Update the mean velocity at a given timestep + void update_umean(); + void update_utau_mean(); ~WallFunction() = default; @@ -27,7 +34,10 @@ private: const amrex::AmrCore& m_mesh; - amrex::Real m_utau{0.0}; + //! LogLaw instance + LogLaw m_log_law; + int m_direction{2}; ///< Direction normal to wall, hardcoded to z + VelPlaneAveragingFine m_pa_vel; }; /** Applies a shear-stress value at the domain boundary @@ -38,15 +48,21 @@ private: class VelWallFunc : public FieldBCIface { public: - VelWallFunc(Field& velocity, const WallFunction& wall_func); + VelWallFunc(Field& velocity, WallFunction& wall_func); void operator()(Field& velocity, const FieldState rho_state) override; static void wall_model( Field& velocity, const FieldState rho_state, const amrex::Real utau); + template + static void wall_model( + Field& velocity, + const FieldState rho_state, + const ShearStressSimple& tau); + private: - const WallFunction& m_wall_func; + WallFunction& m_wall_func; std::string m_wall_shear_stress_type{"constant"}; }; } // namespace amr_wind diff --git a/amr-wind/boundary_conditions/wall_models/WallFunction.cpp b/amr-wind/boundary_conditions/wall_models/WallFunction.cpp index d33a2298a1..342e49b314 100644 --- a/amr-wind/boundary_conditions/wall_models/WallFunction.cpp +++ b/amr-wind/boundary_conditions/wall_models/WallFunction.cpp @@ -1,4 +1,5 @@ #include "amr-wind/boundary_conditions/wall_models/WallFunction.H" +#include "amr-wind/boundary_conditions/wall_models/ShearStressSimple.H" #include "amr-wind/utilities/tensor_ops.H" #include "amr-wind/utilities/trig_ops.H" #include "amr-wind/diffusion/diffusion.H" @@ -11,30 +12,46 @@ namespace amr_wind { -WallFunction::WallFunction(const CFDSim& sim) : m_sim(sim), m_mesh(m_sim.mesh()) +WallFunction::WallFunction(CFDSim& sim) + : m_sim(sim), m_mesh(m_sim.mesh()), m_pa_vel(sim, m_direction) { { amrex::ParmParse pp("BodyForce"); amrex::Vector body_force{{0.0, 0.0, 0.0}}; pp.getarr("magnitude", body_force); - m_utau = std::sqrt(body_force[0]); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - std::abs(body_force[1]) < 1e-16, - "body force in y should be zero for this wall function"); + m_log_law.utau_mean = std::sqrt(std::sqrt( + body_force[0] * body_force[0] + body_force[1] * body_force[1])); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( std::abs(body_force[2]) < 1e-16, "body force in z should be zero for this wall function"); } + { + amrex::ParmParse pp("transport"); + pp.query("viscosity", m_log_law.nu); + } + { + amrex::ParmParse pp("WallFunction"); + // Reference height for log-law + if (pp.contains("log_law_ref_index")) { + pp.get("log_law_ref_index", m_log_law.ref_index); + } + const auto& geom = m_mesh.Geom(0); + m_log_law.zref = + (geom.ProbLo(m_direction) + + (m_log_law.ref_index + 0.5) * geom.CellSize(m_direction)); + } } -VelWallFunc::VelWallFunc(Field& /*unused*/, const WallFunction& wall_func) +VelWallFunc::VelWallFunc(Field& /*unused*/, WallFunction& wall_func) : m_wall_func(wall_func) { amrex::ParmParse pp("WallFunction"); pp.query("wall_shear_stress_type", m_wall_shear_stress_type); m_wall_shear_stress_type = amrex::toLower(m_wall_shear_stress_type); - if (m_wall_shear_stress_type == "constant") { + if (m_wall_shear_stress_type == "constant" || + m_wall_shear_stress_type == "log_law" || + m_wall_shear_stress_type == "schumann") { amrex::Print() << "Shear Stress model: " << m_wall_shear_stress_type << std::endl; } else { @@ -42,6 +59,95 @@ VelWallFunc::VelWallFunc(Field& /*unused*/, const WallFunction& wall_func) } } +template +void VelWallFunc::wall_model( + Field& velocity, const FieldState rho_state, const ShearStressSimple& tau) +{ + BL_PROFILE("amr-wind::VelWallFunc"); + constexpr int idim = 2; + const auto& repo = velocity.repo(); + const auto& density = repo.get_field("density", rho_state); + const auto& viscosity = repo.get_field("velocity_mueff"); + const int nlevels = repo.num_active_levels(); + const auto idx_offset = tau.m_ll.ref_index; + + amrex::Orientation zlo(amrex::Direction::z, amrex::Orientation::low); + amrex::Orientation zhi(amrex::Direction::z, amrex::Orientation::high); + if ((velocity.bc_type()[zlo] != BC::wall_model) && + (velocity.bc_type()[zhi] != BC::wall_model)) { + return; + } + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = repo.mesh().Geom(lev); + const auto& domain = geom.Domain(); + amrex::MFItInfo mfi_info{}; + + const auto& rho_lev = density(lev); + auto& vel_lev = velocity(lev); + auto& vold_lev = velocity.state(FieldState::Old)(lev); + const auto& eta_lev = viscosity(lev); + + if (amrex::Gpu::notInLaunchRegion()) { + mfi_info.SetDynamic(true); + } +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(vel_lev, mfi_info); mfi.isValid(); ++mfi) { + const auto& bx = mfi.validbox(); + auto varr = vel_lev.array(mfi); + auto vold_arr = vold_lev.array(mfi); + auto den = rho_lev.array(mfi); + auto eta = eta_lev.array(mfi); + + if (bx.smallEnd(idim) == domain.smallEnd(idim) && + velocity.bc_type()[zlo] == BC::wall_model) { + amrex::ParallelFor( + amrex::bdryLo(bx, idim), + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real mu = eta(i, j, k); + const amrex::Real uu = + vold_arr(i, j, k + idx_offset, 0); + const amrex::Real vv = + vold_arr(i, j, k + idx_offset, 1); + const amrex::Real wspd = std::sqrt(uu * uu + vv * vv); + // Dirichlet BC + varr(i, j, k - 1, 2) = 0.0; + + // Shear stress BC + varr(i, j, k - 1, 0) = + tau.get_shear(uu, wspd) / mu * den(i, j, k); + varr(i, j, k - 1, 1) = + tau.get_shear(vv, wspd) / mu * den(i, j, k); + }); + } + + if (bx.bigEnd(idim) == domain.bigEnd(idim) && + velocity.bc_type()[zhi] == BC::wall_model) { + amrex::ParallelFor( + amrex::bdryHi(bx, idim), + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real mu = eta(i, j, k - 1); + const amrex::Real uu = + vold_arr(i, j, k - 1 - idx_offset, 0); + const amrex::Real vv = + vold_arr(i, j, k - 1 - idx_offset, 1); + const amrex::Real wspd = std::sqrt(uu * uu + vv * vv); + // Dirichlet BC + varr(i, j, k, 2) = 0.0; + + // Shear stress BC + varr(i, j, k, 0) = + -tau.get_shear(uu, wspd) / mu * den(i, j, k); + varr(i, j, k, 1) = + -tau.get_shear(vv, wspd) / mu * den(i, j, k); + }); + } + } + } +} + void VelWallFunc::wall_model( Field& velocity, const FieldState rho_state, const amrex::Real utau) { @@ -119,7 +225,24 @@ void VelWallFunc::operator()(Field& velocity, const FieldState rho_state) { if (m_wall_shear_stress_type == "constant") { wall_model(velocity, rho_state, m_wall_func.utau()); + } else if (m_wall_shear_stress_type == "log_law") { + m_wall_func.update_umean(); + m_wall_func.update_utau_mean(); + auto tau = SimpleShearLogLaw(m_wall_func.log_law()); + wall_model(velocity, rho_state, tau); + } else if (m_wall_shear_stress_type == "schumann") { + m_wall_func.update_umean(); + auto tau = SimpleShearSchumann(m_wall_func.log_law()); + wall_model(velocity, rho_state, tau); } } +void WallFunction::update_umean() +{ + m_pa_vel(); + m_log_law.wspd_mean = + m_pa_vel.line_hvelmag_average_interpolated(m_log_law.zref); +} + +void WallFunction::update_utau_mean() { m_log_law.update_utau_mean(); } } // namespace amr_wind diff --git a/amr-wind/turbulence/LES/AMDNoTherm.H b/amr-wind/turbulence/LES/AMDNoTherm.H new file mode 100644 index 0000000000..e937bfdbbf --- /dev/null +++ b/amr-wind/turbulence/LES/AMDNoTherm.H @@ -0,0 +1,84 @@ +#ifndef AMDNOTHERM_H +#define AMDNOTHERM_H + +#include +#include +#include +#include "amr-wind/incflo_enums.H" +#include "amr-wind/turbulence/TurbModelBase.H" +#include "amr-wind/core/FieldRepo.H" +#include "amr-wind/fvm/stencils.H" + +namespace amr_wind::turbulence { +/** AMD LES Model + * \ingroup turb_model + */ +template +class AMDNoTherm : public TurbModelBase +{ +public: + static std::string identifier() + { + return "AMDNoTherm-" + Transport::identifier(); + } + + explicit AMDNoTherm(CFDSim& sim); + + //! Model name for debugging purposes + std::string model_name() const override { return "AMDNoTherm"; } + + //! Update the turbulent viscosity field + void update_turbulent_viscosity( + const FieldState fstate, const DiffusionType /*unused*/) override; + + //! Return model coefficients dictionary + TurbulenceModel::CoeffsDictType model_coeffs() const override; + + //! Parse turbulence model coefficients for this model + void parse_model_coeffs() override; + + //! No post advance work for this model + void post_advance_work() override{}; + +private: + //! Poincare coefficient (default value set for 2nd order AMR-wind + //! discretization) + amrex::Real m_C{0.333333333333333}; + const Field& m_vel; + const Field& m_rho; +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real amd_base_muvel( + int i, + int j, + int k, + const amrex::GpuArray dx, // Grid spacing + amrex::Real C, // Poincare const + amrex::Array4 const& gradVel) noexcept +{ + + amrex::Real num_shear = 0; + amrex::Real denom = 0; + for (int ii = 0; ii < AMREX_SPACEDIM; ++ii) { + for (int jj = 0; jj < AMREX_SPACEDIM; ++jj) { + denom = denom + gradVel(i, j, k, ii * AMREX_SPACEDIM + jj) * + gradVel(i, j, k, ii * AMREX_SPACEDIM + jj); + amrex::Real sij = + 0.5 * (gradVel(i, j, k, ii * AMREX_SPACEDIM + jj) + + gradVel(i, j, k, jj * AMREX_SPACEDIM + ii)); + for (int kk = 0; kk < AMREX_SPACEDIM; ++kk) { + amrex::Real dkui = gradVel(i, j, k, ii * AMREX_SPACEDIM + kk); + amrex::Real dkuj = gradVel(i, j, k, jj * AMREX_SPACEDIM + kk); + num_shear = num_shear + dkui * dkuj * dx[kk] * dx[kk] * sij; + } + } + } + denom = std::max(1e-15, denom); + num_shear = -C * num_shear; + + return std::max(1e-15, num_shear / denom); +} + +} // namespace amr_wind::turbulence + +#endif /* AMDNoTherm_H */ diff --git a/amr-wind/turbulence/LES/AMDNoTherm.cpp b/amr-wind/turbulence/LES/AMDNoTherm.cpp new file mode 100644 index 0000000000..1358c16073 --- /dev/null +++ b/amr-wind/turbulence/LES/AMDNoTherm.cpp @@ -0,0 +1,80 @@ +#include +#include + +#include "amr-wind/fvm/gradient.H" +#include "amr-wind/turbulence/LES/AMDNoTherm.H" +#include "amr-wind/turbulence/TurbModelDefs.H" + +#include "AMReX_REAL.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParmParse.H" + +namespace amr_wind { +namespace turbulence { + +template +AMDNoTherm::AMDNoTherm(CFDSim& sim) + : TurbModelBase(sim) + , m_vel(sim.repo().get_field("velocity")) + , m_rho(sim.repo().get_field("density")) +{} + +template +void AMDNoTherm::parse_model_coeffs() +{ + const std::string coeffs_dict = this->model_name() + "_coeffs"; + amrex::ParmParse pp(coeffs_dict); + pp.query("C_poincare", m_C); +} + +template +void AMDNoTherm::update_turbulent_viscosity( + const FieldState fstate, const DiffusionType /*unused*/) +{ + BL_PROFILE( + "amr-wind::" + this->identifier() + "::update_turbulent_viscosity"); + + auto& mu_turb = this->mu_turb(); + auto& repo = mu_turb.repo(); + const auto& vel = m_vel.state(fstate); + const auto& den = m_rho.state(fstate); + auto gradVel = repo.create_scratch_field(AMREX_SPACEDIM * AMREX_SPACEDIM); + fvm::gradient(*gradVel, vel); + const auto& geom_vec = repo.mesh().Geom(); + + const amrex::Real C_poincare = this->m_C; + + const int nlevels = repo.num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = geom_vec[lev]; + + const auto& dx = geom.CellSizeArray(); + for (amrex::MFIter mfi(mu_turb(lev)); mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + const auto& gradVel_arr = (*gradVel)(lev).array(mfi); + const auto& mu_arr = mu_turb(lev).array(mfi); + const auto& rho_arr = den(lev).const_array(mfi); + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real rho = rho_arr(i, j, k); + mu_arr(i, j, k) = + rho * + amd_base_muvel(i, j, k, dx, C_poincare, gradVel_arr); + }); + } + } + + mu_turb.fillpatch(this->m_sim.time().current_time()); +} + +template +TurbulenceModel::CoeffsDictType AMDNoTherm::model_coeffs() const +{ + return TurbulenceModel::CoeffsDictType{{"C_poincare", this->m_C}}; +} + +} // namespace turbulence + +INSTANTIATE_TURBULENCE_MODEL(AMDNoTherm); + +} // namespace amr_wind diff --git a/amr-wind/turbulence/LES/CMakeLists.txt b/amr-wind/turbulence/LES/CMakeLists.txt index 1baced8e00..bae52243a3 100644 --- a/amr-wind/turbulence/LES/CMakeLists.txt +++ b/amr-wind/turbulence/LES/CMakeLists.txt @@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name} PRIVATE Smagorinsky.cpp OneEqKsgs.cpp AMD.cpp + AMDNoTherm.cpp ) diff --git a/docs/sphinx/theory/theory.rst b/docs/sphinx/theory/theory.rst index 47a87083a3..4ebb455dbe 100644 --- a/docs/sphinx/theory/theory.rst +++ b/docs/sphinx/theory/theory.rst @@ -192,9 +192,47 @@ Turbulence Models LES models for subgrid scales ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Smagorinsky model +^^^^^^^^^^^^^^^^^ -AMD model -^^^^^^^^^^^ +Simple eddy viscosity model, the dissipation is calculated using the +resolved strain rate tensor and the grid resolution as + +.. math:: + + \begin{aligned} + \tau_{ij} &= -2 \nu_t \widetilde{S}_{ij} \\ + \nu_t &= C_s^2 \Delta^2 (2 \langle S_{ij} S_{ij} \rangle)^{\frac{1}{2}} + \end{aligned} + + +AMDNoTherm model +^^^^^^^^^^^^^^^^^ +This is the implementation of the base AMD model, useful for flows without a temperature field. + +The eddy viscosity is calculated using an anisotropic derivative with a +different filter width in each direction + +.. math:: + + \begin{aligned} + \hat{\partial}_i &= \sqrt{C} \delta_i \partial_i \textrm{ for } i=1,2,3 \\ + C &= 1/3, \textrm{ Poincare coefficient for } 2^{nd} \textrm{ order gradient} \\ + \delta_i &= \textrm{Filter width along dimension } i \textrm{ for anisotropic grids} + \end{aligned} + +The anisotropic derivative is used to define the eddy viscosity as + +.. math:: + + \begin{aligned} + \tau_{ij} &= -2 \nu_t \widetilde{S}_{ij} \\ + \nu_t &= \frac{- (\hat{\partial}_k \widetilde{u}_i) (\hat{\partial}_k \widetilde{u}_j) \widetilde{S}_{ij}}{ (\partial_l \widetilde{u}_m) (\partial_l \widetilde{u}_m) } + \end{aligned} + + +AMD model (for ABL) +^^^^^^^^^^^^^^^^^^^ The eddy viscosity is calculated using an anisotropic derivative with a different filter width in each direction @@ -316,6 +354,99 @@ There is a simple unit test for both :math:`\nu_t` and :math:`D_e` in ``unit_tests/turbulence/test_turbulence_LES.cpp`` under ``test_AMD_setup_calc``. +Wall models +----------- +The wall models descibed in this section are implemented in ``AMR-wind`` for +running wall-bounded flows (non-ABL cases). + +Log-law wall model +~~~~~~~~~~~~~~~~~~ + +This wall model computes the local :math:`u_\tau` from the velocity at +the first grid cell, and uses this to compute the shear stress, which is +then used as a boundary condition. + +The log law: + +.. math:: u_{\mathrm{mag}} = u_\tau \left(\frac{1}{\kappa}\log\left(\frac{u_\tau z}{\nu}\right) + B\right). \label{eq:loglaw} + +Given a horizontal velocity magnitude +:math:`u_{\mathrm{mag}} = \sqrt{u^2 + v^2}` at +:math:`z = z_{\mathrm{ref}}`, :math:`u_\tau` can be computed using a +non-linear solve to satisfy `[eq:loglaw] <#eq:loglaw>`__. + +In ``AMR-wind`` Newton-Raphson iterations are used with a convergence +criterion of :math:`\lvert u_\tau^{n+1} - u_\tau^n \rvert < 10^{-5}`. +For this, derivative of +:math:`\frac{\partial u_{\mathrm{mag}}}{\partial {u_\tau}}` is used, + +.. math:: \frac{\partial u_{\mathrm{mag}}}{\partial {u_\tau}} = \left(\frac{1}{\kappa}\left(1+\log\left(\frac{u_\tau z_{\mathrm{ref}}}{\nu}\right)\right) + B\right) + +.. math:: u_\tau^{n+1} = u_\tau^{n} - \left(u_\tau^n \left(\frac{1}{\kappa}\log\left(\frac{u_\tau^n z_{\mathrm{ref}}}{\nu}\right) + B\right) - u_{\mathrm{mag}}\right)/\frac{\partial u_{\mathrm{mag}}}{\partial {u_\tau}}. + +Finally, the shear stress is calculated as, + +.. math:: + + \begin{aligned} + \tau_{xz} &= u_\tau^2 \frac{u}{u_\mathrm{mag}} \\ + \tau_{yz} &= u_\tau^2 \frac{v}{u_\mathrm{mag}} + \end{aligned} + +Constant stress model +~~~~~~~~~~~~~~~~~~~~~ + +NOTE: This wall model will be ill-posed unless combined with a Dirichlet +boundary condition on the other wall, :math:`\langle u \rangle` can +drift by a constant otherwise. + +This is a trivial wall model, where the shear stresses are specified as +constants. For a pressure gradient driven channel, + +.. math:: + + \begin{aligned} + u_\tau^2 &= -\frac{\mathrm{d} P}{\mathrm{d} x} \\ + \tau_{xz} &= u_\tau^2 \\ + \tau_{yz} &= 0 + \end{aligned} + +Schumann model +~~~~~~~~~~~~~~ + +NOTE: This wall model will be ill-posed unless combined with a Dirichlet +boundary condition on the other wall, :math:`\langle u \rangle` can +drift by a constant otherwise. + +This model is a modified version of the constant stress model, where the +fluctuations from a reference height :math:`z_\mathrm{ref}` are used to +add fluctuations in the shear stress. + +.. math:: + + \begin{aligned} + u_\tau^2 &= -\frac{\mathrm{d} P}{\mathrm{d} x} \\ + \tau_{xz} &= u_\tau^2 \frac{u}{\langle u_\mathrm{mag} \rangle} \\ + \tau_{yz} &= u_\tau^2 \frac{v}{\langle u_\mathrm{mag} \rangle} + \end{aligned} + +where, :math:`\langle u_\mathrm{mag} \rangle` is the planar average of +:math:`u_{\mathrm{mag}} = \sqrt{u^2 + v^2}` at :math:`z_\mathrm{ref}`. + +Symmetric wall boundary +~~~~~~~~~~~~~~~~~~~~~~~ + +This is a boundary condition to for flows with a symmetry across the +z direction (example: *half-channel* simulations) at the centerline. + +.. math:: + + \begin{aligned} + \tau_{xz} &= 0 \\ + \tau_{yz} &= 0 \\ + w &= 0 + \end{aligned} + Navigating source code ------------------------ diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt index 5fab38a2e7..065884d88c 100644 --- a/unit_tests/CMakeLists.txt +++ b/unit_tests/CMakeLists.txt @@ -14,6 +14,7 @@ add_subdirectory(fvm) add_subdirectory(multiphase) add_subdirectory(ocean_waves) add_subdirectory(projection) +add_subdirectory(boundary_conditions) if(AMR_WIND_ENABLE_MASA) add_subdirectory(mms) diff --git a/unit_tests/boundary_conditions/CMakeLists.txt b/unit_tests/boundary_conditions/CMakeLists.txt new file mode 100644 index 0000000000..7d6460d6f2 --- /dev/null +++ b/unit_tests/boundary_conditions/CMakeLists.txt @@ -0,0 +1,4 @@ +target_sources(${amr_wind_unit_test_exe_name} PRIVATE + # test cases + test_log_law.cpp + ) diff --git a/unit_tests/boundary_conditions/test_log_law.cpp b/unit_tests/boundary_conditions/test_log_law.cpp new file mode 100644 index 0000000000..dd3c159357 --- /dev/null +++ b/unit_tests/boundary_conditions/test_log_law.cpp @@ -0,0 +1,47 @@ +#include "gtest/gtest.h" +#include "aw_test_utils/MeshTest.H" +#include "aw_test_utils/test_utils.H" +#include "amr-wind/boundary_conditions/wall_models/LogLaw.H" + +namespace amr_wind_tests { + +class LogLawTest : public MeshTest +{ +protected: + void populate_parameters() override { MeshTest::populate_parameters(); } + amrex::Real log_law_actual(const amrex::Real utau) const + { + return utau * (std::log(zref * utau / nu) / 0.384 + 4.27); + } + + const amrex::Real zref = 1.0 / 32.0; + // Molecular viscosity + const amrex::Real nu = 1e-4; + // Set test tolerance to convergence test in LogLaw + const amrex::Real tol = 1.0e-5; +}; + +TEST_F(LogLawTest, test_log_law) +{ + populate_parameters(); + amr_wind::LogLaw ll; + ll.zref = zref; + ll.nu = nu; + + amrex::Real utau_expected = 0.1; + amrex::Real wspd = log_law_actual(utau_expected); + auto utau = ll.get_utau(wspd); + EXPECT_NEAR(utau, 0.1, tol); + + utau_expected = 0.05; + EXPECT_NEAR(ll.get_utau(log_law_actual(utau_expected)), utau_expected, tol); + + utau_expected = 2.0; + EXPECT_NEAR(ll.get_utau(log_law_actual(utau_expected)), utau_expected, tol); + + utau_expected = 0.5; + ll.wspd_mean = log_law_actual(utau_expected); + ll.update_utau_mean(); + EXPECT_NEAR(ll.utau_mean, utau_expected, tol); +} +} // namespace amr_wind_tests diff --git a/unit_tests/turbulence/test_turbulence_LES.cpp b/unit_tests/turbulence/test_turbulence_LES.cpp index 4fd1ec9076..03d028665b 100644 --- a/unit_tests/turbulence/test_turbulence_LES.cpp +++ b/unit_tests/turbulence/test_turbulence_LES.cpp @@ -69,6 +69,37 @@ void init_field_amd(amr_wind::Field& fld, amrex::Real scale) } } +void init_field_incomp(amr_wind::Field& fld, amrex::Real scale) +{ + const auto& mesh = fld.repo().mesh(); + const int nlevels = fld.repo().num_active_levels(); + + amrex::Real offset = 0.0; + if (fld.field_location() == amr_wind::FieldLoc::CELL) { + offset = 0.5; + } + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = mesh.Geom(lev).CellSizeArray(); + const auto& problo = mesh.Geom(lev).ProbLoArray(); + + for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) { + auto bx = mfi.growntilebox(); + const auto& farr = fld(lev).array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const amrex::Real x = problo[0] + (i + offset) * dx[0]; + const amrex::Real y = problo[1] + (j + offset) * dx[1]; + const amrex::Real z = problo[2] + (k + offset) * dx[2]; + + farr(i, j, k, 0) = 1 * x * scale; + farr(i, j, k, 1) = -2 * y * scale; + farr(i, j, k, 2) = 1 * z * scale; + }); + } + } +} + void init_field1(amr_wind::Field& fld, amrex::Real tgrad) { const auto& mesh = fld.repo().mesh(); @@ -402,4 +433,71 @@ TEST_F(TurbLESTest, test_AMD_setup_calc) EXPECT_NEAR(ae_min_val, amd_ae_answer, tol); EXPECT_NEAR(ae_max_val, amd_ae_answer, tol); } + +TEST_F(TurbLESTest, test_AMDNoTherm_setup_calc) +{ + // Parser inputs for turbulence model + const amrex::Real C = 0.3; + const amrex::Real rho0 = 1.0; + { + amrex::ParmParse pp("turbulence"); + pp.add("model", (std::string) "AMDNoTherm"); + } + { + amrex::ParmParse pp("AMDNoTherm_coeffs"); + pp.add("C_poincare", C); + } + { + amrex::ParmParse pp("incflo"); + pp.add("density", rho0); + amrex::Vector vvec{8.0, 0.0, 0.0}; + pp.addarr("velocity", vvec); + } + + // Initialize necessary parts of solver + populate_parameters(); + initialize_mesh(); + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + + // Create turbulence model + sim().create_turbulence_model(); + // Get turbulence model + auto& tmodel = sim().turbulence_model(); + + // Get coefficients + auto model_dict = tmodel.model_coeffs(); + + for (const std::pair n : model_dict) { + // Only a single model parameter, Cs + EXPECT_EQ(n.first, "C_poincare"); + EXPECT_EQ(n.second, C); + } + + // Constants for fields + const amrex::Real scale = 2.0; + // Set up velocity field with constant strainrate + auto& vel = sim().repo().get_field("velocity"); + init_field_incomp(vel, scale); + // Set up uniform unity density field + auto& dens = sim().repo().get_field("density"); + dens.setVal(rho0); + + // Update turbulent viscosity directly + tmodel.update_turbulent_viscosity( + amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); + auto& muturb = sim().repo().get_field("mu_turb"); + + // Check values of turbulent viscosity + const auto min_val = utils::field_min(muturb); + const auto max_val = utils::field_max(muturb); + const amrex::Real tol = 1e-12; + + const amrex::Real amd_answer = -C * std::pow(scale, 3) * + (dx * dx - 8 * dy * dy + dz * dz) / + (6 * scale * scale); + EXPECT_NEAR(min_val, amd_answer, tol); + EXPECT_NEAR(max_val, amd_answer, tol); +} } // namespace amr_wind_tests From ef0a19870b15b6025a3e70a62ed4fb896a9e23b5 Mon Sep 17 00:00:00 2001 From: sbidadi9 <96149491+sbidadi9@users.noreply.github.com> Date: Thu, 29 Jun 2023 16:29:24 -0600 Subject: [PATCH 06/16] Implicit treatment of source terms for SST and IDDES (#866) * Implicit treatment of source terms for SST and IDDES * Reference to Menter's paper --- amr-wind/turbulence/RANS/KOmegaSST.cpp | 21 +++++++++++++++++++-- amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp | 18 +++++++++++++++++- 2 files changed, 36 insertions(+), 3 deletions(-) diff --git a/amr-wind/turbulence/RANS/KOmegaSST.cpp b/amr-wind/turbulence/RANS/KOmegaSST.cpp index 5ebbc21fdd..d54e3208c2 100644 --- a/amr-wind/turbulence/RANS/KOmegaSST.cpp +++ b/amr-wind/turbulence/RANS/KOmegaSST.cpp @@ -197,12 +197,13 @@ void KOmegaSST::update_turbulent_viscosity( const amrex::Real diss_amb = beta_star * rho_arr(i, j, k) * sdr_amb * tke_amb; + diss_arr(i, j, k) = -beta_star * rho_arr(i, j, k) * tke_arr(i, j, k) * sdr_arr(i, j, k) + diss_amb; - tke_lhs_arr(i, j, k) = 0.5 * beta_star * rho_arr(i, j, k) * + tke_lhs_arr(i, j, k) = beta_star * rho_arr(i, j, k) * sdr_arr(i, j, k) * deltaT; // For SDR equation: @@ -221,6 +222,8 @@ void KOmegaSST::update_turbulent_viscosity( if (diff_type == DiffusionType::Crank_Nicolson) { + tke_lhs_arr(i, j, k) = 0.5 * tke_lhs_arr(i, j, k); + sdr_src_arr(i, j, k) = production_omega; sdr_diss_arr(i, j, k) = cross_diffusion; @@ -230,8 +233,22 @@ void KOmegaSST::update_turbulent_viscosity( 0.5 * std::abs(cross_diffusion) / (sdr_arr(i, j, k) + 1e-15)) * deltaT; - } else { + } else if (diff_type == DiffusionType::Implicit) { + /* Source term linearization is based on Florian + Menter's (1993) AIAA paper */ + diss_arr(i, j, k) = 0.0; + + sdr_src_arr(i, j, k) = production_omega; + + sdr_diss_arr(i, j, k) = 0.0; + sdr_lhs_arr(i, j, k) = + (2.0 * rho_arr(i, j, k) * beta * sdr_arr(i, j, k) + + std::abs(cross_diffusion) / + (sdr_arr(i, j, k) + 1e-15)) * + deltaT; + + } else { sdr_src_arr(i, j, k) = production_omega + cross_diffusion; diff --git a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp index e7127b3c55..3d863cc4ee 100644 --- a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp +++ b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp @@ -251,7 +251,7 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( tke_arr(i, j, k) / l_iddes + diss_amb; - tke_lhs_arr(i, j, k) = 0.5 * rho_arr(i, j, k) * + tke_lhs_arr(i, j, k) = rho_arr(i, j, k) * std::sqrt(tke_arr(i, j, k)) / l_iddes * deltaT; @@ -273,6 +273,8 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( if (diff_type == DiffusionType::Crank_Nicolson) { + tke_lhs_arr(i, j, k) = 0.5 * tke_lhs_arr(i, j, k); + sdr_src_arr(i, j, k) = production_omega; sdr_diss_arr(i, j, k) = cross_diffusion; @@ -283,6 +285,20 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( (sdr_arr(i, j, k) + 1e-15)) * deltaT; + } else if (diff_type == DiffusionType::Implicit) { + /* Source term linearization is based on Florian + Menter's (1993) AIAA paper */ + diss_arr(i, j, k) = 0.0; + + sdr_src_arr(i, j, k) = production_omega; + + sdr_diss_arr(i, j, k) = 0.0; + + sdr_lhs_arr(i, j, k) = + (2.0 * rho_arr(i, j, k) * beta * sdr_arr(i, j, k) + + std::abs(cross_diffusion) / + (sdr_arr(i, j, k) + 1e-15)) * + deltaT; } else { sdr_src_arr(i, j, k) = production_omega + cross_diffusion; From 05367b31460c0d9af4d36ac5c6582c56f07d25e2 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Wed, 12 Jul 2023 11:53:30 -0400 Subject: [PATCH 07/16] Remove 90 deg Latitude Constraint from Hurricane Forcing (#872) * Removing the constraint that latitude must be 90 degrees for the hurricane forcing. * Small update for formatting. * Ran clang-format... * Putting back what I deleted with clang-format... Nice job, Matt... --------- Co-authored-by: Matt Churchfield --- .../icns/source_terms/HurricaneForcing.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp b/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp index 72cc204108..0df4a05378 100644 --- a/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/HurricaneForcing.cpp @@ -21,14 +21,12 @@ HurricaneForcing::HurricaneForcing(const CFDSim& sim) : m_mesh(sim.mesh()) amrex::ParmParse pp("CoriolisForcing"); amrex::Real rot_time_period = 86400.0; pp.query("rotational_time_period", rot_time_period); - m_coriolis_factor = 2.0 * utils::two_pi() / rot_time_period; - amrex::Print() << "Geostrophic forcing: Coriolis factor = " - << m_coriolis_factor << std::endl; amrex::Real latitude = 90.0; pp.query("latitude", latitude); - AMREX_ALWAYS_ASSERT( - amrex::Math::abs(latitude - 90.0) < - static_cast(vs::DTraits::eps())); + m_coriolis_factor = (2.0 * utils::two_pi() / rot_time_period) * + std::sin(utils::radians(latitude)); + amrex::Print() << "Geostrophic forcing: Coriolis factor = " + << m_coriolis_factor << std::endl; } { From 7439a2bf587008e04c47cd2659fac03c2cfd1833 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Wed, 12 Jul 2023 11:42:12 -0600 Subject: [PATCH 08/16] Actuators: sample mac velocity instead of cell-centered velocity (#868) * new sampling is default for TurbineFastLine and TurbineFastDisk * big improvement shown for ALM * new input parameter available for Actuator, added to documentation * unit test confirms implementation of face velocity sampling --- amr-wind/wind_energy/actuator/Actuator.H | 2 + amr-wind/wind_energy/actuator/Actuator.cpp | 23 ++- .../wind_energy/actuator/ActuatorContainer.H | 12 ++ .../actuator/ActuatorContainer.cpp | 145 ++++++++++++++++++ docs/sphinx/user/inputs_Actuator.rst | 15 ++ .../wind_energy/actuator/test_act_utils.H | 26 +++- .../actuator/test_actuator_sampling.cpp | 125 +++++++++++++++ 7 files changed, 340 insertions(+), 8 deletions(-) diff --git a/amr-wind/wind_energy/actuator/Actuator.H b/amr-wind/wind_energy/actuator/Actuator.H index 7d6c2bda36..1486c55083 100644 --- a/amr-wind/wind_energy/actuator/Actuator.H +++ b/amr-wind/wind_energy/actuator/Actuator.H @@ -75,6 +75,8 @@ private: std::vector> m_actuators; std::unique_ptr m_container; + + bool m_sample_nmhalf{false}; }; } // namespace actuator diff --git a/amr-wind/wind_energy/actuator/Actuator.cpp b/amr-wind/wind_energy/actuator/Actuator.cpp index d8cb7e5698..0ecc66ae2e 100644 --- a/amr-wind/wind_energy/actuator/Actuator.cpp +++ b/amr-wind/wind_energy/actuator/Actuator.cpp @@ -34,6 +34,12 @@ void Actuator::pre_init_actions() std::string type; pp.query("type", type); pp1.query("type", type); + if (type == "TurbineFastLine" || type == "TurbineFastDisk") { + // Only one kind of sampling can be chosen. If OpenFAST is involved + // in the Actuator type, the default behavior is to sample velocity + // at n-1/2 so that forcing is at n+1/2 + m_sample_nmhalf = true; + } AMREX_ALWAYS_ASSERT(!type.empty()); auto obj = ActuatorModel::create(type, m_sim, tname, i); @@ -44,6 +50,9 @@ void Actuator::pre_init_actions() obj->read_inputs(inp); m_actuators.emplace_back(std::move(obj)); } + + // Check if sampling should be modified aside from default behavior + pp.query("sample_vel_nmhalf", m_sample_nmhalf); } void Actuator::post_init_actions() @@ -179,9 +188,17 @@ void Actuator::update_positions() m_container->update_positions(); // Sample velocities at the new locations - const auto& vel = m_sim.repo().get_field("velocity"); - const auto& density = m_sim.repo().get_field("density"); - m_container->sample_fields(vel, density); + if (m_sample_nmhalf) { + const auto& umac = m_sim.repo().get_field("u_mac"); + const auto& vmac = m_sim.repo().get_field("v_mac"); + const auto& wmac = m_sim.repo().get_field("w_mac"); + const auto& density = m_sim.repo().get_field("density"); + m_container->sample_fields(umac, vmac, wmac, density); + } else { + const auto& vel = m_sim.repo().get_field("velocity"); + const auto& density = m_sim.repo().get_field("density"); + m_container->sample_fields(vel, density); + } } /** Provide updated velocities from container to actuator instances diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.H b/amr-wind/wind_energy/actuator/ActuatorContainer.H index 7cf57cf345..030862b2b2 100644 --- a/amr-wind/wind_energy/actuator/ActuatorContainer.H +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.H @@ -78,6 +78,12 @@ public: void sample_fields(const Field& vel, const Field& density); + void sample_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density); + int num_actuator_points() const { return static_cast(m_data.position.size()); @@ -87,6 +93,12 @@ public: void interpolate_fields(const Field& vel, const Field& density); + void interpolate_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density); + void populate_field_buffers(); void initialize_particles(const int total_pts); diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp index 8a448d493f..0b4db03e40 100644 --- a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp @@ -213,6 +213,29 @@ void ActuatorContainer::sample_fields(const Field& vel, const Field& density) m_is_scattered = false; } +void ActuatorContainer::sample_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density) +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::sample_velocities"); + AMREX_ALWAYS_ASSERT(m_container_initialized && m_is_scattered); + + // Sample velocity field + interpolate_fields(umac, vmac, wmac, density); + + // Recall particles to the MPI ranks that contains their corresponding + // turbines + // Redistribute(); + + // Populate the velocity buffer that all actuator instances can access + populate_field_buffers(); + + // Indicate that the particles have been restored to their original MPI rank + m_is_scattered = false; +} + /** Helper method for ActuatorContainer::sample_fields * * Loops over the particle tiles and copies velocity data from particles to the @@ -356,6 +379,128 @@ void ActuatorContainer::interpolate_fields( } } +void ActuatorContainer::interpolate_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density) +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::interpolate_velocities"); + auto* dptr = m_pos_device.data(); + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = m_mesh.Geom(lev); + const auto dx = geom.CellSizeArray(); + const auto dxi = geom.InvCellSizeArray(); + const auto plo = geom.ProbLoArray(); + + for (ParIterType pti(*this, lev); pti.isValid(); ++pti) { + const int np = pti.numParticles(); + auto* pstruct = pti.GetArrayOfStructs()().data(); + const auto uarr = umac(lev).const_array(pti); + const auto varr = vmac(lev).const_array(pti); + const auto warr = wmac(lev).const_array(pti); + const auto darr = density(lev).const_array(pti); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(const int ip) noexcept { + auto& pp = pstruct[ip]; + // Determine offsets within the containing cell + const amrex::Real x = + (pp.pos(0) - plo[0] - 0.5 * dx[0]) * dxi[0]; + const amrex::Real y = + (pp.pos(1) - plo[1] - 0.5 * dx[1]) * dxi[1]; + const amrex::Real z = + (pp.pos(2) - plo[2] - 0.5 * dx[2]) * dxi[2]; + + const amrex::Real xf = (pp.pos(0) - plo[0]) * dxi[0]; + const amrex::Real yf = (pp.pos(1) - plo[1]) * dxi[1]; + const amrex::Real zf = (pp.pos(2) - plo[2]) * dxi[2]; + + // Index of the low corner + const int i = static_cast(std::floor(x)); + const int j = static_cast(std::floor(y)); + const int k = static_cast(std::floor(z)); + + const int i_f = static_cast(std::floor(xf)); + const int j_f = static_cast(std::floor(yf)); + const int k_f = static_cast(std::floor(zf)); + + // Interpolation weights in each direction (linear basis) + const amrex::Real wx_hi = (x - i); + const amrex::Real wy_hi = (y - j); + const amrex::Real wz_hi = (z - k); + + const amrex::Real wx_lo = 1.0 - wx_hi; + const amrex::Real wy_lo = 1.0 - wy_hi; + const amrex::Real wz_lo = 1.0 - wz_hi; + + const amrex::Real wxf_hi = (xf - i_f); + const amrex::Real wyf_hi = (yf - j_f); + const amrex::Real wzf_hi = (zf - k_f); + + const amrex::Real wxf_lo = 1.0 - wxf_hi; + const amrex::Real wyf_lo = 1.0 - wyf_hi; + const amrex::Real wzf_lo = 1.0 - wzf_hi; + + const int iproc = pp.cpu(); + + // u + int ic = 0; + pp.rdata(ic) = + wxf_lo * wy_lo * wz_lo * uarr(i_f, j, k) + + wxf_lo * wy_lo * wz_hi * uarr(i_f, j, k + 1) + + wxf_lo * wy_hi * wz_lo * uarr(i_f, j + 1, k) + + wxf_lo * wy_hi * wz_hi * uarr(i_f, j + 1, k + 1) + + wxf_hi * wy_lo * wz_lo * uarr(i_f + 1, j, k) + + wxf_hi * wy_lo * wz_hi * uarr(i_f + 1, j, k + 1) + + wxf_hi * wy_hi * wz_lo * uarr(i_f + 1, j + 1, k) + + wxf_hi * wy_hi * wz_hi * uarr(i_f + 1, j + 1, k + 1); + + // Reset position vectors so that the particles return back + // to the MPI ranks with the turbines upon redistribution + pp.pos(ic) = dptr[iproc][ic]; + + // v + ic = 1; + pp.rdata(ic) = + wx_lo * wyf_lo * wz_lo * varr(i, j_f, k) + + wx_lo * wyf_lo * wz_hi * varr(i, j_f, k + 1) + + wx_lo * wyf_hi * wz_lo * varr(i, j_f + 1, k) + + wx_lo * wyf_hi * wz_hi * varr(i, j_f + 1, k + 1) + + wx_hi * wyf_lo * wz_lo * varr(i + 1, j_f, k) + + wx_hi * wyf_lo * wz_hi * varr(i + 1, j_f, k + 1) + + wx_hi * wyf_hi * wz_lo * varr(i + 1, j_f + 1, k) + + wx_hi * wyf_hi * wz_hi * varr(i + 1, j_f + 1, k + 1); + pp.pos(ic) = dptr[iproc][ic]; + + // w + ic = 2; + pp.rdata(ic) = + wx_lo * wy_lo * wzf_lo * warr(i, j, k_f) + + wx_lo * wy_lo * wzf_hi * warr(i, j, k_f + 1) + + wx_lo * wy_hi * wzf_lo * warr(i, j + 1, k_f) + + wx_lo * wy_hi * wzf_hi * warr(i, j + 1, k_f + 1) + + wx_hi * wy_lo * wzf_lo * warr(i + 1, j, k_f) + + wx_hi * wy_lo * wzf_hi * warr(i + 1, j, k_f + 1) + + wx_hi * wy_hi * wzf_lo * warr(i + 1, j + 1, k_f) + + wx_hi * wy_hi * wzf_hi * warr(i + 1, j + 1, k_f + 1); + pp.pos(ic) = dptr[iproc][ic]; + + // density + pp.rdata(AMREX_SPACEDIM) = + wx_lo * wy_lo * wz_lo * darr(i, j, k) + + wx_lo * wy_lo * wz_hi * darr(i, j, k + 1) + + wx_lo * wy_hi * wz_lo * darr(i, j + 1, k) + + wx_lo * wy_hi * wz_hi * darr(i, j + 1, k + 1) + + wx_hi * wy_lo * wz_lo * darr(i + 1, j, k) + + wx_hi * wy_lo * wz_hi * darr(i + 1, j, k + 1) + + wx_hi * wy_hi * wz_lo * darr(i + 1, j + 1, k) + + wx_hi * wy_hi * wz_hi * darr(i + 1, j + 1, k + 1); + }); + } + } +} + /** Determine position vector of a point within each MPI rank * * Loops over the boxArray and determines the first patch that belongs to the diff --git a/docs/sphinx/user/inputs_Actuator.rst b/docs/sphinx/user/inputs_Actuator.rst index 6c0d43d140..a14203639b 100644 --- a/docs/sphinx/user/inputs_Actuator.rst +++ b/docs/sphinx/user/inputs_Actuator.rst @@ -29,6 +29,21 @@ turbines as actuator disks and actuator line models. supported are: ``TurbineFastLine``, ``TurbineFastDisk``, and ``FixedWingLine``. +.. input_param:: Actuator.sample_vel_nmhalf + + **type:** Bool, optional + + This option sets the type of velocity sampling used to inform the actuator model. + Setting this variable to true (or `1`) makes the Actuator velocity sampling use + the face-centered velocity field at the time instant of `n-1/2`. For simulations + coupled to OpenFAST, this enables the actuator forces to be implemented at the + time instant of `n+1/2`, which fits best with the underlying numerical model of AMR-Wind + and has been shown to dramatically improve the accuracy of the Actuator Line Model. + This option defaults to true (`1`) when the Actuator type is ``TurbineFastLine`` or + ``TurbineFastDisk``, and it defaults to false (`0`) in all other cases. When the + option is false, the actuator routine samples the cell-centered velocity + at the time instant `n`. + FixedWingLine """"""""""""" diff --git a/unit_tests/wind_energy/actuator/test_act_utils.H b/unit_tests/wind_energy/actuator/test_act_utils.H index 9bfc058908..9a3e41b0a6 100644 --- a/unit_tests/wind_energy/actuator/test_act_utils.H +++ b/unit_tests/wind_energy/actuator/test_act_utils.H @@ -13,9 +13,25 @@ inline void init_field(amr_wind::Field& fld) const int nlevels = fld.repo().num_active_levels(); const int ncomp = fld.num_comp(); - amrex::Real offset = 0.0; + amrex::Real offsetx = 0.0; + amrex::Real offsety = 0.0; + amrex::Real offsetz = 0.0; if (fld.field_location() == amr_wind::FieldLoc::CELL) { - offset = 0.5; + offsetx = 0.5; + offsety = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::XFACE) { + offsety = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::YFACE) { + offsetx = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::ZFACE) { + offsetx = 0.5; + offsety = 0.5; } for (int lev = 0; lev < nlevels; ++lev) { @@ -27,9 +43,9 @@ inline void init_field(amr_wind::Field& fld) const auto& farr = fld(lev).array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const amrex::Real x = problo[0] + (i + offset) * dx[0]; - const amrex::Real y = problo[1] + (j + offset) * dx[1]; - const amrex::Real z = problo[2] + (k + offset) * dx[2]; + const amrex::Real x = problo[0] + (i + offsetx) * dx[0]; + const amrex::Real y = problo[1] + (j + offsety) * dx[1]; + const amrex::Real z = problo[2] + (k + offsetz) * dx[2]; for (int d = 0; d < ncomp; d++) { farr(i, j, k, d) = x + y + z; diff --git a/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp index cad8ef04c8..c45b709ff9 100644 --- a/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp +++ b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp @@ -168,4 +168,129 @@ TEST_F(ActuatorTest, act_container) } } +TEST_F(ActuatorTest, act_container_macvels) +{ + const int nprocs = amrex::ParallelDescriptor::NProcs(); + if (nprocs > 2) { + GTEST_SKIP(); + } + + const int iproc = amrex::ParallelDescriptor::MyProc(); + initialize_mesh(); + auto& umac = + sim().repo().declare_field("u_mac", 1, 3, 1, amr_wind::FieldLoc::XFACE); + auto& vmac = + sim().repo().declare_field("v_mac", 1, 3, 1, amr_wind::FieldLoc::YFACE); + auto& wmac = + sim().repo().declare_field("w_mac", 1, 3, 1, amr_wind::FieldLoc::ZFACE); + auto& density = sim().repo().declare_field("density", 1, 3); + init_field(umac); + init_field(vmac); + init_field(wmac); + density.setVal(1.0); + + // Number of turbines in an MPI rank + const int num_turbines = 2; + // Number of points per turbine + const int num_nodes = 16; + + TestActContainer ac(mesh(), num_turbines); + auto& data = ac.get_data_obj(); + + for (int it = 0; it < num_turbines; ++it) { + data.num_pts[it] = num_nodes; + } + + ac.initialize_container(); + + { + const int lev = 0; + int idx = 0; + const amrex::Real dz = mesh().Geom(lev).CellSize(2); + const amrex::Real ypos = 32.0 * (iproc + 1); + auto& pvec = data.position; + for (int it = 0; it < num_turbines; ++it) { + const amrex::Real xpos = 32.0 * (it + 1); + for (int ni = 0; ni < num_nodes; ++ni) { + const amrex::Real zpos = (ni + 0.5) * dz; + + pvec[idx].x() = xpos; + pvec[idx].y() = ypos; + pvec[idx].z() = zpos; + ++idx; + } + } + ASSERT_EQ(idx, ac.num_actuator_points()); + } + + ac.update_positions(); + + // Check that the update positions scattered the particles to the MPI rank + // that contains the cell enclosing the particle location +#ifndef AMREX_USE_GPU + ASSERT_EQ( + ac.num_actuator_points() * nprocs, ac.NumberOfParticlesAtLevel(0)); +#endif + { + using ParIter = amr_wind::actuator::ActuatorContainer::ParIterType; + int total_particles = 0; + const int lev = 0; + for (ParIter pti(ac, lev); pti.isValid(); ++pti) { + total_particles += pti.numParticles(); + } + + if (iproc == 0) { + ASSERT_EQ(num_turbines * num_nodes * nprocs, total_particles); + } else { + ASSERT_EQ(total_particles, 0); + } + } + + ac.sample_fields(umac, vmac, wmac, density); + ac.Redistribute(); + + // Check to make sure that the velocity sampling gathered the particles back + // to their original MPI ranks +#ifndef AMREX_USE_GPU + ASSERT_EQ( + ac.num_actuator_points() * nprocs, ac.NumberOfParticlesAtLevel(0)); +#endif + { + using ParIter = amr_wind::actuator::ActuatorContainer::ParIterType; + int counter = 0; + int total_particles = 0; + const int lev = 0; + for (ParIter pti(ac, lev); pti.isValid(); ++pti) { + ++counter; + total_particles += pti.numParticles(); + } + + // All particles should have been recalled back to the same patch within + // each MPI rank + ASSERT_EQ(counter, 1); + // Total number of particles should be the same as what this MPI rank + // created + ASSERT_EQ(num_turbines * num_nodes, total_particles); + } + + // Check the interpolated velocity field + { + namespace vs = amr_wind::vs; + constexpr amrex::Real rtol = 1.0e-12; + amrex::Real rerr = 0.0; + const int npts = ac.num_actuator_points(); + const auto& pvec = data.position; + const auto& vvec = data.velocity; + for (int ip = 0; ip < npts; ++ip) { + const auto& pos = pvec[ip]; + const auto& pvel = vvec[ip]; + + const amrex::Real vval = pos.x() + pos.y() + pos.z(); + const vs::Vector vgold{vval, vval, vval}; + rerr += vs::mag_sqr(pvel - vgold); + } + EXPECT_NEAR(rerr, 0.0, rtol); + } +} + } // namespace amr_wind_tests From 9dd9e0e6503e32b8271538186b9dd5c8c5212c5f Mon Sep 17 00:00:00 2001 From: Bruce Perry <53018946+baperry2@users.noreply.github.com> Date: Fri, 14 Jul 2023 12:35:32 -0600 Subject: [PATCH 09/16] update amrex submod (#876) * update amrex submod * clang-tidy fixes for AMREX_SPACEDIM --- amr-wind/core/Field.cpp | 4 +++- amr-wind/overset/TiogaInterface.cpp | 10 +++++----- amr-wind/utilities/sampling/SamplingContainer.H | 2 +- cmake/set_compile_flags.cmake | 4 ++-- submods/amrex | 2 +- 5 files changed, 12 insertions(+), 10 deletions(-) diff --git a/amr-wind/core/Field.cpp b/amr-wind/core/Field.cpp index 17dcc7cc78..0ccf55f492 100644 --- a/amr-wind/core/Field.cpp +++ b/amr-wind/core/Field.cpp @@ -20,7 +20,9 @@ FieldInfo::FieldInfo( , m_ngrow(ngrow) , m_nstates(nstates) , m_floc(floc) - , m_bc_values(AMREX_SPACEDIM * 2, amrex::Vector(ncomp, 0.0)) + , m_bc_values( + static_cast(AMREX_SPACEDIM) * 2, + amrex::Vector(ncomp, 0.0)) , m_bc_values_dview(static_cast(ncomp) * AMREX_SPACEDIM * 2) , m_bcrec(ncomp) , m_bcrec_d(ncomp) diff --git a/amr-wind/overset/TiogaInterface.cpp b/amr-wind/overset/TiogaInterface.cpp index 347bc6af96..486c36d148 100644 --- a/amr-wind/overset/TiogaInterface.cpp +++ b/amr-wind/overset/TiogaInterface.cpp @@ -48,11 +48,11 @@ AMROversetInfo::AMROversetInfo(const int nglobal, const int nlocal) : level(nglobal) , mpi_rank(nglobal) , local_id(nglobal) - , ilow(AMREX_SPACEDIM * nglobal) - , ihigh(AMREX_SPACEDIM * nglobal) - , dims(AMREX_SPACEDIM * nglobal) - , xlo(AMREX_SPACEDIM * nglobal) - , dx(AMREX_SPACEDIM * nglobal) + , ilow(static_cast(AMREX_SPACEDIM) * nglobal) + , ihigh(static_cast(AMREX_SPACEDIM) * nglobal) + , dims(static_cast(AMREX_SPACEDIM) * nglobal) + , xlo(static_cast(AMREX_SPACEDIM) * nglobal) + , dx(static_cast(AMREX_SPACEDIM) * nglobal) , global_idmap(nlocal) , iblank_node(nlocal) , iblank_cell(nlocal) diff --git a/amr-wind/utilities/sampling/SamplingContainer.H b/amr-wind/utilities/sampling/SamplingContainer.H index d7133ac934..d29a500375 100644 --- a/amr-wind/utilities/sampling/SamplingContainer.H +++ b/amr-wind/utilities/sampling/SamplingContainer.H @@ -64,7 +64,7 @@ class SamplingContainer { public: explicit SamplingContainer(amrex::AmrCore& mesh) - : AmrParticleContainer< + : amrex::AmrParticleContainer< SNStructReal, SNStructInt, SNArrayReal, diff --git a/cmake/set_compile_flags.cmake b/cmake/set_compile_flags.cmake index 239c3f3fd9..3f4140f686 100644 --- a/cmake/set_compile_flags.cmake +++ b/cmake/set_compile_flags.cmake @@ -35,10 +35,10 @@ if (CMAKE_CXX_COMPILER_ID STREQUAL "Clang" OR if (${AMR_WIND_USE_INTERNAL_AMREX}) if ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND AMR_WIND_ENABLE_FPE_TRAP_FOR_TESTS) target_compile_options( - amrex PUBLIC $<$:-ffp-exception-behavior=maytrap>) + amrex_3d PUBLIC $<$:-ffp-exception-behavior=maytrap>) endif() target_compile_options( - amrex PUBLIC $<$:-Wno-pass-failed>) + amrex_3d PUBLIC $<$:-Wno-pass-failed>) else() if ((CMAKE_CXX_COMPILER_ID STREQUAL "Clang") AND AMR_WIND_ENABLE_FPE_TRAP_FOR_TESTS) target_compile_options( diff --git a/submods/amrex b/submods/amrex index 2f47d132d7..59a3106293 160000 --- a/submods/amrex +++ b/submods/amrex @@ -1 +1 @@ -Subproject commit 2f47d132d71c49abbeac5e081fe9ce886c621d54 +Subproject commit 59a310629397b69d4d5242ee0eba11ef8718f2aa From c4be30826216bfdd65060420af7ac7185ca9ca97 Mon Sep 17 00:00:00 2001 From: Jon Rood Date: Fri, 14 Jul 2023 13:53:07 -0600 Subject: [PATCH 10/16] Suppress some cppcheck warnings. (#878) --- cmake/amr-wind-utils.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/amr-wind-utils.cmake b/cmake/amr-wind-utils.cmake index 5b2c159b39..cfd55b4cea 100644 --- a/cmake/amr-wind-utils.cmake +++ b/cmake/amr-wind-utils.cmake @@ -108,7 +108,7 @@ macro(init_code_checks) COMMAND ${CMAKE_COMMAND} -E make_directory cppcheck # cppcheck ignores -isystem directories, so we change them to regular -I include directories (with no spaces either) COMMAND sed "s/isystem /I/g" ${CMAKE_BINARY_DIR}/compile_commands.json > cppcheck_compile_commands.json - COMMAND ${CPPCHECK_EXE} --template=gcc --inline-suppr --suppress=unusedFunction --suppress=useStlAlgorithm --std=c++17 --language=c++ --enable=all --project=cppcheck_compile_commands.json -i ${CMAKE_SOURCE_DIR}/submods/amrex/Src -i ${CMAKE_SOURCE_DIR}/submods/AMReX-Hydro -i ${CMAKE_SOURCE_DIR}/submods/googletest --output-file=cppcheck-full-report.txt -j ${NP} + COMMAND ${CPPCHECK_EXE} --template=gcc --inline-suppr --suppress=unusedFunction --suppress=useStlAlgorithm --suppress=missingIncludeSystem --std=c++17 --language=c++ --enable=all --project=cppcheck_compile_commands.json -i ${CMAKE_SOURCE_DIR}/submods/amrex/Src -i ${CMAKE_SOURCE_DIR}/submods/AMReX-Hydro -i ${CMAKE_SOURCE_DIR}/submods/googletest --output-file=cppcheck-full-report.txt -j ${NP} COMMENT "Run cppcheck on project compile_commands.json" BYPRODUCTS cppcheck-full-report.txt WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/cppcheck From 00895fbbd6ab4c22cdcfa05cb99334fe5d0871ca Mon Sep 17 00:00:00 2001 From: Jon Rood Date: Fri, 14 Jul 2023 16:53:38 -0600 Subject: [PATCH 11/16] Fix Intel OneAPI CI. (#877) --- .github/workflows/ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d92e90714c..ec308c570f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -178,11 +178,11 @@ jobs: - name: Prepare SyCL environment run: | export DEBIAN_FRONTEND=noninteractive - sudo wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB - sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB + sudo wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB + sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB echo "deb https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list sudo apt-get update - sudo apt-get install -y --no-install-recommends ninja-build intel-oneapi-dpcpp-cpp-compiler intel-oneapi-mkl-devel + sudo apt-get install -y --no-install-recommends ninja-build intel-oneapi-compiler-dpcpp-cpp intel-oneapi-mkl-devel - name: Configure and build run: | set +e From 1af032a6e42b37c4726d7006e958ab93637e2f13 Mon Sep 17 00:00:00 2001 From: sbidadi9 <96149491+sbidadi9@users.noreply.github.com> Date: Fri, 14 Jul 2023 17:46:29 -0600 Subject: [PATCH 12/16] Adding clipping of SDR production term for IDDES model (#873) Co-authored-by: Jon Rood --- amr-wind/turbulence/RANS/KOmegaSST.H | 2 -- amr-wind/turbulence/RANS/KOmegaSSTIDDES.H | 1 + amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp | 7 +++++-- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/amr-wind/turbulence/RANS/KOmegaSST.H b/amr-wind/turbulence/RANS/KOmegaSST.H index fba60bbca2..e89929fc8c 100644 --- a/amr-wind/turbulence/RANS/KOmegaSST.H +++ b/amr-wind/turbulence/RANS/KOmegaSST.H @@ -94,8 +94,6 @@ protected: amrex::Real m_buoyancy_factor = 0.0; amrex::Real m_sigma_t{0.85}; amrex::Vector m_gravity{{0.0, 0.0, -9.81}}; - - DiffusionType m_diff_type = DiffusionType::Crank_Nicolson; }; } // namespace amr_wind::turbulence diff --git a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.H b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.H index f7ab18105a..2628c68773 100644 --- a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.H +++ b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.H @@ -60,6 +60,7 @@ protected: amrex::Real m_Ct{1.87}; amrex::Real m_Cw{0.15}; amrex::Real m_kappa{0.41}; + amrex::Real m_sdr_prod_clip_factor{10.0}; }; } // namespace amr_wind::turbulence diff --git a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp index 3d863cc4ee..9178369c6e 100644 --- a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp +++ b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp @@ -38,6 +38,7 @@ void KOmegaSSTIDDES::parse_model_coeffs() pp.query("Ct", this->m_Ct); pp.query("Cw", this->m_Cw); pp.query("kappa", this->m_kappa); + pp.query("sdr_prod_clip_factor", this->m_sdr_prod_clip_factor); } template @@ -63,7 +64,8 @@ TurbulenceModel::CoeffsDictType KOmegaSSTIDDES::model_coeffs() const {"Cl", this->m_Cl}, {"Ct", this->m_Ct}, {"Cw", this->m_Cw}, - {"kappa", this->m_kappa}}; + {"kappa", this->m_kappa}, + {"sdr_prod_clip_factor", this->m_sdr_prod_clip_factor}}; } template @@ -90,6 +92,7 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( const amrex::Real Ct = this->m_Ct; const amrex::Real Cw = this->m_Cw; const amrex::Real kappa = this->m_kappa; + const amrex::Real sdr_prod_clip_factor = this->m_sdr_prod_clip_factor; auto& mu_turb = this->mu_turb(); auto lam_mu = (this->m_transport).mu(); @@ -260,7 +263,7 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( rho_arr(i, j, k) * alpha * amrex::min( tmp4 * tmp4, - 10.0 * + sdr_prod_clip_factor * amrex::max(sdr_arr(i, j, k), 0.0) * std::sqrt(tke_arr(i, j, k)) / l_iddes); From caebcd1f2841cf780ea3dcd566f05dbc21bc476f Mon Sep 17 00:00:00 2001 From: Tony Martinez Date: Mon, 17 Jul 2023 15:50:40 -0500 Subject: [PATCH 13/16] Update turbine_fast_ops.H (#879) fix bug with number of points. --- amr-wind/wind_energy/actuator/turbine/fast/turbine_fast_ops.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/actuator/turbine/fast/turbine_fast_ops.H b/amr-wind/wind_energy/actuator/turbine/fast/turbine_fast_ops.H index 56afb90c94..e32e37b60f 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/turbine_fast_ops.H +++ b/amr-wind/wind_energy/actuator/turbine/fast/turbine_fast_ops.H @@ -165,7 +165,7 @@ struct InitDataOp tdata.num_blades = sz_info[0]; // back calculate what the value per blade for number of points in // the openfast data structure - tdata.num_vel_pts_blade = sz_info[2] / tdata.num_blades - 1; + tdata.num_vel_pts_blade = sz_info[2] / tdata.num_blades; data.grid().resize(sz_info[1], sz_info[2]); tdata.chord.resize(sz_info[1]); tdata.num_pts_tower = sz_info[3]; From 487b948016840069454d9b90c212a84da05c64d2 Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Thu, 20 Jul 2023 15:47:05 -0600 Subject: [PATCH 14/16] Fix some warnings (#883) --- amr-wind/CFDSim.cpp | 13 +++++----- amr-wind/boundary_conditions/BCInterface.cpp | 2 ++ amr-wind/core/Field.cpp | 6 +++-- amr-wind/equation_systems/AdvOp_Godunov.H | 2 +- amr-wind/equation_systems/AdvOp_MOL.H | 2 +- .../equation_systems/icns/icns_advection.H | 9 ++++--- .../equation_systems/icns/icns_diffusion.H | 3 +-- .../icns/source_terms/ABLMeanBoussinesq.H | 2 +- .../icns/source_terms/ABLMeanBoussinesq.cpp | 3 ++- .../tke/source_terms/KsgsM84Src.cpp | 2 +- amr-wind/equation_systems/tke/tke_ops.H | 1 + .../bluff_body/bluff_body_ops.cpp | 18 +++++++------- .../immersed_boundary/bluff_body/box_ops.H | 8 +++---- .../bluff_body/cylinder_ops.H | 16 ++++++------- .../immersed_boundary/bluff_body/sphere_ops.H | 10 ++++---- .../relaxation_zones/hos_waves_ops.H | 2 +- .../relaxation_zones/linear_waves_ops.H | 15 ++++++------ .../relaxation_zones/stokes_waves_ops.H | 16 ++++++------- amr-wind/overset/TiogaInterface.cpp | 6 ++--- amr-wind/physics/multiphase/MultiPhase.cpp | 12 +++++----- .../incflo_apply_nodal_projection.cpp | 2 +- amr-wind/setup/init.cpp | 7 +++--- amr-wind/transport_models/TwoPhaseTransport.H | 18 +++++--------- amr-wind/turbulence/LES/AMD.cpp | 5 ++-- amr-wind/turbulence/LES/AMDNoTherm.cpp | 3 ++- amr-wind/turbulence/LES/OneEqKsgs.cpp | 2 +- amr-wind/turbulence/RANS/KOmegaSST.cpp | 5 ++-- amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp | 6 ++--- amr-wind/utilities/SecondMomentAveraging.cpp | 4 ++-- amr-wind/utilities/ThirdMomentAveraging.cpp | 6 ++--- .../utilities/averaging/TimeAveraging.cpp | 2 +- amr-wind/utilities/diagnostics.H | 18 +++++++------- amr-wind/utilities/diagnostics.cpp | 24 +++++++++---------- amr-wind/utilities/sampling/FieldNorms.cpp | 2 +- amr-wind/utilities/sampling/Sampling.cpp | 2 +- amr-wind/wind_energy/ABLFieldInit.cpp | 5 +--- amr-wind/wind_energy/ABLStats.cpp | 2 +- .../wind_energy/actuator/disk/Joukowsky_ops.H | 5 ++-- .../actuator/disk/uniform_ct_ops.H | 5 ++-- .../actuator/turbine/fast/FastIface.cpp | 3 +-- unit_tests/aw_test_utils/MeshTest.cpp | 2 +- unit_tests/core/test_field.cpp | 8 +++---- unit_tests/multiphase/test_momflux.cpp | 14 +++++------ unit_tests/turbulence/test_turbulence_LES.cpp | 8 +++---- unit_tests/utilities/test_free_surface.cpp | 7 +----- .../actuator/test_actuator_joukowsky_disk.cpp | 4 ++-- .../actuator/test_disk_functions.cpp | 2 -- unit_tests/wind_energy/test_abl_bc.cpp | 2 +- unit_tests/wind_energy/test_abl_src.cpp | 4 ++-- 49 files changed, 158 insertions(+), 167 deletions(-) diff --git a/amr-wind/CFDSim.cpp b/amr-wind/CFDSim.cpp index f85ff79e1d..b02e05f78a 100644 --- a/amr-wind/CFDSim.cpp +++ b/amr-wind/CFDSim.cpp @@ -23,18 +23,19 @@ CFDSim::~CFDSim() = default; void CFDSim::create_turbulence_model() { - std::string transport_model = "ConstTransport"; - std::string turbulence_model = "Laminar"; + std::string transport_model_name = "ConstTransport"; + std::string turbulence_model_name = "Laminar"; { amrex::ParmParse pp("transport"); - pp.query("model", transport_model); + pp.query("model", transport_model_name); } { amrex::ParmParse pp("turbulence"); - pp.query("model", turbulence_model); + pp.query("model", turbulence_model_name); } - const std::string identifier = turbulence_model + "-" + transport_model; + const std::string identifier = + turbulence_model_name + "-" + transport_model_name; m_turbulence = turbulence::TurbulenceModel::create(identifier, *this); m_turbulence->parse_model_coeffs(); } @@ -45,7 +46,7 @@ void CFDSim::init_physics() amrex::Vector phys_names; pp.queryarr("physics", phys_names); - for (auto& phy : phys_names) { + for (const auto& phy : phys_names) { m_physics_mgr.create(phy, *this); } } diff --git a/amr-wind/boundary_conditions/BCInterface.cpp b/amr-wind/boundary_conditions/BCInterface.cpp index d0c4acb72a..b17013afce 100644 --- a/amr-wind/boundary_conditions/BCInterface.cpp +++ b/amr-wind/boundary_conditions/BCInterface.cpp @@ -131,6 +131,7 @@ std::pair BCIface::get_dirichlet_udfs() "faces"); } else { inflow_udf = val; + has_inflow_udf = true; } } } @@ -146,6 +147,7 @@ std::pair BCIface::get_dirichlet_udfs() "faces"); } else { wall_udf = val; + has_wall_udf = true; } } } diff --git a/amr-wind/core/Field.cpp b/amr-wind/core/Field.cpp index 0ccf55f492..dd44acad04 100644 --- a/amr-wind/core/Field.cpp +++ b/amr-wind/core/Field.cpp @@ -268,7 +268,7 @@ void Field::fillphysbc(amrex::Real time) noexcept void Field::apply_bc_funcs(const FieldState rho_state) noexcept { BL_ASSERT(m_info->bc_initialized() && m_info->m_bc_copied_to_device); - for (auto& func : m_info->m_bc_func) { + for (const auto& func : m_info->m_bc_func) { (*func)(*this, rho_state); } } @@ -296,8 +296,9 @@ void Field::advance_states() noexcept for (int i = num_time_states() - 1; i > 0; --i) { const auto sold = static_cast(i); const auto snew = static_cast(i - 1); + // cppcheck-suppress constVariableReference auto& old_field = state(sold); - auto& new_field = state(snew); + const auto& new_field = state(snew); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { amrex::MultiFab::Copy( old_field(lev), new_field(lev), 0, 0, num_comp(), num_grow()); @@ -308,6 +309,7 @@ void Field::advance_states() noexcept void Field::copy_state(FieldState to_state, FieldState from_state) noexcept { BL_PROFILE("amr-wind::Field::copy_state"); + // cppcheck-suppress constVariableReference auto& to_field = state(to_state); const auto& from_field = state(from_state); diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 178c221812..2a93b8e6ea 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -84,7 +84,7 @@ struct AdvectionOp< const auto& geom = repo.mesh().Geom(); const auto& src_term = fields.src_term; - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; const auto& dof_field = fields.field.state(fstate); diff --git a/amr-wind/equation_systems/AdvOp_MOL.H b/amr-wind/equation_systems/AdvOp_MOL.H index 1b61ba8b04..256876ef14 100644 --- a/amr-wind/equation_systems/AdvOp_MOL.H +++ b/amr-wind/equation_systems/AdvOp_MOL.H @@ -49,7 +49,7 @@ struct AdvectionOp< const auto& repo = fields.repo; const auto& geom = repo.mesh().Geom(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term.state(fstate); const auto& dof_field = fields.field.state(fstate); const auto& den = density.state(fstate); diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index bf7e434380..d89e9c4dcf 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -314,11 +314,11 @@ struct AdvectionOp void operator()(const FieldState fstate, const amrex::Real dt) { - auto& repo = fields.repo; + const auto& repo = fields.repo; const auto& geom = repo.mesh().Geom(); const auto& src_term = fields.src_term; - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; const auto& dof_field = fields.field.state(fstate); @@ -550,8 +550,7 @@ struct AdvectionOp const amrex::Real /*time*/) { - // cppcheck-suppress constVariable - auto& repo = fields.repo; + const auto& repo = fields.repo; auto& dof_field = fields.field.state(fstate); // computation of velocity on faces requires @@ -579,7 +578,7 @@ struct AdvectionOp const auto& repo = fields.repo; const auto& geom = repo.mesh().Geom(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term.state(fstate); const auto& dof_field = fields.field.state(fstate); const auto& rho = repo.get_field("density").state(fstate); diff --git a/amr-wind/equation_systems/icns/icns_diffusion.H b/amr-wind/equation_systems/icns/icns_diffusion.H index c4872ad310..eb9f2daa72 100644 --- a/amr-wind/equation_systems/icns/icns_diffusion.H +++ b/amr-wind/equation_systems/icns/icns_diffusion.H @@ -443,8 +443,7 @@ public: { const FieldState fstate = FieldState::New; auto& repo = m_pdefields.repo; - // cppcheck-suppress constVariable - auto& field = m_pdefields.field; + const auto& field = m_pdefields.field; const auto& density = m_density.state(fstate); const int nlevels = repo.num_active_levels(); const int ndim = field.num_comp(); diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H index 328b3f83f5..d479cd9965 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H @@ -53,7 +53,7 @@ private: bool m_const_profile{false}; //! Read a temperature profile from a text file - void read_temperature_profile(std::string profile_file_name); + void read_temperature_profile(const std::string& profile_file_name); }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp index b15ad07aa7..390e37b705 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp @@ -126,7 +126,8 @@ void ABLMeanBoussinesq::mean_temperature_update(const FieldPlaneAveraging& tavg) tavg.line_average().end(), m_theta_vals.begin()); } -void ABLMeanBoussinesq::read_temperature_profile(std::string profile_file_name) +void ABLMeanBoussinesq::read_temperature_profile( + const std::string& profile_file_name) { amrex::Vector theta_ht, theta_vals; diff --git a/amr-wind/equation_systems/tke/source_terms/KsgsM84Src.cpp b/amr-wind/equation_systems/tke/source_terms/KsgsM84Src.cpp index 5e27b189c8..5855218366 100644 --- a/amr-wind/equation_systems/tke/source_terms/KsgsM84Src.cpp +++ b/amr-wind/equation_systems/tke/source_terms/KsgsM84Src.cpp @@ -36,7 +36,7 @@ void KsgsM84Src::operator()( const amrex::Real Ceps = this->m_Ceps; const amrex::Real CepsGround = this->m_CepsGround; - auto& repo = (this->m_tke).repo(); + const auto& repo = (this->m_tke).repo(); const auto geom = repo.mesh().Geom(lev); const amrex::Real dx = geom.CellSize()[0]; diff --git a/amr-wind/equation_systems/tke/tke_ops.H b/amr-wind/equation_systems/tke/tke_ops.H index 35b10bf318..11e32aef0c 100644 --- a/amr-wind/equation_systems/tke/tke_ops.H +++ b/amr-wind/equation_systems/tke/tke_ops.H @@ -36,6 +36,7 @@ struct PostSolveOp void operator()(const amrex::Real time) { + // cppcheck-suppress constVariableReference auto& field = m_fields.field; const auto& repo = field.repo(); const int nlevels = repo.num_active_levels(); diff --git a/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp b/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp index b2ef98916a..9988921584 100644 --- a/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp +++ b/amr-wind/immersed_boundary/bluff_body/bluff_body_ops.cpp @@ -31,7 +31,7 @@ void apply_mms_vel(CFDSim& sim) const int nlevels = sim.repo().num_active_levels(); const auto& levelset = sim.repo().get_field("ib_levelset"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& velocity = sim.repo().get_field("velocity"); auto& m_conv_taylor_green = sim.physics_manager().get(); @@ -50,8 +50,8 @@ void apply_mms_vel(CFDSim& sim) for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) { const auto& bx = mfi.growntilebox(); - auto phi = levelset(lev).array(mfi); - auto varr = velocity(lev).array(mfi); + const auto& phi = levelset(lev).const_array(mfi); + const auto& varr = velocity(lev).array(mfi); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; @@ -77,7 +77,7 @@ void apply_dirichlet_vel(CFDSim& sim, const amrex::Vector& vel_bc) { const int nlevels = sim.repo().num_active_levels(); auto& geom = sim.mesh().Geom(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& velocity = sim.repo().get_field("velocity"); auto& levelset = sim.repo().get_field("ib_levelset"); levelset.fillpatch(sim.time().current_time()); @@ -94,13 +94,13 @@ void apply_dirichlet_vel(CFDSim& sim, const amrex::Vector& vel_bc) for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) { const auto& bx = mfi.tilebox(); - auto varr = velocity(lev).array(mfi); - auto phi_arr = levelset(lev).array(mfi); + const auto& varr = velocity(lev).array(mfi); + const auto& phi_arr = levelset(lev).const_array(mfi); auto norm_arr = normal(lev).array(mfi); - amrex::Real velx = vel_bc[0]; - amrex::Real vely = vel_bc[1]; - amrex::Real velz = vel_bc[2]; + const amrex::Real velx = vel_bc[0]; + const amrex::Real vely = vel_bc[1]; + const amrex::Real velz = vel_bc[2]; amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/amr-wind/immersed_boundary/bluff_body/box_ops.H b/amr-wind/immersed_boundary/bluff_body/box_ops.H index 11bbaa1095..4c2175f7a0 100644 --- a/amr-wind/immersed_boundary/bluff_body/box_ops.H +++ b/amr-wind/immersed_boundary/bluff_body/box_ops.H @@ -48,9 +48,9 @@ struct InitDataOp const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& mask_node = sim.repo().get_int_field("mask_node"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& levelset = sim.repo().get_field("ib_levelset"); auto nlevels = sim.repo().num_active_levels(); @@ -61,8 +61,8 @@ struct InitDataOp const auto& dx = geom[lev].CellSizeArray(); for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) { const auto& bx = mfi.growntilebox(); - auto epsilon_node = mask_node(lev).array(mfi); - auto phi = levelset(lev).array(mfi); + const auto& epsilon_node = mask_node(lev).array(mfi); + const auto& phi = levelset(lev).array(mfi); const amrex::Real x0 = wdata.center_loc[0]; const amrex::Real y0 = wdata.center_loc[1]; diff --git a/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H b/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H index 82139ce601..2002f661b6 100644 --- a/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H +++ b/amr-wind/immersed_boundary/bluff_body/cylinder_ops.H @@ -46,9 +46,9 @@ struct InitDataOp const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& mask_node = sim.repo().get_int_field("mask_node"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& levelset = sim.repo().get_field("ib_levelset"); auto nlevels = sim.repo().num_active_levels(); @@ -59,7 +59,7 @@ struct InitDataOp const auto& dx = geom[lev].CellSizeArray(); for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) { - auto phi = levelset(lev).array(mfi); + const auto& phi = levelset(lev).array(mfi); const amrex::Real x0 = wdata.center_loc[0]; const amrex::Real y0 = wdata.center_loc[1]; @@ -69,23 +69,23 @@ struct InitDataOp bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real y = problo[1] + (j + 0.5) * dx[1]; - amrex::Real phi_glob = phi(i, j, k); - amrex::Real r = std::sqrt( + const amrex::Real phi_glob = phi(i, j, k); + const amrex::Real r = std::sqrt( (x - x0) * (x - x0) + (y - y0) * (y - y0)); - amrex::Real phi_loc = r - R; + const amrex::Real phi_loc = r - R; phi(i, j, k) = std::min(phi_loc, phi_glob); }); const auto& nbx = mfi.nodaltilebox(); - auto epsilon_node = mask_node(lev).array(mfi); + const auto& epsilon_node = mask_node(lev).array(mfi); amrex::ParallelFor( nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real x = problo[0] + i * dx[0]; const amrex::Real y = problo[1] + j * dx[1]; - amrex::Real r = std::sqrt( + const amrex::Real r = std::sqrt( (x - x0) * (x - x0) + (y - y0) * (y - y0)); if (r <= R) { diff --git a/amr-wind/immersed_boundary/bluff_body/sphere_ops.H b/amr-wind/immersed_boundary/bluff_body/sphere_ops.H index 630a5704cd..ec8fc94709 100644 --- a/amr-wind/immersed_boundary/bluff_body/sphere_ops.H +++ b/amr-wind/immersed_boundary/bluff_body/sphere_ops.H @@ -44,9 +44,9 @@ struct InitDataOp const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& mask_node = sim.repo().get_int_field("mask_node"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& levelset = sim.repo().get_field("ib_levelset"); auto nlevels = sim.repo().num_active_levels(); @@ -58,7 +58,7 @@ struct InitDataOp for (amrex::MFIter mfi(levelset(lev)); mfi.isValid(); ++mfi) { const auto& bx = mfi.growntilebox(); - auto phi = levelset(lev).array(mfi); + const auto& phi = levelset(lev).array(mfi); const amrex::Real x0 = wdata.center_loc[0]; const amrex::Real y0 = wdata.center_loc[1]; @@ -74,12 +74,12 @@ struct InitDataOp (x - x0) * (x - x0) + (y - y0) * (y - y0) + (z - z0) * (z - z0)); - amrex::Real phi_loc = r - R; + const amrex::Real phi_loc = r - R; phi(i, j, k) = std::min(phi_loc, phi_glob); }); const auto& nbx = mfi.nodaltilebox(); - auto epsilon_node = mask_node(lev).array(mfi); + const auto& epsilon_node = mask_node(lev).array(mfi); amrex::ParallelFor( nbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real x = problo[0] + i * dx[0]; diff --git a/amr-wind/ocean_waves/relaxation_zones/hos_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/hos_waves_ops.H index 105a12787e..408333129c 100644 --- a/amr-wind/ocean_waves/relaxation_zones/hos_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/hos_waves_ops.H @@ -118,7 +118,7 @@ void ReadHOSFileLev( } void StoreHOSDataLoop( - HOSWaves::MetaType& wdata, + const HOSWaves::MetaType& wdata, amrex::Array4 const& phi, amrex::Array4 const& vel, const amrex::Real* dev_eta_data, diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index adb2eb23ec..eb79ebe0b5 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -35,16 +35,17 @@ struct InitDataOp auto& sim = data.sim(); + // cppcheck-suppress constVariableReference auto& m_levelset = sim.repo().get_field("levelset"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& dx = geom.CellSizeArray(); for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) { - auto phi = m_levelset(level).array(mfi); - auto vel = m_velocity(level).array(mfi); + const auto& phi = m_levelset(level).array(mfi); + const auto& vel = m_velocity(level).array(mfi); const amrex::Real zsl = wdata.zsl; const auto& gbx3 = mfi.growntilebox(3); @@ -108,9 +109,9 @@ struct UpdateRelaxZonesOp auto& sim = data.sim(); const auto& time = sim.time().new_time(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); @@ -121,8 +122,8 @@ struct UpdateRelaxZonesOp const auto& dx = geom[lev].CellSizeArray(); for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) { - auto phi = m_ow_levelset(lev).array(mfi); - auto vel = m_ow_velocity(lev).array(mfi); + const auto& phi = m_ow_levelset(lev).array(mfi); + const auto& vel = m_ow_velocity(lev).array(mfi); const amrex::Real waveheight = wdata.wave_height; const amrex::Real wavelength = wdata.wave_length; diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index dbe5bea220..3fd83b9e5b 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -36,16 +36,16 @@ struct InitDataOp const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_levelset = sim.repo().get_field("levelset"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& dx = geom.CellSizeArray(); for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) { - auto phi = m_levelset(level).array(mfi); - auto vel = velocity(level).array(mfi); + const auto& phi = m_levelset(level).array(mfi); + const auto& vel = velocity(level).array(mfi); const amrex::Real zsl = wdata.zsl; const auto& gbx3 = mfi.growntilebox(3); @@ -71,9 +71,9 @@ struct UpdateRelaxZonesOp auto& sim = data.sim(); const auto& time = sim.time().new_time(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); @@ -84,8 +84,8 @@ struct UpdateRelaxZonesOp const auto& dx = geom[lev].CellSizeArray(); for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) { - auto phi = m_ow_levelset(lev).array(mfi); - auto vel = m_ow_velocity(lev).array(mfi); + const auto& phi = m_ow_levelset(lev).array(mfi); + const auto& vel = m_ow_velocity(lev).array(mfi); const amrex::Real waveheight = wdata.wave_height; const amrex::Real wavelength = wdata.wave_length; diff --git a/amr-wind/overset/TiogaInterface.cpp b/amr-wind/overset/TiogaInterface.cpp index 486c36d148..59e42d414d 100644 --- a/amr-wind/overset/TiogaInterface.cpp +++ b/amr-wind/overset/TiogaInterface.cpp @@ -103,7 +103,7 @@ void TiogaInterface::post_regrid_actions() void TiogaInterface::pre_overset_conn_work() { - auto& repo = m_sim.repo(); + const auto& repo = m_sim.repo(); const int num_ghost = m_sim.pde_manager().num_ghost_state(); m_iblank_cell_host = repo.create_int_scratch_field_on_host( "iblank_cell_host", 1, num_ghost, FieldLoc::CELL); @@ -122,7 +122,7 @@ void TiogaInterface::pre_overset_conn_work() void TiogaInterface::post_overset_conn_work() { - auto& repo = m_sim.repo(); + const auto& repo = m_sim.repo(); const int nlevels = repo.num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { htod_memcpy(m_iblank_cell(lev), (*m_iblank_cell_host)(lev), 0, 0, 1); @@ -389,7 +389,7 @@ void TiogaInterface::amr_to_tioga_mesh() void TiogaInterface::amr_to_tioga_iblank() { BL_PROFILE("amr-wind::TiogaInterface::amr_to_tioga_iblank"); - auto& mesh = m_sim.mesh(); + const auto& mesh = m_sim.mesh(); const int nlevels = mesh.finestLevel() + 1; // Reset local patch counter diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp index 5dbc09aeb0..1dc99e2a77 100644 --- a/amr-wind/physics/multiphase/MultiPhase.cpp +++ b/amr-wind/physics/multiphase/MultiPhase.cpp @@ -366,11 +366,11 @@ void MultiPhase::calculate_advected_facedensity() amrex::Real c_r2 = m_rho2; // Get advected vof terms at each face - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& advalpha_x = m_sim.repo().get_field("advalpha_x"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& advalpha_y = m_sim.repo().get_field("advalpha_y"); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& advalpha_z = m_sim.repo().get_field("advalpha_z"); for (int lev = 0; lev < nlevels; ++lev) { @@ -382,9 +382,9 @@ void MultiPhase::calculate_advected_facedensity() const auto& ybx = amrex::surroundingNodes(bx, 1); const auto& zbx = amrex::surroundingNodes(bx, 2); - auto aa_x = advalpha_x(lev).array(mfi); - auto aa_y = advalpha_y(lev).array(mfi); - auto aa_z = advalpha_z(lev).array(mfi); + const auto& aa_x = advalpha_x(lev).array(mfi); + const auto& aa_y = advalpha_y(lev).array(mfi); + const auto& aa_z = advalpha_z(lev).array(mfi); amrex::ParallelFor( bxg1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp index 586ae26863..b7c1a9c5a2 100644 --- a/amr-wind/projection/incflo_apply_nodal_projection.cpp +++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp @@ -414,7 +414,7 @@ void incflo::ApplyProjection( // Determine if reference pressure should be added back if (m_repo.field_exists("reference_pressure") && time != 0.0) { - auto& p0 = m_repo.get_field("reference_pressure"); + const auto& p0 = m_repo.get_field("reference_pressure"); for (int lev = 0; lev <= finest_level; lev++) { amrex::MultiFab::Add( pressure(lev), p0(lev), 0, 0, 1, pressure.num_grow()[0]); diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp index 1181470147..2482a98baa 100644 --- a/amr-wind/setup/init.cpp +++ b/amr-wind/setup/init.cpp @@ -128,11 +128,12 @@ void incflo::InitialIterations() // Add mean pressure back if available if (m_repo.field_exists("reference_pressure")) { - auto& pressure = m_repo.get_field("p"); - auto& p0 = m_repo.get_field("reference_pressure"); + // cppcheck-suppress constVariableReference + auto& press = m_repo.get_field("p"); + const auto& p0 = m_repo.get_field("reference_pressure"); for (int lev = 0; lev <= finest_level; lev++) { amrex::MultiFab::Add( - pressure(lev), p0(lev), 0, 0, 1, pressure.num_grow()[0]); + press(lev), p0(lev), 0, 0, 1, press.num_grow()[0]); } } amrex::Print() << "Completed initial pressure iterations" << std::endl diff --git a/amr-wind/transport_models/TwoPhaseTransport.H b/amr-wind/transport_models/TwoPhaseTransport.H index 9aed3ea906..52b580fa30 100644 --- a/amr-wind/transport_models/TwoPhaseTransport.H +++ b/amr-wind/transport_models/TwoPhaseTransport.H @@ -65,16 +65,13 @@ public: if (m_ifacetype == InterfaceCapturingMethod::VOF) { - // cppcheck-suppress constVariable - auto& vof = m_repo.get_field("vof"); + const auto& vof = m_repo.get_field("vof"); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { for (amrex::MFIter mfi((*mu)(lev)); mfi.isValid(); ++mfi) { const auto& vbx = mfi.growntilebox(); - const amrex::Array4& volfrac = - vof(lev).array(mfi); - const amrex::Array4& visc = - (*mu)(lev).array(mfi); + const auto& volfrac = vof(lev).const_array(mfi); + const auto& visc = (*mu)(lev).array(mfi); const amrex::Real mu1 = m_mu1; const amrex::Real mu2 = m_mu2; amrex::ParallelFor( @@ -88,18 +85,15 @@ public: } else if (m_ifacetype == InterfaceCapturingMethod::LS) { - // cppcheck-suppress constVariable - auto& levelset = m_repo.get_field("levelset"); + const auto& levelset = m_repo.get_field("levelset"); const auto& geom = m_repo.mesh().Geom(); for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) { for (amrex::MFIter mfi((*mu)(lev)); mfi.isValid(); ++mfi) { const auto& vbx = mfi.growntilebox(); const auto& dx = geom[lev].CellSizeArray(); - const amrex::Array4& visc = - (*mu)(lev).array(mfi); - const amrex::Array4& phi = - levelset(lev).array(mfi); + const auto& visc = (*mu)(lev).array(mfi); + const auto& phi = levelset(lev).const_array(mfi); const amrex::Real eps = std::cbrt(2. * dx[0] * dx[1] * dx[2]); const amrex::Real mu1 = m_mu1; diff --git a/amr-wind/turbulence/LES/AMD.cpp b/amr-wind/turbulence/LES/AMD.cpp index 7629252fb6..18f8edd39d 100644 --- a/amr-wind/turbulence/LES/AMD.cpp +++ b/amr-wind/turbulence/LES/AMD.cpp @@ -11,6 +11,7 @@ namespace amr_wind { namespace turbulence { template +// cppcheck-suppress uninitMemberVar AMD::AMD(CFDSim& sim) : TurbModelBase(sim) , m_vel(sim.repo().get_field("velocity")) @@ -46,7 +47,7 @@ void AMD::update_turbulent_viscosity( "amr-wind::" + this->identifier() + "::update_turbulent_viscosity"); auto& mu_turb = this->mu_turb(); - auto& repo = mu_turb.repo(); + const auto& repo = mu_turb.repo(); const auto& vel = m_vel.state(fstate); const auto& temp = m_temperature.state(fstate); const auto& den = m_rho.state(fstate); @@ -227,7 +228,7 @@ void AMD::update_alphaeff(Field& alphaeff) BL_PROFILE("amr-wind::" + this->identifier() + "::update_alphaeff"); - auto& repo = alphaeff.repo(); + const auto& repo = alphaeff.repo(); const auto& geom_vec = repo.mesh().Geom(); const amrex::Real C_poincare = this->m_C; namespace stencil = amr_wind::fvm::stencil; diff --git a/amr-wind/turbulence/LES/AMDNoTherm.cpp b/amr-wind/turbulence/LES/AMDNoTherm.cpp index 1358c16073..310aac57f3 100644 --- a/amr-wind/turbulence/LES/AMDNoTherm.cpp +++ b/amr-wind/turbulence/LES/AMDNoTherm.cpp @@ -13,6 +13,7 @@ namespace amr_wind { namespace turbulence { template +// cppcheck-suppress uninitMemberVar AMDNoTherm::AMDNoTherm(CFDSim& sim) : TurbModelBase(sim) , m_vel(sim.repo().get_field("velocity")) @@ -35,7 +36,7 @@ void AMDNoTherm::update_turbulent_viscosity( "amr-wind::" + this->identifier() + "::update_turbulent_viscosity"); auto& mu_turb = this->mu_turb(); - auto& repo = mu_turb.repo(); + const auto& repo = mu_turb.repo(); const auto& vel = m_vel.state(fstate); const auto& den = m_rho.state(fstate); auto gradVel = repo.create_scratch_field(AMREX_SPACEDIM * AMREX_SPACEDIM); diff --git a/amr-wind/turbulence/LES/OneEqKsgs.cpp b/amr-wind/turbulence/LES/OneEqKsgs.cpp index 4660c5ff34..864f02eb34 100644 --- a/amr-wind/turbulence/LES/OneEqKsgs.cpp +++ b/amr-wind/turbulence/LES/OneEqKsgs.cpp @@ -233,7 +233,7 @@ void OneEqKsgsM84::post_advance_work() // Update sdr field based on sfs ke auto& tke = *(this->m_tke); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& sdr = *(this->m_sdr); const amrex::Real Ce = this->m_Ce; diff --git a/amr-wind/turbulence/RANS/KOmegaSST.cpp b/amr-wind/turbulence/RANS/KOmegaSST.cpp index d54e3208c2..bd107705d0 100644 --- a/amr-wind/turbulence/RANS/KOmegaSST.cpp +++ b/amr-wind/turbulence/RANS/KOmegaSST.cpp @@ -86,11 +86,10 @@ void KOmegaSST::update_turbulent_viscosity( const auto& den = this->m_rho.state(fstate); const auto& tke = (*this->m_tke).state(fstate); const auto& sdr = (*this->m_sdr).state(fstate); - // cppcheck-suppress constVariable - auto& repo = mu_turb.repo(); + const auto& repo = mu_turb.repo(); auto& tke_lhs = (this->m_sim).repo().get_field("tke_lhs_src_term"); tke_lhs.setVal(0.0); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& sdr_lhs = (this->m_sim).repo().get_field("sdr_lhs_src_term"); auto gradK = (this->m_sim.repo()).create_scratch_field(3, 0); diff --git a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp index 9178369c6e..1d3f2cd6e3 100644 --- a/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp +++ b/amr-wind/turbulence/RANS/KOmegaSSTIDDES.cpp @@ -19,6 +19,7 @@ template KOmegaSSTIDDES::~KOmegaSSTIDDES() = default; template +// cppcheck-suppress uninitMemberVar KOmegaSSTIDDES::KOmegaSSTIDDES(CFDSim& sim) : KOmegaSST(sim) , m_rans_ind(sim.repo().declare_field("rans_indicator", 1, 1, 1)) @@ -99,12 +100,11 @@ void KOmegaSSTIDDES::update_turbulent_viscosity( const auto& den = this->m_rho.state(fstate); const auto& tke = (*this->m_tke).state(fstate); const auto& sdr = (*this->m_sdr).state(fstate); - // cppcheck-suppress constVariable - auto& repo = mu_turb.repo(); + const auto& repo = mu_turb.repo(); const auto& geom_vec = repo.mesh().Geom(); auto& tke_lhs = (this->m_sim).repo().get_field("tke_lhs_src_term"); tke_lhs.setVal(0.0); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& sdr_lhs = (this->m_sim).repo().get_field("sdr_lhs_src_term"); auto gradK = (this->m_sim.repo()).create_scratch_field(3, 0); diff --git a/amr-wind/utilities/SecondMomentAveraging.cpp b/amr-wind/utilities/SecondMomentAveraging.cpp index 334674d86e..53d6a9ea08 100644 --- a/amr-wind/utilities/SecondMomentAveraging.cpp +++ b/amr-wind/utilities/SecondMomentAveraging.cpp @@ -137,8 +137,8 @@ void SecondMomentAveraging::compute_average( m_plane_average2.line_average().data(), m_plane_average2.line_average().size()); - amrex::Real* line_avg1 = lavg1.data(); - amrex::Real* line_avg2 = lavg2.data(); + const auto* line_avg1 = lavg1.data(); + const auto* line_avg2 = lavg2.data(); amrex::Real denom = 1.0 / (amrex::Real)m_plane_average1.ncell_plane(); diff --git a/amr-wind/utilities/ThirdMomentAveraging.cpp b/amr-wind/utilities/ThirdMomentAveraging.cpp index e31b11d489..e7cc421bb0 100644 --- a/amr-wind/utilities/ThirdMomentAveraging.cpp +++ b/amr-wind/utilities/ThirdMomentAveraging.cpp @@ -158,9 +158,9 @@ void ThirdMomentAveraging::compute_average( m_plane_average3.line_average().data(), m_plane_average3.line_average().size()); - amrex::Real* line_avg1 = lavg1.data(); - amrex::Real* line_avg2 = lavg2.data(); - amrex::Real* line_avg3 = lavg3.data(); + const auto* line_avg1 = lavg1.data(); + const auto* line_avg2 = lavg2.data(); + const auto* line_avg3 = lavg3.data(); amrex::Real denom = 1.0 / (amrex::Real)m_plane_average1.ncell_plane(); diff --git a/amr-wind/utilities/averaging/TimeAveraging.cpp b/amr-wind/utilities/averaging/TimeAveraging.cpp index 55e2e723e1..4c1d35c753 100644 --- a/amr-wind/utilities/averaging/TimeAveraging.cpp +++ b/amr-wind/utilities/averaging/TimeAveraging.cpp @@ -87,7 +87,7 @@ void TimeAveraging::post_advance_work() } const amrex::Real elapsed_time = (cur_time - m_start_time); - for (auto& avg : m_averages) { + for (const auto& avg : m_averages) { (*avg)(time, m_filter, elapsed_time); } } diff --git a/amr-wind/utilities/diagnostics.H b/amr-wind/utilities/diagnostics.H index 91f4b6d1e6..7605db0121 100644 --- a/amr-wind/utilities/diagnostics.H +++ b/amr-wind/utilities/diagnostics.H @@ -22,10 +22,10 @@ amrex::Real get_vel_min( const int vdir); amrex::Real get_vel_loc( - amrex::MultiFab& vel, - amrex::iMultiFab& level_mask, - int vdir, - int ldir, + const amrex::MultiFab& vel, + const amrex::iMultiFab& level_mask, + const int vdir, + const int ldir, amrex::Real vel_max, const amrex::GpuArray problo, const amrex::GpuArray dx); @@ -47,10 +47,10 @@ amrex::Real get_macvel_min( const int vdir); amrex::Real get_macvel_loc( - amrex::MultiFab& macvel, - amrex::iMultiFab& level_mask, - int vdir, - int ldir, + const amrex::MultiFab& macvel, + const amrex::iMultiFab& level_mask, + const int vdir, + const int ldir, amrex::Real vel_max, const amrex::GpuArray problo, const amrex::GpuArray dx); @@ -63,4 +63,4 @@ amrex::Array PrintMaxMACVelLocations( } // namespace amr_wind::diagnostics -#endif \ No newline at end of file +#endif diff --git a/amr-wind/utilities/diagnostics.cpp b/amr-wind/utilities/diagnostics.cpp index 5f2cc159e6..5cc4271602 100644 --- a/amr-wind/utilities/diagnostics.cpp +++ b/amr-wind/utilities/diagnostics.cpp @@ -41,10 +41,10 @@ amrex::Real amr_wind::diagnostics::get_vel_min( } amrex::Real amr_wind::diagnostics::get_vel_loc( - amrex::MultiFab& vel, - amrex::iMultiFab& level_mask, - int vdir, - int ldir, + const amrex::MultiFab& vel, + const amrex::iMultiFab& level_mask, + const int vdir, + const int ldir, amrex::Real vel_max, const amrex::GpuArray problo, const amrex::GpuArray dx) @@ -114,10 +114,10 @@ amrex::Real amr_wind::diagnostics::get_macvel_min( } amrex::Real amr_wind::diagnostics::get_macvel_loc( - amrex::MultiFab& macvel, - amrex::iMultiFab& level_mask, - int vdir, - int ldir, + const amrex::MultiFab& macvel, + const amrex::iMultiFab& level_mask, + const int vdir, + const int ldir, amrex::Real mvel_max, const amrex::GpuArray problo, const amrex::GpuArray dx) @@ -154,7 +154,7 @@ amrex::Array amr_wind::diagnostics::PrintMaxVelLocations( BL_PROFILE("amr-wind::diagnostics::PrintMaxVelLocations"); // Get fields - auto& vel = repo.get_field("velocity"); + const auto& vel = repo.get_field("velocity"); const int finest_level = repo.num_active_levels() - 1; // Get infinity norm of velocities @@ -335,9 +335,9 @@ amrex::Array amr_wind::diagnostics::PrintMaxMACVelLocations( BL_PROFILE("amr-wind::diagnostics::PrintMaxMACVelLocations"); // Get fields - auto& u_mac = repo.get_field("u_mac"); - auto& v_mac = repo.get_field("v_mac"); - auto& w_mac = repo.get_field("w_mac"); + const auto& u_mac = repo.get_field("u_mac"); + const auto& v_mac = repo.get_field("v_mac"); + const auto& w_mac = repo.get_field("w_mac"); const int finest_level = repo.num_active_levels() - 1; // Get infinity norm of mac velocities diff --git a/amr-wind/utilities/sampling/FieldNorms.cpp b/amr-wind/utilities/sampling/FieldNorms.cpp index 8798c28576..bd58dc540d 100644 --- a/amr-wind/utilities/sampling/FieldNorms.cpp +++ b/amr-wind/utilities/sampling/FieldNorms.cpp @@ -23,7 +23,7 @@ void FieldNorms::initialize() pp.query("output_frequency", m_out_freq); } - auto& io_mng = m_sim.io_manager(); + const auto& io_mng = m_sim.io_manager(); for (const auto& fld : io_mng.plot_fields()) { ioutils::add_var_names(m_var_names, fld->name(), fld->num_comp()); } diff --git a/amr-wind/utilities/sampling/Sampling.cpp b/amr-wind/utilities/sampling/Sampling.cpp index 36281e9133..40b745fd8d 100644 --- a/amr-wind/utilities/sampling/Sampling.cpp +++ b/amr-wind/utilities/sampling/Sampling.cpp @@ -52,7 +52,7 @@ void Sampling::initialize() // Load different probe types, default probe type is line int idx = 0; m_total_particles = 0; - for (auto& lbl : labels) { + for (const auto& lbl : labels) { const std::string key = m_label + "." + lbl; amrex::ParmParse pp1(key); std::string stype = "LineSampler"; diff --git a/amr-wind/wind_energy/ABLFieldInit.cpp b/amr-wind/wind_energy/ABLFieldInit.cpp index cd696d9d61..47d82173d8 100644 --- a/amr-wind/wind_energy/ABLFieldInit.cpp +++ b/amr-wind/wind_energy/ABLFieldInit.cpp @@ -160,10 +160,7 @@ void ABLFieldInit::operator()( } void ABLFieldInit::perturb_temperature( - const int lev, - const amrex::Geometry& geom, - // cppcheck-suppress constParameter - Field& temperature) const + const int lev, const amrex::Geometry& geom, Field& temperature) const { /** Perturbations for the temperature field * diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index 3640cf2a89..b75d5bf68f 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -102,7 +102,7 @@ void ABLStats::calc_sfs_stress_avgs( BL_PROFILE("amr-wind::ABLStats::calc_sfs_stress_avgs"); - auto& repo = m_sim.repo(); + const auto& repo = m_sim.repo(); const auto& m_vel = repo.get_field("velocity"); auto gradVel = repo.create_scratch_field(9); diff --git a/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H b/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H index a3bc010396..81a03c0f19 100644 --- a/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H +++ b/amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H @@ -328,7 +328,7 @@ struct ProcessOutputsOp { private: // cppcheck-suppress uninitMemberVarPrivate - Joukowsky::DataType& m_data; + const Joukowsky::DataType& m_data; //! Path to the output directory (specified by Actuator physics class) std::string m_out_dir; @@ -340,7 +340,8 @@ private: public: // cppcheck-suppress constParameter - explicit ProcessOutputsOp(Joukowsky::DataType& data) + explicit ProcessOutputsOp( + const Joukowsky::DataType& data) : m_data(data) {} void operator()(Joukowsky::DataType& /*data*/) {} diff --git a/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H b/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H index d6616e0718..72ac4d7bdf 100644 --- a/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H +++ b/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H @@ -127,7 +127,7 @@ struct ProcessOutputsOp { private: // cppcheck-suppress uninitMemberVarPrivate - UniformCt::DataType& m_data; + const UniformCt::DataType& m_data; //! Path to the output directory (specified by Actuator physics class) std::string m_out_dir; @@ -139,7 +139,8 @@ private: public: // cppcheck-suppress constParameter - explicit ProcessOutputsOp(UniformCt::DataType& data) + explicit ProcessOutputsOp( + const UniformCt::DataType& data) : m_data(data) {} void operator()(UniformCt::DataType& /*unused*/) {} diff --git a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp index 353078b0ca..4c4df585b9 100644 --- a/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp +++ b/amr-wind/wind_energy/actuator/turbine/fast/FastIface.cpp @@ -253,7 +253,7 @@ void FastIface::fast_init_turbine(FastTurbine& fi) } } -// cppcheck-suppress constParameter +// cppcheck-suppress constParameterReference // NOLINTNEXTLINE(readability-convert-member-functions-to-static) void FastIface::fast_replay_turbine(FastTurbine& fi) { @@ -301,7 +301,6 @@ void FastIface::fast_replay_turbine(FastTurbine& fi) #endif } -// cppcheck-suppress constParameter // NOLINTNEXTLINE(readability-convert-member-functions-to-static) void FastIface::fast_restart_turbine(FastTurbine& fi) { diff --git a/unit_tests/aw_test_utils/MeshTest.cpp b/unit_tests/aw_test_utils/MeshTest.cpp index be13909f79..66555afdae 100644 --- a/unit_tests/aw_test_utils/MeshTest.cpp +++ b/unit_tests/aw_test_utils/MeshTest.cpp @@ -34,7 +34,7 @@ void MeshTest::reset_prob_domain() pp.getarr("is_periodic", periodic, 0, AMREX_SPACEDIM); amrex::RealBox rb(problo.data(), probhi.data()); - amrex::Geometry* gg = amrex::AMReX::top()->getDefaultGeometry(); + const amrex::Geometry* gg = amrex::AMReX::top()->getDefaultGeometry(); if (gg != nullptr) { amrex::Geometry::ResetDefaultProbDomain(rb); diff --git a/unit_tests/core/test_field.cpp b/unit_tests/core/test_field.cpp index 2eef72185c..7beaf67f2d 100644 --- a/unit_tests/core/test_field.cpp +++ b/unit_tests/core/test_field.cpp @@ -318,7 +318,7 @@ TEST_F(FieldRepoTest, int_scratch_fields) populate_parameters(); create_mesh_instance(); - auto& frepo = mesh().field_repo(); + const auto& frepo = mesh().field_repo(); // Check that int scratch field creation is disallowed before mesh is created #if !(defined(AMREX_USE_MPI) && defined(__APPLE__)) EXPECT_THROW( @@ -326,10 +326,8 @@ TEST_F(FieldRepoTest, int_scratch_fields) #endif initialize_mesh(); - // cppcheck-suppress constVariable auto ibcell_host = frepo.create_int_scratch_field_on_host( "iblank_cell_host", 1, 0, amr_wind::FieldLoc::CELL); - // cppcheck-suppress constVariable auto ibnode_host = frepo.create_int_scratch_field_on_host( "iblank_node_host", 1, 0, amr_wind::FieldLoc::NODE); @@ -351,9 +349,9 @@ TEST_F(FieldRepoTest, int_fields) initialize_mesh(); auto& frepo = mesh().field_repo(); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& ibcell = frepo.declare_int_field("iblank_cell", 1, 0, 1); - // cppcheck-suppress constVariable + // cppcheck-suppress constVariableReference auto& ibnode = frepo.declare_int_field( "iblank_node", 1, 0, 1, amr_wind::FieldLoc::NODE); diff --git a/unit_tests/multiphase/test_momflux.cpp b/unit_tests/multiphase/test_momflux.cpp index a098a1a2d0..5c1f3320a4 100644 --- a/unit_tests/multiphase/test_momflux.cpp +++ b/unit_tests/multiphase/test_momflux.cpp @@ -204,7 +204,7 @@ class MassMomFluxOpTest : public MeshTest const auto& advrho_f = repo.get_field("advalpha_" + sdir); // Get convective term - auto& conv_term = + const auto& conv_term = mom_eqn.fields().conv_term.state(amr_wind::FieldState::New); // Base level @@ -214,12 +214,12 @@ class MassMomFluxOpTest : public MeshTest const auto& problo = geom[lev].ProbLoArray(); for (amrex::MFIter mfi(vof(lev)); mfi.isValid(); ++mfi) { - const auto& um = umac(lev).array(mfi); - const auto& vm = vmac(lev).array(mfi); - const auto& wm = wmac(lev).array(mfi); - const auto& rf = advrho_f(lev).array(mfi); - const auto& vel = velocity(lev).array(mfi); - const auto& dqdt = conv_term(lev).array(mfi); + const auto& um = umac(lev).const_array(mfi); + const auto& vm = vmac(lev).const_array(mfi); + const auto& wm = wmac(lev).const_array(mfi); + const auto& rf = advrho_f(lev).const_array(mfi); + const auto& vel = velocity(lev).const_array(mfi); + const auto& dqdt = conv_term(lev).const_array(mfi); // Small mesh, loop in serial for check for (int i = 0; i < 2; ++i) { diff --git a/unit_tests/turbulence/test_turbulence_LES.cpp b/unit_tests/turbulence/test_turbulence_LES.cpp index 03d028665b..ef1cee932f 100644 --- a/unit_tests/turbulence/test_turbulence_LES.cpp +++ b/unit_tests/turbulence/test_turbulence_LES.cpp @@ -209,7 +209,7 @@ TEST_F(TurbLESTest, test_smag_setup_calc) // Update turbulent viscosity directly tmodel.update_turbulent_viscosity( amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); - auto& muturb = sim().repo().get_field("mu_turb"); + const auto& muturb = sim().repo().get_field("mu_turb"); // Check values of turbulent viscosity auto min_val = utils::field_min(muturb); @@ -321,7 +321,7 @@ TEST_F(TurbLESTest, test_1eqKsgs_setup_calc) // Update turbulent viscosity directly tmodel.update_turbulent_viscosity( amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); - auto& muturb = sim().repo().get_field("mu_turb"); + const auto& muturb = sim().repo().get_field("mu_turb"); // Check values of turbulent viscosity const auto min_val = utils::field_min(muturb); @@ -408,7 +408,7 @@ TEST_F(TurbLESTest, test_AMD_setup_calc) // Update turbulent viscosity directly tmodel.update_turbulent_viscosity( amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); - auto& muturb = sim().repo().get_field("mu_turb"); + const auto& muturb = sim().repo().get_field("mu_turb"); // Check values of turbulent viscosity const auto min_val = utils::field_min(muturb); @@ -487,7 +487,7 @@ TEST_F(TurbLESTest, test_AMDNoTherm_setup_calc) // Update turbulent viscosity directly tmodel.update_turbulent_viscosity( amr_wind::FieldState::New, DiffusionType::Crank_Nicolson); - auto& muturb = sim().repo().get_field("mu_turb"); + const auto& muturb = sim().repo().get_field("mu_turb"); // Check values of turbulent viscosity const auto min_val = utils::field_min(muturb); diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index d85c560414..aef9912540 100644 --- a/unit_tests/utilities/test_free_surface.cpp +++ b/unit_tests/utilities/test_free_surface.cpp @@ -137,10 +137,7 @@ void init_vof_slope( class FSRefineMesh : public AmrTestMesh { public: - FSRefineMesh() - : m_sim(*this) - , m_repo(m_sim.repo()) - , m_mesh_refiner(new amr_wind::RefineCriteriaManager(m_sim)) + FSRefineMesh() : m_mesh_refiner(new amr_wind::RefineCriteriaManager(m_sim)) {} amr_wind::FieldRepo& field_repo() { return m_repo; } amr_wind::CFDSim& sim() { return m_sim; } @@ -193,8 +190,6 @@ class FSRefineMesh : public AmrTestMesh } private: - amr_wind::CFDSim m_sim; - amr_wind::FieldRepo& m_repo; std::unique_ptr m_mesh_refiner; }; diff --git a/unit_tests/wind_energy/actuator/test_actuator_joukowsky_disk.cpp b/unit_tests/wind_energy/actuator/test_actuator_joukowsky_disk.cpp index 23f9f064e0..d12293c861 100644 --- a/unit_tests/wind_energy/actuator/test_actuator_joukowsky_disk.cpp +++ b/unit_tests/wind_energy/actuator/test_actuator_joukowsky_disk.cpp @@ -163,7 +163,7 @@ struct ProcessOutputsOp<::amr_wind_tests::Joukowsky, ActSrcDisk> { ProcessOutputsOp<::amr_wind_tests::Joukowsky, ActSrcDisk>( ::amr_wind_tests::Joukowsky::DataType& /**/) - {} // cppcheck-suppress missingReturn + {} void operator()(::amr_wind_tests::Joukowsky::DataType& /*data*/) {} void read_io_options(const utils::ActParser& /**/) {} void prepare_outputs(const std::string& /**/) {} @@ -205,4 +205,4 @@ TEST_F(ActJoukowskyTest, execution) act.pre_init_actions(); act.post_init_actions(); } -} // namespace amr_wind_tests \ No newline at end of file +} // namespace amr_wind_tests diff --git a/unit_tests/wind_energy/actuator/test_disk_functions.cpp b/unit_tests/wind_energy/actuator/test_disk_functions.cpp index cd75fcbfa3..dca60c2dd4 100644 --- a/unit_tests/wind_energy/actuator/test_disk_functions.cpp +++ b/unit_tests/wind_energy/actuator/test_disk_functions.cpp @@ -39,7 +39,6 @@ class TestAreaComputer : public ::testing::TestWithParam amrex::Real area; }; -// cppcheck-suppress uninitDerivedMemberVar TEST_P(TestAreaComputer, area_matches) { const auto params = GetParam(); @@ -53,7 +52,6 @@ TEST_P(TestAreaComputer, area_matches) EXPECT_NEAR(area, area_computed, 1.e-12 * area); } -// cppcheck-suppress uninitDerivedMemberVar TEST_P(TestAreaComputer, weight_sums_to_one) { const auto params = GetParam(); diff --git a/unit_tests/wind_energy/test_abl_bc.cpp b/unit_tests/wind_energy/test_abl_bc.cpp index 811d1a4c8a..a702065612 100644 --- a/unit_tests/wind_energy/test_abl_bc.cpp +++ b/unit_tests/wind_energy/test_abl_bc.cpp @@ -143,7 +143,7 @@ TEST_F(ABLMeshTest, abl_wall_model) icns_eq.compute_mueff(amr_wind::FieldState::Old); // Check test setup by verifying mu - auto& viscosity = sim().repo().get_field("velocity_mueff"); + const auto& viscosity = sim().repo().get_field("velocity_mueff"); EXPECT_NEAR(mu, utils::field_max(viscosity), tol); EXPECT_NEAR(mu, utils::field_min(viscosity), tol); diff --git a/unit_tests/wind_energy/test_abl_src.cpp b/unit_tests/wind_energy/test_abl_src.cpp index 9a4a4377b4..ba4ccbaec6 100644 --- a/unit_tests/wind_energy/test_abl_src.cpp +++ b/unit_tests/wind_energy/test_abl_src.cpp @@ -364,7 +364,7 @@ TEST_F(ABLMeshTest, coriolis_const_vel) // Velocity in x-direction test { - amrex::Real golds[AMREX_SPACEDIM] = { + const amrex::Real golds[AMREX_SPACEDIM] = { 0.0, -corfac * latfac * vel_comp, corfac * latfac * vel_comp}; vel.setVal(0.0); src_term.setVal(0.0); @@ -388,7 +388,7 @@ TEST_F(ABLMeshTest, coriolis_const_vel) // Velocity in y-direction test { - amrex::Real golds[AMREX_SPACEDIM] = { + const amrex::Real golds[AMREX_SPACEDIM] = { corfac * latfac * vel_comp, 0.0, 0.0}; vel.setVal(0.0); src_term.setVal(0.0); From f701e87727b191f14ce167c89f5d0ec7fcc0dcfe Mon Sep 17 00:00:00 2001 From: "Marc T. Henry de Frahan" Date: Fri, 21 Jul 2023 10:14:17 -0600 Subject: [PATCH 15/16] Remove cppcheck suppressions (#884) --- amr-wind/equation_systems/PDEOps.H | 1 - amr-wind/main.cpp | 17 +++++++---------- .../wind_energy/actuator/disk/uniform_ct_ops.H | 1 - 3 files changed, 7 insertions(+), 12 deletions(-) diff --git a/amr-wind/equation_systems/PDEOps.H b/amr-wind/equation_systems/PDEOps.H index 005f9f9a13..f845bb14e3 100644 --- a/amr-wind/equation_systems/PDEOps.H +++ b/amr-wind/equation_systems/PDEOps.H @@ -83,7 +83,6 @@ struct SrcTermOpBase for (auto& src_name : src_terms) { // Prefer to use emplace_back here - // cppcheck-suppress useStlAlgorithm sources.emplace_back(PDE::SrcTerm::create(src_name, sim)); } } diff --git a/amr-wind/main.cpp b/amr-wind/main.cpp index c4f1855e55..6e70fafdea 100644 --- a/amr-wind/main.cpp +++ b/amr-wind/main.cpp @@ -17,16 +17,13 @@ int main(int argc, char* argv[]) return 1; } - // cppcheck-suppress knownConditionTrueFalse - if (argc >= 2) { - // Look for "-h" or "--help" flag and print usage - for (auto i = 1; i < argc; i++) { - const std::string param(argv[i]); - if ((param == "--help") || (param == "-h")) { - amr_wind::io::print_banner(MPI_COMM_WORLD, std::cout); - amr_wind::io::print_usage(MPI_COMM_WORLD, std::cout); - return 0; - } + // Look for "-h" or "--help" flag and print usage + for (auto i = 1; i < argc; i++) { + const std::string param(argv[i]); + if ((param == "--help") || (param == "-h")) { + amr_wind::io::print_banner(MPI_COMM_WORLD, std::cout); + amr_wind::io::print_usage(MPI_COMM_WORLD, std::cout); + return 0; } } diff --git a/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H b/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H index 72ac4d7bdf..5de30a0760 100644 --- a/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H +++ b/amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H @@ -138,7 +138,6 @@ private: int m_out_freq{10}; public: - // cppcheck-suppress constParameter explicit ProcessOutputsOp( const UniformCt::DataType& data) : m_data(data) From 761b26dc2a542b4c7c8ba47d6ba8482c3dfdfe8b Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Mon, 24 Jul 2023 10:16:28 -0600 Subject: [PATCH 16/16] Time intervals for plt and chkpt output (#881) * ability to output plt and chkpt files based on time intervals (not just based on timestep intervals) * option to enforce dt (to line up with output intervals) is available as input bool --------- Co-authored-by: Jon Rood --- amr-wind/core/SimTime.H | 46 +++++++- amr-wind/core/SimTime.cpp | 112 ++++++++++++++++++- unit_tests/core/test_simtime.cpp | 186 +++++++++++++++++++++++++++++++ 3 files changed, 336 insertions(+), 8 deletions(-) diff --git a/amr-wind/core/SimTime.H b/amr-wind/core/SimTime.H index 0674ed3ac0..8a1a7604d4 100644 --- a/amr-wind/core/SimTime.H +++ b/amr-wind/core/SimTime.H @@ -4,6 +4,21 @@ #include "AMReX_REAL.H" #include "AMReX_GpuQualifiers.H" #include "AMReX_Extension.H" +#include +#include + +namespace { +AMREX_FORCE_INLINE amrex::Real get_enforced_dt_for_output( + const amrex::Real dt, + const amrex::Real cur_time, + const amrex::Real interval, + const amrex::Real tol) +{ + // Tolerance for being slightly under is relative to requested interval + return std::min( + dt, (std::floor(cur_time / interval + tol) + 1) * interval - cur_time); +} +} // namespace namespace amr_wind { @@ -158,6 +173,9 @@ private: //! Initial reduction in timestep size for startup amrex::Real m_init_shrink{0.1}; + //! Maximum growth of dt between timesteps + amrex::Real m_dt_growth{0.1}; + //! Counter for the number of timesteps since start of simulation int m_time_index{0}; @@ -176,12 +194,36 @@ private: //! Maximum timesteps for simulation int m_stop_time_index{-1}; - //! Time interval for plot file output + //! Time step interval for plot file output int m_plt_interval{-1}; - //! Time interval for writing checkpoint/restart files + //! Time interval for plot file output + amrex::Real m_plt_t_interval{-1.0}; + + //! Bool for if plt time interval should be forced + bool m_force_plt_dt{false}; + + //! Relative (to dt) tolerance for plot time interval output + amrex::Real m_plt_t_tol{1e-8}; + + //! Relative (to plt_t_interval) tolerance for enforcing dt + amrex::Real m_force_plt_tol{1e-3}; + + //! Time step interval for writing checkpoint/restart files int m_chkpt_interval{-1}; + //! Time interval for writing checkpoint/restart files + amrex::Real m_chkpt_t_interval{-1.0}; + + //! Bool for if checkpoint time interval should be forced + bool m_force_chkpt_dt{false}; + + //! Relative (to dt) tolerance for checkpoint time interval output + amrex::Real m_chkpt_t_tol{1e-8}; + + //! Relative (to chkpt_t_interval) tolerance for enforcing dt + amrex::Real m_force_chkpt_tol{1e-3}; + //! Time interval for regridding int m_regrid_interval{-1}; diff --git a/amr-wind/core/SimTime.cpp b/amr-wind/core/SimTime.cpp index d4fbf1f01b..3e0947313f 100644 --- a/amr-wind/core/SimTime.cpp +++ b/amr-wind/core/SimTime.cpp @@ -19,21 +19,75 @@ void SimTime::parse_parameters() pp.query("fixed_dt", m_fixed_dt); pp.query("initial_dt", m_initial_dt); pp.query("init_shrink", m_init_shrink); + pp.query("max_dt_growth", m_dt_growth); pp.query("cfl", m_max_cfl); pp.query("verbose", m_verbose); pp.query("regrid_interval", m_regrid_interval); pp.query("plot_interval", m_plt_interval); + pp.query("plot_time_interval", m_plt_t_interval); + pp.query("enforce_plot_time_dt", m_force_plt_dt); pp.query("checkpoint_interval", m_chkpt_interval); + pp.query("checkpoint_time_interval", m_chkpt_t_interval); + pp.query("enforce_checkpoint_time_dt", m_force_chkpt_dt); pp.query("regrid_start", m_regrid_start_index); pp.query("plot_start", m_plt_start_index); pp.query("checkpoint_start", m_chkpt_start_index); pp.query("use_force_cfl", m_use_force_cfl); + // Tolerances + pp.query("plot_time_interval_reltol", m_plt_t_tol); + pp.query("enforce_plot_dt_reltol", m_force_plt_tol); + pp.query("checkpoint_time_interval_reltol", m_chkpt_t_tol); + pp.query("enforce_checkpoint_dt_reltol", m_force_chkpt_tol); + if (m_fixed_dt > 0.0) { m_dt[0] = m_fixed_dt; } else { m_adaptive = true; } + + if (m_plt_interval > 0 && m_plt_t_interval > 0.0) { + amrex::Abort( + "plot_interval and plot_time_interval are both specified. " + "timestep- and time-based plotting should not be used together; " + "please only specify one."); + } + + if (m_chkpt_interval > 0 && m_chkpt_t_interval > 0.0) { + amrex::Abort( + "checkpoint_interval and checkpoint_time_interval are both " + "specified. timestep- and time-based checkpointing should not be " + "used together; please only specify one."); + } + + if (m_plt_t_interval <= 0.0 && m_force_plt_dt) { + amrex::Abort( + "enforce_plot_time_dt is true, but no plot time interval has been " + "provided."); + } + + if (m_chkpt_t_interval <= 0.0 && m_force_chkpt_dt) { + amrex::Abort( + "enforce_checkpoint_time_dt is true, but no checkpoint time " + "interval has been provided."); + } + + if (!m_adaptive && (m_force_plt_dt || m_force_chkpt_dt)) { + amrex::Abort( + "an output time interval has been specified to be enforced upon " + "dt, but dt is not adaptive."); + } + + if (m_force_plt_dt && m_force_chkpt_dt) { + amrex::Print() + << "WARNING: Time intervals will be enforced upon dt for both " + "plotfiles and checkpoint files."; + amrex::Print() + << " -- If these time intervals are different and not factors of " + "each other, inefficient behavior may result, depending on " + "tolerances: dt may be shortened often and outputs may occur in " + "consecutive timesteps."; + } } bool SimTime::new_timestep() @@ -110,9 +164,9 @@ void SimTime::set_current_cfl( dt_new *= m_init_shrink; } - // Limit timestep growth to 10% per timestep + // Limit timestep growth to % per timestep if (m_dt[0] > 0.0) { - dt_new = amrex::min(dt_new, 1.1 * m_dt[0]); + dt_new = amrex::min(dt_new, (1.0 + m_dt_growth) * m_dt[0]); } // Don't overshoot stop time @@ -123,6 +177,25 @@ void SimTime::set_current_cfl( if (m_adaptive) { m_dt[0] = dt_new; + // Shorten timestep to hit output frequency exactly + if (m_chkpt_t_interval > 0.0 && m_force_chkpt_dt) { + // Shorten dt if going to overshoot next output time + m_dt[0] = get_enforced_dt_for_output( + m_dt[0], m_cur_time, m_chkpt_interval, m_force_chkpt_tol); + // how it works: the floor operator gets the index of the last + // output, with a tolerance proportional to the current dt. + // adding 1 and multiplying by the time interval finds the next + // expected output time. if the time increment between the current + // time and the next expected output is less than the current dt, + // then the dt becomes this increment to get to the next output + // time. + } + if (m_plt_t_interval > 0.0 && m_force_plt_dt) { + // Shorten dt if going to overshoot next output time + m_dt[0] = get_enforced_dt_for_output( + m_dt[0], m_cur_time, m_plt_t_interval, m_force_plt_tol); + } + if (m_is_init && m_initial_dt > 0.0) { m_dt[0] = amrex::min(dt_new, m_initial_dt); } @@ -200,16 +273,43 @@ bool SimTime::do_regrid() const bool SimTime::write_plot_file() const { + // If dt is enforced, allow smallest tolerance to be in effect. This avoids + // unintentionally plotting in consecutive timesteps because of shortened dt + amrex::Real tol = m_plt_t_tol * m_dt[0]; + tol = + (m_force_chkpt_dt + ? std::min(tol, m_force_chkpt_tol * m_chkpt_t_interval) + : tol); + tol = + (m_force_plt_dt ? std::min(tol, m_force_plt_tol * m_plt_t_interval) + : tol); return ( - (m_plt_interval > 0) && - ((m_time_index - m_plt_start_index) % m_plt_interval == 0)); + ((m_plt_interval > 0) && + ((m_time_index - m_plt_start_index) % m_plt_interval == 0)) || + (m_plt_t_interval > 0.0 && + ((m_new_time + tol) / m_plt_t_interval - + std::floor((m_new_time + tol) / m_plt_t_interval) < + m_dt[0] / m_plt_t_interval))); } bool SimTime::write_checkpoint() const { + // If dt is enforced, use smallest tolerance + amrex::Real tol = m_plt_t_tol * m_dt[0]; + tol = + (m_force_chkpt_dt + ? std::min(tol, m_force_chkpt_tol * m_chkpt_t_interval) + : tol); + tol = + (m_force_plt_dt ? std::min(tol, m_force_plt_tol * m_plt_t_interval) + : tol); return ( - (m_chkpt_interval > 0) && - ((m_time_index - m_chkpt_start_index) % m_chkpt_interval == 0)); + ((m_chkpt_interval > 0) && + ((m_time_index - m_chkpt_start_index) % m_chkpt_interval == 0)) || + (m_chkpt_t_interval > 0.0 && + ((m_new_time + tol) / m_chkpt_t_interval - + std::floor((m_new_time + tol) / m_chkpt_t_interval) < + m_dt[0] / m_chkpt_t_interval))); } bool SimTime::write_last_plot_file() const diff --git a/unit_tests/core/test_simtime.cpp b/unit_tests/core/test_simtime.cpp index f70e25d165..ba18d8165c 100644 --- a/unit_tests/core/test_simtime.cpp +++ b/unit_tests/core/test_simtime.cpp @@ -83,6 +83,7 @@ TEST_F(SimTimeTest, time_loop) if (time.do_regrid()) { ++regrid_counter; } + std::cout << time.new_time() << std::endl; } EXPECT_EQ(counter, 5); EXPECT_EQ(plot_counter, 5); @@ -130,4 +131,189 @@ TEST_F(SimTimeTest, fixed_dt_loop) EXPECT_FALSE(time.write_last_plot_file()); } +TEST_F(SimTimeTest, plt_chk_timeinterval_loop) +{ + build_simtime_params(); + { + amrex::ParmParse pp("time"); + pp.add("regrid_interval", -1); + pp.add("checkpoint_interval", -1); + pp.add("plot_interval", -1); + pp.add("checkpoint_time_interval", 4.0); + pp.add("plot_time_interval", 1.0); + pp.add("fixed_dt", 0.3); + pp.add("stop_time", 5.0); + pp.add("max_step", 100); + } + amr_wind::SimTime time; + time.parse_parameters(); + + int counter = 0; + int plot_counter = 0; + int plot_step_sum = 0; + int chkpt_counter = 0; + int chkpt_step_sum = 0; + while (time.new_timestep()) { + time.set_current_cfl(0.45 / 0.3, 0.0, 0.0); + ++counter; + + if (time.write_plot_file()) { + ++plot_counter; + plot_step_sum += counter; + } + if (time.write_checkpoint()) { + ++chkpt_counter; + chkpt_step_sum += counter; + } + } + EXPECT_EQ(plot_counter, 5); + EXPECT_EQ(plot_step_sum, 4 + 7 + 10 + 14 + 17); + EXPECT_EQ(chkpt_counter, 1); + EXPECT_EQ(chkpt_step_sum, 14); +} + +TEST_F(SimTimeTest, enforce_dt_out) +{ + // Should not change if already correct + amrex::Real result = get_enforced_dt_for_output(0.1, 3.9, 2.0, 1e-3); + EXPECT_NEAR(result, 0.1, 1e-12); + + // Should not change if short of interval + result = get_enforced_dt_for_output(0.1, 3.9 - 2.0 * 5e-4, 2.0, 1e-3); + EXPECT_NEAR(result, 0.1, 1e-12); + + // Should not change if starting near interval + result = get_enforced_dt_for_output(0.1, 4.0, 2.0, 1e-3); + EXPECT_NEAR(result, 0.1, 1e-12); + result = get_enforced_dt_for_output(0.1, 4.0 - 2.0 * 0.99e-3, 2.0, 1e-3); + EXPECT_NEAR(result, 0.1, 1e-12); + // Past the tolerance, will change + result = get_enforced_dt_for_output(0.1, 4.0 - 2.0 * 1.01e-3, 2.0, 1e-3); + EXPECT_LT(result, 0.1); + + // Shortens dt if overlapping with intervals + result = get_enforced_dt_for_output(0.1, 3.95, 2.0, 1e-3); + EXPECT_NEAR(result, 0.05, 1e-12); +} + +TEST_F(SimTimeTest, enforce_timeinterval) +{ + build_simtime_params(); + { + amrex::ParmParse pp("time"); + pp.add("regrid_interval", -1); + pp.add("checkpoint_interval", -1); + pp.add("plot_interval", -1); + + pp.add("plot_time_interval", 0.5); + pp.add("enforce_plot_time_dt", true); + + // Default values for tolerances + pp.add("stop_time", 1.0); + pp.add("max_step", 10); + } + amr_wind::SimTime time; + time.parse_parameters(); + + int counter = 0; + int plot_counter = 0; + amrex::Real plot_time_sum = 0.0; + int plot_step_sum = 0; + while (time.new_timestep()) { + time.set_current_cfl(0.45 / 0.4, 0.0, 0.0); + ++counter; + if (time.write_plot_file()) { + ++plot_counter; + plot_time_sum += time.new_time(); + plot_step_sum += counter; + } + } + EXPECT_EQ(plot_counter, 2); + EXPECT_NEAR(plot_time_sum, 1.5, 1e-8); + EXPECT_EQ(plot_step_sum, 2 + 6); +} + +TEST_F(SimTimeTest, enforce_timeinterval_bigtimetol) +{ + build_simtime_params(); + { + amrex::ParmParse pp("time"); + pp.add("regrid_interval", -1); + pp.add("checkpoint_interval", -1); + pp.add("plot_interval", -1); + + pp.add("plot_time_interval", 0.5); + pp.add("enforce_plot_time_dt", true); + + // Choose values that could give issues + pp.add("enforce_plot_dt_reltol", 1e-8); + pp.add("plot_time_interval_reltol", 1e0); + + pp.add("stop_time", 1.0); + pp.add("max_step", 10); + } + amr_wind::SimTime time; + time.parse_parameters(); + + int counter = 0; + int plot_counter = 0; + amrex::Real plot_time_sum = 0.0; + int plot_step_sum = 0; + while (time.new_timestep()) { + time.set_current_cfl(0.45 / 0.4, 0.0, 0.0); + ++counter; + if (time.write_plot_file()) { + ++plot_counter; + plot_time_sum += time.new_time(); + plot_step_sum += counter; + } + } + EXPECT_EQ(plot_counter, 2); + EXPECT_NEAR(plot_time_sum, 1.5, 1e-8); + EXPECT_EQ(plot_step_sum, 2 + 6); +} + +TEST_F(SimTimeTest, enforce_timeinterval_bigdttol) +{ + build_simtime_params(); + { + amrex::ParmParse pp("time"); + pp.add("regrid_interval", -1); + pp.add("checkpoint_interval", -1); + pp.add("plot_interval", -1); + + pp.add("plot_time_interval", 0.5); + pp.add("enforce_plot_time_dt", true); + + // Weak enforcement of plot time interval on dt + pp.add("enforce_plot_dt_reltol", 1e0); + pp.add("plot_time_interval_reltol", 1e-8); + + pp.add("stop_time", 1.0); + pp.add("max_step", 10); + } + amr_wind::SimTime time; + time.parse_parameters(); + + int counter = 0; + int plot_counter = 0; + amrex::Real plot_time_sum = 0.0; + int plot_step_sum = 0; + while (time.new_timestep()) { + time.set_current_cfl(0.45 / 0.4, 0.0, 0.0); + ++counter; + if (time.write_plot_file()) { + ++plot_counter; + plot_time_sum += time.new_time(); + plot_step_sum += counter; + } + } + // With big dt tolerance, dt never gets shortened except at the end + // (reaching t = 1.0). Ordinary plot time interval tolerance ensures that + // plot files are still written at the first step after interval is passed. + EXPECT_EQ(plot_counter, 2); + EXPECT_NEAR(plot_time_sum, 0.8 + 1.0, 1e-8); + EXPECT_EQ(plot_step_sum, 2 + 3); +} + } // namespace amr_wind_tests