radar: add in-situ averaging version

docs/source/user/nalu_run/nalu_inp.rst

@@ -1586,6 +1586,22 @@ Data probes
    lines_per_cone_circle      Required. Number of rays around the cone circumference
    ========================== ===================================================================
 
+.. inpfile:: dataprobes.lidar_specifications.radar_cone_filter
+
+   Implements a few options for filtering the spherical cap of a cone. `truncated_normal{n}` rules with `n=1,2,3`
+   weight the filtering based on truncated normal distribution, with with the circle of the cone being 1,2, or 3
+   sigma away. This means that the sampling is more weighted toward the center of the cone with higher `n`. `radau`
+   has weight function = 1.
+
+   ========================== ===================================================================
+   Parameter                  Description
+   ========================== ===================================================================
+   cone_angle                 Required. cone half angle in degrees centered on radar_specifications.axis
+   quadrature_type            Required. Type of quadrature. `radau` or `truncated_normal1`, `truncated_normal2`, `truncated_normal3`.
+   radau_points               Optional. If `radau` quadrature is used, number of integration points
+   lines_per_cone_circle      Required. Number of rays around the cone circumference
+   ========================== ===================================================================
+
 .. inpfile:: dataprobes.lidar_specifications.misc
 
    The user may also set a number of parameters corresponding to the hardware

include/wind_energy/LidarPatterns.h

@@ -130,7 +130,7 @@ class RadarSegmentGenerator final : public SegmentGenerator
 public:
   void load(const YAML::Node& node) final;
   Segment generate(double t) const final;
-  void set_axis(vs::Vector axis);
+  vs::Vector center() const { return {center_[0], center_[1], center_[2]}; }
 
 private:
   enum class phase { FORWARD, FORWARD_PAUSE, REVERSE, REVERSE_PAUSE };

include/wind_energy/SyntheticLidar.h

@@ -22,6 +22,12 @@
 namespace sierra {
 namespace nalu {
 
+struct RadarFilter
+{
+  std::vector<vs::Vector> rays;
+  std::vector<double> weights;
+};
+
 class LidarLineOfSite
 {
 public:
@@ -40,6 +46,12 @@ class LidarLineOfSite
     const std::string& coordinates_name,
     double dtratio);
 
+  void output_cone_filtered(
+    const stk::mesh::BulkData& bulk,
+    const stk::mesh::Selector& active,
+    const std::string& coordinates_name,
+    double dtratio);
+
 private:
   enum class Output { NETCDF, TEXT, DATAPROBE } output_type_{Output::NETCDF};
   enum class Predictor {
@@ -58,6 +70,11 @@ class LidarLineOfSite
     double time,
     const std::vector<std::array<double, 3>>& x,
     const std::vector<std::array<double, 3>>& u);
+  void output_txt_filtered(
+    double time,
+    const std::vector<std::array<double, 3>>& x,
+    const std::vector<double>& u,
+    int n);
   std::map<std::string, int> ncVarIDs_;
 
   mutable double lidar_time_{0};
@@ -75,12 +92,23 @@ class LidarLineOfSite
   bool warn_on_missing_{false};
   bool reuse_search_data_{true};
   bool always_output_{false};
+
+  RadarFilter radar_data_;
 };
 
 namespace details {
 std::vector<vs::Vector>
 make_radar_grid(double phi, int nphi, int ntheta, vs::Vector axis);
-}
+
+RadarFilter parse_radar_filter(const YAML::Node& node);
+
+std::pair<std::vector<vs::Vector>, std::vector<double>>
+spherical_cap_radau(double gammav, int ntheta, int nphi);
+
+std::pair<std::vector<vs::Vector>, std::vector<double>>
+spherical_cap_truncated_normal(double gammav, int ntheta, int nsigma);
+
+} // namespace details
 
 class LidarLOS
 {
@@ -97,6 +125,7 @@ class LidarLOS
 
 private:
   std::vector<LidarLineOfSite> lidars_;
+  std::vector<LidarLineOfSite> radars_;
 };
 
 } // namespace nalu

reg_tests/test_files/ablNeutralEdgeNoSlip/ablNeutralEdgeNoSlip.yaml

@@ -348,6 +348,26 @@ realms:
               beam_length: 12500 # in m, maximum possible length of beam from the center
               elevation_angles: [0,1,2,3] # elevation angles to scan after reset
 
+          - name: radar/radar-beam-filter
+            type: radar
+            frequency: 4
+            points_along_line: 5
+            always_output: yes
+            output: text
+            radar_cone_filter:
+              cone_angle: 0.25 #half angle of cone, in degrees
+              quadrature_type: radau # type of quadrature. `radau` or `trunc_normal{1,2,3}`
+              radau_points: 9 # number of points for radau -> 8 circles
+              lines_per_cone_circle: 6 # number of lines around the cone circle
+            radar_specifications:
+              bbox: [0., 0., 0., 5000., 5000., 1000.] # where to cut the radar line(m)
+              center: [-10000, 2500, 100] #ideally outside of the box
+              angular_speed: 1 #degrees to sweep in a second
+              axis: [1, 0, 0] #main axis, cone grid is set around this
+              sweep_angle: 10 #angle of the sweep
+              reset_time:  0 # time to pause after a sweep
+              beam_length: 12500 # in m, maximum possible length of beam from the center
+              elevation_angles: [0,1,2,3] # elevation angles to scan after reset
 
     # This defines the ABL forcing to drive the winds to 8 m/s from
     # 245 degrees (southwest) at 90 m above the surface in a planar

