diff --git a/amr-wind/boundary_conditions/wall_models/WallFunction.H b/amr-wind/boundary_conditions/wall_models/WallFunction.H index 3ae240acb4..c195d2a475 100644 --- a/amr-wind/boundary_conditions/wall_models/WallFunction.H +++ b/amr-wind/boundary_conditions/wall_models/WallFunction.H @@ -3,7 +3,7 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/core/FieldBCOps.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" #include "amr-wind/boundary_conditions/wall_models/LogLaw.H" #include "amr-wind/boundary_conditions/wall_models/MOSD.H" @@ -40,7 +40,7 @@ private: //! LogLaw instance LogLaw m_log_law; int m_direction{2}; ///< Direction normal to wall, hardcoded to z - VelPlaneAveragingFine m_pa_vel; + VelPlaneAveraging m_pa_vel; MOSD m_mosd; }; diff --git a/amr-wind/boundary_conditions/wall_models/WallFunction.cpp b/amr-wind/boundary_conditions/wall_models/WallFunction.cpp index 24527c64eb..2d10426a48 100644 --- a/amr-wind/boundary_conditions/wall_models/WallFunction.cpp +++ b/amr-wind/boundary_conditions/wall_models/WallFunction.cpp @@ -16,7 +16,7 @@ using namespace amrex::literals; namespace amr_wind { WallFunction::WallFunction(CFDSim& sim) - : m_sim(sim), m_mesh(m_sim.mesh()), m_pa_vel(sim, m_direction) + : m_sim(sim), m_mesh(m_sim.mesh()), m_pa_vel(sim, m_direction, -1) { amrex::Real mu; amrex::Real rho{1.0_rt}; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.H b/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.H index c96c1951eb..db7aed8974 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.H @@ -7,7 +7,7 @@ #include "amr-wind/core/FieldRepo.H" #include "amr-wind/equation_systems/icns/MomentumSource.H" #include "amr-wind/core/SimTime.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" #include "amr-wind/utilities/ncutils/nc_interface.H" #include "amr-wind/wind_energy/ABLMesoscaleForcing.H" #include "amr-wind/wind_energy/ABLMesoscaleInput.H" @@ -41,13 +41,13 @@ public: void mean_velocity_init(const ABLMesoscaleInput& ncfile); void mean_velocity_init( - const VelPlaneAveragingFine& vavg, const ABLMesoscaleInput& ncfile); + const VelPlaneAveraging& vavg, const ABLMesoscaleInput& ncfile); void mean_velocity_heights(std::unique_ptr const& ncfile); void mean_velocity_heights( - const VelPlaneAveragingFine& vavg, + const VelPlaneAveraging& vavg, std::unique_ptr const& ncfile); amrex::Vector& mom_u_error() { return m_err_U; } diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.cpp b/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.cpp index c5d43a16fc..abf0a7285d 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLMesoForcingMom.cpp @@ -55,7 +55,7 @@ void ABLMesoForcingMom::mean_velocity_init(const ABLMesoscaleInput& ncfile) } void ABLMesoForcingMom::mean_velocity_init( - const VelPlaneAveragingFine& vavg, const ABLMesoscaleInput& ncfile) + const VelPlaneAveraging& vavg, const ABLMesoscaleInput& ncfile) { const int num_meso_ht = ncfile.nheights(); m_nht = vavg.ncell_line(); @@ -122,7 +122,7 @@ void ABLMesoForcingMom::mean_velocity_heights( } void ABLMesoForcingMom::mean_velocity_heights( - const VelPlaneAveragingFine& vavg, + const VelPlaneAveraging& vavg, std::unique_ptr const& ncfile) { if (m_forcing_scheme.empty()) { diff --git a/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.H b/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.H index 417029bf85..5cb6fde6f1 100644 --- a/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.H +++ b/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.H @@ -5,7 +5,7 @@ #include "amr-wind/core/FieldRepo.H" #include "amr-wind/equation_systems/temperature/TemperatureSource.H" #include "amr-wind/core/SimTime.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" #include "amr-wind/utilities/ncutils/nc_interface.H" #include "amr-wind/wind_energy/ABLMesoscaleForcing.H" #include "amr-wind/wind_energy/ABLMesoscaleInput.H" @@ -34,13 +34,13 @@ public: void mean_temperature_init(const ABLMesoscaleInput& ncfile); void mean_temperature_init( - const FieldPlaneAveragingFine& tavg, const ABLMesoscaleInput& ncfile); + const FieldPlaneAveraging& tavg, const ABLMesoscaleInput& ncfile); amrex::Real mean_temperature_heights(std::unique_ptr const& ncfile); amrex::Real mean_temperature_heights( - const FieldPlaneAveragingFine& tavg, + const FieldPlaneAveraging& tavg, std::unique_ptr const& ncfile); amrex::Vector& theta_error() { return m_err_Theta; } diff --git a/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.cpp b/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.cpp index 6fee9d4f23..cd5ee71e4b 100644 --- a/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/ABLMesoForcingTemp.cpp @@ -51,7 +51,7 @@ void ABLMesoForcingTemp::mean_temperature_init(const ABLMesoscaleInput& ncfile) ncfile.meso_heights().end(), m_meso_ht.begin()); } void ABLMesoForcingTemp::mean_temperature_init( - const FieldPlaneAveragingFine& tavg, const ABLMesoscaleInput& ncfile) + const FieldPlaneAveraging& tavg, const ABLMesoscaleInput& ncfile) { const int num_meso_ht = ncfile.nheights(); m_nht = tavg.ncell_line(); @@ -116,7 +116,7 @@ amrex::Real ABLMesoForcingTemp::mean_temperature_heights( } amrex::Real ABLMesoForcingTemp::mean_temperature_heights( - const FieldPlaneAveragingFine& tavg, + const FieldPlaneAveraging& tavg, std::unique_ptr const& ncfile) { amrex::Real currtime; diff --git a/amr-wind/turbulence/LES/AMD.H b/amr-wind/turbulence/LES/AMD.H index 98ac915329..83885b2b18 100644 --- a/amr-wind/turbulence/LES/AMD.H +++ b/amr-wind/turbulence/LES/AMD.H @@ -4,7 +4,7 @@ #include #include #include "amr-wind/turbulence/TurbModelBase.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" #include "amr-wind/core/FieldRepo.H" #include "amr-wind/fvm/stencils.H" #include "amr-wind/utilities/DirectionSelector.H" @@ -59,7 +59,8 @@ private: const Field& m_vel; const Field& m_temperature; const Field& m_rho; - FieldPlaneAveragingFine m_pa_temp; + + FieldPlaneAveraging m_pa_temp; amrex::Vector m_gravity{0.0_rt, 0.0_rt, -9.81_rt}; }; diff --git a/amr-wind/turbulence/LES/AMD.cpp b/amr-wind/turbulence/LES/AMD.cpp index 5ac9b96e03..69e438b8da 100644 --- a/amr-wind/turbulence/LES/AMD.cpp +++ b/amr-wind/turbulence/LES/AMD.cpp @@ -23,7 +23,7 @@ AMD::AMD(CFDSim& sim) , m_vel(sim.repo().get_field("velocity")) , m_temperature(sim.repo().get_field("temperature")) , m_rho(sim.repo().get_field("density")) - , m_pa_temp(m_temperature, sim.time(), m_normal_dir, true) + , m_pa_temp(m_temperature, sim.time(), m_normal_dir, -1, true) { { amrex::ParmParse pp("incflo"); diff --git a/amr-wind/utilities/CMakeLists.txt b/amr-wind/utilities/CMakeLists.txt index 6310c6db40..8cf2b60caa 100644 --- a/amr-wind/utilities/CMakeLists.txt +++ b/amr-wind/utilities/CMakeLists.txt @@ -9,7 +9,6 @@ target_sources(${amr_wind_lib_name} io_utils.cpp IOManager.cpp FieldPlaneAveraging.cpp - FieldPlaneAveragingFine.cpp SecondMomentAveraging.cpp ThirdMomentAveraging.cpp diff --git a/amr-wind/utilities/FieldPlaneAveraging.H b/amr-wind/utilities/FieldPlaneAveraging.H index 3b4e25f0b4..d42bd8e0d5 100644 --- a/amr-wind/utilities/FieldPlaneAveraging.H +++ b/amr-wind/utilities/FieldPlaneAveraging.H @@ -33,6 +33,8 @@ public: * \param field_in [in] Field to be averaged * \param time [in] Time instance to determine output frequencies * \param axis_in [in] Direction along which planes are computed + * \param max_level [in] Maximum mesh level to consider when calculating + * averages * \param compute_deriv [in] Should the derivative of the averages be * computed */ @@ -40,12 +42,15 @@ public: const FType& field_in, const amr_wind::SimTime& time, int axis_in, + int max_level, bool compute_deriv = false); virtual ~FPlaneAveraging() = default; virtual void operator()(); + void convert_x_to_ind(amrex::Real x, int& ind, amrex::Real& c) const; + /** evaluate line average at specific location for any average component */ [[nodiscard]] amrex::Real line_average_interpolated(amrex::Real x, int comp) const; @@ -72,7 +77,7 @@ public: [[nodiscard]] amrex::Real xlo() const { return m_xlo; }; [[nodiscard]] int axis() const { return m_axis; }; - [[nodiscard]] int level() const { return m_level; }; + [[nodiscard]] int level() const { return m_max_level; }; [[nodiscard]] int ncomp() const { return m_ncomp; }; [[nodiscard]] int ncell_plane() const { return m_ncell_plane; }; [[nodiscard]] int ncell_line() const { return m_ncell_line; }; @@ -110,27 +115,26 @@ protected: of each cell along a line*/ amrex::Real m_dx; /** mesh spacing in axis direction*/ - amrex::Real m_xlo; /** bottom of domain in axis direction */ + amrex::Real m_xlo; /** bottom of line */ + amrex::Real m_xhi; /** top of line */ int m_ncell_plane; /** number of cells in plane */ int m_ncell_line; /** number of cells along line */ - int m_precision = 4; /** precision for line plot text file */ - const int m_level = - 0; /** level for plane averaging for now fixed at level=0 */ + int m_precision = 6; /** precision for line plot text file */ int m_last_updated_index = -1; /** keep track of the last time index that the operator was called */ const FType& m_field; const SimTime& m_time; const int m_axis; + int m_max_level; const bool m_comp_deriv; public: // public for GPU /** fill line storage with averages */ template - void - compute_averages(const IndexSelector& idxOp, const amrex::MultiFab& mfab); + void compute_averages(const IndexSelector& idxOp); /** fill derivatives of line storage */ void compute_line_derivatives(); @@ -145,7 +149,7 @@ using ScratchFieldPlaneAveraging = FPlaneAveraging; class VelPlaneAveraging : public FieldPlaneAveraging { public: - VelPlaneAveraging(CFDSim& sim, int axis_in); + VelPlaneAveraging(CFDSim& sim, int axis_in, int max_level); ~VelPlaneAveraging() override = default; @@ -155,17 +159,19 @@ private: amrex::Vector m_line_hvelmag_average; /** line storage for the average horizontal velocity magnitude */ - //! line storage for the derivative of average horizontal velocity magnitude - amrex::Vector m_line_hvelmag_deriv; + + amrex::Vector + m_line_Su_average; /** line storage for the average horizontal + velocity magnitude time x-velocity */ + + amrex::Vector + m_line_Sv_average; /** line storage for the average horizontal + velocity magnitude time y-velocity */ public: // public for GPU /** fill line storage with horizontal velocity magnitude averages */ template - void compute_hvelmag_averages( - const IndexSelector& idx_op, - int h1_idx, - int h2_idx, - const amrex::MultiFab& mfab); + void compute_hvelmag_averages(const IndexSelector& idxOp); /** return vector containing horizontal velocity magnitude average */ const amrex::Vector& line_hvelmag_average() @@ -177,24 +183,12 @@ public: // public for GPU * magnitude */ [[nodiscard]] amrex::Real line_hvelmag_average_interpolated(amrex::Real x) const; - /** evaluate line average at specific cell for horizontal velocity magnitude - */ - [[nodiscard]] amrex::Real line_hvelmag_average_cell(int ind) const; - - /** compute derivatives of horizontal velocity magnitude */ - void compute_line_hvelmag_derivatives(); - /** evaluate line average derivative at specific location for horizontal - velocity magnitude */ - [[nodiscard]] amrex::Real - line_hvelmag_derivative_interpolated(amrex::Real x) const; - /** evaluate derivative of a line average at specific cell horizontal - velocity magnitude */ - [[nodiscard]] amrex::Real - line_hvelmag_derivative_of_average_cell(int ind) const; - - void output_line_average_ascii( - const std::string& filename, int step, amrex::Real time); - void output_line_average_ascii(int step, amrex::Real time); + /** evaluate line haverage at specific location for horizontal velocity + * magnitude times x-velocity */ + [[nodiscard]] amrex::Real line_su_average_interpolated(amrex::Real x) const; + /** evaluate line haverage at specific location for horizontal velocity + * magnitude times y-velocity */ + [[nodiscard]] amrex::Real line_sv_average_interpolated(amrex::Real x) const; }; } // namespace amr_wind diff --git a/amr-wind/utilities/FieldPlaneAveraging.cpp b/amr-wind/utilities/FieldPlaneAveraging.cpp index f8be4dfa7d..be9d1d4289 100644 --- a/amr-wind/utilities/FieldPlaneAveraging.cpp +++ b/amr-wind/utilities/FieldPlaneAveraging.cpp @@ -1,5 +1,7 @@ #include "amr-wind/utilities/FieldPlaneAveraging.H" - +#include "amr-wind/utilities/constants.H" +#include "AMReX_iMultiFab.H" +#include "AMReX_MultiFabUtil.H" #include #include "AMReX_REAL.H" @@ -12,32 +14,49 @@ FPlaneAveraging::FPlaneAveraging( const FType& field_in, const amr_wind::SimTime& time, int axis_in, + int max_level, bool compute_deriv) : m_field(field_in) , m_time(time) , m_axis(axis_in) + , m_max_level(max_level) , m_comp_deriv(compute_deriv) { AMREX_ALWAYS_ASSERT(m_axis >= 0 && m_axis < AMREX_SPACEDIM); auto geom = m_field.repo().mesh().Geom(); - // level=0 is default, could later make this an input. - // Might only makes sense for fully covered levels - m_xlo = geom[m_level].ProbLo(m_axis); - m_dx = geom[m_level].CellSize(m_axis); - m_ncomp = m_field.num_comp(); + // beginning and end of line, for now assuming line is the length of the + // entire domain + m_xlo = geom[0].ProbLo(m_axis); + m_xhi = geom[0].ProbHi(m_axis); + + const int finestLevel = + m_max_level < 0 ? m_field.repo().mesh().maxLevel() : m_max_level; + const amrex::IntVect dom_lo_vec(geom[finestLevel].Domain().smallEnd()); + const amrex::IntVect dom_hi_vec(geom[finestLevel].Domain().bigEnd()); + const int dom_lo = dom_lo_vec[m_axis]; + const int dom_hi = dom_hi_vec[m_axis]; + + AMREX_ALWAYS_ASSERT(dom_lo == 0); + const auto& mesh = m_field.repo().mesh(); + int dom_hi2 = geom[0].Domain().bigEnd()[m_axis] + 1; + for (int i = 0; i < finestLevel; ++i) { + dom_hi2 *= mesh.refRatio(i)[m_axis]; + } + AMREX_ALWAYS_ASSERT(dom_hi + 1 == dom_hi2); + + // TODO: make an input maybe? + m_ncell_line = dom_hi - dom_lo + 1; - const amrex::Box& domain = geom[m_level].Domain(); - const amrex::IntVect dom_lo(domain.loVect()); - const amrex::IntVect dom_hi(domain.hiVect()); + m_dx = (m_xhi - m_xlo) / static_cast(m_ncell_line); - m_ncell_line = dom_hi[m_axis] - dom_lo[m_axis] + 1; + m_ncomp = m_field.num_comp(); // count number of cells in plane m_ncell_plane = 1; for (int i = 0; i < AMREX_SPACEDIM; ++i) { if (i != m_axis) { - m_ncell_plane *= (dom_hi[i] - dom_lo[i] + 1); + m_ncell_plane *= (dom_hi_vec[i] - dom_lo_vec[i] + 1); } } @@ -53,6 +72,27 @@ FPlaneAveraging::FPlaneAveraging( } } +template +void FPlaneAveraging::convert_x_to_ind( + amrex::Real x, int& ind, amrex::Real& c) const +{ + c = 0.0_rt; + ind = 0; + + if (x > m_xlo + 0.5_rt * m_dx) { + ind = static_cast(floor((x - m_xlo) / m_dx - 0.5_rt)); + const amrex::Real x1 = m_xlo + (ind + 0.5_rt) * m_dx; + c = (x - x1) / m_dx; + } + + if (ind + 1 >= m_ncell_line) { + ind = m_ncell_line - 2; + c = 1.0_rt; + } + + AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line); +} + template void FPlaneAveraging::output_line_average_ascii( const std::string& filename, int step, amrex::Real time) @@ -113,21 +153,9 @@ FPlaneAveraging::line_average_interpolated(amrex::Real x, int comp) const AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp); - amrex::Real c = 0.0_rt; - int ind = 0; - - if (x > m_xlo + (0.5_rt * m_dx)) { - ind = static_cast(std::floor(((x - m_xlo) / m_dx) - 0.5_rt)); - const amrex::Real x1 = m_xlo + ((ind + 0.5_rt) * m_dx); - c = (x - x1) / m_dx; - } - - if (ind + 1 >= m_ncell_line) { - ind = m_ncell_line - 2; - c = 1.0_rt; - } - - AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line); + int ind; + amrex::Real c; + convert_x_to_ind(x, ind, c); return (m_line_average[(m_ncomp * ind) + comp] * (1.0_rt - c)) + (m_line_average[(m_ncomp * (ind + 1)) + comp] * c); @@ -157,10 +185,238 @@ amrex::Real FPlaneAveraging::line_average_cell(int ind, int comp) const return m_line_average[(m_ncomp * ind) + comp]; } +template +void FPlaneAveraging::operator()() +{ + BL_PROFILE("amr-wind::FPlaneAveraging::operator"); + + m_last_updated_index = m_time.time_index(); + + std::fill(m_line_average.begin(), m_line_average.end(), 0.0_rt); + + switch (m_axis) { + case 0: + compute_averages(XDir()); + break; + case 1: + compute_averages(YDir()); + break; + case 2: + compute_averages(ZDir()); + break; + default: + amrex::Abort("axis must be equal to 0, 1, or 2"); + break; + } + + if (m_comp_deriv) { + compute_line_derivatives(); + } +} + +template +template +void FPlaneAveraging::compute_averages(const IndexSelector& idxOp) +{ + BL_PROFILE("amr-wind::PlaneAveraging::compute_averages"); + + amrex::AsyncArray lavg( + m_line_average.data(), m_line_average.size()); + + amrex::Real* line_avg = lavg.data(); + const int num_comps = m_ncomp; + const int num_cells = m_ncell_line; + const amrex::Real line_dx = m_dx; + const amrex::Real xlo = m_xlo; + + auto g0 = m_field.repo().mesh().Geom(0); + const amrex::Real denom = + (amrex::Real)m_ncell_line / + ((g0.ProbHi(0) - g0.ProbLo(0)) * (g0.ProbHi(1) - g0.ProbLo(1)) * + (g0.ProbHi(2) - g0.ProbLo(2))); + + const bool periodic_dir = g0.periodicity().isPeriodic(m_axis); + + const auto& mesh = m_field.repo().mesh(); + const int finestLevel = m_max_level < 0 ? mesh.finestLevel() : m_max_level; + const auto dir = m_axis; + const bool no_ghost = (m_field.num_grow()[dir] == 0); + + for (int lev = 0; lev <= finestLevel; ++lev) { + + const auto& geom = m_field.repo().mesh().Geom(lev); + const amrex::Real dx = geom.CellSize()[dir]; + const amrex::Real dy = geom.CellSize()[idxOp.odir1]; + const amrex::Real dz = geom.CellSize()[idxOp.odir2]; + + const amrex::Real problo_x = geom.ProbLo(dir); + const amrex::Real probhi_x = geom.ProbHi(dir); + + amrex::iMultiFab level_mask; + if (lev < finestLevel) { + level_mask = makeFineMask( + mesh.boxArray(lev), mesh.DistributionMap(lev), + mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& mfab = m_field(lev); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); + ++mfi) { + amrex::Box bx = mfi.tilebox(); + + auto fab_arr = mfab.const_array(mfi); + auto mask_arr = level_mask.const_array(mfi); + + amrex::Box pbx = + perpendicular_box(bx, amrex::IntVect{0, 0, 0}); + + amrex::ParallelFor( + amrex::Gpu::KernelInfo().setReduction(true), pbx, + [=] AMREX_GPU_DEVICE( + int p_i, int p_j, int p_k, + amrex::Gpu::Handler const& handler) noexcept { + // Loop over the direction perpendicular to the plane. + // This reduces the atomic pressure on the destination + // arrays. + + amrex::Box lbx = parallel_box( + bx, amrex::IntVect{p_i, p_j, p_k}); + + for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { + for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { + for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); + ++i) { + + // cell coordinates + const amrex::Real cell_xlo = + xlo + idxOp(i, j, k) * dx; + const amrex::Real cell_xhi = cell_xlo + dx; + + // line indices + const int line_ind_lo = amrex::min( + amrex::max( + static_cast( + (cell_xlo - xlo) / line_dx), + 0), + num_cells - 1); + const int line_ind_hi = amrex::min( + amrex::max( + static_cast( + (cell_xhi - + amr_wind::constants::TIGHT_TOL - + xlo) / + line_dx), + 0), + num_cells - 1); + + AMREX_ASSERT(line_ind_lo >= 0); + AMREX_ASSERT(line_ind_hi >= 0); + AMREX_ASSERT(line_ind_lo < num_cells); + AMREX_ASSERT(line_ind_hi < num_cells); + + for (int ind = line_ind_lo; ind <= line_ind_hi; + ++ind) { + + // line coordinates + const amrex::Real line_xlo = + xlo + ind * line_dx; + const amrex::Real line_xhi = + line_xlo + line_dx; + + amrex::Real deltax; + + if (line_xlo <= cell_xlo) { + deltax = line_xhi - cell_xlo; + } else if (line_xhi >= cell_xhi) { + deltax = cell_xhi - line_xlo; + } else { + deltax = line_dx; + } + deltax = amrex::min(deltax, dx); + const amrex::Real vol = deltax * dy * dz; + + // Calculate location of target + const auto x_targ = + 0.5_rt * (line_xlo + line_xhi); + // Calculate location of cell center + const amrex::IntVect iv{i, j, k}; + const auto idx = iv[dir]; + const auto x_cell = + problo_x + + (static_cast(idx) + + 0.5_rt) * + dx; + // Get location of neighboring cell centers + auto x_up = x_cell + dx; + auto x_down = x_cell - dx; + // Bound locations by domain limits + if (!periodic_dir) { + x_up = amrex::min(probhi_x, x_up); + x_down = amrex::max(problo_x, x_down); + } + // Pick indices of closest neighbor + // Bound indices in case of no ghost cells + auto iv_nb = iv; + auto x_nb = x_cell; + if (std::abs(x_up - x_targ) < + std::abs(x_down - x_targ)) { + x_nb = x_up; + iv_nb[dir] += 1; + if (no_ghost) { + iv_nb[dir] = amrex::min( + iv_nb[dir], bx.bigEnd(dir)); + } + } else { + x_nb = x_down; + iv_nb[dir] -= 1; + if (no_ghost) { + iv_nb[dir] = amrex::max( + iv_nb[dir], bx.smallEnd(dir)); + } + } + // Interpolate to target location using + // closest neighbor + // (will do nothing if already at cell + // center) + for (int n = 0; n < num_comps; ++n) { + const auto f_interp = + fab_arr(iv, n) + + (fab_arr(iv_nb, n) - + fab_arr(iv, n)) * + ((x_targ - x_cell) / + (x_nb - x_cell)); + amrex::Gpu::deviceReduceSum( + &line_avg[num_comps * ind + n], + mask_arr(i, j, k) * f_interp * vol * + denom, + handler); + } + } + } + } + } + }); + } + } + + lavg.copyToHost(m_line_average.data(), m_line_average.size()); + amrex::ParallelDescriptor::ReduceRealSum( + m_line_average.data(), m_line_average.size()); +} + template void FPlaneAveraging::compute_line_derivatives() { - BL_PROFILE("amr-wind::PlaneAveraging::compute_line_derivatives"); + BL_PROFILE("amr-wind::FPlaneAveraging::compute_line_derivatives"); for (int i = 0; i < m_ncell_line; ++i) { for (int n = 0; n < m_ncomp; ++n) { m_line_deriv[(m_ncomp * i) + n] = @@ -173,7 +429,7 @@ template amrex::Real FPlaneAveraging::line_derivative_of_average_cell(int ind, int comp) const { - BL_PROFILE("amr-wind::PlaneAveraging::line_derivative_of_average_cell"); + BL_PROFILE("amr-wind::FPlaneAveraging::line_derivative_of_average_cell"); AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp); AMREX_ALWAYS_ASSERT(ind >= 0 && ind < m_ncell_line); @@ -202,134 +458,28 @@ template amrex::Real FPlaneAveraging::line_derivative_interpolated( amrex::Real x, int comp) const { - BL_PROFILE("amr-wind::PlaneAveraging::line_derivative_interpolated"); + BL_PROFILE("amr-wind::FPlaneAveraging::line_derivative_interpolated"); AMREX_ALWAYS_ASSERT(comp >= 0 && comp < m_ncomp); - amrex::Real c = 0.0_rt; - int ind = 0; - - if (x > m_xlo + (0.5_rt * m_dx)) { - ind = static_cast(std::floor(((x - m_xlo) / m_dx) - 0.5_rt)); - const amrex::Real x1 = m_xlo + ((ind + 0.5_rt) * m_dx); - c = (x - x1) / m_dx; - } - - if (ind + 1 >= m_ncell_line) { - ind = m_ncell_line - 2; - c = 1.0_rt; - } - - AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line); + int ind; + amrex::Real c; + convert_x_to_ind(x, ind, c); return (m_line_deriv[(m_ncomp * ind) + comp] * (1.0_rt - c)) + (m_line_deriv[(m_ncomp * (ind + 1)) + comp] * c); } -template -void FPlaneAveraging::operator()() -{ - BL_PROFILE("amr-wind::FPlaneAveraging::operator"); - - m_last_updated_index = m_time.time_index(); - - std::ranges::fill(m_line_average, 0.0_rt); - - switch (m_axis) { - case 0: - compute_averages(XDir(), m_field(m_level)); - break; - case 1: - compute_averages(YDir(), m_field(m_level)); - break; - case 2: - compute_averages(ZDir(), m_field(m_level)); - break; - default: - amrex::Abort("axis must be equal to 0, 1, or 2"); - break; - } - - if (m_comp_deriv) { - compute_line_derivatives(); - } -} - -template -template -void FPlaneAveraging::compute_averages( - const IndexSelector& idxOp, const amrex::MultiFab& mfab) -{ - BL_PROFILE("amr-wind::PlaneAveraging::compute_averages"); - - const amrex::Real denom = 1.0_rt / (amrex::Real)m_ncell_plane; - - amrex::AsyncArray lavg( - m_line_average.data(), m_line_average.size()); - - amrex::Real* line_avg = lavg.data(); - const int captured_ncomp = m_ncomp; - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - amrex::Box bx = mfi.tilebox(); - - auto fab_arr = mfab.const_array(mfi); - - amrex::Box pbx = - perpendicular_box(bx, amrex::IntVect{0, 0, 0}); - - amrex::ParallelFor( - amrex::Gpu::KernelInfo().setReduction(true), pbx, - [=] AMREX_GPU_DEVICE( - int p_i, int p_j, int p_k, amrex::Gpu::Handler const& handler) { - // Loop over the direction perpendicular to the plane. - // This reduces the atomic pressure on the destination arrays. - - amrex::Box lbx = parallel_box( - bx, amrex::IntVect{p_i, p_j, p_k}); - - for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { - for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { - for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); ++i) { - - const int ind = idxOp(i, j, k); - - for (int n = 0; n < captured_ncomp; ++n) { - amrex::Gpu::deviceReduceSum( - &line_avg[(captured_ncomp * ind) + n], - fab_arr(i, j, k, n) * denom, handler); - } - } - } - } - }); - } - - lavg.copyToHost(m_line_average.data(), m_line_average.size()); - amrex::ParallelDescriptor::ReduceRealSum( - m_line_average.data(), m_line_average.size()); - - // fixme later remove denom above and replace with this - // std::for_each(m_line_average.begin(), m_line_average.end(), - // [denom](amrex::Real &el){el *= denom; }); -} - template class FPlaneAveraging; template class FPlaneAveraging; -// NOLINTBEGIN(clang-analyzer-security.ArrayBound) -VelPlaneAveraging::VelPlaneAveraging(CFDSim& sim, int axis_in) +VelPlaneAveraging::VelPlaneAveraging(CFDSim& sim, int axis_in, int max_level) : FieldPlaneAveraging( - sim.repo().get_field("velocity"), sim.time(), axis_in, true) + sim.repo().get_field("velocity"), sim.time(), axis_in, max_level) { m_line_hvelmag_average.resize(m_ncell_line, 0.0_rt); - if (m_comp_deriv) { - m_line_hvelmag_deriv.resize(m_ncell_line, 0.0_rt); - } + m_line_Su_average.resize(m_ncell_line, 0.0_rt); + m_line_Sv_average.resize(m_ncell_line, 0.0_rt); } // NOLINTEND(clang-analyzer-security.ArrayBound) @@ -338,209 +488,227 @@ void VelPlaneAveraging::operator()() BL_PROFILE("amr-wind::VelPlaneAveraging::operator"); + // velocity averages FieldPlaneAveraging::operator()(); - std::ranges::fill(m_line_hvelmag_average, 0.0_rt); + std::fill( + m_line_hvelmag_average.begin(), m_line_hvelmag_average.end(), 0.0_rt); + std::fill(m_line_Su_average.begin(), m_line_Su_average.end(), 0.0_rt); + std::fill(m_line_Sv_average.begin(), m_line_Sv_average.end(), 0.0_rt); switch (m_axis) { case 0: - compute_hvelmag_averages(XDir(), 1, 2, m_field(m_level)); + compute_hvelmag_averages(XDir()); break; case 1: - compute_hvelmag_averages(YDir(), 0, 2, m_field(m_level)); + compute_hvelmag_averages(YDir()); break; case 2: - compute_hvelmag_averages(ZDir(), 0, 1, m_field(m_level)); + compute_hvelmag_averages(ZDir()); break; default: amrex::Abort("axis must be equal to 0, 1, or 2"); break; } - - if (m_comp_deriv) { - compute_line_hvelmag_derivatives(); - } } template -void VelPlaneAveraging::compute_hvelmag_averages( - const IndexSelector& idx_op, - const int h1_idx, - const int h2_idx, - const amrex::MultiFab& mfab) +void VelPlaneAveraging::compute_hvelmag_averages(const IndexSelector& idxOp) { - BL_PROFILE("amr-wind::VelPlaneAveraging::compute_hvelmag_averages"); - const amrex::Real denom = 1.0_rt / (amrex::Real)m_ncell_plane; + BL_PROFILE("amr-wind::VelPlaneAveraging::compute_hvelmag_averages"); - amrex::AsyncArray lavg( + amrex::AsyncArray lavg_vm( m_line_hvelmag_average.data(), m_line_hvelmag_average.size()); - amrex::Real* line_avg = lavg.data(); + amrex::AsyncArray lavg_Su( + m_line_Su_average.data(), m_line_Su_average.size()); + amrex::AsyncArray lavg_Sv( + m_line_Sv_average.data(), m_line_Sv_average.size()); + + amrex::Real* line_avg_vm = lavg_vm.data(); + amrex::Real* line_avg_Su = lavg_Su.data(); + amrex::Real* line_avg_Sv = lavg_Sv.data(); + + const int num_cells = m_ncell_line; + const amrex::Real line_dx = m_dx; + const amrex::Real xlo = m_xlo; + + auto g0 = m_field.repo().mesh().Geom(0); + const amrex::Real denom = + (amrex::Real)m_ncell_line / + ((g0.ProbHi(0) - g0.ProbLo(0)) * (g0.ProbHi(1) - g0.ProbLo(1)) * + (g0.ProbHi(2) - g0.ProbLo(2))); + + const auto& mesh = m_field.repo().mesh(); + const int finestLevel = m_max_level < 0 ? mesh.finestLevel() : m_max_level; + + for (int lev = 0; lev <= finestLevel; ++lev) { + + const auto& geom = m_field.repo().mesh().Geom(lev); + const amrex::Real dx = geom.CellSize()[m_axis]; + const amrex::Real dy = geom.CellSize()[idxOp.odir1]; + const amrex::Real dz = geom.CellSize()[idxOp.odir2]; + + amrex::iMultiFab level_mask; + if (lev < finestLevel) { + level_mask = makeFineMask( + mesh.boxArray(lev), mesh.DistributionMap(lev), + mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& mfab = m_field(lev); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - amrex::Box bx = mfi.tilebox(); - - auto fab_arr = mfab.const_array(mfi); - - amrex::Box pbx = - perpendicular_box(bx, amrex::IntVect{0, 0, 0}); - - amrex::ParallelFor( - amrex::Gpu::KernelInfo().setReduction(true), pbx, - [=] AMREX_GPU_DEVICE( - int p_i, int p_j, int p_k, amrex::Gpu::Handler const& handler) { - // Loop over the direction perpendicular to the plane. - // This reduces the atomic pressure on the destination arrays. - - amrex::Box lbx = parallel_box( - bx, amrex::IntVect{p_i, p_j, p_k}); - - for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { - for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { - for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); ++i) { - - const int ind = idx_op(i, j, k); - const amrex::Real hvelmag = std::sqrt( - (fab_arr(i, j, k, h1_idx) * - fab_arr(i, j, k, h1_idx)) + - (fab_arr(i, j, k, h2_idx) * - fab_arr(i, j, k, h2_idx))); - amrex::Gpu::deviceReduceSum( - &line_avg[ind], hvelmag * denom, handler); + for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); + ++mfi) { + amrex::Box bx = mfi.tilebox(); + + auto fab_arr = mfab.const_array(mfi); + auto mask_arr = level_mask.const_array(mfi); + + amrex::Box pbx = + perpendicular_box(bx, amrex::IntVect{0, 0, 0}); + + amrex::ParallelFor( + amrex::Gpu::KernelInfo().setReduction(true), pbx, + [=] AMREX_GPU_DEVICE( + int p_i, int p_j, int p_k, + amrex::Gpu::Handler const& handler) noexcept { + // Loop over the direction perpendicular to the plane. + // This reduces the atomic pressure on the destination + // arrays. + + amrex::Box lbx = parallel_box( + bx, amrex::IntVect{p_i, p_j, p_k}); + + for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { + for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { + for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); + ++i) { + + // cell coordinates + const amrex::Real cell_xlo = + xlo + idxOp(i, j, k) * dx; + const amrex::Real cell_xhi = cell_xlo + dx; + + // line indices + const int line_ind_lo = amrex::min( + amrex::max( + static_cast( + (cell_xlo - xlo) / line_dx), + 0), + num_cells - 1); + const int line_ind_hi = amrex::min( + amrex::max( + static_cast( + (cell_xhi - xlo) / line_dx), + 0), + num_cells - 1); + + AMREX_ASSERT(line_ind_lo >= 0); + AMREX_ASSERT(line_ind_hi >= 0); + AMREX_ASSERT(line_ind_lo < num_cells); + AMREX_ASSERT(line_ind_hi < num_cells); + + for (int ind = line_ind_lo; ind <= line_ind_hi; + ++ind) { + + // line coordinates + const amrex::Real line_xlo = + xlo + ind * line_dx; + const amrex::Real line_xhi = + line_xlo + line_dx; + + amrex::Real deltax; + + if (line_xlo <= cell_xlo) { + deltax = line_xhi - cell_xlo; + } else if (line_xhi >= cell_xhi) { + deltax = cell_xhi - line_xlo; + } else { + deltax = line_dx; + } + + deltax = amrex::min(deltax, dx); + const amrex::Real vol = + mask_arr(i, j, k) * deltax * dy * dz; + + const amrex::Real hvelmag = std::sqrt( + fab_arr(i, j, k, idxOp.odir1) * + fab_arr(i, j, k, idxOp.odir1) + + fab_arr(i, j, k, idxOp.odir2) * + fab_arr(i, j, k, idxOp.odir2)); + const amrex::Real Su = + hvelmag * fab_arr(i, j, k, idxOp.odir1); + const amrex::Real Sv = + hvelmag * fab_arr(i, j, k, idxOp.odir2); + + amrex::Gpu::deviceReduceSum( + &line_avg_vm[ind], + hvelmag * vol * denom, handler); + amrex::Gpu::deviceReduceSum( + &line_avg_Su[ind], Su * vol * denom, + handler); + amrex::Gpu::deviceReduceSum( + &line_avg_Sv[ind], Sv * vol * denom, + handler); + } + } } } - } - }); + }); + } } - lavg.copyToHost( + lavg_vm.copyToHost( m_line_hvelmag_average.data(), m_line_hvelmag_average.size()); + lavg_Su.copyToHost(m_line_Su_average.data(), m_line_Su_average.size()); + lavg_Sv.copyToHost(m_line_Sv_average.data(), m_line_Sv_average.size()); amrex::ParallelDescriptor::ReduceRealSum( m_line_hvelmag_average.data(), static_cast(m_line_hvelmag_average.size())); -} - -void VelPlaneAveraging::compute_line_hvelmag_derivatives() -{ - BL_PROFILE("amr-wind::VelPlaneAveraging::compute_line_hvelmag_derivatives"); - for (int i = 0; i < m_ncell_line; ++i) { - m_line_hvelmag_deriv[i] = line_hvelmag_derivative_of_average_cell(i); - } -} - -amrex::Real -VelPlaneAveraging::line_hvelmag_derivative_of_average_cell(int ind) const -{ - BL_PROFILE("amr-wind::VelPlaneAveraging::line_derivative_of_average_cell"); - - AMREX_ALWAYS_ASSERT(ind >= 0 && ind < m_ncell_line); - - amrex::Real dudx; - - if (ind == 0) { - dudx = - (m_line_hvelmag_average[(ind + 1)] - m_line_hvelmag_average[ind]) / - m_dx; - } else if (ind == m_ncell_line - 1) { - dudx = (m_line_hvelmag_average[ind] - m_line_hvelmag_average[ind - 1]) / - m_dx; - } else { - dudx = 0.5_rt * - (m_line_hvelmag_average[ind + 1] - - m_line_hvelmag_average[ind - 1]) / - m_dx; - } - - return dudx; + amrex::ParallelDescriptor::ReduceRealSum( + m_line_Su_average.data(), static_cast(m_line_Su_average.size())); + amrex::ParallelDescriptor::ReduceRealSum( + m_line_Sv_average.data(), static_cast(m_line_Sv_average.size())); } amrex::Real VelPlaneAveraging::line_hvelmag_average_interpolated(amrex::Real x) const { - - BL_PROFILE("amr-wind::PlaneAveraging::line_average_interpolated"); - - amrex::Real c = 0.0_rt; - int ind = 0; - - if (x > m_xlo + (0.5_rt * m_dx)) { - ind = static_cast(std::floor(((x - m_xlo) / m_dx) - 0.5_rt)); - const amrex::Real x1 = m_xlo + ((ind + 0.5_rt) * m_dx); - c = (x - x1) / m_dx; - } - - if (ind + 1 >= m_ncell_line) { - ind = m_ncell_line - 2; - c = 1.0_rt; - } - - AMREX_ALWAYS_ASSERT(ind >= 0 && ind + 1 < m_ncell_line); + int ind; + amrex::Real c; + convert_x_to_ind(x, ind, c); return (m_line_hvelmag_average[ind] * (1.0_rt - c)) + (m_line_hvelmag_average[ind + 1] * c); } -amrex::Real VelPlaneAveraging::line_hvelmag_average_cell(int ind) const +amrex::Real VelPlaneAveraging::line_su_average_interpolated(amrex::Real x) const { - BL_PROFILE("amr-wind::VelPlaneAveraging::line_hvelmag_average_cell"); - - AMREX_ALWAYS_ASSERT(ind >= 0 && ind < m_ncell_line); + int ind; + amrex::Real c; + convert_x_to_ind(x, ind, c); - return m_line_average[ind]; + return m_line_Su_average[ind] * (1.0_rt - c) + + m_line_Su_average[ind + 1] * c; } -void VelPlaneAveraging::output_line_average_ascii( - const std::string& filename, int step, amrex::Real time) +amrex::Real VelPlaneAveraging::line_sv_average_interpolated(amrex::Real x) const { - BL_PROFILE("amr-wind::VelPlaneAveraging::output_line_average_ascii"); - - if (step != m_last_updated_index) { - operator()(); - } - - if (!amrex::ParallelDescriptor::IOProcessor()) { - return; - } - - std::ofstream outfile; - outfile.precision(m_precision); - - if (step == 1) { - // make new file - outfile.open(filename.c_str(), std::ios_base::out); - outfile << "#ncell,ncomp" << '\n'; - outfile << m_ncell_line << ", " << m_ncomp + 4 << '\n'; - outfile << "#step,time,z"; - for (int i = 0; i < m_ncomp; ++i) { - outfile << ",<" + m_field.name() + std::to_string(i) + ">"; - } - outfile << ", "; - outfile << '\n'; - } else { - // append file - outfile.open(filename.c_str(), std::ios_base::out | std::ios_base::app); - } + int ind; + amrex::Real c; + convert_x_to_ind(x, ind, c); - for (int i = 0; i < m_ncell_line; ++i) { - outfile << step << ", " << std::scientific << time << ", " - << m_line_xcentroid[i]; - for (int n = 0; n < m_ncomp; ++n) { - outfile << ", " << std::scientific - << m_line_average[(m_ncomp * i) + n]; - } - outfile << ", " << std::scientific << m_line_hvelmag_average[i]; - outfile << '\n'; - } -} - -void VelPlaneAveraging::output_line_average_ascii(int step, amrex::Real time) -{ - const std::string filename = "plane_average_" + m_field.name() + ".txt"; - output_line_average_ascii(filename, step, time); + return m_line_Sv_average[ind] * (1.0_rt - c) + + m_line_Sv_average[ind + 1] * c; } } // namespace amr_wind diff --git a/amr-wind/utilities/FieldPlaneAveragingFine.H b/amr-wind/utilities/FieldPlaneAveragingFine.H deleted file mode 100644 index cea8596122..0000000000 --- a/amr-wind/utilities/FieldPlaneAveragingFine.H +++ /dev/null @@ -1,183 +0,0 @@ -#ifndef FIELDPLANEAVERAGINGFINE_H -#define FIELDPLANEAVERAGINGFINE_H - -#include "amr-wind/utilities/DirectionSelector.H" -#include "amr-wind/CFDSim.H" -#include "amr-wind/core/Field.H" -#include "amr-wind/core/SimTime.H" - -/** - * \defgroup statistics Field statistics - * Field statistics - * - * This group contains utilities for performing turbulence averaging and - * outputting statistics during wind simulations. - * - * \ingroup utilities - */ - -namespace amr_wind { - -/** Output average of a field on planes normal to a given direction - * \ingroup statistics we_abl - * - * The user can choose a direction (x, y, or z), the default value is the - * z-direction. The field is then averaged on planes at the cell-centers at - * level 0 (coarsest level) along the specified direction. - */ - -template -class FPlaneAveragingFine -{ -public: - /** - * \param field_in [in] Field to be averaged - * \param time [in] Time instance to determine output frequencies - * \param axis_in [in] Direction along which planes are computed - * \param compute_deriv [in] Should the derivative of the averages be - * computed - */ - FPlaneAveragingFine( - const FType& field_in, - const amr_wind::SimTime& time, - int axis_in, - bool compute_deriv = false); - - virtual ~FPlaneAveragingFine() = default; - - virtual void operator()(); - - void convert_x_to_ind(amrex::Real x, int& ind, amrex::Real& c) const; - - /** evaluate line average at specific location for any average component */ - [[nodiscard]] amrex::Real - line_average_interpolated(amrex::Real x, int comp) const; - /** evaluate line average at specific cell for any average component */ - [[nodiscard]] amrex::Real line_average_cell(int ind, int comp) const; - - /** evaluate line average derivative at specific location for any average - component */ - [[nodiscard]] amrex::Real - line_derivative_interpolated(amrex::Real x, int comp) const; - /** evaluate derivative of a line average at specific cell for any component - */ - [[nodiscard]] amrex::Real - line_derivative_of_average_cell(int ind, int comp) const; - - void output_line_average_ascii( - const std::string& filename, int step, amrex::Real time); - void output_line_average_ascii(int step, amrex::Real time); - - /** change precision of text file output */ - void set_precision(int p) { m_precision = p; }; - - [[nodiscard]] amrex::Real dx() const { return m_dx; }; - [[nodiscard]] amrex::Real xlo() const { return m_xlo; }; - - [[nodiscard]] int axis() const { return m_axis; }; - [[nodiscard]] int ncomp() const { return m_ncomp; }; - [[nodiscard]] int ncell_line() const { return m_ncell_line; }; - [[nodiscard]] int last_updated_index() const - { - return m_last_updated_index; - }; - - [[nodiscard]] const amrex::Vector& line_average() const - { - return m_line_average; - }; - [[nodiscard]] const amrex::Vector& line_deriv() const - { - return m_line_deriv; - }; - void line_average(int comp, amrex::Vector& l_vec); - [[nodiscard]] const amrex::Vector& line_centroids() const - { - return m_line_xcentroid; - }; - - [[nodiscard]] const FType& field() const { return m_field; }; - -protected: - int m_ncomp; /** number of average components */ - - amrex::Vector - m_line_average; /** line storage for the average velocity and tracer - variables */ - //! line storage for the derivative of average field - amrex::Vector m_line_deriv; - - amrex::Vector m_line_xcentroid; /** line storage for centroids - of each cell along a line*/ - - amrex::Real m_dx; /** mesh spacing in axis direction*/ - amrex::Real m_xlo; /** bottom of line */ - amrex::Real m_xhi; /** top of line */ - - int m_ncell_line; /** number of cells along line */ - - int m_precision = 6; /** precision for line plot text file */ - int m_last_updated_index = -1; /** keep track of the last time index that - the operator was called */ - - const FType& m_field; - const SimTime& m_time; - const int m_axis; - const bool m_comp_deriv; - -public: // public for GPU - /** fill line storage with averages */ - template - void compute_averages(const IndexSelector& idxOp); - - /** fill derivatives of line storage */ - void compute_line_derivatives(); -}; - -using FieldPlaneAveragingFine = FPlaneAveragingFine; -using ScratchFieldPlaneAveragingFine = FPlaneAveragingFine; - -/** Specific application of FieldPlaneAveraging to the velocity field - * \ingroup statistics we_abl - */ -class VelPlaneAveragingFine : public FieldPlaneAveragingFine -{ -public: - VelPlaneAveragingFine(CFDSim& sim, int axis_in); - - ~VelPlaneAveragingFine() override = default; - - void operator()() override; - -private: - amrex::Vector - m_line_hvelmag_average; /** line storage for the average horizontal - velocity magnitude */ - amrex::Vector - m_line_Su_average; /** line storage for the average horizontal - velocity magnitude time x-velocity */ - - amrex::Vector - m_line_Sv_average; /** line storage for the average horizontal - velocity magnitude time y-velocity */ - -public: // public for GPU - /** fill line storage with horizontal velocity magnitude averages */ - template - void compute_hvelmag_averages(const IndexSelector& idxOp); - - /** evaluate line haverage at specific location for horizontal velocity - * magnitude */ - [[nodiscard]] amrex::Real - line_hvelmag_average_interpolated(amrex::Real x) const; - /** evaluate line haverage at specific location for horizontal velocity - * magnitude times x-velocity */ - [[nodiscard]] amrex::Real line_su_average_interpolated(amrex::Real x) const; - /** evaluate line haverage at specific location for horizontal velocity - * magnitude times y-velocity */ - [[nodiscard]] amrex::Real line_sv_average_interpolated(amrex::Real x) const; -}; - -} // namespace amr_wind - -#endif /* FieldPlaneAveragingFine_H */ diff --git a/amr-wind/utilities/FieldPlaneAveragingFine.cpp b/amr-wind/utilities/FieldPlaneAveragingFine.cpp deleted file mode 100644 index 1faa6cc49f..0000000000 --- a/amr-wind/utilities/FieldPlaneAveragingFine.cpp +++ /dev/null @@ -1,694 +0,0 @@ -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" -#include "amr-wind/utilities/constants.H" -#include "AMReX_iMultiFab.H" -#include "AMReX_MultiFabUtil.H" -#include "AMReX_REAL.H" -#include - -using namespace amrex::literals; - -namespace amr_wind { - -template -FPlaneAveragingFine::FPlaneAveragingFine( - const FType& field_in, - const amr_wind::SimTime& time, - int axis_in, - bool compute_deriv) - : m_field(field_in) - , m_time(time) - , m_axis(axis_in) - , m_comp_deriv(compute_deriv) -{ - AMREX_ALWAYS_ASSERT((m_axis >= 0) && (m_axis < AMREX_SPACEDIM)); - auto geom = m_field.repo().mesh().Geom(); - - // beginning and end of line, for now assuming line is the length of the - // entire domain - m_xlo = geom[0].ProbLo(m_axis); - m_xhi = geom[0].ProbHi(m_axis); - - const int finestLevel = m_field.repo().mesh().maxLevel(); - const int dom_lo = geom[finestLevel].Domain().smallEnd()[m_axis]; - const int dom_hi = geom[finestLevel].Domain().bigEnd()[m_axis]; - - AMREX_ALWAYS_ASSERT(dom_lo == 0); - const auto& mesh = m_field.repo().mesh(); - int dom_hi2 = geom[0].Domain().bigEnd()[m_axis] + 1; - for (int i = 0; i < finestLevel; ++i) { - dom_hi2 *= mesh.refRatio(i)[m_axis]; - } - AMREX_ALWAYS_ASSERT(dom_hi + 1 == dom_hi2); - - // TODO: make an input maybe? - m_ncell_line = dom_hi - dom_lo + 1; - - m_dx = (m_xhi - m_xlo) / static_cast(m_ncell_line); - - m_ncomp = m_field.num_comp(); - - m_line_average.resize(static_cast(m_ncell_line) * m_ncomp, 0.0_rt); - if (m_comp_deriv) { - m_line_deriv.resize( - static_cast(m_ncell_line) * m_ncomp, 0.0_rt); - } - m_line_xcentroid.resize(m_ncell_line); - - for (int i = 0; i < m_ncell_line; ++i) { - m_line_xcentroid[i] = m_xlo + ((i + 0.5_rt) * m_dx); - } -} - -template -void FPlaneAveragingFine::convert_x_to_ind( - amrex::Real x, int& ind, amrex::Real& c) const -{ - c = 0.0_rt; - ind = 0; - - if (x > m_xlo + (0.5_rt * m_dx)) { - ind = static_cast(std::floor(((x - m_xlo) / m_dx) - 0.5_rt)); - const amrex::Real x1 = m_xlo + ((ind + 0.5_rt) * m_dx); - c = (x - x1) / m_dx; - } - - if (ind + 1 >= m_ncell_line) { - ind = m_ncell_line - 2; - c = 1.0_rt; - } - - AMREX_ALWAYS_ASSERT((ind >= 0) && ((ind + 1) < m_ncell_line)); -} - -template -void FPlaneAveragingFine::output_line_average_ascii( - const std::string& filename, int step, amrex::Real time) -{ - BL_PROFILE("amr-wind::FPlaneAveragingFine::output_line_average_ascii"); - - if (step != m_last_updated_index) { - operator()(); - } - - if (!amrex::ParallelDescriptor::IOProcessor()) { - return; - } - - std::ofstream outfile; - outfile.precision(m_precision); - - if (step == 1) { - // make new file - outfile.open(filename.c_str(), std::ios_base::out); - outfile << "#ncell,ncomp" << '\n'; - outfile << m_ncell_line << ", " << m_ncomp + 3 << '\n'; - outfile << "#step,time,z"; - for (int i = 0; i < m_ncomp; ++i) { - outfile << ",<" + m_field.name() + std::to_string(i) + ">"; - } - outfile << '\n'; - } else { - // append file - outfile.open(filename.c_str(), std::ios_base::out | std::ios_base::app); - } - - for (int i = 0; i < m_ncell_line; ++i) { - outfile << step << ", " << std::scientific << time << ", " - << m_line_xcentroid[i]; - for (int n = 0; n < m_ncomp; ++n) { - outfile << ", " << std::scientific - << m_line_average[(m_ncomp * i) + n]; - } - outfile << '\n'; - } -} - -template -void FPlaneAveragingFine::output_line_average_ascii( - int step, amrex::Real time) -{ - const std::string filename = "plane_average_" + m_field.name() + ".txt"; - output_line_average_ascii(filename, step, time); -} - -template -amrex::Real FPlaneAveragingFine::line_average_interpolated( - amrex::Real x, int comp) const -{ - BL_PROFILE("amr-wind::PlaneAveragingFine::line_average_interpolated"); - - AMREX_ALWAYS_ASSERT((comp >= 0) && (comp < m_ncomp)); - - int ind; - amrex::Real c; - convert_x_to_ind(x, ind, c); - - return (m_line_average[(m_ncomp * ind) + comp] * (1.0_rt - c)) + - (m_line_average[(m_ncomp * (ind + 1)) + comp] * c); -} - -template -void FPlaneAveragingFine::line_average( - int comp, amrex::Vector& l_vec) -{ - BL_PROFILE("amr-wind::PlaneAveragingFine::line_average"); - - AMREX_ALWAYS_ASSERT((comp >= 0) && (comp < m_ncomp)); - - for (int i = 0; i < m_ncell_line; i++) { - l_vec[i] = m_line_average[(m_ncomp * i) + comp]; - } -} - -template -amrex::Real -FPlaneAveragingFine::line_average_cell(int ind, int comp) const -{ - BL_PROFILE("amr-wind::PlaneAveragingFine::line_average_cell"); - - AMREX_ALWAYS_ASSERT((comp >= 0) && (comp < m_ncomp)); - AMREX_ALWAYS_ASSERT((ind >= 0) && (ind < m_ncell_line)); - - return m_line_average[(m_ncomp * ind) + comp]; -} - -template -void FPlaneAveragingFine::operator()() -{ - BL_PROFILE("amr-wind::FPlaneAveragingFine::operator"); - - m_last_updated_index = m_time.time_index(); - - std::ranges::fill(m_line_average, 0.0_rt); - - switch (m_axis) { - case 0: - compute_averages(XDir()); - break; - case 1: - compute_averages(YDir()); - break; - case 2: - compute_averages(ZDir()); - break; - default: - amrex::Abort("axis must be equal to 0, 1, or 2"); - break; - } - - if (m_comp_deriv) { - compute_line_derivatives(); - } -} - -template -template -void FPlaneAveragingFine::compute_averages(const IndexSelector& idxOp) -{ - BL_PROFILE("amr-wind::PlaneAveragingFine::compute_averages"); - - amrex::AsyncArray lavg( - m_line_average.data(), m_line_average.size()); - - amrex::Real* line_avg = lavg.data(); - const int num_comps = m_ncomp; - const int num_cells = m_ncell_line; - const amrex::Real line_dx = m_dx; - const amrex::Real xlo = m_xlo; - - auto g0 = m_field.repo().mesh().Geom(0); - const amrex::Real denom = - (amrex::Real)m_ncell_line / - ((g0.ProbHi(0) - g0.ProbLo(0)) * (g0.ProbHi(1) - g0.ProbLo(1)) * - (g0.ProbHi(2) - g0.ProbLo(2))); - - const bool periodic_dir = g0.periodicity().isPeriodic(m_axis); - - const auto& mesh = m_field.repo().mesh(); - const int finestLevel = mesh.finestLevel(); - - for (int lev = 0; lev <= finestLevel; ++lev) { - - const auto& geom = m_field.repo().mesh().Geom(lev); - const amrex::Real dx = geom.CellSize()[m_axis]; - const amrex::Real dy = geom.CellSize()[idxOp.odir1]; - const amrex::Real dz = geom.CellSize()[idxOp.odir2]; - - const amrex::Real problo_x = geom.ProbLo(m_axis); - const amrex::Real probhi_x = geom.ProbHi(m_axis); - - const auto dir = m_axis; - - amrex::iMultiFab level_mask; - if (lev < finestLevel) { - level_mask = makeFineMask( - mesh.boxArray(lev), mesh.DistributionMap(lev), - mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); - } else { - level_mask.define( - mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, - amrex::MFInfo()); - level_mask.setVal(1); - } - - const auto& mfab = m_field(lev); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - amrex::Box bx = mfi.tilebox(); - - auto fab_arr = mfab.const_array(mfi); - auto mask_arr = level_mask.const_array(mfi); - - amrex::Box pbx = - perpendicular_box(bx, amrex::IntVect{0, 0, 0}); - - amrex::ParallelFor( - amrex::Gpu::KernelInfo().setReduction(true), pbx, - [=] AMREX_GPU_DEVICE( - int p_i, int p_j, int p_k, - amrex::Gpu::Handler const& handler) { - // Loop over the direction perpendicular to the plane. - // This reduces the atomic pressure on the destination - // arrays. - - amrex::Box lbx = parallel_box( - bx, amrex::IntVect{p_i, p_j, p_k}); - - for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { - for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { - for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); - ++i) { - - // cell coordinates - const amrex::Real cell_xlo = - xlo + (idxOp(i, j, k) * dx); - const amrex::Real cell_xhi = cell_xlo + dx; - - // line indices - const int line_ind_lo = amrex::min( - amrex::max( - static_cast( - (cell_xlo - xlo) / line_dx), - 0), - num_cells - 1); - const int line_ind_hi = amrex::min( - amrex::max( - static_cast( - (cell_xhi - - amr_wind::constants::TIGHT_TOL - - xlo) / - line_dx), - 0), - num_cells - 1); - - AMREX_ASSERT(line_ind_lo >= 0); - AMREX_ASSERT(line_ind_hi >= 0); - AMREX_ASSERT(line_ind_lo < num_cells); - AMREX_ASSERT(line_ind_hi < num_cells); - - for (int ind = line_ind_lo; ind <= line_ind_hi; - ++ind) { - - // line coordinates - const amrex::Real line_xlo = - xlo + (ind * line_dx); - const amrex::Real line_xhi = - line_xlo + line_dx; - - amrex::Real deltax; - - if (line_xlo <= cell_xlo) { - deltax = line_xhi - cell_xlo; - } else if (line_xhi >= cell_xhi) { - deltax = cell_xhi - line_xlo; - } else { - deltax = line_dx; - } - deltax = amrex::min(deltax, dx); - const amrex::Real vol = deltax * dy * dz; - - // Calculate location of target - const auto x_targ = - 0.5_rt * (line_xlo + line_xhi); - // Calculate location of cell center - const amrex::IntVect iv{i, j, k}; - const auto idx = iv[dir]; - const auto x_cell = - problo_x + ((idx + 0.5_rt) * dx); - // Get location of neighboring cell centers - auto x_up = x_cell + dx; - auto x_down = x_cell - dx; - // Bound locations by domain limits - if (!periodic_dir) { - x_up = amrex::min( - probhi_x, x_up); - x_down = amrex::max( - problo_x, x_down); - } - // Pick indices of closest neighbor - auto iv_nb = iv; - auto x_nb = x_cell; - if (std::abs(x_up - x_targ) < - std::abs(x_down - x_targ)) { - x_nb = x_up; - iv_nb[dir] += 1; - } else { - x_nb = x_down; - iv_nb[dir] -= 1; - } - // Interpolate to target location using - // closest neighbor - // (will do nothing if already at cell - // center) - for (int n = 0; n < num_comps; ++n) { - const auto f_interp = - fab_arr(iv, n) + - ((fab_arr(iv_nb, n) - - fab_arr(iv, n)) * - ((x_targ - x_cell) / - (x_nb - x_cell))); - amrex::Gpu::deviceReduceSum( - &line_avg[(num_comps * ind) + n], - mask_arr(i, j, k) * f_interp * vol * - denom, - handler); - } - } - } - } - } - }); - } - } - - lavg.copyToHost(m_line_average.data(), m_line_average.size()); - amrex::ParallelDescriptor::ReduceRealSum( - m_line_average.data(), m_line_average.size()); -} - -template -void FPlaneAveragingFine::compute_line_derivatives() -{ - BL_PROFILE("amr-wind::FPlaneAveragingFine::compute_line_derivatives"); - for (int i = 0; i < m_ncell_line; ++i) { - for (int n = 0; n < m_ncomp; ++n) { - m_line_deriv[(m_ncomp * i) + n] = - line_derivative_of_average_cell(i, n); - } - } -} - -template -amrex::Real FPlaneAveragingFine::line_derivative_of_average_cell( - int ind, int comp) const -{ - BL_PROFILE( - "amr-wind::FPlaneAveragingFine::line_derivative_of_average_cell"); - - AMREX_ALWAYS_ASSERT((comp >= 0) && (comp < m_ncomp)); - AMREX_ALWAYS_ASSERT((ind >= 0) && (ind < m_ncell_line)); - - amrex::Real dudx; - - if (ind == 0) { - dudx = (m_line_average[(m_ncomp * (ind + 1)) + comp] - - m_line_average[(m_ncomp * ind) + comp]) / - m_dx; - } else if (ind == m_ncell_line - 1) { - dudx = (m_line_average[(m_ncomp * (ind)) + comp] - - m_line_average[(m_ncomp * (ind - 1)) + comp]) / - m_dx; - } else { - dudx = 0.5_rt * - (m_line_average[(m_ncomp * (ind + 1)) + comp] - - m_line_average[(m_ncomp * (ind - 1)) + comp]) / - m_dx; - } - - return dudx; -} - -template -amrex::Real FPlaneAveragingFine::line_derivative_interpolated( - amrex::Real x, int comp) const -{ - BL_PROFILE("amr-wind::FPlaneAveragingFine::line_derivative_interpolated"); - - AMREX_ALWAYS_ASSERT((comp >= 0) && (comp < m_ncomp)); - - int ind; - amrex::Real c; - convert_x_to_ind(x, ind, c); - - return (m_line_deriv[(m_ncomp * ind) + comp] * (1.0_rt - c)) + - (m_line_deriv[(m_ncomp * (ind + 1)) + comp] * c); -} - -template class FPlaneAveragingFine; -template class FPlaneAveragingFine; - -// NOLINTBEGIN(clang-analyzer-security.ArrayBound) -VelPlaneAveragingFine::VelPlaneAveragingFine(CFDSim& sim, int axis_in) - : FieldPlaneAveragingFine( - sim.repo().get_field("velocity"), sim.time(), axis_in) -{ - m_line_hvelmag_average.resize(m_ncell_line, 0.0_rt); - m_line_Su_average.resize(m_ncell_line, 0.0_rt); - m_line_Sv_average.resize(m_ncell_line, 0.0_rt); -} -// NOLINTEND(clang-analyzer-security.ArrayBound) - -void VelPlaneAveragingFine::operator()() -{ - - BL_PROFILE("amr-wind::VelPlaneAveraging::operator"); - - // velocity averages - FieldPlaneAveragingFine::operator()(); - - std::ranges::fill(m_line_hvelmag_average, 0.0_rt); - std::ranges::fill(m_line_Su_average, 0.0_rt); - std::ranges::fill(m_line_Sv_average, 0.0_rt); - - switch (m_axis) { - case 0: - compute_hvelmag_averages(XDir()); - break; - case 1: - compute_hvelmag_averages(YDir()); - break; - case 2: - compute_hvelmag_averages(ZDir()); - break; - default: - amrex::Abort("axis must be equal to 0, 1, or 2"); - break; - } -} - -template -void VelPlaneAveragingFine::compute_hvelmag_averages(const IndexSelector& idxOp) -{ - - BL_PROFILE("amr-wind::VelPlaneAveragingFine::compute_hvelmag_averages"); - - amrex::AsyncArray lavg_vm( - m_line_hvelmag_average.data(), m_line_hvelmag_average.size()); - amrex::AsyncArray lavg_Su( - m_line_Su_average.data(), m_line_Su_average.size()); - amrex::AsyncArray lavg_Sv( - m_line_Sv_average.data(), m_line_Sv_average.size()); - - amrex::Real* line_avg_vm = lavg_vm.data(); - amrex::Real* line_avg_Su = lavg_Su.data(); - amrex::Real* line_avg_Sv = lavg_Sv.data(); - - const int num_cells = m_ncell_line; - const amrex::Real line_dx = m_dx; - const amrex::Real xlo = m_xlo; - - auto g0 = m_field.repo().mesh().Geom(0); - const amrex::Real denom = - (amrex::Real)m_ncell_line / - ((g0.ProbHi(0) - g0.ProbLo(0)) * (g0.ProbHi(1) - g0.ProbLo(1)) * - (g0.ProbHi(2) - g0.ProbLo(2))); - - const auto& mesh = m_field.repo().mesh(); - const int finestLevel = mesh.finestLevel(); - - for (int lev = 0; lev <= finestLevel; ++lev) { - - const auto& geom = m_field.repo().mesh().Geom(lev); - const amrex::Real dx = geom.CellSize()[m_axis]; - const amrex::Real dy = geom.CellSize()[idxOp.odir1]; - const amrex::Real dz = geom.CellSize()[idxOp.odir2]; - - amrex::iMultiFab level_mask; - if (lev < finestLevel) { - level_mask = makeFineMask( - mesh.boxArray(lev), mesh.DistributionMap(lev), - mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); - } else { - level_mask.define( - mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, - amrex::MFInfo()); - level_mask.setVal(1); - } - - const auto& mfab = m_field(lev); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (amrex::MFIter mfi(mfab, amrex::TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - amrex::Box bx = mfi.tilebox(); - - auto fab_arr = mfab.const_array(mfi); - auto mask_arr = level_mask.const_array(mfi); - - amrex::Box pbx = - perpendicular_box(bx, amrex::IntVect{0, 0, 0}); - - amrex::ParallelFor( - amrex::Gpu::KernelInfo().setReduction(true), pbx, - [=] AMREX_GPU_DEVICE( - int p_i, int p_j, int p_k, - amrex::Gpu::Handler const& handler) { - // Loop over the direction perpendicular to the plane. - // This reduces the atomic pressure on the destination - // arrays. - - amrex::Box lbx = parallel_box( - bx, amrex::IntVect{p_i, p_j, p_k}); - - for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { - for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { - for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); - ++i) { - - // cell coordinates - const amrex::Real cell_xlo = - xlo + (idxOp(i, j, k) * dx); - const amrex::Real cell_xhi = cell_xlo + dx; - - // line indices - const int line_ind_lo = amrex::min( - amrex::max( - static_cast( - (cell_xlo - xlo) / line_dx), - 0), - num_cells - 1); - const int line_ind_hi = amrex::min( - amrex::max( - static_cast( - (cell_xhi - xlo) / line_dx), - 0), - num_cells - 1); - - AMREX_ASSERT(line_ind_lo >= 0); - AMREX_ASSERT(line_ind_hi >= 0); - AMREX_ASSERT(line_ind_lo < num_cells); - AMREX_ASSERT(line_ind_hi < num_cells); - - for (int ind = line_ind_lo; ind <= line_ind_hi; - ++ind) { - - // line coordinates - const amrex::Real line_xlo = - xlo + (ind * line_dx); - const amrex::Real line_xhi = - line_xlo + line_dx; - - amrex::Real deltax; - - if (line_xlo <= cell_xlo) { - deltax = line_xhi - cell_xlo; - } else if (line_xhi >= cell_xhi) { - deltax = cell_xhi - line_xlo; - } else { - deltax = line_dx; - } - - deltax = amrex::min(deltax, dx); - const amrex::Real vol = - mask_arr(i, j, k) * deltax * dy * dz; - - const amrex::Real hvelmag = std::sqrt( - (fab_arr(i, j, k, idxOp.odir1) * - fab_arr(i, j, k, idxOp.odir1)) + - (fab_arr(i, j, k, idxOp.odir2) * - fab_arr(i, j, k, idxOp.odir2))); - const amrex::Real Su = - hvelmag * fab_arr(i, j, k, idxOp.odir1); - const amrex::Real Sv = - hvelmag * fab_arr(i, j, k, idxOp.odir2); - - amrex::Gpu::deviceReduceSum( - &line_avg_vm[ind], - hvelmag * vol * denom, handler); - amrex::Gpu::deviceReduceSum( - &line_avg_Su[ind], Su * vol * denom, - handler); - amrex::Gpu::deviceReduceSum( - &line_avg_Sv[ind], Sv * vol * denom, - handler); - } - } - } - } - }); - } - } - - lavg_vm.copyToHost( - m_line_hvelmag_average.data(), m_line_hvelmag_average.size()); - lavg_Su.copyToHost(m_line_Su_average.data(), m_line_Su_average.size()); - lavg_Sv.copyToHost(m_line_Sv_average.data(), m_line_Sv_average.size()); - amrex::ParallelDescriptor::ReduceRealSum( - m_line_hvelmag_average.data(), - static_cast(m_line_hvelmag_average.size())); - amrex::ParallelDescriptor::ReduceRealSum( - m_line_Su_average.data(), static_cast(m_line_Su_average.size())); - amrex::ParallelDescriptor::ReduceRealSum( - m_line_Sv_average.data(), static_cast(m_line_Sv_average.size())); -} - -amrex::Real -VelPlaneAveragingFine::line_hvelmag_average_interpolated(amrex::Real x) const -{ - int ind; - amrex::Real c; - convert_x_to_ind(x, ind, c); - - return (m_line_hvelmag_average[ind] * (1.0_rt - c)) + - (m_line_hvelmag_average[ind + 1] * c); -} - -amrex::Real -VelPlaneAveragingFine::line_su_average_interpolated(amrex::Real x) const -{ - int ind; - amrex::Real c; - convert_x_to_ind(x, ind, c); - - return (m_line_Su_average[ind] * (1.0_rt - c)) + - (m_line_Su_average[ind + 1] * c); -} - -amrex::Real -VelPlaneAveragingFine::line_sv_average_interpolated(amrex::Real x) const -{ - int ind; - amrex::Real c; - convert_x_to_ind(x, ind, c); - - return (m_line_Sv_average[ind] * (1.0_rt - c)) + - (m_line_Sv_average[ind + 1] * c); -} - -} // namespace amr_wind diff --git a/amr-wind/utilities/SecondMomentAveraging.H b/amr-wind/utilities/SecondMomentAveraging.H index 5d85d1d422..692b7887b5 100644 --- a/amr-wind/utilities/SecondMomentAveraging.H +++ b/amr-wind/utilities/SecondMomentAveraging.H @@ -71,10 +71,7 @@ private: public: // public for GPU /** fill line storage with averages */ template - void compute_average( - const IndexSelector& idxOp, - const amrex::MultiFab& mfab1, - const amrex::MultiFab& mfab2); + void compute_average(const IndexSelector& idxOp); }; } // namespace amr_wind diff --git a/amr-wind/utilities/SecondMomentAveraging.cpp b/amr-wind/utilities/SecondMomentAveraging.cpp index 2bcd73a164..cfb493adc8 100644 --- a/amr-wind/utilities/SecondMomentAveraging.cpp +++ b/amr-wind/utilities/SecondMomentAveraging.cpp @@ -1,6 +1,9 @@ #include #include "SecondMomentAveraging.H" +#include "amr-wind/utilities/constants.H" +#include "AMReX_iMultiFab.H" +#include "AMReX_MultiFabUtil.H" #include "AMReX_REAL.H" using namespace amrex::literals; @@ -101,20 +104,15 @@ void SecondMomentAveraging::operator()() std::ranges::fill(m_second_moments_line, 0.0_rt); - const auto& field1 = m_plane_average1.field(); - const auto& field2 = m_plane_average2.field(); - - const int level = m_plane_average1.level(); - switch (m_plane_average1.axis()) { case 0: - compute_average(XDir(), field1(level), field2(level)); + compute_average(XDir()); break; case 1: - compute_average(YDir(), field1(level), field2(level)); + compute_average(YDir()); break; case 2: - compute_average(ZDir(), field1(level), field2(level)); + compute_average(ZDir()); break; default: amrex::Abort("axis must be equal to 0, 1, or 2"); @@ -123,10 +121,7 @@ void SecondMomentAveraging::operator()() } template -void SecondMomentAveraging::compute_average( - const IndexSelector& idxOp, - const amrex::MultiFab& mfab1, - const amrex::MultiFab& mfab2) +void SecondMomentAveraging::compute_average(const IndexSelector& idxOp) { BL_PROFILE("amr-wind::SecondMomentAveraging::compute_average"); @@ -150,56 +145,203 @@ void SecondMomentAveraging::compute_average( const int ncomp1 = m_plane_average1.ncomp(); const int ncomp2 = m_plane_average2.ncomp(); const int nmoments = m_num_moments; + const int num_cells = m_plane_average1.ncell_line(); + const amrex::Real line_dx = m_plane_average1.dx(); + + const auto max_lev = m_plane_average1.level(); + const auto& mesh = m_plane_average1.field().repo().mesh(); + const int finestLevel = max_lev < 0 ? mesh.finestLevel() : max_lev; + const auto dir = m_plane_average1.axis(); + const bool no_ghost = + (amrex::min( + m_plane_average1.field().num_grow()[dir], + m_plane_average2.field().num_grow()[dir]) == 0); + const bool periodic_dir = mesh.Geom(0).periodicity().isPeriodic(dir); + + for (int lev = 0; lev <= finestLevel; ++lev) { + + const auto& geom = mesh.Geom(lev); + const amrex::Real dx = geom.CellSize()[dir]; + const amrex::Real dy = geom.CellSize()[idxOp.odir1]; + const amrex::Real dz = geom.CellSize()[idxOp.odir2]; + + const amrex::Real problo_x = geom.ProbLo(dir); + const amrex::Real probhi_x = geom.ProbHi(dir); + + amrex::iMultiFab level_mask; + if (lev < finestLevel) { + level_mask = makeFineMask( + mesh.boxArray(lev), mesh.DistributionMap(lev), + mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& mfab1 = m_plane_average1.field()(lev); + const auto& mfab2 = m_plane_average2.field()(lev); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (amrex::MFIter mfi(mfab1, amrex::TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - amrex::Box bx = mfi.tilebox(); - - auto mfab_arr1 = mfab1.const_array(mfi); - auto mfab_arr2 = mfab2.const_array(mfi); - - amrex::Box pbx = - perpendicular_box(bx, amrex::IntVect{0, 0, 0}); - - amrex::ParallelFor( - amrex::Gpu::KernelInfo().setReduction(true), pbx, - [=] AMREX_GPU_DEVICE( - int p_i, int p_j, int p_k, amrex::Gpu::Handler const& handler) { - // Loop over the direction perpendicular to the plane. - // This reduces the atomic pressure on the destination arrays. - - amrex::Box lbx = parallel_box( - bx, amrex::IntVect{p_i, p_j, p_k}); - - for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { - for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { - for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); ++i) { - - const int ind = idxOp(i, j, k); - - int nf = 0; - for (int m = 0; m < ncomp1; ++m) { - const amrex::Real up1 = - mfab_arr1(i, j, k, m) - - line_avg1[(ncomp1 * ind) + m]; - for (int n = 0; n < ncomp2; ++n) { - const amrex::Real up2 = - mfab_arr2(i, j, k, n) - - line_avg2[(ncomp2 * ind) + n]; - - amrex::Gpu::deviceReduceSum( - &line_fluc[(nmoments * ind) + nf], - up1 * up2 * denom, handler); - ++nf; + for (amrex::MFIter mfi(mfab1, amrex::TilingIfNotGPU()); mfi.isValid(); + ++mfi) { + amrex::Box bx = mfi.tilebox(); + + auto mfab_arr1 = mfab1.const_array(mfi); + auto mfab_arr2 = mfab2.const_array(mfi); + auto mask_arr = level_mask.const_array(mfi); + + amrex::Box pbx = + perpendicular_box(bx, amrex::IntVect{0, 0, 0}); + + amrex::ParallelFor( + amrex::Gpu::KernelInfo().setReduction(true), pbx, + [=] AMREX_GPU_DEVICE( + int p_i, int p_j, int p_k, + amrex::Gpu::Handler const& handler) { + // Loop over the direction perpendicular to the plane. + // This reduces the atomic pressure on the destination + // arrays. + + amrex::Box lbx = parallel_box( + bx, amrex::IntVect{p_i, p_j, p_k}); + + for (int k = lbx.smallEnd(2); k <= lbx.bigEnd(2); ++k) { + for (int j = lbx.smallEnd(1); j <= lbx.bigEnd(1); ++j) { + for (int i = lbx.smallEnd(0); i <= lbx.bigEnd(0); + ++i) { + + // cell coordinates + const amrex::Real cell_xlo = + problo_x + idxOp(i, j, k) * dx; + const amrex::Real cell_xhi = cell_xlo + dx; + + // line indices + const int line_ind_lo = amrex::min( + amrex::max( + static_cast( + (cell_xlo - problo_x) / line_dx), + 0), + num_cells - 1); + const int line_ind_hi = amrex::min( + amrex::max( + static_cast( + (cell_xhi - + amr_wind::constants::TIGHT_TOL - + problo_x) / + line_dx), + 0), + num_cells - 1); + + AMREX_ASSERT(line_ind_lo >= 0); + AMREX_ASSERT(line_ind_hi >= 0); + AMREX_ASSERT(line_ind_lo < num_cells); + AMREX_ASSERT(line_ind_hi < num_cells); + + for (int ind = line_ind_lo; ind <= line_ind_hi; + ++ind) { + + // line coordinates + const amrex::Real line_xlo = + problo_x + ind * line_dx; + const amrex::Real line_xhi = + line_xlo + line_dx; + + amrex::Real deltax; + + if (line_xlo <= cell_xlo) { + deltax = line_xhi - cell_xlo; + } else if (line_xhi >= cell_xhi) { + deltax = cell_xhi - line_xlo; + } else { + deltax = line_dx; + } + deltax = amrex::min(deltax, dx); + const amrex::Real vol = deltax * dy * dz; + + // Calculate location of target + const auto x_targ = + 0.5_rt * (line_xlo + line_xhi); + // Calculate location of cell center + const amrex::IntVect iv{i, j, k}; + const auto idx = iv[dir]; + const auto x_cell = + problo_x + + (static_cast(idx) + + 0.5_rt) * + dx; + // Get location of neighboring cell centers + auto x_up = x_cell + dx; + auto x_down = x_cell - dx; + // Bound locations by domain limits + if (!periodic_dir) { + x_up = amrex::min(probhi_x, x_up); + x_down = amrex::max(problo_x, x_down); + } + // Pick indices of closest neighbor + // Bound indices in case of no ghost cells + auto iv_nb = iv; + auto x_nb = x_cell; + if (std::abs(x_up - x_targ) < + std::abs(x_down - x_targ)) { + x_nb = x_up; + iv_nb[dir] += 1; + if (no_ghost) { + iv_nb[dir] = amrex::min( + iv_nb[dir], bx.bigEnd(dir)); + } + } else { + x_nb = x_down; + iv_nb[dir] -= 1; + if (no_ghost) { + iv_nb[dir] = amrex::max( + iv_nb[dir], bx.smallEnd(dir)); + } + } + // Interpolate to target location using + // closest neighbor + // (will do nothing if already at cell + // center) + int nf = 0; + for (int m = 0; m < ncomp1; ++m) { + const auto arr1_interp = + mfab_arr1(iv, m) + + (mfab_arr1(iv_nb, m) - + mfab_arr1(iv, m)) * + ((x_targ - x_cell) / + (x_nb - x_cell)); + const amrex::Real up1 = + arr1_interp - + line_avg1[(ncomp1 * ind) + m]; + for (int n = 0; n < ncomp2; ++n) { + const auto arr2_interp = + mfab_arr2(iv, n) + + (mfab_arr2(iv_nb, n) - + mfab_arr2(iv, n)) * + ((x_targ - x_cell) / + (x_nb - x_cell)); + const amrex::Real up2 = + arr2_interp - + line_avg2[(ncomp2 * ind) + n]; + + amrex::Gpu::deviceReduceSum( + &line_fluc + [(nmoments * ind) + nf], + mask_arr(i, j, k) * up1 * up2 * + vol * denom, + handler); + ++nf; + } + } } } } } - } - }); + }); + } } lfluc.copyToHost( diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index f8f0cb88ca..6b61327ad3 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -37,13 +37,16 @@ ABL::ABL(CFDSim& sim) { std::string statistics_mode = "precursor"; int dir = 2; + // Use current behavior for regression (!!CHANGE for final PR!!) + int stats_max_level = 0; amrex::ParmParse pp("ABL"); pp.query("enable_hybrid_rl_mode", m_hybrid_rl); pp.query("initial_sdr_value", m_init_sdr); pp.query("normal_direction", dir); pp.query("statistics_mode", statistics_mode); - m_stats = - ABLStatsBase::create(statistics_mode, sim, m_abl_wall_func, dir); + pp.query("stats_max_level", stats_max_level); + m_stats = ABLStatsBase::create( + statistics_mode, sim, m_abl_wall_func, dir, stats_max_level); // Check for file input m_file_input = pp.contains("initial_condition_input_file"); } diff --git a/amr-wind/wind_energy/ABLStats.H b/amr-wind/wind_energy/ABLStats.H index 920837de7f..9081fe13aa 100644 --- a/amr-wind/wind_energy/ABLStats.H +++ b/amr-wind/wind_energy/ABLStats.H @@ -5,7 +5,6 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/transport_models/TransportModel.H" #include "amr-wind/utilities/FieldPlaneAveraging.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" #include "amr-wind/utilities/SecondMomentAveraging.H" #include "amr-wind/utilities/ThirdMomentAveraging.H" #include "amr-wind/utilities/PostProcessing.H" @@ -41,7 +40,10 @@ public: static std::string identifier() { return "precursor"; } ABLStats( - CFDSim& /*sim*/, const ABLWallFunction& /*abl_wall_func*/, int dir); + CFDSim& /*sim*/, + const ABLWallFunction& /*abl_wall_func*/, + int dir, + int max_level); ~ABLStats() override; @@ -60,7 +62,7 @@ public: void compute_zi(); //! Return vel plane averaging instance - const VelPlaneAveragingFine& vel_profile() const override + const VelPlaneAveraging& vel_profile() const override { return m_pa_vel_fine; }; @@ -78,7 +80,7 @@ public: } //! Return instance that handles temperature statistics - const FieldPlaneAveragingFine& theta_profile_fine() const override + const FieldPlaneAveraging& theta_profile_fine() const override { return m_pa_temp_fine; } @@ -145,8 +147,8 @@ private: VelPlaneAveraging m_pa_vel; FieldPlaneAveraging m_pa_temp; - VelPlaneAveragingFine m_pa_vel_fine; - FieldPlaneAveragingFine m_pa_temp_fine; + VelPlaneAveraging m_pa_vel_fine; + FieldPlaneAveraging m_pa_temp_fine; FieldPlaneAveraging m_pa_mueff; SecondMomentAveraging m_pa_tt; SecondMomentAveraging m_pa_tu; @@ -169,6 +171,13 @@ private: #endif std::string m_ascii_file_name; +#ifdef AMR_WIND_USE_NETCDF + //! Maximum level to consider in planar averages: -1 means no limit + int m_max_level{-1}; + // - when netcdf is not used, max level is still used but it is only needed + // - at initialization, does not need to be stored +#endif + //! Frequency of data sampling and output int m_out_freq{100}; diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index 8b58974d0d..eda826daf8 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -22,20 +22,26 @@ using namespace amrex::literals; namespace amr_wind { ABLStats::ABLStats( - CFDSim& sim, const ABLWallFunction& abl_wall_func, const int dir) + CFDSim& sim, + const ABLWallFunction& abl_wall_func, + const int dir, + const int max_level) : m_sim(sim) , m_abl_wall_func(abl_wall_func) , m_temperature(sim.repo().get_field("temperature")) , m_mueff(sim.pde_manager().icns().fields().mueff) - , m_pa_vel(sim, dir) - , m_pa_temp(m_temperature, sim.time(), dir) - , m_pa_vel_fine(sim, dir) - , m_pa_temp_fine(m_temperature, sim.time(), dir) - , m_pa_mueff(m_mueff, sim.time(), dir) + , m_pa_vel(sim, dir, max_level) + , m_pa_temp(m_temperature, sim.time(), dir, max_level) + , m_pa_vel_fine(sim, dir, -1) + , m_pa_temp_fine(m_temperature, sim.time(), dir, -1) + , m_pa_mueff(m_mueff, sim.time(), dir, max_level) , m_pa_tt(m_pa_temp, m_pa_temp) , m_pa_tu(m_pa_vel, m_pa_temp) , m_pa_uu(m_pa_vel, m_pa_vel) , m_pa_uuu(m_pa_vel, m_pa_vel, m_pa_vel) +#ifdef AMR_WIND_USE_NETCDF + , m_max_level(max_level) +#endif {} ABLStats::~ABLStats() = default; @@ -588,15 +594,17 @@ void ABLStats::write_netcdf() auto sfs_stress = m_sim.repo().create_scratch_field("sfs_stress", 3); auto t_sfs_stress = m_sim.repo().create_scratch_field("tsfs_stress", 3); calc_sfs_stress_avgs(*sfs_stress, *t_sfs_stress); - ScratchFieldPlaneAveraging pa_sfs(*sfs_stress, m_sim.time(), m_normal_dir); + ScratchFieldPlaneAveraging pa_sfs( + *sfs_stress, m_sim.time(), m_normal_dir, m_max_level); pa_sfs(); ScratchFieldPlaneAveraging pa_tsfs( - *t_sfs_stress, m_sim.time(), m_normal_dir); + *t_sfs_stress, m_sim.time(), m_normal_dir, m_max_level); pa_tsfs(); if (m_sim.repo().field_exists("tke")) { auto& m_ksgs = m_sim.repo().get_field("tke"); - FieldPlaneAveraging pa_ksgs(m_ksgs, m_sim.time(), m_normal_dir); + FieldPlaneAveraging pa_ksgs( + m_ksgs, m_sim.time(), m_normal_dir, m_max_level); pa_ksgs(); } @@ -757,28 +765,28 @@ void ABLStats::write_netcdf() m_sim.time().delta_t()); { FieldPlaneAveraging pa_tke_buoy_prod( - tke_dissip, m_sim.time(), m_normal_dir); + tke_dissip, m_sim.time(), m_normal_dir, m_max_level); pa_tke_buoy_prod(); auto var = grp.var("tke_buoy"); var.put(pa_tke_buoy_prod.line_average().data(), start, count); } { FieldPlaneAveraging pa_tke_shear_prod( - tke_shear_prod, m_sim.time(), m_normal_dir); + tke_shear_prod, m_sim.time(), m_normal_dir, m_max_level); pa_tke_shear_prod(); auto var = grp.var("tke_shear"); var.put(pa_tke_shear_prod.line_average().data(), start, count); } { FieldPlaneAveraging pa_tke_dissip( - tke_dissip, m_sim.time(), m_normal_dir); + tke_dissip, m_sim.time(), m_normal_dir, m_max_level); pa_tke_dissip(); auto var = grp.var("tke_dissip"); var.put(pa_tke_dissip.line_average().data(), start, count); } { ScratchFieldPlaneAveraging pa_tke_diff( - *tke_diffusion, m_sim.time(), m_normal_dir); + *tke_diffusion, m_sim.time(), m_normal_dir, m_max_level); pa_tke_diff(); auto var = grp.var("tke_diff"); var.put(pa_tke_diff.line_average().data(), start, count); diff --git a/amr-wind/wind_energy/ABLStatsBase.H b/amr-wind/wind_energy/ABLStatsBase.H index e9d0b6048b..f3df42390c 100644 --- a/amr-wind/wind_energy/ABLStatsBase.H +++ b/amr-wind/wind_energy/ABLStatsBase.H @@ -21,7 +21,7 @@ enum class ABLStatsMode : std::uint8_t { }; class ABLStatsBase - : public Factory + : public Factory { public: static std::string base_identifier() { return "ABLStatsBase"; } @@ -36,11 +36,11 @@ public: //! Interpolating object for vertical velocity profile [[nodiscard]] virtual const VelPlaneAveraging& vel_profile_coarse() const = 0; - [[nodiscard]] virtual const VelPlaneAveragingFine& vel_profile() const = 0; + [[nodiscard]] virtual const VelPlaneAveraging& vel_profile() const = 0; //! Interpolating object for vertical temperature profile [[nodiscard]] virtual const FieldPlaneAveraging& theta_profile() const = 0; - [[nodiscard]] virtual const FieldPlaneAveragingFine& + [[nodiscard]] virtual const FieldPlaneAveraging& theta_profile_fine() const = 0; //! Perform initialization actions after the mesh has been created diff --git a/amr-wind/wind_energy/ABLWallFunction.H b/amr-wind/wind_energy/ABLWallFunction.H index bfa18e1fc2..d49a8fefc4 100644 --- a/amr-wind/wind_energy/ABLWallFunction.H +++ b/amr-wind/wind_energy/ABLWallFunction.H @@ -3,7 +3,6 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/utilities/FieldPlaneAveraging.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" #include "amr-wind/core/FieldBCOps.H" #include "amr-wind/wind_energy/MOData.H" #include "amr-wind/utilities/constants.H" @@ -38,8 +37,8 @@ public: void init_log_law_height(); //! Update the mean velocity at a given timestep - void update_umean( - const VelPlaneAveragingFine& vpa, const FieldPlaneAveragingFine& tpa); + void + update_umean(const VelPlaneAveraging& vpa, const FieldPlaneAveraging& tpa); void update_tflux(amrex::Real tflux); diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 86a8caf06b..ae76f803ad 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -150,7 +150,7 @@ void ABLWallFunction::init_log_law_height() } void ABLWallFunction::update_umean( - const VelPlaneAveragingFine& vpa, const FieldPlaneAveragingFine& tpa) + const VelPlaneAveraging& vpa, const FieldPlaneAveraging& tpa) { const auto& time = m_sim.time(); diff --git a/unit_tests/utilities/test_field_plane_averaging.cpp b/unit_tests/utilities/test_field_plane_averaging.cpp index 7b41b6dc9b..788860e5a5 100644 --- a/unit_tests/utilities/test_field_plane_averaging.cpp +++ b/unit_tests/utilities/test_field_plane_averaging.cpp @@ -45,7 +45,7 @@ TEST_F(FieldPlaneAveragingTest, test_constant) // test the average of a constant is the same constant for (int dir = 0; dir < 3; ++dir) { - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); amrex::Real z = 0.5_rt * (problo[dir] + probhi[dir]); @@ -126,7 +126,7 @@ TEST_F(FieldPlaneAveragingTest, test_linear) add_linear(dir, u0, mesh().Geom(0), bx, vel); }); - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); constexpr int n = 20; @@ -241,7 +241,7 @@ void FieldPlaneAveragingTest::test_dir(int dir) add_periodic(dir, a, mesh().Geom(lev), bx, vel); }); - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); amrex::Real x = 0.5_rt * (problo[dir] + probhi[dir]); diff --git a/unit_tests/utilities/test_field_plane_averaging_fine.cpp b/unit_tests/utilities/test_field_plane_averaging_fine.cpp index df237a6881..460518b898 100644 --- a/unit_tests/utilities/test_field_plane_averaging_fine.cpp +++ b/unit_tests/utilities/test_field_plane_averaging_fine.cpp @@ -10,7 +10,7 @@ #include "AMReX_Vector.H" #include "AMReX_REAL.H" -#include "amr-wind/utilities/FieldPlaneAveragingFine.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" using namespace amrex::literals; @@ -127,8 +127,8 @@ TEST_F(FieldPlaneAveragingFineTest, test_linear_fine_only) constexpr int dir = 2; init_field_linear(velocityf, u0, dir); - amr_wind::FieldPlaneAveragingFine pa_fine( - velocityf, sim().time(), dir, true); + amr_wind::FieldPlaneAveraging pa_fine( + velocityf, sim().time(), dir, -1, true); pa_fine(); constexpr int n = 20; @@ -180,8 +180,8 @@ TEST_F(FieldPlaneAveragingFineTest, test_linear) constexpr int dir = 2; init_field_linear(velocityf, u0, dir); - amr_wind::FieldPlaneAveragingFine pa_fine( - velocityf, sim().time(), dir, true); + amr_wind::FieldPlaneAveraging pa_fine( + velocityf, sim().time(), dir, -1, true); pa_fine(); constexpr int n = 20; diff --git a/unit_tests/utilities/test_plane_averaging.cpp b/unit_tests/utilities/test_plane_averaging.cpp index 9fead03aa0..13cdd84b4e 100644 --- a/unit_tests/utilities/test_plane_averaging.cpp +++ b/unit_tests/utilities/test_plane_averaging.cpp @@ -52,7 +52,7 @@ TEST_F(PlaneAveragingTest, test_constant) // test the average of a constant is the same constant for (int dir = 0; dir < 3; ++dir) { - amr_wind::VelPlaneAveraging pa(sim(), dir); + amr_wind::VelPlaneAveraging pa(sim(), dir, 0); pa(); amrex::Real z = 0.5_rt * (problo[dir] + probhi[dir]); @@ -133,7 +133,7 @@ TEST_F(PlaneAveragingTest, test_linear) add_linear(dir, u0, mesh().Geom(0), bx, vel); }); - amr_wind::VelPlaneAveraging pa(sim(), dir); + amr_wind::VelPlaneAveraging pa(sim(), dir, 0); pa(); constexpr int n = 20; @@ -251,7 +251,7 @@ void PlaneAveragingTest::test_dir(int dir) add_periodic(dir, a, mesh().Geom(lev), bx, vel); }); - amr_wind::VelPlaneAveraging pa(sim(), dir); + amr_wind::VelPlaneAveraging pa(sim(), dir, 0); pa(); amrex::Real x = 0.5_rt * (problo[dir] + probhi[dir]); diff --git a/unit_tests/utilities/test_second_moment.cpp b/unit_tests/utilities/test_second_moment.cpp index 379c647e61..49cc620dc9 100644 --- a/unit_tests/utilities/test_second_moment.cpp +++ b/unit_tests/utilities/test_second_moment.cpp @@ -48,7 +48,7 @@ TEST_F(SecondMomentAveragingTest, test_constant) // test the average of a constant is the same constant for (int dir = 0; dir < 3; ++dir) { - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); amr_wind::SecondMomentAveraging uu(pa, pa); @@ -120,7 +120,7 @@ TEST_F(SecondMomentAveragingTest, test_linear) add_linear(dir, u0, mesh().Geom(0), bx, vel); }); - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); amr_wind::SecondMomentAveraging uu(pa, pa); uu(); @@ -232,7 +232,7 @@ void SecondMomentAveragingTest::test_dir(int dir) add_periodic(a, mesh().Geom(lev), bx, vel); }); - amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir); + amr_wind::FieldPlaneAveraging pa(velocityf, sim().time(), dir, 0); pa(); amr_wind::SecondMomentAveraging uu(pa, pa); diff --git a/unit_tests/wind_energy/test_abl_stats.cpp b/unit_tests/wind_energy/test_abl_stats.cpp index 05b5966a4f..fc1d515dae 100644 --- a/unit_tests/wind_energy/test_abl_stats.cpp +++ b/unit_tests/wind_energy/test_abl_stats.cpp @@ -142,7 +142,7 @@ TEST_F(ABLMeshTest, stats_tke_diffusion) // Initialize ABL Stats amr_wind::ABLWallFunction wall_func(sim()); - amr_wind::ABLStats stats(sim(), wall_func, 2); + amr_wind::ABLStats stats(sim(), wall_func, 2, 0); // Calculate diffusion term stats.calc_tke_diffusion(*diff, buoy, shear, dissip, dt); @@ -211,7 +211,7 @@ TEST_F(ABLMeshTest, stats_energy_budget) // Initialize ABL Stats amr_wind::ABLWallFunction wall_func(sim()); - amr_wind::ABLStats stats(sim(), wall_func, 2); + amr_wind::ABLStats stats(sim(), wall_func, 2, 0); // Set initial tke value and advance states init_field1(tke);