ADM Turbine Level Search

include/aero/actuator/ActuatorBulk.h

@@ -48,6 +48,7 @@ struct ActuatorMeta
   stk::search::SearchMethod searchMethod_;
   ActScalarIntDv numPointsTurbine_;
   bool useFLLC_ = false;
+  bool turbineLevelSearch_ = false; //Adding turbine level search option to meta. 
   ActVectorDblDv epsilonChord_;
   ActVectorDblDv epsilon_;
   ActFixScalarBool entityFLLC_;
@@ -67,6 +68,16 @@ struct ActuatorBulk
 
   void stk_search_act_pnts(
     const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk);
+
+  //Generic stk search
+  void stk_search(
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine = false);
+
+  //do I need this to be virtual function so the generic stk_search can call it, even though 
+  //the turbine search only makes sense in the ActuatorBulkFast context as implemented? 
+  virtual void stk_turbine_search(
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine = false); 
+
   void zero_source_terms(stk::mesh::BulkData& stkBulk);
   void parallel_sum_source_term(stk::mesh::BulkData& stkBulk);
   void compute_offsets(const ActuatorMeta& actMeta);

include/aero/actuator/ActuatorBulkFAST.h

@@ -55,13 +55,18 @@ struct ActuatorBulkFAST : public ActuatorBulk
   void init_epsilon(const ActuatorMetaFAST& actMeta);
   bool is_tstep_ratio_admissable(
     const double fastTimeStep, const double naluTimeStep);
+
+  // This is placed in ActuatorBulkFAST instead of ActuatorBulk because hublocations are needed
+  void stk_turbine_search (
+    const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk,bool onlyFine = false) override; 
 
   virtual ~ActuatorBulkFAST();
 
   ActFixVectorDbl turbineThrust_;
   ActFixVectorDbl turbineTorque_;
   ActFixVectorDbl hubLocations_;
   ActFixVectorDbl hubOrientation_;
+  ActFixVectorDbl turbineSearchRadius_; //need vector for turbine search...this will be a different size than searchRadius_ in Actuatorbulk.h
 
   ActTensorDblDv orientationTensor_;
 
@@ -70,7 +75,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/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;
@@ -68,17 +69,25 @@ struct GenericLoopOverCoarseSearchResults
     innerLoopFunctor_.preloop();
   }
 
+  // see ActuatorExecutorFASTSngp.C line 58
   void operator()(int index) const
   {
+    // TODO (GOPAL): Index gives us actuator id and element id from coarse search 
+    // coarse search associates actuator points with elements 
+    //
+    // element associates a group of points
+    // 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 +102,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 +113,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,6 +127,8 @@ struct GenericLoopOverCoarseSearchResults
       // anything else that is required should be stashed on the functor
       // during functor construction i.e. ActuatorBulk, flags, ActuatorMeta,
       // etc.
+      //
+      // pointID helps look up data from openfast
       innerLoopFunctor_(pointId, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
     }
   }

include/aero/actuator/ActuatorGenericTurbineSearchFunctor.h

