ADM Turbine Level Search

include/aero/actuator/ActuatorBulk.h

@@ -65,8 +65,12 @@ struct ActuatorBulk
   ActuatorBulk(const ActuatorMeta& actMeta);
   virtual ~ActuatorBulk() {}
 
+  virtual void stk_search(
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine = false);
+
   void stk_search_act_pnts(
     const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk);
+
   void zero_source_terms(stk::mesh::BulkData& stkBulk);
   void parallel_sum_source_term(stk::mesh::BulkData& stkBulk);
   void compute_offsets(const ActuatorMeta& actMeta);
@@ -98,6 +102,7 @@ struct ActuatorBulk
   ActFixElemIds elemContainingPoint_;
 
   const int localTurbineId_;
+  bool singlePointCoarseSearch_ = false; 
 };
 
 } // namespace nalu

include/aero/actuator/ActuatorBulkFAST.h

@@ -55,13 +55,20 @@ struct ActuatorBulkFAST : public ActuatorBulk
   void init_epsilon(const ActuatorMetaFAST& actMeta);
   bool is_tstep_ratio_admissable(
     const double fastTimeStep, const double naluTimeStep);
+
+  void stk_turbine_search (
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk,bool onlyFine = false); 
+
+  void stk_search(
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine = false) override;
 
   virtual ~ActuatorBulkFAST();
 
   ActFixVectorDbl turbineThrust_;
   ActFixVectorDbl turbineTorque_;
   ActFixVectorDbl hubLocations_;
   ActFixVectorDbl hubOrientation_;
+  ActFixScalarDbl turbineSearchRadius_; //need vector for turbine search...this will be a different size than searchRadius_ in Actuatorbulk.h
 
   ActTensorDblDv orientationTensor_;
 
@@ -70,7 +77,7 @@ struct ActuatorBulkFAST : public ActuatorBulk
   ActDualViewHelper<ActuatorMemSpace> dvHelper_;
 };
 
-// helper functions to
+// helper function to
 // squash calls to std::cout from OpenFAST
 inline void
 squash_fast_output(std::function<void()> func)

include/aero/actuator/ActuatorFunctors.h

@@ -10,6 +10,7 @@
 #ifndef ACTUATORFUNCTORS_H_
 #define ACTUATORFUNCTORS_H_
 
+#include <aero/actuator/ActuatorGenericTurbineSearchFunctor.h>
 #include <aero/actuator/ActuatorGenericSearchFunctor.h>
 #include <aero/actuator/ActuatorBulk.h>
 #include <FieldTypeDef.h>

include/aero/actuator/ActuatorFunctorsFAST.h

