diff --git a/include/aero/actuator/ActuatorBulk.h b/include/aero/actuator/ActuatorBulk.h index af1d146086..1d812cccbf 100644 --- a/include/aero/actuator/ActuatorBulk.h +++ b/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 diff --git a/include/aero/actuator/ActuatorBulkFAST.h b/include/aero/actuator/ActuatorBulkFAST.h index 28301bc433..64c37b5afb 100644 --- a/include/aero/actuator/ActuatorBulkFAST.h +++ b/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 dvHelper_; }; -// helper functions to +// helper function to // squash calls to std::cout from OpenFAST inline void squash_fast_output(std::function func) diff --git a/include/aero/actuator/ActuatorGenericSearchFunctor.h b/include/aero/actuator/ActuatorGenericSearchFunctor.h index 408a3fa5a6..8af9793add 100644 --- a/include/aero/actuator/ActuatorGenericSearchFunctor.h +++ b/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]); + } + } } } diff --git a/src/aero/actuator/ActuatorBulk.C b/src/aero/actuator/ActuatorBulk.C index 69b4c80283..1babfd5944 100644 --- a/src/aero/actuator/ActuatorBulk.C +++ b/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(); - auto radius = searchRadius_.template view(); + auto radius = searchRadius_.template view(); 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) { diff --git a/src/aero/actuator/ActuatorBulkFAST.C b/src/aero/actuator/ActuatorBulkFAST.C index 02c2915db7..aeb03f91fd 100644 --- a/src/aero/actuator/ActuatorBulkFAST.C +++ b/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 ActuatorBulkFAST::local_range_policy() { diff --git a/src/aero/actuator/ActuatorExecutorsFASTNgp.C b/src/aero/actuator/ActuatorExecutorsFASTNgp.C index ef6b4dc917..f11f417379 100644 --- a/src/aero/actuator/ActuatorExecutorsFASTNgp.C +++ b/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(); diff --git a/src/aero/actuator/ActuatorFunctors.C b/src/aero/actuator/ActuatorFunctors.C index ca55f00032..e5a31463cf 100644 --- a/src/aero/actuator/ActuatorFunctors.C +++ b/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; + } } } diff --git a/src/aero/actuator/ActuatorModel.C b/src/aero/actuator/ActuatorModel.C index 056b26d1ea..1ec9267390 100644 --- a/src/aero/actuator/ActuatorModel.C +++ b/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 }