@@ -0,0 +1,158 @@
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+// Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC
+// (NTESS), National Renewable Energy Laboratory, University of Texas Austin,
+// Northwest Research Associates. Under the terms of Contract DE-NA0003525
+// with NTESS, the U.S. Government retains certain rights in this software.
+//
+// This software is released under the BSD 3-clause license. See LICENSE file
+// for more details.
+//
+#ifndef ACTUATORGENERICSEARCHFUNCTOR_H_
+#define ACTUATORGENERICSEARCHFUNCTOR_H_
+
+#include <aero/actuator/UtilitiesActuator.h>
+#include <aero/actuator/ActuatorTypes.h>
+#include <stk_mesh/base/BulkData.hpp>
+#include <FieldTypeDef.h>
+
+namespace sierra {
+namespace nalu {
+
+
+//TODO: Should this be inherited somehow from GenericLoopOverCoarseSearchResults functor?
+//TODO: the actuator meta is needed below. Do we have access to it?
+template <typename ActuatorBulk, typename functor>
+  //coarse search actuatorbulk.c L96
+struct GenericLoopOverCoarseTurbineSearchResults
+{
+  using execution_space = ActuatorFixedExecutionSpace;
+
+  // ctor if functor only requires actBulk for constructor
+  GenericLoopOverCoarseTurbineSearchResults(
+    ActuatorBulk& actBulk, stk::mesh::BulkData& stkBulk)
+    : actBulk_(actBulk),
+      stkBulk_(stkBulk),
+      coordinates_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "coordinates")),
+      actuatorSource_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "actuator_source")),
+      dualNodalVolume_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "dual_nodal_volume")),
+      innerLoopFunctor_(actBulk)
+  {
+    actBulk_.coarseSearchElemIds_.sync_host();
+    actBulk_.coarseSearchPointIds_.sync_host();
+    innerLoopFunctor_.preloop();
+  }
+
+  // ctor for functor constructor taking multiple args
+  GenericLoopOverCoarseTurbineSearchResults(
+    ActuatorBulk& actBulk,
+    stk::mesh::BulkData& stkBulk,
+    functor innerLoopFunctor)
+    : actBulk_(actBulk),
+      stkBulk_(stkBulk),
+      coordinates_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "coordinates")),
+      actuatorSource_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "actuator_source")),
+      dualNodalVolume_(stkBulk_.mesh_meta_data().template get_field<double>(
+        stk::topology::NODE_RANK, "dual_nodal_volume")),
+      innerLoopFunctor_(innerLoopFunctor)
+  {
+    actBulk_.coarseSearchElemIds_.sync_host();
+    actBulk_.coarseSearchPointIds_.sync_host();
+    innerLoopFunctor_.preloop();
+  }
+
+  // see ActuatorExecutorFASTSngp.C line 58
+  // TODO: should the index be loop over actuator points or loop over turbine points. I think it would be
+  // more efficient to loop over all turbine points and then over each actuator point so index should be 
+  // turbine points. 
+  //
+  // might be able to replace operator with individual cases
+  void operator()(int index) const
+  {
+    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();
+
+    // just allocate for largest expected size (hex27)
+    STK_ThrowAssert(numIp <= 216);
+    STK_ThrowAssert(numNodes <= 27);
+
+    double scvip[216];
+    double elemcoords[27 * 3];
+    sierra::nalu::SharedMemView<double*> scvIp(&scvip[0], 216);
+    sierra::nalu::SharedMemView<double**> elemCoords(&elemcoords[0], 27, 3);
+
+    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]);
+      for (int j = 0; j < 3; j++) {
+        elemCoords(i, j) = coords[j];
+      }
+    }
+
+    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];
+      const double* nodeCoords = stk::mesh::field_data(*coordinates_, node);
+      const double dual_vol = *stk::mesh::field_data(*dualNodalVolume_, node);
+      double* sourceTerm = stk::mesh::field_data(*actuatorSource_, node);
+
+      //TODO: Need to get access to actMeta for this to work
+      //This loop should also be parallelize I would think.. 
+      //Option 1:
+      //use kokkoss size function to get numPointsTotal from the size of pointCentroids. Cannot use meta here
+      for (int actPtInd = 0; actPtInd < actBulk_.pointCentroids_.size(); actPtInd ++){
+        //auto pointCoords = Kokkos::subview(actBulk_.pointCentroid_.view_host(), actPtInd, Kokkos::ALL);
+      //Option 2:
+      //Kokkos::parallel_for("GenericActPtLoop", HostRangePolicy(0, actMeta.numPointsTotal_), [&](int actPtInd) {
+      //auto point = pointCentroid_.h_view(actPtInd)  don't need point just index
+      //TODO: Is actPtInd the right index?
+      //TODO: How can we avoid all the zero integration values if actPtId is not associated with element box? 
+      //      Or is this OK because we spread the force to the entire disk
+      innerLoopFunctor_(actPtInd, nodeCoords, sourceTerm, dual_vol, scvIp[nIp]);
+      //
+    }
+  }
+
+  ActuatorBulk& actBulk_;
+  stk::mesh::BulkData& stkBulk_;
+  VectorFieldType* coordinates_;
+  VectorFieldType* actuatorSource_;
+  ScalarFieldType* dualNodalVolume_;
+  functor innerLoopFunctor_;
+};
+
+} // namespace nalu
+} // namespace sierra
+
+#endif /* ACTUATORGENERICSEARCHFUNCTOR_H_ */

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,27 @@ ActuatorBulk::stk_search_act_pnts(
   actuator_utils::reduce_view_on_host(localParallelRedundancy_);
 }
 
+void stk_turbine_search(
+  const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk, bool onlyFine /*= false*/)
+{
+    STK_ThrowErrorMsg("Turbine Search Requires ActuatorBulkFAST data");
+}
+
+void
+ActuatorBulk::stk_search(
+  const ActuatorMeta& actMeta, stk::mesh::BulkData& stkBulk,bool onlyFine /*= false*/)
+{
+    if (actMeta.turbineLevelSearch_){
+      // perform turbine level search and cache to the bulk data
+      stk_turbine_search(actMeta.get(), stkBulk,onlyFine);
+    }
+    else{
+      //TODO: Does it make sense for actuator point search to have onlyFine option?
+      stk_search_act_pnts(actMeta.get(), stkBulk);
+    }
+}
+
+
 void
 ActuatorBulk::zero_source_terms(stk::mesh::BulkData& stkBulk)
 {