diff --git a/amr-wind/utilities/index_operations.H b/amr-wind/utilities/index_operations.H index 49182fa36f..0b357fe664 100644 --- a/amr-wind/utilities/index_operations.H +++ b/amr-wind/utilities/index_operations.H @@ -4,6 +4,7 @@ #include #include #include +#include #include namespace amr_wind::utils { @@ -36,6 +37,21 @@ AMREX_FORCE_INLINE T perpendicular_idx(const int normal) return T{-1, -1}; } +//! Check if a point is inside a box +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool contains( + const amrex::Box& box, + const amrex::RealVect& pos, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& dxinv) +{ + const amrex::IntVect iv(AMREX_D_DECL( + static_cast(amrex::Math::floor((pos[0] - prob_lo[0]) * dxinv[0])), + static_cast(amrex::Math::floor((pos[1] - prob_lo[1]) * dxinv[1])), + static_cast( + amrex::Math::floor((pos[2] - prob_lo[2]) * dxinv[2])))); + return box.contains(iv); +} + /** Convert a bounding box into amrex::Box index space at a given level * * \param rbx Bounding box as defined in global domain coordinates diff --git a/amr-wind/utilities/sampling/DTUSpinnerSampler.H b/amr-wind/utilities/sampling/DTUSpinnerSampler.H index 3f03c82c9c..cbe2e60538 100644 --- a/amr-wind/utilities/sampling/DTUSpinnerSampler.H +++ b/amr-wind/utilities/sampling/DTUSpinnerSampler.H @@ -41,7 +41,13 @@ public: * * */ - void sampling_locations(SampleLocType& /*locs*/) const override; + void sampling_locations(SampleLocType& /*sample_locs**/) const override; + + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; static vs::Vector generate_lidar_pattern( PrismParameters InnerPrism, PrismParameters OuterPrism, double time); diff --git a/amr-wind/utilities/sampling/DTUSpinnerSampler.cpp b/amr-wind/utilities/sampling/DTUSpinnerSampler.cpp index cf87ab5679..471fad0272 100644 --- a/amr-wind/utilities/sampling/DTUSpinnerSampler.cpp +++ b/amr-wind/utilities/sampling/DTUSpinnerSampler.cpp @@ -2,6 +2,7 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/utilities/tensor_ops.H" #include "amr-wind/utilities/linear_interpolation.H" +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" #ifdef AMR_WIND_USE_OPENFAST @@ -132,22 +133,27 @@ vs::Vector DTUSpinnerSampler::generate_lidar_pattern( reflection_2, sampling_utils::reflect(reflection_1, axis)); } -// - -void DTUSpinnerSampler::sampling_locations(SampleLocType& locs) const +void DTUSpinnerSampler::sampling_locations(SampleLocType& sample_locs) const { + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); - // The total number of points at this time step - long n_samples = m_beam_points * m_ntotal; + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); - // Resize to number of points in line times number of sampling times - if (locs.size() < n_samples) { - locs.resize(n_samples); - } + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} +void DTUSpinnerSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); const amrex::Real ndiv = amrex::max(m_beam_points - 1, 1); amrex::Array dx; - // Loop per subsampling for (int k = 0; k < m_ntotal; ++k) { @@ -159,9 +165,15 @@ void DTUSpinnerSampler::sampling_locations(SampleLocType& locs) const } for (int i = 0; i < m_beam_points; ++i) { - for (int d = 0; d < AMREX_SPACEDIM; ++d) { - locs[i + k * m_beam_points][d] = - m_start[d + offset] + i * dx[d]; + const amrex::RealVect loc = {AMREX_D_DECL( + m_start[0 + offset] + i * dx[0], + m_start[1 + offset] + i * dx[1], + m_start[2 + offset] + i * dx[2])}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, i + k * m_beam_points); + break; + } } } } @@ -512,8 +524,8 @@ void DTUSpinnerSampler::output_netcdf_data( std::vector count{1, 0, AMREX_SPACEDIM}; std::vector starti{nt, 0}; std::vector counti{1, 0}; - SamplerBase::SampleLocType locs; - sampling_locations(locs); + SampleLocType sample_locs; + sampling_locations(sample_locs); auto xyz = grp.var("points"); auto xp = grp.var("points_x"); @@ -522,7 +534,8 @@ void DTUSpinnerSampler::output_netcdf_data( count[1] = num_points(); counti[1] = num_points(); - xyz.put(locs[0].data(), start, count); + const auto& locs = sample_locs.locations(); + xyz.put(locs[0].begin(), start, count); auto n_samples = m_beam_points * m_ntotal; diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.H b/amr-wind/utilities/sampling/FreeSurfaceSampler.H index 83d5056f42..6a0e72080f 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.H +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.H @@ -23,10 +23,17 @@ public: void check_bounds() override; //! Populate and return a vector of probe locations to be sampled - void sampling_locations(SampleLocType& /*locs*/) const override; - void output_locations(SampleLocType& locs) const override + void sampling_locations(SampleLocType& /*sample_locs**/) const override; + + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + + void output_locations(SampleLocType& sample_locs) const override { - return sampling_locations(locs); + return sampling_locations(sample_locs); } //! Find heights associated with 2D sample locations diff --git a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp index a560eec23f..8da73bd617 100644 --- a/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp +++ b/amr-wind/utilities/sampling/FreeSurfaceSampler.cpp @@ -3,6 +3,7 @@ #include #include #include "amr-wind/equation_systems/vof/volume_fractions.H" +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" @@ -365,22 +366,39 @@ void FreeSurfaceSampler::check_bounds() } } -void FreeSurfaceSampler::sampling_locations(SampleLocType& locs) const +void FreeSurfaceSampler::sampling_locations(SampleLocType& sample_locs) const { - locs.resize(num_output_points()); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void FreeSurfaceSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); int idx = 0; + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); for (int j = 0; j < m_npts_dir[1]; ++j) { for (int i = 0; i < m_npts_dir[0]; ++i) { - // Initialize output values to 0.0 for (int ni = 0; ni < m_ninst; ++ni) { - // Grid direction 1 - locs[idx * m_ninst + ni][m_gc0] = m_grid_locs[idx][0]; - // Grid direction 2 - locs[idx * m_ninst + ni][m_gc1] = m_grid_locs[idx][1]; - // Output direction - locs[idx * m_ninst + ni][m_coorddir] = - m_out[idx * m_ninst + ni]; + amrex::RealVect loc; + loc[m_gc0] = m_grid_locs[idx][0]; + loc[m_gc1] = m_grid_locs[idx][1]; + loc[m_coorddir] = m_out[idx * m_ninst + ni]; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, idx * m_ninst + ni); + break; + } + } } ++idx; } @@ -727,11 +745,12 @@ void FreeSurfaceSampler::output_netcdf_data( // Write the coordinates every time std::vector start{nt, 0, 0}; std::vector count{1, 0, AMREX_SPACEDIM}; - SamplerBase::SampleLocType locs; - sampling_locations(locs); + SampleLocType sample_locs; + sampling_locations(sample_locs); auto xyz = grp.var("points"); count[1] = num_output_points(); - xyz.put(locs[0].data(), start, count); + const auto& locs = sample_locs.locations(); + xyz.put(locs[0].begin(), start, count); } #else void FreeSurfaceSampler::define_netcdf_metadata( diff --git a/amr-wind/utilities/sampling/LidarSampler.cpp b/amr-wind/utilities/sampling/LidarSampler.cpp index 83dbd6a665..81970406d9 100644 --- a/amr-wind/utilities/sampling/LidarSampler.cpp +++ b/amr-wind/utilities/sampling/LidarSampler.cpp @@ -109,11 +109,12 @@ void LidarSampler::output_netcdf_data( // Write the coordinates every time std::vector start{nt, 0, 0}; std::vector count{1, 0, AMREX_SPACEDIM}; - SamplerBase::SampleLocType locs; - sampling_locations(locs); + SampleLocType sample_locs; + sampling_locations(sample_locs); auto xyz = grp.var("points"); count[1] = num_points(); - xyz.put(locs[0].data(), start, count); + const auto& locs = sample_locs.locations(); + xyz.put(locs[0].begin(), start, count); } #else void LidarSampler::define_netcdf_metadata( diff --git a/amr-wind/utilities/sampling/LineSampler.H b/amr-wind/utilities/sampling/LineSampler.H index bc0029f3b2..eb651903c1 100644 --- a/amr-wind/utilities/sampling/LineSampler.H +++ b/amr-wind/utilities/sampling/LineSampler.H @@ -31,11 +31,17 @@ public: void check_bounds() override; //! Populate and return a vector of probe locations to be sampled - void sampling_locations(SampleLocType& /*locs*/) const override; + void sampling_locations(SampleLocType& /*sample_locs**/) const override; - void output_locations(SampleLocType& locs) const override + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + + void output_locations(SampleLocType& sample_locs) const override { - return sampling_locations(locs); + return sampling_locations(sample_locs); } void diff --git a/amr-wind/utilities/sampling/LineSampler.cpp b/amr-wind/utilities/sampling/LineSampler.cpp index 15cfbcfe3f..9b39086e4e 100644 --- a/amr-wind/utilities/sampling/LineSampler.cpp +++ b/amr-wind/utilities/sampling/LineSampler.cpp @@ -1,6 +1,7 @@ #include "amr-wind/utilities/sampling/LineSampler.H" #include "amr-wind/CFDSim.H" #include "amr-wind/utilities/tensor_ops.H" +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" @@ -53,10 +54,25 @@ void LineSampler::check_bounds() } } -void LineSampler::sampling_locations(SampleLocType& locs) const +void LineSampler::sampling_locations(SampleLocType& sample_locs) const { - locs.resize(m_npts); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void LineSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); const amrex::Real ndiv = amrex::max(m_npts - 1, 1); amrex::Array dx; @@ -65,8 +81,14 @@ void LineSampler::sampling_locations(SampleLocType& locs) const } for (int i = 0; i < m_npts; ++i) { - for (int d = 0; d < AMREX_SPACEDIM; ++d) { - locs[i][d] = m_start[d] + i * dx[d]; + const amrex::RealVect loc = {AMREX_D_DECL( + m_start[0] + i * dx[0], m_start[1] + i * dx[1], + m_start[2] + i * dx[2])}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, i); + break; + } } } } diff --git a/amr-wind/utilities/sampling/PlaneSampler.H b/amr-wind/utilities/sampling/PlaneSampler.H index 3c2fa6db72..d61f5c1b53 100644 --- a/amr-wind/utilities/sampling/PlaneSampler.H +++ b/amr-wind/utilities/sampling/PlaneSampler.H @@ -42,11 +42,17 @@ public: void check_bounds() override; //! Populate and return a vector of probe locations to be sampled - void sampling_locations(SampleLocType& /*locs*/) const override; + void sampling_locations(SampleLocType& /*sample_locs**/) const override; - void output_locations(SampleLocType& locs) const override + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + + void output_locations(SampleLocType& sample_locs) const override { - return sampling_locations(locs); + return sampling_locations(sample_locs); } void diff --git a/amr-wind/utilities/sampling/PlaneSampler.cpp b/amr-wind/utilities/sampling/PlaneSampler.cpp index 7c16e53ace..5098cca096 100644 --- a/amr-wind/utilities/sampling/PlaneSampler.cpp +++ b/amr-wind/utilities/sampling/PlaneSampler.cpp @@ -2,7 +2,7 @@ #include "amr-wind/utilities/sampling/PlaneSampler.H" #include "amr-wind/CFDSim.H" - +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" namespace amr_wind::sampling { @@ -75,9 +75,21 @@ void PlaneSampler::check_bounds() } } -void PlaneSampler::sampling_locations(SampleLocType& locs) const +void PlaneSampler::sampling_locations(SampleLocType& sample_locs) const { - locs.resize(m_npts); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void PlaneSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); amrex::Array dx; amrex::Array dy; @@ -89,12 +101,22 @@ void PlaneSampler::sampling_locations(SampleLocType& locs) const int idx = 0; const int nplanes = static_cast(m_poffsets.size()); + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); for (int k = 0; k < nplanes; ++k) { for (int j = 0; j < m_npts_dir[1]; ++j) { for (int i = 0; i < m_npts_dir[0]; ++i) { + amrex::RealVect loc; for (int d = 0; d < AMREX_SPACEDIM; ++d) { - locs[idx][d] = m_origin[d] + dx[d] * i + dy[d] * j + - m_poffsets[k] * m_offset_vector[d]; + loc[d] = m_origin[d] + dx[d] * i + dy[d] * j + + m_poffsets[k] * m_offset_vector[d]; + } + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, idx); + break; + } } ++idx; } diff --git a/amr-wind/utilities/sampling/ProbeSampler.H b/amr-wind/utilities/sampling/ProbeSampler.H index de59322241..eedd12a3fc 100644 --- a/amr-wind/utilities/sampling/ProbeSampler.H +++ b/amr-wind/utilities/sampling/ProbeSampler.H @@ -28,11 +28,18 @@ public: //! Check and fix the bounds of the sampler so the probes are in the domain void check_bounds() override; - void sampling_locations(SampleLocType& /*locs*/) const override; + //! Populate the vector with coordinates of the sampling locations + void sampling_locations(SampleLocType& /*sample_locs**/) const override; - void output_locations(SampleLocType& locs) const override + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + + void output_locations(SampleLocType& sample_locs) const override { - return sampling_locations(locs); + return sampling_locations(sample_locs); } void diff --git a/amr-wind/utilities/sampling/ProbeSampler.cpp b/amr-wind/utilities/sampling/ProbeSampler.cpp index dfe7aed42e..dceb277c6e 100644 --- a/amr-wind/utilities/sampling/ProbeSampler.cpp +++ b/amr-wind/utilities/sampling/ProbeSampler.cpp @@ -1,6 +1,6 @@ #include "amr-wind/utilities/sampling/ProbeSampler.H" #include "amr-wind/CFDSim.H" - +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" namespace amr_wind::sampling { @@ -22,10 +22,11 @@ void ProbeSampler::initialize(const std::string& key) ifh >> m_npts; ifh.ignore(std::numeric_limits::max(), '\n'); - m_probes.resize(m_npts); for (int i = 0; i < m_npts; ++i) { - ifh >> m_probes[i][0] >> m_probes[i][1] >> m_probes[i][2]; + amrex::RealVect loc; + ifh >> loc[0] >> loc[1] >> loc[2]; ifh.ignore(std::numeric_limits::max(), '\n'); + m_probes.push_back(loc, i); } check_bounds(); } @@ -35,16 +36,17 @@ void ProbeSampler::check_bounds() const int lev = 0; const auto* prob_lo = m_sim.mesh().Geom(lev).ProbLo(); const auto* prob_hi = m_sim.mesh().Geom(lev).ProbHi(); + auto& probe_locs = m_probes.locations(); bool all_ok = true; for (int i = 0; i < m_npts; ++i) { for (int d = 0; d < AMREX_SPACEDIM; ++d) { - if (m_probes[i][d] <= prob_lo[d]) { + if (probe_locs[i][d] <= prob_lo[d]) { all_ok = false; - m_probes[i][d] = prob_lo[d] + bounds_tol; + probe_locs[i][d] = prob_lo[d] + bounds_tol; } - if (m_probes[i][d] >= prob_hi[d]) { + if (probe_locs[i][d] >= prob_hi[d]) { all_ok = false; - m_probes[i][d] = prob_hi[d] - bounds_tol; + probe_locs[i][d] = prob_hi[d] - bounds_tol; } } } @@ -55,12 +57,36 @@ void ProbeSampler::check_bounds() } } -void ProbeSampler::sampling_locations(SampleLocType& locs) const +void ProbeSampler::sampling_locations(SampleLocType& sample_locs) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void ProbeSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const { - locs.resize(m_npts); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const auto& probe_locs = m_probes.locations(); + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); for (int i = 0; i < m_npts; ++i) { for (int d = 0; d < AMREX_SPACEDIM; ++d) { - locs[i][d] = m_probes[i][d]; + const amrex::RealVect loc = {AMREX_D_DECL( + probe_locs[i][0], probe_locs[i][1], probe_locs[i][2])}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, i); + break; + } + } } } } diff --git a/amr-wind/utilities/sampling/RadarSampler.H b/amr-wind/utilities/sampling/RadarSampler.H index 2bf8ace73a..5383dc2217 100644 --- a/amr-wind/utilities/sampling/RadarSampler.H +++ b/amr-wind/utilities/sampling/RadarSampler.H @@ -47,12 +47,19 @@ public: double determine_current_sweep_angle() const; //! Populate and return a vector of probe locations to be sampled - void sampling_locations(SampleLocType& /*locs*/) const override; + void sampling_locations(SampleLocType& /*sample_locs**/) const override; + + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + bool update_sampling_locations() override; - void cone_axis_locations(SampleLocType& /*locs*/) const; - void output_locations(SampleLocType& locs) const override + void cone_axis_locations(SampleLocType& /*sample_locs**/) const; + void output_locations(SampleLocType& sample_locs) const override { - return cone_axis_locations(locs); + return cone_axis_locations(sample_locs); } void post_sample_actions() override; @@ -140,10 +147,10 @@ protected: amrex::Vector m_start; amrex::Vector m_end; - SamplerBase::SampleLocType m_initial_cone; - SamplerBase::SampleLocType m_current_cones; - SamplerBase::SampleLocType m_prior_cones; - SamplerBase::SampleLocType m_sample_cones; + SampleLocType::LocType m_initial_cone; + SampleLocType::LocType m_current_cones; + SampleLocType::LocType m_prior_cones; + SampleLocType::LocType m_sample_cones; LosRotType m_los_proj; LosUnitType m_los_unit; @@ -193,7 +200,7 @@ protected: long m_num_points_prior; long m_num_output_points_prior; - double m_fill_val{-99999.99}; + amrex::Real m_fill_val{-99999.99}; }; } // namespace amr_wind::sampling diff --git a/amr-wind/utilities/sampling/RadarSampler.cpp b/amr-wind/utilities/sampling/RadarSampler.cpp index 4f98bdc579..4502581bc3 100644 --- a/amr-wind/utilities/sampling/RadarSampler.cpp +++ b/amr-wind/utilities/sampling/RadarSampler.cpp @@ -2,11 +2,9 @@ #include "amr-wind/CFDSim.H" #include "amr-wind/utilities/tensor_ops.H" #include "amr-wind/utilities/trig_ops.H" - +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" -#define PI 3.14159265 - namespace amr_wind::sampling { RadarSampler::RadarSampler(const CFDSim& sim) : m_sim(sim) {} @@ -219,8 +217,8 @@ bool RadarSampler::update_sampling_locations() << "-------------------------" << std::endl; } - long conetipbegin = m_cone_size - 1 - num_points_quad(); - long conetipend = m_cone_size; + const amrex::Long conetipbegin = m_cone_size - 1 - num_points_quad(); + const amrex::Long conetipend = m_cone_size; // Loop for oversampling for (int k = 0; k < m_ntotal; k++) { @@ -251,7 +249,7 @@ bool RadarSampler::update_sampling_locations() vertical_ref, sweep_angle, radar_ref)); vs::Vector elevation_axis(-vertical_ref ^ swept_axis); - long cq_idx = k * num_points_quad(); + amrex::Long cq_idx = k * num_points_quad(); // Add rotated cone to current cones for (int i = 0; i < m_cone_size; ++i) { @@ -263,7 +261,7 @@ bool RadarSampler::update_sampling_locations() vertical_ref, sweep_angle, temp_point)); vs::Vector rotated_point(sampling_utils::rotate_euler_vector( elevation_axis, elevation_angle, swept_point)); - long point_index = i + k * m_cone_size; + const amrex::Long point_index = i + k * m_cone_size; m_current_cones[point_index][0] = rotated_point[0] + m_start[0]; m_current_cones[point_index][1] = rotated_point[1] + m_start[1]; m_current_cones[point_index][2] = rotated_point[2] + m_start[2]; @@ -284,10 +282,10 @@ bool RadarSampler::update_sampling_locations() } else { // This cone falls outside time bounds // For this timestep - long cq_idx = k * num_points_quad(); + amrex::Long cq_idx = k * num_points_quad(); for (int i = 0; i < m_cone_size; ++i) { - long point_index = i + k * m_cone_size; + const amrex::Long point_index = i + k * m_cone_size; m_current_cones[point_index][0] = m_fill_val; m_current_cones[point_index][1] = m_fill_val; m_current_cones[point_index][2] = m_fill_val; @@ -324,7 +322,7 @@ void RadarSampler::new_cone() utils::radians(m_cone_angle), m_ntheta, sampling_utils::NormalRule::HALFPOWER, m_rays, m_weights); - long nquad = num_points_quad(); + const amrex::Long nquad = num_points_quad(); if (m_debug_print) { for (int j = 0; j < nquad; ++j) { @@ -351,7 +349,7 @@ void RadarSampler::new_cone() for (int n = 0; n < m_npts; ++n) { for (int j = 0; j < nquad; ++j) { - long pt_idx = j + n * nquad; + const amrex::Long pt_idx = j + n * nquad; const auto radius = dx[2] * n; m_initial_cone[pt_idx][0] = radius * m_rays[j][0]; m_initial_cone[pt_idx][1] = radius * m_rays[j][1]; @@ -360,7 +358,7 @@ void RadarSampler::new_cone() } vs::Vector tc_axis = canon_vector ^ init_axis; - amrex::Real tc_angle = std::acos(init_axis & canon_vector) * 180.0 / PI; + amrex::Real tc_angle = std::acos(init_axis & canon_vector) * 180.0 / M_PI; // Initial cone points along input-file-given axis and is full size for (int i = 0; i < num_points_cone(); ++i) { @@ -388,8 +386,8 @@ void RadarSampler::calc_lineofsight_velocity( for (int k = 0; k < m_ntotal; k++) { for (int i = 0; i < m_cone_size; ++i) { - long p_idx = i + k * m_cone_size; - long cq_idx = i % num_points_quad(); + const amrex::Long p_idx = i + k * m_cone_size; + const amrex::Long cq_idx = i % num_points_quad(); vs::Vector temp_vel( velocity_raw[p_idx][0], velocity_raw[p_idx][1], velocity_raw[p_idx][2]); @@ -399,7 +397,7 @@ void RadarSampler::calc_lineofsight_velocity( } for (int k = 0; k < m_ntotal; k++) { - long a_start = k * num_points_axis(); + const amrex::Long a_start = k * num_points_axis(); std::vector temp_vals( los_temp.begin() + k * num_points_cone(), los_temp.begin() + (k + 1) * num_points_cone()); @@ -419,13 +417,13 @@ std::vector RadarSampler::modify_sample_data( // there are m_ntotal steps (cones) based on sampling rate AMREX_ALWAYS_ASSERT(static_cast(sample_data.size()) == num_points()); - const long n_cones = m_ntotal; + const amrex::Long n_cones = m_ntotal; std::vector mod_data(num_output_points()); for (int ic = 0; ic < n_cones; ++ic) { - long c_start = ic * num_points_cone(); - long c_end = (ic + 1) * num_points_cone(); - long a_start = ic * num_points_axis(); + const amrex::Long c_start = ic * num_points_cone(); + const amrex::Long c_end = (ic + 1) * num_points_cone(); + const amrex::Long a_start = ic * num_points_axis(); // Send a single cone to be line averaged std::vector temp_vals( @@ -440,7 +438,7 @@ void RadarSampler::line_average( const std::vector& weights, const std::vector& values, std::vector& reduced, - long offset) + const amrex::Long offset) { const int nline = static_cast(values.size() / weights.size()); const int nquad = static_cast(weights.size()); @@ -460,38 +458,67 @@ void RadarSampler::line_average( } } -void RadarSampler::sampling_locations(SampleLocType& locs) const +void RadarSampler::sampling_locations(SampleLocType& sample_locs) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void RadarSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const { - locs.resize(num_points()); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); for (int i = 0; i < num_points_scan(); ++i) { - locs[i][0] = (m_radar_iter > 0) ? m_prior_cones[i][0] : m_fill_val; - locs[i][1] = (m_radar_iter > 0) ? m_prior_cones[i][1] : m_fill_val; - locs[i][2] = (m_radar_iter > 0) ? m_prior_cones[i][2] : m_fill_val; + const amrex::RealVect loc = {AMREX_D_DECL( + (m_radar_iter > 0) ? m_prior_cones[i][0] : m_fill_val, + (m_radar_iter > 0) ? m_prior_cones[i][1] : m_fill_val, + (m_radar_iter > 0) ? m_prior_cones[i][2] : m_fill_val)}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, i); + break; + } + } } for (int i = 0; i < num_points_scan(); ++i) { - long ioff = i + num_points_scan(); - locs[ioff][0] = m_current_cones[i][0]; - locs[ioff][1] = m_current_cones[i][1]; - locs[ioff][2] = m_current_cones[i][2]; + const amrex::Long ioff = i + num_points_scan(); + const amrex::RealVect loc = {AMREX_D_DECL( + m_current_cones[i][0], m_current_cones[i][1], + m_current_cones[i][2])}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, ioff); + break; + } + } } } // TODO: Create test for this calculation void RadarSampler::cone_axis_locations(SampleLocType& axis_locs) const { - axis_locs.resize(num_output_points()); - for (int k = 0; k < m_ntotal; k++) { for (int i = 0; i < m_npts; ++i) { - int pi = i + k * m_npts; - long ci = i * num_points_quad() + k * num_points_cone(); - axis_locs[pi][0] = m_prior_cones[ci][0]; - axis_locs[pi][1] = m_prior_cones[ci][1]; - axis_locs[pi][2] = m_prior_cones[ci][2]; + const amrex::Long pi = i + k * m_npts; + const amrex::Long ci = + i * num_points_quad() + k * num_points_cone(); + const amrex::RealVect loc = {AMREX_D_DECL( + m_prior_cones[ci][0], m_prior_cones[ci][1], + m_prior_cones[ci][2])}; + axis_locs.push_back(loc, pi); } } + AMREX_ALWAYS_ASSERT(axis_locs.locations().size() == num_output_points()); } void RadarSampler::post_sample_actions() @@ -537,21 +564,23 @@ void RadarSampler::output_netcdf_data( std::vector start{nt, 0, 0}; std::vector count{1, 0, AMREX_SPACEDIM}; - SamplerBase::SampleLocType locs; - cone_axis_locations(locs); + SampleLocType sample_locs; + cone_axis_locations(sample_locs); count[1] = num_output_points(); auto xyz = grp.var("points"); - xyz.put(locs[0].data(), start, count); + const auto& locs = sample_locs.locations(); + xyz.put(locs[0].begin(), start, count); if (m_output_cone_points) { std::vector cone_start{nt, 0, 0}; std::vector cone_count{1, 0, AMREX_SPACEDIM}; - SamplerBase::SampleLocType conelocs; + SampleLocType conelocs; sampling_locations(conelocs); cone_count[1] = num_points(); auto cone_xyz = grp.var("conepoints"); - cone_xyz.put(conelocs[0].data(), cone_start, cone_count); + const auto& clocs = sample_locs.locations(); + cone_xyz.put(clocs[0].begin(), cone_start, cone_count); } } diff --git a/amr-wind/utilities/sampling/SamplerBase.H b/amr-wind/utilities/sampling/SamplerBase.H index fe108473ca..23e5f21d4f 100644 --- a/amr-wind/utilities/sampling/SamplerBase.H +++ b/amr-wind/utilities/sampling/SamplerBase.H @@ -1,6 +1,7 @@ #ifndef SAMPLERBASE_H #define SAMPLERBASE_H +#include #include "amr-wind/core/Factory.H" #include "amr-wind/utilities/ncutils/nc_interface.H" @@ -10,6 +11,39 @@ class CFDSim; namespace sampling { +struct SampleLocType +{ + using LocType = amrex::Vector; + using IdType = amrex::Vector; + + void push_back(const amrex::RealVect& loc, const amrex::Long id) + { + m_locations.push_back(loc); + m_ids.push_back(id); + } + + const LocType& locations() const + { + AMREX_ALWAYS_ASSERT(m_locations.size() == m_ids.size()); + return m_locations; + } + + LocType& locations() + { + AMREX_ALWAYS_ASSERT(m_locations.size() == m_ids.size()); + return m_locations; + } + + const IdType& ids() const + { + AMREX_ALWAYS_ASSERT(m_locations.size() == m_ids.size()); + return m_ids; + } + + LocType m_locations; + IdType m_ids; +}; + /** Abstract representation of data probes to sample flow data * \ingroup sampling * @@ -20,9 +54,6 @@ namespace sampling { class SamplerBase : public Factory { public: - using SampleLocType = - amrex::Vector>; - static std::string base_identifier() { return "SamplerBase"; } static constexpr amrex::Real bounds_tol = @@ -55,6 +86,11 @@ public: //! Populate the vector with coordinates of the sampling locations virtual void sampling_locations(SampleLocType&) const = 0; + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + virtual void sampling_locations( + SampleLocType&, const amrex::Vector&) const = 0; + //! Check and fix the bounds of the sampler so the probes are in the domain virtual void check_bounds() = 0; diff --git a/amr-wind/utilities/sampling/Sampling.cpp b/amr-wind/utilities/sampling/Sampling.cpp index a15d181aef..3d6e5c5d90 100644 --- a/amr-wind/utilities/sampling/Sampling.cpp +++ b/amr-wind/utilities/sampling/Sampling.cpp @@ -470,14 +470,15 @@ void Sampling::prepare_netcdf_file() { const std::vector start{0, 0}; std::vector count{0, AMREX_SPACEDIM}; - SamplerBase::SampleLocType locs; for (const auto& obj : m_samplers) { auto grp = ncf.group(obj->label()); obj->populate_netcdf_metadata(grp); - obj->output_locations(locs); + SampleLocType sample_locs; + obj->output_locations(sample_locs); auto xyz = grp.var("coordinates"); count[0] = obj->num_output_points(); - xyz.put(locs[0].data(), start, count); + const auto& locs = sample_locs.locations(); + xyz.put(locs[0].begin(), start, count); } } diff --git a/amr-wind/utilities/sampling/SamplingContainer.cpp b/amr-wind/utilities/sampling/SamplingContainer.cpp index 7fc7e127c4..560db791b5 100644 --- a/amr-wind/utilities/sampling/SamplingContainer.cpp +++ b/amr-wind/utilities/sampling/SamplingContainer.cpp @@ -30,16 +30,8 @@ void SamplingContainer::initialize_particles( { BL_PROFILE("amr-wind::SamplingContainer::initialize"); - // We will assign all particles to the first box in level 0 and let - // redistribute scatter it to the appropriate rank and box. const int lev = 0; - const int iproc = amrex::ParallelDescriptor::MyProc(); - const int owner = ParticleDistributionMap(lev)[0]; - - // Let only the MPI rank owning the first box do the work - if (owner != iproc) { - return; - } + const auto iproc = amrex::ParallelDescriptor::MyProc(); long num_particles = 0; for (const auto& probes : samplers) { @@ -47,45 +39,179 @@ void SamplingContainer::initialize_particles( } m_total_particles = num_particles; - const int grid_id = 0; - const int tile_id = 0; - // Setup has already processed runtime array information, safe to do - // GetParticles here. - auto& ptile = GetParticles(lev)[std::make_pair(grid_id, tile_id)]; - ptile.resize(num_particles); - - int pidx = 0; - const int nextid = static_cast(ParticleType::NextID()); - auto* pstruct = ptile.GetArrayOfStructs()().data(); - SamplerBase::SampleLocType locs; - for (const auto& probe : samplers) { - probe->sampling_locations(locs); + const int nprobes = static_cast(samplers.size()); + amrex::iMultiFab particle_counts( + ParticleBoxArray(lev), ParticleDistributionMap(lev), nprobes, 0, + amrex::MFInfo()); + amrex::iMultiFab offsets( + ParticleBoxArray(lev), ParticleDistributionMap(lev), nprobes, 0, + amrex::MFInfo()); + particle_counts.setVal(0); + offsets.setVal(0); + + const auto& dxinv = m_mesh.Geom(lev).InvCellSizeArray(); + const auto& plo = m_mesh.Geom(lev).ProbLoArray(); + for (int iprobe = 0; iprobe < nprobes; iprobe++) { + const auto& probe = samplers[iprobe]; + SampleLocType sample_locs; + probe->sampling_locations(sample_locs); + const auto& locs = sample_locs.locations(); const int npts = static_cast(locs.size()); - const auto probe_id = probe->id(); - amrex::Gpu::DeviceVector> - dlocs(npts); + amrex::Gpu::DeviceVector dlocs(npts); + amrex::Gpu::copy( amrex::Gpu::hostToDevice, locs.begin(), locs.end(), dlocs.begin()); const auto* dpos = dlocs.data(); - amrex::ParallelFor(npts, [=] AMREX_GPU_DEVICE(const int ip) noexcept { - const auto uid = pidx + ip; - auto& pp = pstruct[uid]; - pp.id() = nextid + uid; - pp.cpu() = iproc; + // don't use openmp + for (amrex::MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const amrex::Box& box = mfi.tilebox(); + + // count the number of particles in each cell + const auto& np_arr = particle_counts[mfi].array(); + amrex::ParallelFor( + dlocs.size(), [=] AMREX_GPU_DEVICE(long ip) noexcept { + const amrex::RealVect pos( + AMREX_D_DECL(dpos[ip][0], dpos[ip][1], dpos[ip][2])); - for (int n = 0; n < AMREX_SPACEDIM; ++n) { - pp.pos(n) = dpos[ip][n]; - } - pp.idata(IIx::uid) = uid; - pp.idata(IIx::sid) = probe_id; - pp.idata(IIx::nid) = ip; - }); - amrex::Gpu::streamSynchronize(); - pidx += npts; + const amrex::IntVect div(AMREX_D_DECL( + static_cast( + amrex::Math::floor((pos[0] - plo[0]) * dxinv[0])), + static_cast( + amrex::Math::floor((pos[1] - plo[1]) * dxinv[1])), + static_cast( + amrex::Math::floor((pos[2] - plo[2]) * dxinv[2])))); + if (box.contains(div)) { + amrex::Gpu::Atomic::AddNoRet(&np_arr(div, iprobe), 1); + } + }); + } + AMREX_ALWAYS_ASSERT(particle_counts.sum(iprobe) == npts); + } + + // compute the offsets and size the particle tiles + // don't use openmp + for (amrex::MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const amrex::Box& box = mfi.tilebox(); + const auto ncells = static_cast(box.numPts()); + + const auto& np_arr = particle_counts[mfi].const_array(); + int* p_offsets = offsets[mfi].dataPtr(); + const auto np_box = amrex::Scan::PrefixSum( + ncells, + [=] AMREX_GPU_DEVICE(int i) -> int { + const auto iv = box.atOffset(i); + int total_np = 0; + for (int iprobe = 0; iprobe < nprobes; iprobe++) { + total_np += np_arr(iv, iprobe); + } + return total_np; + }, + [=] AMREX_GPU_DEVICE(int i, const int& xi) { p_offsets[i] = xi; }, + amrex::Scan::Type::exclusive); + + const auto& offsets_arr = offsets[mfi].array(); + amrex::ParallelFor( + box, [=] AMREX_GPU_DEVICE( + int i, int j, int AMREX_D_PICK(, , k)) noexcept { + const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); + for (int iprobe = 1; iprobe < nprobes; iprobe++) { + offsets_arr(iv, iprobe) = + offsets_arr(iv, iprobe - 1) + np_arr(iv, iprobe - 1); + } + }); + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + auto& ptile = GetParticles(lev)[std::make_pair(grid_id, tile_id)]; + ptile.resize(np_box); + + AMREX_ASSERT( + np_box == + particle_counts[mfi].sum(box, 0, nprobes)); } - AMREX_ALWAYS_ASSERT(pidx == num_particles); + for (int iprobe = 0; iprobe < nprobes; iprobe++) { + const auto& probe = samplers[iprobe]; + const auto probe_id = probe->id(); + SampleLocType sample_locs; + probe->sampling_locations(sample_locs); + const auto& locs = sample_locs.locations(); + const int npts = static_cast(locs.size()); + amrex::Gpu::DeviceVector dlocs(npts); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, locs.begin(), locs.end(), dlocs.begin()); + const auto* p_dlocs = dlocs.data(); + const auto& ids = sample_locs.ids(); + amrex::Gpu::DeviceVector dids(npts); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, ids.begin(), ids.end(), dids.begin()); + const auto* p_dids = dids.data(); + + // don't use openmp + for (amrex::MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { + const amrex::Box& box = mfi.tilebox(); + const auto& np_arr = particle_counts[mfi].const_array(); + const auto& offset_arr = offsets[mfi].const_array(); + + const int grid_id = mfi.index(); + const int tile_id = mfi.LocalTileIndex(); + auto& ptile = GetParticles(lev)[std::make_pair(grid_id, tile_id)]; + if (ptile.numParticles() == 0) { + continue; + } + + const int np_probe = + particle_counts[mfi].sum(box, iprobe, 1); + if (np_probe == 0) { + continue; + } + + const amrex::Long nextid = ParticleType::NextID(); + ParticleType::NextID(nextid + np_probe); + + auto* pstruct = ptile.GetArrayOfStructs()().data(); + amrex::ParallelFor( + box, [=] AMREX_GPU_DEVICE( + int i, int j, int AMREX_D_PICK(, , k)) noexcept { + const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); + + if (np_arr(iv, iprobe) > 0) { + const int start = offset_arr(iv, iprobe); + int n = start; + for (int ip = 0; ip < npts; ip++) { + const amrex::RealVect loc(AMREX_D_DECL( + p_dlocs[ip][0], p_dlocs[ip][1], + p_dlocs[ip][2])); + + const amrex::IntVect div(AMREX_D_DECL( + static_cast(amrex::Math::floor( + (loc[0] - plo[0]) * dxinv[0])), + static_cast(amrex::Math::floor( + (loc[1] - plo[1]) * dxinv[1])), + static_cast(amrex::Math::floor( + (loc[2] - plo[2]) * dxinv[2])))); + if (div == iv) { + auto& pp = pstruct[n]; + pp.id() = nextid + n; + pp.cpu() = iproc; + for (int idim = 0; idim < AMREX_SPACEDIM; + ++idim) { + pp.pos(idim) = loc[idim]; + } + pp.idata(IIx::uid) = ip + npts * iprobe; + pp.idata(IIx::sid) = probe_id; + pp.idata(IIx::nid) = + static_cast(p_dids[ip]); + n++; + } + } + AMREX_ALWAYS_ASSERT(np_arr(iv, iprobe) == (n - start)); + } + }); + amrex::Gpu::streamSynchronize(); + } + } } void SamplingContainer::interpolate_derived_fields( diff --git a/amr-wind/utilities/sampling/VolumeSampler.H b/amr-wind/utilities/sampling/VolumeSampler.H index a531ce7a46..4738c296ae 100644 --- a/amr-wind/utilities/sampling/VolumeSampler.H +++ b/amr-wind/utilities/sampling/VolumeSampler.H @@ -32,11 +32,17 @@ public: void check_bounds() override; //! Populate and return a vector of probe locations to be sampled - void sampling_locations(SampleLocType& /*locs*/) const override; + void sampling_locations(SampleLocType& /*sample_locs**/) const override; - void output_locations(SampleLocType& locs) const override + //! Populate the vector with coordinates of the sampling locations inside + //! boxes + void sampling_locations( + SampleLocType& /*sample_locs**/, + const amrex::Vector& /*boxes*/) const override; + + void output_locations(SampleLocType& sample_locs) const override { - return sampling_locations(locs); + return sampling_locations(sample_locs); } void diff --git a/amr-wind/utilities/sampling/VolumeSampler.cpp b/amr-wind/utilities/sampling/VolumeSampler.cpp index e03130f578..3d38a23f5b 100644 --- a/amr-wind/utilities/sampling/VolumeSampler.cpp +++ b/amr-wind/utilities/sampling/VolumeSampler.cpp @@ -2,7 +2,7 @@ #include "amr-wind/utilities/sampling/VolumeSampler.H" #include "amr-wind/CFDSim.H" - +#include "amr-wind/utilities/index_operations.H" #include "AMReX_ParmParse.H" namespace amr_wind::sampling { @@ -66,10 +66,25 @@ void VolumeSampler::check_bounds() } } -void VolumeSampler::sampling_locations(SampleLocType& locs) const +void VolumeSampler::sampling_locations(SampleLocType& sample_locs) const +{ + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + + const int lev = 0; + const auto domain = m_sim.mesh().Geom(lev).Domain(); + sampling_locations(sample_locs, {domain}); + + AMREX_ALWAYS_ASSERT(sample_locs.locations().size() == num_points()); +} + +void VolumeSampler::sampling_locations( + SampleLocType& sample_locs, const amrex::Vector& boxes) const { - locs.resize(m_npts); + AMREX_ALWAYS_ASSERT(sample_locs.locations().empty()); + const int lev = 0; + const auto& dxinv = m_sim.mesh().Geom(lev).InvCellSizeArray(); + const auto& plo = m_sim.mesh().Geom(lev).ProbLoArray(); const amrex::Array dx = { (m_hi[0] - m_lo[0]) / m_npts_dir[0], (m_hi[1] - m_lo[1]) / m_npts_dir[1], @@ -79,9 +94,15 @@ void VolumeSampler::sampling_locations(SampleLocType& locs) const for (int k = 0; k < m_npts_dir[2]; ++k) { for (int j = 0; j < m_npts_dir[1]; ++j) { for (int i = 0; i < m_npts_dir[0]; ++i) { - locs[idx][0] = m_lo[0] + dx[0] * i; - locs[idx][1] = m_lo[1] + dx[1] * j; - locs[idx][2] = m_lo[2] + dx[2] * k; + const amrex::RealVect loc = {AMREX_D_DECL( + m_lo[0] + dx[0] * i, m_lo[1] + dx[1] * j, + m_lo[2] + dx[2] * k)}; + for (const auto& box : boxes) { + if (utils::contains(box, loc, plo, dxinv)) { + sample_locs.push_back(loc, idx); + break; + } + } ++idx; } } diff --git a/test/test_files/abl_sampling/abl_sampling.inp b/test/test_files/abl_sampling/abl_sampling.inp index 04f6105de9..02af083c13 100644 --- a/test/test_files/abl_sampling/abl_sampling.inp +++ b/test/test_files/abl_sampling/abl_sampling.inp @@ -54,6 +54,7 @@ ABL.surface_roughness_z0 = 0.15 #.......................................# amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel amr.max_level = 0 # Max AMR level in hierarchy +amr.max_grid_size = 8 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # GEOMETRY # @@ -83,11 +84,15 @@ sampling.output_format = native sampling.fields = velocity sampling.int_fields = mask_cell sampling.derived_fields = mag_vorticity "components(velocity,0,1)" "grad(temperature)" "grad(density)" -sampling.labels = volume1 +sampling.labels = volume1 volume2 sampling.volume1.type = VolumeSampler sampling.volume1.lo = 200.0 200.0 200.0 sampling.volume1.hi = 500.0 500.0 500.0 sampling.volume1.num_points = 30 10 10 +sampling.volume2.type = VolumeSampler +sampling.volume2.lo = 205.0 200.0 200.0 +sampling.volume2.hi = 505.0 500.0 500.0 +sampling.volume2.num_points = 30 10 10 sampling2.output_frequency = 1 sampling2.output_format = native diff --git a/tools/fcompare_particles.py b/tools/fcompare_particles.py index fa8c44aeb6..b41c4e2190 100644 --- a/tools/fcompare_particles.py +++ b/tools/fcompare_particles.py @@ -5,7 +5,6 @@ import amrex_particle from amrex_particle import AmrexParticleFile - def main(): """Compare particle files""" parser = argparse.ArgumentParser(description="A comparison tool for particles") @@ -30,9 +29,11 @@ def main(): assert Path(args.f0).is_dir() assert Path(args.f1).is_dir() - cols_to_drop = ["uid", "set_id", "probe_id"] - p0df = AmrexParticleFile(Path(args.f0) / "particles")().drop(cols_to_drop, axis=1) - p1df = AmrexParticleFile(Path(args.f1) / "particles")().drop(cols_to_drop, axis=1) + p0df = AmrexParticleFile(Path(args.f0) / "particles")() + p1df = AmrexParticleFile(Path(args.f1) / "particles")() + assert p0df.shape == p1df.shape + p0df.sort_values(by=["uid"], inplace=True, kind="stable", ignore_index=True) + p1df.sort_values(by=["uid"], inplace=True, kind="stable", ignore_index=True) adiff = np.sqrt(np.square(p0df - p1df).sum(axis=0)) rdiff = np.sqrt(np.square(p0df - p1df).sum(axis=0)) / np.sqrt( diff --git a/unit_tests/utilities/test_free_surface.cpp b/unit_tests/utilities/test_free_surface.cpp index d8b7cbc680..98af291bb9 100644 --- a/unit_tests/utilities/test_free_surface.cpp +++ b/unit_tests/utilities/test_free_surface.cpp @@ -298,13 +298,14 @@ int FreeSurfaceImpl::check_sloc(const std::string& op) { // Get number of points and sampling locations array auto npts_tot = num_points(); - amrex::Vector> locs; - output_locations(locs); + amr_wind::sampling::SampleLocType sample_locs; + output_locations(sample_locs); // Get locations from other functions auto gridlocs = grid_locations(); auto out = heights(); // Loop through grid points and check output int icheck = 0; + const auto& locs = sample_locs.locations(); for (int n = 0; n < npts_tot; ++n) { if (op == "=") { EXPECT_EQ(locs[n][0], gridlocs[n][0]); diff --git a/unit_tests/utilities/test_sampling.cpp b/unit_tests/utilities/test_sampling.cpp index dbc453a910..641d76ad99 100644 --- a/unit_tests/utilities/test_sampling.cpp +++ b/unit_tests/utilities/test_sampling.cpp @@ -280,10 +280,10 @@ TEST_F(SamplingTest, plane_sampler) amr_wind::sampling::PlaneSampler plane(sim()); plane.initialize("plane"); - amr_wind::sampling::PlaneSampler::SampleLocType locs; - plane.sampling_locations(locs); + amr_wind::sampling::SampleLocType sample_locs; + plane.sampling_locations(sample_locs); - ASSERT_EQ(locs.size(), 3 * 3 * 2); + ASSERT_EQ(sample_locs.locations().size(), 3 * 3 * 2); } TEST_F(SamplingTest, volume_sampler) @@ -296,10 +296,10 @@ TEST_F(SamplingTest, volume_sampler) amr_wind::sampling::VolumeSampler volume(sim()); volume.initialize("volume"); - amr_wind::sampling::VolumeSampler::SampleLocType locs; - volume.sampling_locations(locs); + amr_wind::sampling::SampleLocType sample_locs; + volume.sampling_locations(sample_locs); - ASSERT_EQ(locs.size(), 3 * 5 * 5); + ASSERT_EQ(sample_locs.locations().size(), 3 * 5 * 5); } TEST_F(SamplingTest, spinner_sampler) @@ -327,10 +327,10 @@ TEST_F(SamplingTest, spinner_sampler) amr_wind::sampling::DTUSpinnerSampler spinner(sim()); spinner.initialize("spinner"); - amr_wind::sampling::DTUSpinnerSampler::SampleLocType locs; - spinner.sampling_locations(locs); + amr_wind::sampling::SampleLocType sample_locs; + spinner.sampling_locations(sample_locs); - ASSERT_EQ(locs.size(), 21600); + ASSERT_EQ(sample_locs.locations().size(), 21600); } TEST_F(SamplingTest, radar_sampler) @@ -359,10 +359,10 @@ TEST_F(SamplingTest, radar_sampler) amr_wind::sampling::RadarSampler radar(sim()); radar.initialize("radar"); - amr_wind::sampling::RadarSampler::SampleLocType locs; - radar.sampling_locations(locs); + amr_wind::sampling::SampleLocType sample_locs; + radar.sampling_locations(sample_locs); - ASSERT_EQ(locs.size(), 193536); + ASSERT_EQ(sample_locs.locations().size(), 193536); } TEST_F(SamplingTest, sampling_utils)