@@ -38,7 +38,7 @@ inline void
 RunActFastUpdatePoints(ActuatorBulkFAST& actBulk)
 {
   Kokkos::deep_copy(actBulk.pointCentroid_.view_host(), 0.0);
-  actBulk.pointCentroid_.modify_host();
+  actBulk.pointCentroid_.modify_host(); //actuator point locations in space 
   Kokkos::parallel_for(
     "ActFastUpdatePoints", actBulk.local_range_policy(),
     ActFastUpdatePoints(actBulk));

include/aero/actuator/ActuatorGenericSearchFunctor.h

@@ -26,6 +26,7 @@ namespace sierra {
 namespace nalu {
 
 template <typename ActuatorBulk, typename functor>
+  //coarse search actuatorbulk.c L96
 struct GenericLoopOverCoarseSearchResults
 {
   using execution_space = ActuatorFixedExecutionSpace;
@@ -46,6 +47,7 @@ struct GenericLoopOverCoarseSearchResults
     actBulk_.coarseSearchElemIds_.sync_host();
     actBulk_.coarseSearchPointIds_.sync_host();
     innerLoopFunctor_.preloop();
+    //innerLoopExtent_ = actBulk_.singlePointCoarseSearch ? 1 : actBulk_.pointCentroid_.extent(0)
   }
 
   // ctor for functor constructor taking multiple args
@@ -66,19 +68,24 @@ struct GenericLoopOverCoarseSearchResults
     actBulk_.coarseSearchElemIds_.sync_host();
     actBulk_.coarseSearchPointIds_.sync_host();
     innerLoopFunctor_.preloop();
+    //innerLoopExtent_ = actBulk_.singlePointCoarseSearch ? 1 : actBulk_.pointCentroid_.extent(0)
   }
 
+  // see ActuatorExecutorFASTSngp.C line 58
   void operator()(int index) const
   {
+    // properties of elements are controlled by master element
     auto pointId = actBulk_.coarseSearchPointIds_.h_view(index);
     auto elemId = actBulk_.coarseSearchElemIds_.h_view(index);
 
+    // get element topology 
     const stk::mesh::Entity elem =
       stkBulk_.get_entity(stk::topology::ELEMENT_RANK, elemId);
     const stk::topology& elemTopo = stkBulk_.bucket(elem).topology();
     MasterElement* meSCV =
       MasterElementRepo::get_volume_master_element_on_host(elemTopo);
 
+    // element number of nodes and integration points
     const unsigned numNodes = stkBulk_.num_nodes(elem);
     const int numIp = meSCV->num_integration_points();
 
@@ -93,6 +100,7 @@ struct GenericLoopOverCoarseSearchResults
 
     stk::mesh::Entity const* elem_nod_rels = stkBulk_.begin_nodes(elem);
 
+    // get element coordinates
     for (unsigned i = 0; i < numNodes; i++) {
       const double* coords =
         stk::mesh::field_data(*coordinates_, elem_nod_rels[i]);
@@ -103,8 +111,10 @@ struct GenericLoopOverCoarseSearchResults
 
     meSCV->determinant(elemCoords, scvIp);
 
+    // relationship of element nodes to integration points
     const auto* ipNodeMap = meSCV->ipNodeMap();
 
+    // loop over integration points
     for (int nIp = 0; nIp < numIp; nIp++) {
       const int nodeIndex = ipNodeMap[nIp];
       stk::mesh::Entity node = elem_nod_rels[nodeIndex];
@@ -115,7 +125,20 @@ struct GenericLoopOverCoarseSearchResults
       // anything else that is required should be stashed on the functor
       // during functor construction i.e. ActuatorBulk, flags, ActuatorMeta,
       // etc.
-      innerLoopFunctor_(pointId, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
+      //
+      //
+      // for (int actPtInd = 0; actPtInd < innerLoopExtent_; actPtInd ++){
+      //    innerLoopFunctor_(actPtInd, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
+      //  }
+      //
+      if (actBulk_.singlePointCoarseSearch_) { 
+        innerLoopFunctor_(pointId, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
+      } else {
+        for (int actPtInd = 0; actPtInd < actBulk_.pointCentroid_.extent(0); actPtInd ++){
+          innerLoopFunctor_(actPtInd, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
+        }
+      }
+
     }
   }
 
@@ -125,6 +148,7 @@ struct GenericLoopOverCoarseSearchResults
   VectorFieldType* actuatorSource_;
   ScalarFieldType* dualNodalVolume_;
   functor innerLoopFunctor_;
+  //const size_t innerLoopExtent_;
 };
 
 } // namespace nalu

src/aero/actuator/ActuatorBulk.C

@@ -88,8 +88,8 @@ ActuatorBulk::stk_search_act_pnts(
   const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk)
 {
   auto points = pointCentroid_.template view<ActuatorFixedMemSpace>();
-  auto radius = searchRadius_.template view<ActuatorFixedMemSpace>();
 
+  auto radius = searchRadius_.template view<ActuatorFixedMemSpace>();
   auto boundSpheres = CreateBoundingSpheres(points, radius);
   auto elemBoxes = CreateElementBoxes(stkBulk, actMeta.searchTargetNames_);
 
@@ -105,6 +105,13 @@ ActuatorBulk::stk_search_act_pnts(
   actuator_utils::reduce_view_on_host(localParallelRedundancy_);
 }
 
+void ActuatorBulk::stk_search(
+  const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk,bool onlyFine /*= false*/)
+{
+      stk_search_act_pnts(actMeta, stkBulk);
+}
+
+
 void
 ActuatorBulk::zero_source_terms(stk::mesh::BulkData& stkBulk)
 {

src/aero/actuator/ActuatorBulkFAST.C

@@ -43,13 +43,16 @@ ActuatorMetaFAST::is_disk()
   return actuatorType_ == ActuatorType::ActDiskFASTNGP;
 }
 
+
+
 ActuatorBulkFAST::ActuatorBulkFAST(
   const ActuatorMetaFAST& actMeta, double naluTimeStep)
   : ActuatorBulk(actMeta),
     turbineThrust_("turbineThrust", actMeta.numberOfActuators_),
     turbineTorque_("turbineTorque", actMeta.numberOfActuators_),
     hubLocations_("hubLocations", actMeta.numberOfActuators_),
     hubOrientation_("hubOrientations", actMeta.numberOfActuators_),
+    turbineSearchRadius_("hubLocations", actMeta.numberOfActuators_),
     orientationTensor_(
       "orientationTensor",
       actMeta.isotropicGaussian_ ? 0 : actMeta.numPointsTotal_),
@@ -242,6 +245,65 @@ ActuatorBulkFAST::init_epsilon(const ActuatorMetaFAST& actMeta)
   searchRadius_.sync_host();
 }
 
+void
+ActuatorBulkFAST::stk_search(
+  const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine /* = false */)
+{
+    if (singlePointCoarseSearch_){
+      //TODO: Does it make sense for actuator point search to have onlyFine option?
+      stk_search_act_pnts(actMeta, stkBulk);
+    }
+    else{
+      // perform turbine level search and cache to the bulk data
+      stk_turbine_search(actMeta, stkBulk,onlyFine);
+    }
+}
+
+
+void
+ActuatorBulkFAST::stk_turbine_search(
+  const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine /* = false */)
+{
+  auto points = hubLocations_;
+  if (!onlyFine) {
+  // Loop over all turbines to initialize the search radius 
+  for (int iTurb = 0; iTurb < openFast_.get_nTurbinesGlob(); ++iTurb) {
+    // if my process contains the turbine
+    if (NaluEnv::self().parallel_rank() == openFast_.get_procNo(iTurb)) {
+      auto hubLoc = Kokkos::subview(hubLocations_, iTurb, Kokkos::ALL);
+      double turbineRadius = 0.0;
+      // Approximate turbine radius to define search radius
+      //
+      /* const int nbfp = openFast_.get_numForcePtsBlade(iTurb); */
+      /* Point bladeTip = actuator_utils::get_fast_point(openFast_, iTurb, fast::BLADE, nbfp, 0); */
+      /* for (int i = 0; i < 3; ++i) { */
+      /*   turbineRadius += std::pow(bladeTip[i] - hubLoc[i], 2.0); */
+      /* } */
+      //use the hub height as a surrogate for turbineRadius
+      turbineRadius = hubLoc[2];  //TODO: Is z 1 or 2?
+      turbineSearchRadius_(iTurb) = 1.25 * turbineRadius * std::sqrt(2);  //TODO: Could switch to bounding boxes here instead
+
+    }
+  }
+  //actuator_utils::reduce_view_on_host(turbineSearchRadius_.view_host());
+  //turbineSearchRadius_.sync_host();
+  auto radius = turbineSearchRadius_;
+  auto boundSpheres = CreateBoundingSpheres(points,radius);
+  auto elemBoxes = CreateElementBoxes(stkBulk, actMeta.searchTargetNames_);
+
+  // the coarse search now associates element boxes with turbines
+  ExecuteCoarseSearch(
+    boundSpheres, elemBoxes, coarseSearchPointIds_, coarseSearchElemIds_,
+    actMeta.searchMethod_);
+  }
+
+  ExecuteFineSearch(
+    stkBulk, coarseSearchPointIds_, coarseSearchElemIds_, points,
+    elemContainingPoint_, localCoords_, pointIsLocal_,
+    localParallelRedundancy_);
+
+}
+
 Kokkos::RangePolicy<ActuatorFixedExecutionSpace>
 ActuatorBulkFAST::local_range_policy()
 {

src/aero/actuator/ActuatorExecutorsFASTNgp.C

@@ -25,33 +25,44 @@ ActuatorLineFastNGP::ActuatorLineFastNGP(
 void
 ActuatorLineFastNGP::operator()()
 {
-  actBulk_.zero_source_terms(stkBulk_);
+  //Zero the (body-force) actuator source term 
+  actBulk_.zero_source_terms(stkBulk_); 
 
-  // set range policy to only operating over points owned by local fast turbine
+  //set range policy to only operating over points owned by local fast turbine
   auto fastRangePolicy = actBulk_.local_range_policy();
 
+  //Interpolate velocity to actuator points. 
   RunInterpActuatorVel(actBulk_, stkBulk_);
 
+  // Add FLLC correction to velocity field. 
   apply_fllc(actBulk_);
 
   Kokkos::parallel_for(
     "assignFastVelActuatorNgpFAST", fastRangePolicy,
     ActFastAssignVel(actBulk_));
 
+  // get relative velocity from openFAST
   ActFastCacheRelativeVelocities(actBulk_);
 
+  // Compute filtered lifting line correction 
   compute_fllc();
 
+  // Send interpolated velocities at actuator points to openFAST
   actBulk_.interpolate_velocities_to_fast();
 
+  // Get actuator point centroids 
   RunActFastUpdatePoints(actBulk_);
 
-  actBulk_.stk_search_act_pnts(actMeta_, stkBulk_);
+  // Execute fine and coarse search given point centroids 
+  actBulk_.stk_search(actMeta_, stkBulk_); 
 
+  // call openfast and step
   actBulk_.step_fast();
 
+  // compute the force from openfast at actuator points
   RunActFastComputeForce(actBulk_);
 
+  // Loop over all coarse points
   const int localSizeCoarseSearch =
     actBulk_.coarseSearchElemIds_.view_host().extent_int(0);
 
@@ -113,7 +124,7 @@ ActuatorDiskFastNGP::operator()()
 
     actBulk_.update_ADM_points(actMeta_);
 
-    actBulk_.stk_search_act_pnts(actMeta_, stkBulk_);
+    actBulk_.stk_search(actMeta_, stkBulk_,true);
   }
 
   actBulk_.step_fast();

src/aero/actuator/ActuatorFunctors.C

@@ -91,15 +91,21 @@ SpreadForceInnerLoop::operator()(
   actuator_utils::compute_distance(
     3, nodeCoords, pointCoords.data(), &distance[0]);
 
-  const double gauss =
-    actuator_utils::Gaussian_projection(3, &distance[0], epsilon.data());
-
-  for (int j = 0; j < 3; j++) {
-    projectedForce[j] = gauss * pointForce(j);
-  }
+  //Check distance between actuator point and element centroid. Only needed if singlePointCoarseSearch_==False
+  //auto epsilonRadius = 
+  //  Kokkos::subview(actBulk_.searchRadius_.view_host(), pointId, Kokkos::ALL);
+  auto epsilonRadius = actBulk_.searchRadius_.h_view(pointId);
+  if (std::sqrt(distance[0]*distance[0] + distance[1]*distance[1] + distance[2]*distance[2]) < epsilonRadius) {
+    const double gauss =
+      actuator_utils::Gaussian_projection(3, &distance[0], epsilon.data());
+
+    for (int j = 0; j < 3; j++) {
+      projectedForce[j] = gauss * pointForce(j);
+    }
 
-  for (int j = 0; j < 3; j++) {
-    sourceTerm[j] += projectedForce[j] * scvIp / dual_vol;
+    for (int j = 0; j < 3; j++) {
+      sourceTerm[j] += projectedForce[j] * scvIp / dual_vol;
+    }
   }
 }
 

src/aero/actuator/ActuatorModel.C

@@ -145,8 +145,8 @@ ActuatorModel::init(stk::mesh::BulkData& stkBulk)
     break;
 #endif
 #else
-    // perform search for actline and actdisk
-    actBulk_->stk_search_act_pnts(*actMeta_.get(), stkBulk);
+      //perform stk_search (coarse + fine search)
+      actBulk_->stk_search(*actMeta_.get(), stkBulk);
     break;
 #endif
   }