Add line of site sampling along a Spinner Lidar rosette pattern (#48)

include/DataProbePostProcessing.h

@@ -14,10 +14,14 @@
 #include <string>
 #include <vector>
 #include <utility>
+#include <memory>
 
 // stk_mesh/base/fem
 #include <stk_mesh/base/Selector.hpp>
 
+#include <stk_io/StkMeshIoBroker.hpp>
+
+
 // stk forwards
 namespace stk {
   namespace mesh {
@@ -85,6 +89,8 @@ class DataProbePostProcessing
   void load(
     const YAML::Node & node);
 
+  void add_external_data_probe_spec_info(DataProbeSpecInfo* dpsInfo);
+
   // setup part creation and nodal field registration (before populate_mesh())
   void setup();
 
@@ -106,11 +112,16 @@ class DataProbePostProcessing
   // create the transfer and hold the vector in the DataProbePostProcessing class
   void create_transfer();
 
+  // optionally create an exodus database
+  void create_exodus();
+
   // populate nodal field and output norms (if appropriate)
   void execute();
 
   // output to a file
-  void provide_output(const double currentTime);
+  void provide_output_txt(const double currentTime);
+  void provide_output_exodus(const double currentTime);
+
 
   // provide the inactive selector
   stk::mesh::Selector &get_inactive_selector();
@@ -138,6 +149,13 @@ class DataProbePostProcessing
 
   // hold the transfers
   Transfers *transfers_;
+
+private:
+  std::unique_ptr<stk::io::StkMeshIoBroker> io;
+
+  bool useExo_{false};
+  std::string exoName_;
+  size_t fileIndex_;
 };
 
 } // namespace nalu

include/master_element/TensorOps.h

@@ -23,6 +23,10 @@ namespace nalu{
     return (1.0e6*std::numeric_limits<double>::min());
   }
 
+  inline constexpr double small_positive_value() {
+    return 1.0e2*std::numeric_limits<double>::epsilon();
+  }
+
   template <typename ScalarType>
   KOKKOS_FORCEINLINE_FUNCTION ScalarType vecnorm_sq2(const ScalarType* x) {
     return (x[0] * x[0] + x[1] * x[1]);
@@ -33,6 +37,14 @@ namespace nalu{
     return (x[0] * x[0] + x[1] * x[1] + x[2] * x[2]);
   }
 
+  template <typename ScalarType>
+  void normalize_vec3(ScalarType* x) {
+    const ScalarType invmag = 1.0 / stk::math::sqrt(vecnorm_sq3(x));
+    x[0] *= invmag;
+    x[1] *= invmag;
+    x[2] *= invmag;
+  }
+
   template <typename ScalarType>
   ScalarType ddot(const ScalarType* u, const ScalarType* v, int n)
   {

include/user_functions/CMakeLists.txt

@@ -8,6 +8,7 @@ add_sources(GlobalHeaderList
    ConvectingTaylorVortexPressureAuxFunction.h
    ConvectingTaylorVortexVelocityAuxFunction.h
    FlowPastCylinderTempAuxFunction.h
+   GaussianWakeVelocityAuxFunction.h
    KovasznayPressureAuxFunction.h
    KovasznayVelocityAuxFunction.h
    LinearRampMeshDisplacementAuxFunction.h

include/wind_energy/SyntheticLidar.h

@@ -0,0 +1,92 @@
+
+#ifndef SyntheticLidar_H
+#define SyntheticLidar_H
+
+#include <DataProbePostProcessing.h>
+
+#include <memory>
+#include <array>
+
+namespace sierra {
+namespace nalu {
+
+struct Segment
+{
+  Segment() = default;
+  Segment(std::array<double, 3> tip, std::array<double,3> tail)
+  : tip_{tip}, tail_{tail} {};
+
+  std::array<double, 3> tip_{{}};
+  std::array<double, 3> tail_{{}};
+};
+
+struct PrismParameters{
+  PrismParameters() = default;
+  PrismParameters(double theta0, double rot, double azimuth)
+  : theta0_{theta0}, rot_{rot}, azimuth_{azimuth}
+  {};
+
+  double theta0_{0}; // rad
+  double rot_{0}; // rad / s
+  double azimuth_{0}; // rad
+};
+
+class SpinnerLidarSegmentGenerator
+{
+public:
+
+  SpinnerLidarSegmentGenerator() = default;
+
+  SpinnerLidarSegmentGenerator(
+    PrismParameters inner,
+    PrismParameters outer,
+    double in_beamLength);
+
+  void load(const YAML::Node& node);
+
+  Segment generate_path_segment(double time) const;
+
+  void set_inner_prism(PrismParameters innerPrism) { innerPrism_ = innerPrism; }
+  void set_inner_prism(double theta0, double rot, double azi) { innerPrism_ = {theta0, rot, azi}; }
+  void set_outer_prism(PrismParameters outerPrism) { outerPrism_ = outerPrism; }
+  void set_outer_prism(double theta0, double rot, double azi) { outerPrism_ = {theta0, rot, azi}; }
+  void set_lidar_center(Coordinates lidarCenter) { lidarCenter_ = {{lidarCenter.x_, lidarCenter.y_, lidarCenter.z_}}; }
+  void set_laser_axis(Coordinates laserAxis) { laserAxis_ = {{laserAxis.x_, laserAxis.y_, laserAxis.z_}}; }
+  void set_ground_normal(Coordinates gNormal) { groundNormal_ = {{gNormal.x_, gNormal.y_, gNormal.z_}}; }
+  void set_beam_length(double beamLength) { beamLength_ = beamLength; }
+
+private:
+  PrismParameters innerPrism_;
+  PrismParameters outerPrism_;
+  double beamLength_{1.0};
+
+  std::array<double, 3> lidarCenter_{{0,0,0}};
+  std::array<double, 3> laserAxis_{{1,0,0}};
+  std::array<double, 3> groundNormal_{{0,0,1}};
+};
+
+class LidarLineOfSite
+{
+public:
+  LidarLineOfSite() = default;
+  std::unique_ptr<DataProbeSpecInfo> determine_line_of_site_info(const YAML::Node& node);
+private:
+  void load(const YAML::Node& node);
+
+  SpinnerLidarSegmentGenerator segGen;
+
+  double scanTime_{2};
+  int nsamples_{984};
+  int npoints_{100};
+  std::vector<std::string> fromTargetNames_;
+
+  std::string name_{"lidar_line"};
+
+};
+
+
+
+}
+}
+
+#endif

reg_tests/test_files/ablNeutralEdge/ablNeutralEdge.i

@@ -88,7 +88,7 @@ realms:
           mass_fraction: 1.0
 
       specifications:
- 
+
         # Density here was computed such that P_ref = rho_ref*(R/mw)*300K
         - name: density
           type: constant

reg_tests/test_files/ablUnstableEdge/ablUnstableEdge.i

@@ -184,6 +184,30 @@ realms:
         - source_term_parameters:
             momentum: [0.000135, 0.0, 0.0]
 
+    data_probes:
+
+      output_frequency: 1
+      search_method: stk_kdtree
+      search_tolerance: 1.0e-3
+      search_expansion_factor: 2.0
+
+      lidar_specifications:
+        from_target_part: [Unspecified-2-HEX]
+        inner_prism_initial_theta: 90
+        inner_prism_rotation_rate: 3.5
+        inner_prism_azimuth: 15.2
+        outer_prism_initial_theta: 90
+        outer_prism_rotation_rate: 6.5
+        outer_prism_azimuth: 15.2
+        scan_time: 2 #seconds
+        number_of_samples: 984
+        points_along_line: 100
+        center: [500,500,100]
+        beam_length: 1.0
+        axis: [1,0,0]
+        ground_direction: [0,0,1]
+
+
     output:
       output_data_base_name: abl_1km_cube.e
       output_frequency: 5