ADM Turbine Level Search

include/aero/actuator/ActuatorBulk.h

@@ -65,8 +65,14 @@ 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 +104,7 @@ struct ActuatorBulk
   ActFixElemIds elemContainingPoint_;
 
   const int localTurbineId_;
+  bool singlePointCoarseSearch_ = false;
 };
 
 } // namespace nalu

include/aero/actuator/ActuatorBulkFAST.h

@@ -56,12 +56,23 @@ struct ActuatorBulkFAST : public ActuatorBulk
   bool is_tstep_ratio_admissable(
     const double fastTimeStep, const double naluTimeStep);
 
+  void stk_search_collective_act_pnts(
+    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_;
 
   ActTensorDblDv orientationTensor_;
 
@@ -70,7 +81,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/ActuatorGenericSearchFunctor.h

@@ -68,17 +68,21 @@ struct GenericLoopOverCoarseSearchResults
     innerLoopFunctor_.preloop();
   }
 
+  // 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 +97,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 +108,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 +122,16 @@ 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]);
+      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]);
+        }
+      }
     }
   }
 

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,15 @@ 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

@@ -50,6 +50,7 @@ ActuatorBulkFAST::ActuatorBulkFAST(
     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 +243,69 @@ 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_search_collective_act_pnts(actMeta, stkBulk, onlyFine);
+  }
+}
+
+void
+ActuatorBulkFAST::stk_search_collective_act_pnts(
+  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()()
 {
+  // 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
   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);
-  }
+  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
   }