Remove sampling spike

amr-wind/utilities/index_operations.H

@@ -4,6 +4,7 @@
 #include <AMReX_Vector.H>
 #include <AMReX_Box.H>
 #include <AMReX_RealBox.H>
+#include <AMReX_RealVect.H>
 #include <AMReX_Geometry.H>
 
 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<amrex::Real, AMREX_SPACEDIM>& prob_lo,
+    const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dxinv)
+{
+    const amrex::IntVect iv(AMREX_D_DECL(
+        static_cast<int>(amrex::Math::floor((pos[0] - prob_lo[0]) * dxinv[0])),
+        static_cast<int>(amrex::Math::floor((pos[1] - prob_lo[1]) * dxinv[1])),
+        static_cast<int>(
+            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

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<amrex::Box>& /*boxes*/) const override;
 
     static vs::Vector generate_lidar_pattern(
         PrismParameters InnerPrism, PrismParameters OuterPrism, double time);

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<amrex::Box>& 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<amrex::Real, AMREX_SPACEDIM> 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<size_t> count{1, 0, AMREX_SPACEDIM};
     std::vector<size_t> starti{nt, 0};
     std::vector<size_t> 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;
 

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<amrex::Box>& /*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

amr-wind/utilities/sampling/FreeSurfaceSampler.cpp

@@ -3,6 +3,7 @@
 #include <AMReX_MultiFabUtil.H>
 #include <utility>
 #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<amrex::Box>& 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<size_t> start{nt, 0, 0};
     std::vector<size_t> 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(

amr-wind/utilities/sampling/LidarSampler.cpp

@@ -109,11 +109,12 @@ void LidarSampler::output_netcdf_data(
     // Write the coordinates every time
     std::vector<size_t> start{nt, 0, 0};
     std::vector<size_t> 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(

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<amrex::Box>& /*boxes*/) const override;
+
+    void output_locations(SampleLocType& sample_locs) const override
     {
-        return sampling_locations(locs);
+        return sampling_locations(sample_locs);
     }
 
     void

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<amrex::Box>& 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<amrex::Real, AMREX_SPACEDIM> 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;
+            }
         }
     }
 }

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<amrex::Box>& /*boxes*/) const override;
+
+    void output_locations(SampleLocType& sample_locs) const override
     {
-        return sampling_locations(locs);
+        return sampling_locations(sample_locs);
     }
 
     void

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<amrex::Box>& boxes) const
+{
+    AMREX_ALWAYS_ASSERT(sample_locs.locations().empty());
 
     amrex::Array<amrex::Real, AMREX_SPACEDIM> dx;
     amrex::Array<amrex::Real, AMREX_SPACEDIM> dy;
@@ -89,12 +101,22 @@ void PlaneSampler::sampling_locations(SampleLocType& locs) const
 
     int idx = 0;
     const int nplanes = static_cast<int>(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;
             }