src/wind_energy/LidarPatterns.C

@@ -174,23 +174,27 @@ ScanningLidarSegmentGenerator::load(const YAML::Node& node)
   normalize_vec3(axis_.data());
 
   double sweep_angle_in_degrees = 20;
-  get_if_present(node, "sweep_angle", sweep_angle_in_degrees);
+  get_if_present(
+    node, "sweep_angle", sweep_angle_in_degrees, sweep_angle_in_degrees);
   ThrowRequireMsg(sweep_angle_in_degrees > 0, "Sweep angle must be positive");
   sweep_angle_ = convert::degrees_to_radians(sweep_angle_in_degrees);
 
   double step_delta_angle_in_degrees = 1;
-  get_if_present(node, "step_delta_angle", step_delta_angle_in_degrees);
+  get_if_present(
+    node, "step_delta_angle", step_delta_angle_in_degrees,
+    step_delta_angle_in_degrees);
   step_delta_angle_ = convert::degrees_to_radians(step_delta_angle_in_degrees);
 
   ThrowRequireMsg(step_delta_angle_ > 0, "step delta angle must be positive");
   ThrowRequireMsg(
     step_delta_angle_ <= sweep_angle_,
     "step delta angle must be less than full sweep");
 
-  get_if_present(node, "stare_time", stare_time_);
+  get_if_present(node, "stare_time", stare_time_, stare_time_);
   ThrowRequireMsg(stare_time_ > 0, "stare time must be positive");
 
-  get_if_present(node, "reset_time_delta", reset_time_delta_);
+  get_if_present(
+    node, "reset_time_delta", reset_time_delta_, reset_time_delta_);
   ThrowRequireMsg(
     reset_time_delta_ >= 0, "reset time delta must be semi-positive");
 
@@ -312,22 +316,26 @@ SpinnerLidarSegmentGenerator::load(const YAML::Node& node)
   normalize_vec3(laserAxis_.data());
 
   double innerPrismTheta0 = 90;
-  get_if_present(node, "inner_prism_initial_theta", innerPrismTheta0);
+  get_if_present(
+    node, "inner_prism_initial_theta", innerPrismTheta0, innerPrismTheta0);
 
   double innerPrismRot = 3.5;
-  get_if_present(node, "inner_prism_rotation_rate", innerPrismRot);
+  get_if_present(
+    node, "inner_prism_rotation_rate", innerPrismRot, innerPrismRot);
 
   double innerPrismAzi = 15.2;
-  get_if_present(node, "inner_prism_azimuth", innerPrismAzi);
+  get_if_present(node, "inner_prism_azimuth", innerPrismAzi, innerPrismAzi);
 
   double outerPrismTheta0 = 90;
-  get_if_present(node, "outer_prism_initial_theta", outerPrismTheta0);
+  get_if_present(
+    node, "outer_prism_initial_theta", outerPrismTheta0, outerPrismTheta0);
 
   double outerPrismRot = 6.5;
-  get_if_present(node, "outer_prism_rotation_rate", outerPrismRot);
+  get_if_present(
+    node, "outer_prism_rotation_rate", outerPrismRot, outerPrismRot);
 
   double outerPrismAzi = 15.2;
-  get_if_present(node, "outer_prism_azimuth", outerPrismAzi);
+  get_if_present(node, "outer_prism_azimuth", outerPrismAzi, outerPrismAzi);
 
   innerPrism_ = {
     convert::degrees_to_radians(innerPrismTheta0),
@@ -538,12 +546,13 @@ RadarSegmentGenerator::load(const YAML::Node& node)
   normalize_vec3(axis_.data());
 
   double sweep_angle_in_degrees = 20;
-  get_if_present(node, "sweep_angle", sweep_angle_in_degrees);
-  ThrowRequireMsg(sweep_angle_in_degrees > 0, "Sweep angle must be positive");
+  get_if_present(
+    node, "sweep_angle", sweep_angle_in_degrees, sweep_angle_in_degrees);
+  ThrowRequireMsg(sweep_angle_in_degrees >= 0, "Sweep angle must be semipositive");
   sweep_angle_ = convert::degrees_to_radians(sweep_angle_in_degrees);
 
   double angular_speed = 30; // deg/s
-  get_if_present(node, "angular_speed", angular_speed);
+  get_if_present(node, "angular_speed", angular_speed, angular_speed);
   angular_speed_ = convert::degrees_to_radians(angular_speed);
 
   beam_length_ = 50e3; // m
@@ -566,7 +575,8 @@ RadarSegmentGenerator::load(const YAML::Node& node)
     }
   }
 
-  get_if_present(node, "reset_time_delta", reset_time_delta_);
+  get_if_present(
+    node, "reset_time_delta", reset_time_delta_, reset_time_delta_);
   ThrowRequireMsg(
     reset_time_delta_ >= 0, "reset time delta must be semi-positive");