diff --git a/amr-wind/core/gpu_utils.H b/amr-wind/core/gpu_utils.H new file mode 100644 index 0000000000..af5799316a --- /dev/null +++ b/amr-wind/core/gpu_utils.H @@ -0,0 +1,21 @@ +#ifndef GPU_UTILS_H +#define GPU_UTILS_H + +#include "AMReX_Gpu.H" + +namespace amr_wind { +namespace gpu { + +template +inline amrex::Gpu::DeviceVector device_view(const amrex::Vector& hview) +{ + amrex::Gpu::DeviceVector dview(hview.size()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, hview.begin(), hview.end(), dview.begin()); + return dview; +} + +} // namespace gpu +} // namespace amr_wind + +#endif /* GPU_UTILS_H */ diff --git a/amr-wind/wind_energy/CMakeLists.txt b/amr-wind/wind_energy/CMakeLists.txt index 6f347ee2ee..8b80176ec4 100644 --- a/amr-wind/wind_energy/CMakeLists.txt +++ b/amr-wind/wind_energy/CMakeLists.txt @@ -9,3 +9,5 @@ target_sources(${amr_wind_lib_name} ABLBoundaryPlane.cpp MOData.cpp ) + +add_subdirectory(actuator) diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.H b/amr-wind/wind_energy/actuator/ActuatorContainer.H new file mode 100644 index 0000000000..eee590992f --- /dev/null +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.H @@ -0,0 +1,103 @@ +#ifndef ACTUATORCONTAINER_H +#define ACTUATORCONTAINER_H + +#include "amr-wind/core/vs/vector_space.H" + +#include "AMReX_AmrParticles.H" + +namespace amr_wind { + +class Field; + +namespace actuator { + +/** Cloud of actuator points + * + * Holds the position vector and velocity information for a cloud of actuator + * points belonging to all the turbines within a given MPI rank. + */ +struct ActuatorCloud +{ + //! Number of nodes per turbine on this MPI rank + amrex::Vector num_pts; + + //! Position vectors of all actuator nodes belonging to this MPI rank + amrex::Vector position; + + //! Velocity vectors of all actuator nodes belonging to this MPI rank + amrex::Vector velocity; + + //! Total number of turbines located on this MPI rank + int num_objects; + + explicit ActuatorCloud(const int nobjects); +}; + +//! Number or real entries in Array of Structs (AOS) +static constexpr int NumPStructReal = AMREX_SPACEDIM; +//! Number of integer entries in Array of Structs (AOS) +static constexpr int NumPStructInt = 0; +//! Number of real entries in Struct of Arrays (SOA) +static constexpr int NumPArrayReal = 0; +//! Number of int entries in Struct of Arrays (SOA) +static constexpr int NumPArrayInt = 0; + +/** Specialization of AmrParticleContainer for sampling velocities + */ +class ActuatorContainer + : public amrex::AmrParticleContainer< + NumPStructReal, + NumPStructInt, + NumPArrayReal, + NumPArrayInt> +{ +public: + explicit ActuatorContainer(amrex::AmrCore& mesh, const int num_objects); + + void post_regrid_actions(); + + void initialize_container(); + + void update_positions(); + + void sample_velocities(const Field& vel); + + int num_actuator_points() const { return m_data.position.size(); } + + // public for CUDA, not safe for general access + + void interpolate_velocities(const Field& vel); + + void populate_vel_buffer(); + +protected: + void compute_local_coordinates(); + + // Accessor to allow unit testing + ActuatorCloud& point_data() { return m_data; } + +private: + amrex::AmrCore& m_mesh; + + // Object that holds the position and velocity information + ActuatorCloud m_data; + + //! Coordinates of points that are contained within boxes assigned to each + //! MPI rank + amrex::Vector m_proc_pos; + + //! Device view of the process position vectors + amrex::Gpu::DeviceVector m_pos_device; + + //! Flag indicating whether memory has allocated for all data structures + bool m_container_initialized{false}; + + //! Flag indicating whether the particles are scattered throughout the + //! domain, or if they have been recalled to the original MPI rank + bool m_is_scattered{false}; +}; + +} // namespace actuator +} // namespace amr_wind + +#endif /* ACTUATORCONTAINER_H */ diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp new file mode 100644 index 0000000000..4687bdc223 --- /dev/null +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp @@ -0,0 +1,323 @@ +#include "amr-wind/wind_energy/actuator/ActuatorContainer.H" +#include "amr-wind/core/gpu_utils.H" +#include "amr-wind/core/Field.H" + +#include "AMReX_Scan.H" + +#include + +namespace amr_wind { +namespace actuator { + +ActuatorCloud::ActuatorCloud(const int nobjects) + : num_pts(nobjects, 0), num_objects(nobjects) +{} + +ActuatorContainer::ActuatorContainer( + amrex::AmrCore& mesh, const int num_objects) + : amrex::AmrParticleContainer< + NumPStructReal, + NumPStructInt, + NumPArrayReal, + NumPArrayInt>(&mesh) + , m_mesh(mesh) + , m_data(num_objects) + , m_proc_pos(amrex::ParallelDescriptor::NProcs(), vs::Vector::zero()) + , m_pos_device(amrex::ParallelDescriptor::NProcs(), vs::Vector::zero()) +{} + +/** Allocate memory and initialize the particles within the container + * + * This method is only called once during the simulation. It allocates the + * arrays for holding the position vector and velocity data on host memory and + * also initializes corresponding particles within the first available particle + * tile within the container. It is expected that the actuator manager instance + * has already populated the number of points per turbine before invoking this + * method. + */ +void ActuatorContainer::initialize_container() +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::initialize_container"); + + compute_local_coordinates(); + + // Initialize global data arrays + auto total_pts = + std::accumulate(m_data.num_pts.begin(), m_data.num_pts.end(), 0); + m_data.position.resize(total_pts); + m_data.velocity.resize(total_pts); + + // Initialize particle container data structures. + // + // We assign all particles into the first available container within this + // MPI rank and let redistribute take care of scattering the particles into + // the respective MPI rank containing the cell enclosing this particle. The + // actual redistribution happens after position vectors are updated in + // update_positions. + + // query the particle ID from the container. We should always be starting + // from 1. + const auto id_start = ParticleType::NextID(); + AMREX_ALWAYS_ASSERT(id_start == 1u); + const int iproc = amrex::ParallelDescriptor::MyProc(); + + // Flag indicating if a tile was found where all particles were deposited. + bool assigned = false; + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; ((lev < nlevels) && !assigned); ++lev) { + for (auto mfi = MakeMFIter(lev); (mfi.isValid() && !assigned); ++mfi) { + auto& ptile = GetParticles( + lev)[std::make_pair(mfi.index(), mfi.LocalTileIndex())]; + ptile.resize(total_pts); + auto* pstruct = ptile.GetArrayOfStructs()().data(); + + amrex::ParallelFor( + total_pts, [=] AMREX_GPU_DEVICE(const int ip) noexcept { + auto& pp = pstruct[ip]; + + pp.id() = id_start + ip; + pp.cpu() = iproc; + }); + assigned = true; + } + } + + // Indicate that we have initialized the containers and remaining methods + // are safe to use + m_container_initialized = true; + m_is_scattered = false; +} + +/** Update position vectors of the particles within a container based on the + * data provided by actuator instances. + * + * This method assumes that the particles have been restored to their + * originating MPI ranks after any scattering. This happens during + * initialization (after AcutatorContainer::initialize_particles) and after + * ActuatorContainer::sample_velocities. + * + * After executing this method, the particles have been scattered across the + * domain such that each particle is contained within a tile that has the cell + * enclosing the particle. + * + */ +void ActuatorContainer::update_positions() +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::update_positions"); + AMREX_ALWAYS_ASSERT(m_container_initialized && !m_is_scattered); + + const auto dpos = gpu::device_view(m_data.position); + const auto dptr = dpos.data(); + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; lev < nlevels; ++lev) { + for (ParIterType pti(*this, lev); pti.isValid(); ++pti) { + const int np = pti.numParticles(); + auto* pstruct = pti.GetArrayOfStructs()().data(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(const int ip) noexcept { + auto& pp = pstruct[ip]; + const auto idx = pp.id() - 1; + + auto& pvec = dptr[idx]; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + pp.pos(n) = pvec[n]; + } + }); + } + } + + // Scatter particles to appropriate MPI ranks + Redistribute(); + + // Indicate that it is safe to sample velocities + m_is_scattered = true; +} + +/** Interpolate the velocity field using a trilinear interpolation + * + * This method performs three tasks: + * - Sample the velocity field and interpolate onto particle location + * - Restore the particles back to their original MPI rank + * - Copy data from particles into the buffer used by actuator instances + * + * After invocation if this method, the particles are back in their original + * rank and the position vectors can be updated safely. + */ +void ActuatorContainer::sample_velocities(const Field& vel) +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::sample_velocities"); + AMREX_ALWAYS_ASSERT(m_container_initialized && m_is_scattered); + + // Sample velocity field + interpolate_velocities(vel); + + // Recall particles to the MPI ranks that contains their corresponding + // turbines + Redistribute(); + + // Populate the velocity buffer that all actuator instances can access + populate_vel_buffer(); + + // Indicate that the particles have been restored to their original MPI rank + m_is_scattered = false; +} + +/** Helper method for ActuatorContainer::sample_velocities + * + * Loops over the particle tiles and copies velocity data from particles to the + * velocity array. + */ +void ActuatorContainer::populate_vel_buffer() +{ + amrex::Gpu::DeviceVector vel(m_data.velocity.size()); + auto* varr = vel.data(); + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; lev < nlevels; ++lev) { + for (ParIterType pti(*this, lev); pti.isValid(); ++pti) { + const int np = pti.numParticles(); + auto* pstruct = pti.GetArrayOfStructs()().data(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(const int ip) noexcept { + auto& pp = pstruct[ip]; + const auto idx = pp.id() - 1; + + auto& vvel = varr[idx]; + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + vvel[n] = pp.rdata(n); + } + }); + } + } + + amrex::Gpu::copy( + amrex::Gpu::deviceToHost, vel.begin(), vel.end(), + m_data.velocity.begin()); +} + +/** Helper method for ActuatorContainer::sample_velocities + * + * Performs a trilinear interpolation of the velocity field to particle + * locations. It also updates the particle locations such that the next + * Redistribute call restores the particles back to their original MPI rank + * where they were created. + */ +void ActuatorContainer::interpolate_velocities(const Field& vel) +{ + auto* dptr = m_pos_device.data(); + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; lev < nlevels; ++lev) { + const auto& geom = m_mesh.Geom(lev); + const auto dx = geom.CellSizeArray(); + const auto dxi = geom.InvCellSizeArray(); + const auto plo = geom.ProbLoArray(); + + for (ParIterType pti(*this, lev); pti.isValid(); ++pti) { + const int np = pti.numParticles(); + auto* pstruct = pti.GetArrayOfStructs()().data(); + const auto varr = vel(lev).const_array(pti); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(const int ip) noexcept { + auto& pp = pstruct[ip]; + // Determine offsets within the containing cell + const amrex::Real x = + (pp.pos(0) - plo[0] - 0.5 * dx[0]) * dxi[0]; + const amrex::Real y = + (pp.pos(1) - plo[1] - 0.5 * dx[1]) * dxi[1]; + const amrex::Real z = + (pp.pos(2) - plo[2] - 0.5 * dx[2]) * dxi[2]; + + // Index of the low corner + const int i = static_cast(amrex::Math::floor(x)); + const int j = static_cast(amrex::Math::floor(y)); + const int k = static_cast(amrex::Math::floor(z)); + + // Interpolation weights in each direction (linear basis) + const amrex::Real wx_hi = (x - i); + const amrex::Real wy_hi = (y - j); + const amrex::Real wz_hi = (z - k); + + const amrex::Real wx_lo = 1.0 - wx_hi; + const amrex::Real wy_lo = 1.0 - wy_hi; + const amrex::Real wz_lo = 1.0 - wz_hi; + + const int iproc = pp.cpu(); + + for (int ic = 0; ic < AMREX_SPACEDIM; ++ic) { + pp.rdata(ic) = + wx_lo * wy_lo * wz_lo * varr(i, j, k, ic) + + wx_lo * wy_lo * wz_hi * varr(i, j, k + 1, ic) + + wx_lo * wy_hi * wz_lo * varr(i, j + 1, k, ic) + + wx_lo * wy_hi * wz_hi * varr(i, j + 1, k + 1, ic) + + wx_hi * wy_lo * wz_lo * varr(i + 1, j, k, ic) + + wx_hi * wy_lo * wz_hi * varr(i + 1, j, k + 1, ic) + + wx_hi * wy_hi * wz_lo * varr(i + 1, j + 1, k, ic) + + wx_hi * wy_hi * wz_hi * varr(i + 1, j + 1, k + 1, ic); + + // Reset position vectors so that the particles return back + // to the MPI ranks with the turbines upon redistribution + pp.pos(ic) = dptr[iproc][ic]; + } + }); + } + } +} + +/** Determine position vector of a point within each MPI rank + * + * Loops over the boxArray and determines the first patch that belongs to the + * current MPI rank. Uses that to generate a position vector that is known to + * exist in a given rank. Setting the position vectors of the particles to this + * know location will recall particles belonging to this rank during + * Redistribute. + */ +void ActuatorContainer::compute_local_coordinates() +{ + BL_PROFILE( + "amr-wind::actuator::ActuatorContainer::compute_local_coordinates"); + const int nprocs = amrex::ParallelDescriptor::NProcs(); + const int iproc = amrex::ParallelDescriptor::MyProc(); + + // Reset position vectors to zero (required for parallel reduce sum) + m_proc_pos.assign(nprocs, vs::Vector::zero()); + + // Flag indicating whether a point within the domain belonging to this rank + // has been found + bool assigned = false; + const int nlevels = m_mesh.finestLevel() + 1; + for (int lev = 0; ((lev < nlevels) && !assigned); ++lev) { + const auto& ba = m_mesh.boxArray(lev); + const auto& dm = m_mesh.DistributionMap(lev); + + const int nbx = dm.size(); + for (int i = 0; (i < nbx) && !assigned; ++i) { + if (dm[i] != iproc) continue; + + const auto& geom = m_mesh.Geom(lev); + const auto& dx = geom.CellSize(); + const auto& problo = geom.ProbLo(); + const auto& bx = ba[i]; + const int* lo = bx.loVect(); + + auto& pvec = m_proc_pos[iproc]; + pvec.x() = problo[0] + (lo[0] + 0.5) * dx[0]; + pvec.y() = problo[1] + (lo[1] + 0.5) * dx[1]; + pvec.z() = problo[2] + (lo[2] + 0.5) * dx[2]; + + // Indicate that we have found a point and it is safe to exit the + // loop + assigned = true; + } + } + + // Share position vectors with every process + amrex::ParallelDescriptor::ReduceRealSum( + &(m_proc_pos[0].x()), m_proc_pos.size() * vs::Vector::ncomp); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_proc_pos.begin(), m_proc_pos.end(), + m_pos_device.begin()); +} + +void ActuatorContainer::post_regrid_actions() { compute_local_coordinates(); } + +} // namespace actuator +} // namespace amr_wind diff --git a/amr-wind/wind_energy/actuator/ActuatorModel.H b/amr-wind/wind_energy/actuator/ActuatorModel.H new file mode 100644 index 0000000000..829941aa7f --- /dev/null +++ b/amr-wind/wind_energy/actuator/ActuatorModel.H @@ -0,0 +1,96 @@ +#ifndef ACTUATORMODEL_H +#define ACTUATORMODEL_H + +#include "amr-wind/core/Factory.H" +#include "amr-wind/wind_energy/actuator/actuator_types.H" +#include "amr-wind/wind_energy/actuator/actuator_ops.H" + +namespace amr_wind { + +class CFDSim; + +namespace actuator { + +/** Abstract representation of an actuator source in the flowfield + * + * This class is intended to represent a single object (e.g., turbine, fixed + * wing) that exists in the computational domain. A collection of source + * objects (e.g., turbines in a wind farm) is modeled using the + * amr_wind::Actuator physics class that holds a collection of instances of + * ActuatorModel and acts on them. + */ +class ActuatorModel + : public Factory +{ +public: + static std::string base_identifier() { return "ActuatorModel"; } + + virtual ~ActuatorModel() = default; + + virtual const std::string& label() const = 0; + + virtual int id() const = 0; + + virtual const ActInfo& info() const = 0; + + virtual void pre_init_actions() = 0; + + virtual void post_init_actions() = 0; +}; + +/** Concrete implementation of the ActuatorModel for different actuator types + * + * \tparam ActTrait An actuator type trait that defines the behavior of the + * different actions (initialization, velocity sampling, force computation, and + * momentum source term computations) that are executed by the ActuatorModel + * class through the amr_wind::Actuator physics driver class. + * + * \tparam SrcTrait A source type trait that defines how the forcing term is + * computed, e.g., if the turbine rotor is represented as a line or a disk. For + * most cases, like the fixed wing case, this will just be line, therefore, it + * is set as default for convenience. + */ +template +class ActModel : public ActuatorModel::Register> +{ +private: + //! Instance that holds all data related to a particular actuator type + ActDataHolder m_data; + +public: + static std::string identifier() + { + return ActTrait::identifier() + SrcTrait::identifier(); + } + + /** + * \param sim CFD simulation environment instance + * \param label Unique string identifier for this instance + * \param id Unique integer tag for this instance + */ + ActModel(CFDSim& sim, const std::string& label, const int id) + : m_data(sim, label, id) + {} + + //! Return the unique identifier (name) for this instance + const std::string& label() const override { return m_data.m_info.label; } + + //! Return the unique tag (integer ID) for this instance + int id() const override { return m_data.m_info.id; } + + //! Return the core info object for this actuator instance + const ActInfo& info() const override { return m_data.m_info; } + + void pre_init_actions() override { ops::read_inputs(m_data); } + + void post_init_actions() override + { + ops::init_data_structures(m_data); + ops::determine_influenced_procs(m_data); + } +}; + +} // namespace actuator +} // namespace amr_wind + +#endif /* ACTUATORMODEL_H */ diff --git a/amr-wind/wind_energy/actuator/CMakeLists.txt b/amr-wind/wind_energy/actuator/CMakeLists.txt new file mode 100644 index 0000000000..15678a6dd2 --- /dev/null +++ b/amr-wind/wind_energy/actuator/CMakeLists.txt @@ -0,0 +1,6 @@ +target_sources(${amr_wind_lib_name} + PRIVATE + + actuator_utils.cpp + ActuatorContainer.cpp + ) diff --git a/amr-wind/wind_energy/actuator/actuator_ops.H b/amr-wind/wind_energy/actuator/actuator_ops.H new file mode 100644 index 0000000000..d56cd7e145 --- /dev/null +++ b/amr-wind/wind_energy/actuator/actuator_ops.H @@ -0,0 +1,49 @@ +#ifndef ACT_OPS_H +#define ACT_OPS_H + +namespace amr_wind { +namespace actuator { +namespace ops { + +/** Read user inputs from an input file for a given actuator object + * + * This function is called within Physics::pre_init_actions() and, therefore, + * does not have access to mesh data. Actions that require mesh information for + * initialization should be deferred until ops::init_data_structures call. + * + * \tparam T An actuator traits type + * \param Data object for the specific actuator instance + */ +template +void read_inputs(typename T::DataType&); + +/** Perform one-time initialization of the actuator object + * + * This function is responsible for allocating memory ensuring that the + * actuator instance has been initialized and ready to be used by the Actuator + * physics instance. This function is called once during + * Actuator::post_init_actions + * + * \tparam T An actuator traits type + * \param Data object for the specific actuator instance + */ +template +void init_data_structures(typename T::DataType&); + +/** Determine the list of processes (MPI ranks) that is expected to be + * influenced by this actuator source. + * + * The list of processes is used to determine where velocity sampling data is + * too be gathered, as well as ensuring that all actuator data is properly + * mapped. + */ +template +void determine_influenced_procs(typename T::DataType&); + +} // namespace ops +} // namespace actuator +} // namespace amr_wind + +#include "amr-wind/wind_energy/actuator/actuator_opsI.H" + +#endif /* ACT_OPS_H */ diff --git a/amr-wind/wind_energy/actuator/actuator_opsI.H b/amr-wind/wind_energy/actuator/actuator_opsI.H new file mode 100644 index 0000000000..0b4830f36f --- /dev/null +++ b/amr-wind/wind_energy/actuator/actuator_opsI.H @@ -0,0 +1,31 @@ +#ifndef ACTUATOR_OPSI_H +#define ACTUATOR_OPSI_H + +#include "amr-wind/wind_energy/actuator/actuator_utils.H" +#include "amr-wind/CFDSim.H" + +#include + +namespace amr_wind { +namespace actuator { +namespace ops { + +template +void determine_influenced_procs(typename T::DataType& data) +{ + auto& info = data.m_info; + info.procs = + utils::determine_influenced_procs(data.m_sim.mesh(), info.bound_box); + + info.root_proc = *std::min_element(info.procs.begin(), info.procs.end()); + + const int iproc = amrex::ParallelDescriptor::MyProc(); + auto in_proc = info.procs.find(iproc); + info.actuator_in_proc = (in_proc != info.procs.end()); +} + +} // namespace ops +} // namespace actuator +} // namespace amr_wind + +#endif /* ACTUATOR_OPSI_H */ diff --git a/amr-wind/wind_energy/actuator/actuator_types.H b/amr-wind/wind_energy/actuator/actuator_types.H new file mode 100644 index 0000000000..ed1a00753a --- /dev/null +++ b/amr-wind/wind_energy/actuator/actuator_types.H @@ -0,0 +1,161 @@ +#ifndef ACTUATOR_TYPES_H +#define ACTUATOR_TYPES_H + +#include "amr-wind/core/vs/vector_space.H" + +#include "AMReX_Gpu.H" +#include "AMReX_RealBox.H" + +#include + +namespace amr_wind { + +class CFDSim; + +namespace actuator { + +/** Abstract representation of an actuator type + */ +struct ActuatorType +{}; + +/** Abstract representation of a source type + */ +struct ActSrcType +{}; + +/** Actuator line representation + */ +struct ActSrcLine : ActSrcType +{ + static std::string identifier() { return "Line"; } + + //! Flag indicating if this is a line type + static constexpr bool is_line = true; + + //! Flag indicating if this is a disk type + static constexpr bool is_disk = false; +}; + +/** Acutator disk representation + */ +struct ActSrcDisk : ActSrcType +{ + static std::string identifier() { return "Disk"; } + + static constexpr bool is_line = true; + static constexpr bool is_disk = false; +}; + +using VecList = amrex::Vector; +using TensorList = amrex::Vector; +using DeviceVecList = amrex::Gpu::DeviceVector; +using DeviceTensorList = amrex::Gpu::DeviceVector; + +/** Actuator data at each node for an actuator component + */ +struct ActGrid +{ + //! Position vectors of the actuator forcing points on the grid + VecList pos; + + //! Force vector at the forcing points + VecList force; + + //! Gaussian smearing parameter + VecList epsilon; + + //! Transformation matrix at the forcing points + TensorList orientation; + + //! Position vectors for the points where velocity is sampled + VecList vel_pos; + + //! Velocity vector at the sampled locations + VecList vel; + + /** Helper method to resize the data arrays defined on the grid + * + * \params Number of nodes that contain forcing data + * \params Number of nodes where velocity field is sampled + */ + void resize(int num_force_pts, int num_vel_pts) + { + pos.resize(num_force_pts); + force.resize(num_force_pts); + epsilon.resize(num_force_pts); + orientation.resize(num_force_pts); + vel_pos.resize(num_vel_pts); + vel.resize(num_vel_pts); + } + + /** Convenience function to resize both force/velocity data to same + * locations + */ + inline void resize(int num_pts) { resize(num_pts, num_pts); } +}; + +/** Basic information that is common to all actuator types + */ +struct ActInfo +{ + //! Unique string identifier for the actuator body (wing, turbine) + std::string label; + + //! Set of MPI ranks where this actuator source is active + std::set procs; + + //! Bounding box used to determine processes where this turbine has + //! influence + amrex::RealBox bound_box; + + //! Unique integer identifier for the turbine + int id{-1}; + + //! Root process where this turbine is active + int root_proc{-1}; + + //! Flag indicating whether this actuator component is active in the current + //! MPI rank + bool actuator_in_proc{false}; + + ActInfo(const std::string& label_in, const int id_in) + : label(label_in), id(id_in) + {} +}; + +/** Abstract representation of data holder for specific actuator types + * + */ +template +struct ActDataHolder +{ + //! Instance of the CFD simulation environment + CFDSim& m_sim; + + //! Basic information about this actuator component in relation to the + //! simulation (usually ActInfo) + typename ActTrait::InfoType m_info; + + //! Nodal data on the actuator grid defined for this component (usually + //! ActGrid) + typename ActTrait::GridType m_grid; + + //! Additional data necessary for a given actuator type + typename ActTrait::MetaType m_meta; + + /** Initialize the data structures + * + * \param sim Instance of the CFD simulation environment + * \param label Unique string identifier for this actuator component + * \param id Unique global integer identifier for this actuator component + */ + ActDataHolder(CFDSim& sim, const std::string& label, const int id) + : m_sim(sim), m_info(label, id), m_grid(), m_meta() + {} +}; + +} // namespace actuator +} // namespace amr_wind + +#endif /* ACTUATOR_TYPES_H */ diff --git a/amr-wind/wind_energy/actuator/actuator_utils.H b/amr-wind/wind_energy/actuator/actuator_utils.H new file mode 100644 index 0000000000..dc1b9e4123 --- /dev/null +++ b/amr-wind/wind_energy/actuator/actuator_utils.H @@ -0,0 +1,29 @@ +#ifndef ACUTATOR_UTILS_H +#define ACUTATOR_UTILS_H + +#include "AMReX_AmrCore.H" + +#include + +namespace amr_wind { +namespace actuator { +namespace utils { + +/** Return a set of process IDs (MPI ranks) that contain AMR boxes that interact + * with a given actuator body. + * + * The region of influence of the actuator body is determined by checking for + * intersections with a bounding box ``rbox``. + * + * \param mesh AMReX mesh instance + * \param rbox The bounding box that defines the region of influence of a + * turbine + */ +std::set determine_influenced_procs( + const amrex::AmrCore& mesh, const amrex::RealBox& rbox); + +} // namespace utils +} // namespace actuator +} // namespace amr_wind + +#endif /* ACUTATOR_UTILS_H */ diff --git a/amr-wind/wind_energy/actuator/actuator_utils.cpp b/amr-wind/wind_energy/actuator/actuator_utils.cpp new file mode 100644 index 0000000000..a6fd59df05 --- /dev/null +++ b/amr-wind/wind_energy/actuator/actuator_utils.cpp @@ -0,0 +1,66 @@ +#include "amr-wind/wind_energy/actuator/actuator_utils.H" + +namespace amr_wind { +namespace actuator { +namespace utils { + +namespace { + +/** Convert a bounding box into amrex::Box index space at a given level + * + * \param rbx Bounding box as defined in global domain coordinates + * \param geom AMReX geometry information for a given level + * \return The Box instance that defines the index space equivalent to bounding + * boxt + */ +amrex::Box +realbox_to_box(const amrex::RealBox& rbx, const amrex::Geometry& geom) +{ + const auto* problo = geom.ProbLo(); + const auto* probhi = geom.ProbHi(); + const auto* dxi = geom.InvCellSize(); + + amrex::IntVect lo, hi; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + amrex::Real bbox_min = amrex::max(rbx.lo()[i], problo[i]); + amrex::Real bbox_max = amrex::min(rbx.hi()[i], probhi[i]); + + amrex::Real rlo = amrex::Math::floor((bbox_min - problo[i]) * dxi[i]); + amrex::Real rhi = amrex::Math::ceil((bbox_max - problo[i]) * dxi[i]); + + lo[i] = static_cast(rlo); + hi[i] = static_cast(rhi); + } + + return amrex::Box{lo, hi}; +} + +} // namespace + +std::set determine_influenced_procs( + const amrex::AmrCore& mesh, const amrex::RealBox& rbx) +{ + std::set procs; + const int finest_level = mesh.finestLevel(); + const int nlevels = mesh.finestLevel() + 1; + auto bx = realbox_to_box(rbx, mesh.Geom(0)); + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& ba = mesh.boxArray(lev); + const auto& dm = mesh.DistributionMap(lev); + + // Get all possible intersections at this level + const auto& isects = ba.intersections(bx); + + // Extract the processor ranks + for (const auto& is : isects) procs.insert(dm[is.first]); + + if (lev < finest_level) bx = bx.refine(mesh.refRatio(lev)); + } + + return procs; +} + +} // namespace utils +} // namespace actuator +} // namespace amr_wind diff --git a/unit_tests/wind_energy/CMakeLists.txt b/unit_tests/wind_energy/CMakeLists.txt index aea07bc8bf..e5517c04d3 100644 --- a/unit_tests/wind_energy/CMakeLists.txt +++ b/unit_tests/wind_energy/CMakeLists.txt @@ -6,3 +6,5 @@ target_sources(${amr_wind_unit_test_exe_name} PRIVATE test_abl_init.cpp test_abl_src.cpp ) + +add_subdirectory(actuator) diff --git a/unit_tests/wind_energy/actuator/CMakeLists.txt b/unit_tests/wind_energy/actuator/CMakeLists.txt new file mode 100644 index 0000000000..fdd3b439e0 --- /dev/null +++ b/unit_tests/wind_energy/actuator/CMakeLists.txt @@ -0,0 +1,4 @@ +target_sources(${amr_wind_unit_test_exe_name} PRIVATE + test_actuator_sampling.cpp + test_actuator_flat_plate.cpp + ) diff --git a/unit_tests/wind_energy/actuator/test_actuator_flat_plate.cpp b/unit_tests/wind_energy/actuator/test_actuator_flat_plate.cpp new file mode 100644 index 0000000000..3eecafbb7d --- /dev/null +++ b/unit_tests/wind_energy/actuator/test_actuator_flat_plate.cpp @@ -0,0 +1,152 @@ +#include "aw_test_utils/MeshTest.H" + +#include "amr-wind/wind_energy/actuator/ActuatorContainer.H" +#include "amr-wind/wind_energy/actuator/ActuatorModel.H" +#include "amr-wind/wind_energy/actuator/actuator_types.H" +#include "amr-wind/core/gpu_utils.H" +#include "amr-wind/core/vs/vector_space.H" + +namespace amr_wind_tests { +namespace { +class ActFlatPlateTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + + { + amrex::ParmParse pp("amr"); + amrex::Vector ncell{{32, 32, 32}}; + pp.add("max_level", 0); + pp.add("max_grid_size", 16); + pp.addarr("n_cell", ncell); + } + { + amrex::ParmParse pp("geometry"); + amrex::Vector problo{{0.0, 0.0, 0.0}}; + amrex::Vector probhi{{32.0, 32.0, 32.0}}; + + pp.addarr("prob_lo", problo); + pp.addarr("prob_hi", probhi); + } + } +}; + +namespace act = amr_wind::actuator; +namespace vs = amr_wind::vs; + +struct FlatPlateData +{ + // Number of point along the wing + int num_pts{11}; + // Pitch angle (i.e., angle of attack) of wing + amrex::Real pitch{6.0}; + // User input for gaussian smearing parameter + vs::Vector eps_inp{2.0, 2.0, 2.0}; + // Normal to the wing (used to compute chord direction) + vs::Vector normal{0.0, 0.0, 1.0}; + + vs::Vector begin{0.0, 0.0, 0.0}; + vs::Vector end{0.0, 0.0, 0.0}; +}; + +struct FlatPlate : public act::ActuatorType +{ + using InfoType = act::ActInfo; + using GridType = act::ActGrid; + using MetaType = FlatPlateData; + using DataType = act::ActDataHolder; + + static const std::string identifier() { return "FlatPlate"; } +}; + +} // namespace + +} // namespace amr_wind_tests + +namespace amr_wind { +namespace actuator { + +namespace ops { + +template <> +void read_inputs<::amr_wind_tests::FlatPlate>( + ::amr_wind_tests::FlatPlate::DataType& data) +{ + auto& info = data.m_info; + auto& fpobj = data.m_meta; + + fpobj.begin = {16.0, 12.0, 16.0}; + fpobj.end = {16.0, 20.0, 16.0}; + + // Create a bounding box extending to max search radius in each direction + info.bound_box = amrex::RealBox(8.0, 4.0, 8.0, 24.0, 28.0, 24.0); +} + +template <> +void init_data_structures<::amr_wind_tests::FlatPlate>( + ::amr_wind_tests::FlatPlate::DataType& data) +{ + auto& grid = data.m_grid; + auto& obj = data.m_meta; + + int npts = obj.num_pts; + grid.resize(npts); + + // Wing span + auto wspan = obj.end - obj.begin; + // Compute chord/flow direction as a cross-product + auto chord = (wspan ^ obj.normal); + // Set up global to local transformation matrix + auto tmat = vs::Tensor(chord.unit(), wspan.unit(), obj.normal.unit()); + + // Equal spacing along span + auto dx = (1.0 / static_cast(npts - 1)) * wspan; + + for (int i = 0; i < npts; ++i) { + grid.pos[i] = obj.begin + static_cast(i) * dx; + grid.epsilon[i] = obj.eps_inp; + grid.orientation[i] = tmat; + } + + // Initialize remaining data + grid.force.assign(npts, vs::Vector::zero()); + grid.vel_pos.assign(grid.pos.begin(), grid.pos.end()); + grid.vel.assign(npts, vs::Vector::zero()); + + { + // Check that the wing coordinates were initialized correctly + auto wing_len = vs::mag(grid.pos.back() - grid.pos.front()); + EXPECT_NEAR(wing_len, 8.0, 1.0e-12); + + // Check that the orientation tensor was initialized correctly + EXPECT_NEAR(chord.unit() & vs::Vector::ihat(), 1.0, 1.0e-12); + EXPECT_NEAR(vs::mag_sqr(tmat), 3.0, 1.0e-12); + } +} + +} // namespace ops + +template class ::amr_wind::actuator:: + ActModel<::amr_wind_tests::FlatPlate, ::amr_wind::actuator::ActSrcLine>; +} // namespace actuator +} // namespace amr_wind + +namespace amr_wind_tests { + +TEST_F(ActFlatPlateTest, test_init) +{ + initialize_mesh(); + + ::amr_wind::actuator::ActModel flat_plate( + sim(), "flat_plate1", 0); + flat_plate.pre_init_actions(); + flat_plate.post_init_actions(); + + const auto& info = flat_plate.info(); + EXPECT_EQ(info.root_proc, 0); + EXPECT_EQ(info.procs.size(), amrex::ParallelDescriptor::NProcs()); +} + +} // namespace amr_wind_tests diff --git a/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp new file mode 100644 index 0000000000..b3510c3716 --- /dev/null +++ b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp @@ -0,0 +1,199 @@ +#include "aw_test_utils/MeshTest.H" + +#include "amr-wind/wind_energy/actuator/ActuatorContainer.H" +#include "amr-wind/core/gpu_utils.H" +#include "amr-wind/core/vs/vector_space.H" + +#include + +namespace amr_wind_tests { +namespace { + +// Utility function to populate the velocity field used for tests +void init_field(amr_wind::Field& fld) +{ + const auto& mesh = fld.repo().mesh(); + const int nlevels = fld.repo().num_active_levels(); + const int ncomp = fld.num_comp(); + + amrex::Real offset = 0.0; + if (fld.field_location() == amr_wind::FieldLoc::CELL) offset = 0.5; + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = mesh.Geom(lev).CellSizeArray(); + const auto& problo = mesh.Geom(lev).ProbLoArray(); + + for (amrex::MFIter mfi(fld(lev)); mfi.isValid(); ++mfi) { + auto bx = mfi.growntilebox(); + const auto& farr = fld(lev).array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + const amrex::Real x = problo[0] + (i + offset) * dx[0]; + const amrex::Real y = problo[1] + (j + offset) * dx[1]; + const amrex::Real z = problo[2] + (k + offset) * dx[2]; + + for (int d = 0; d < ncomp; d++) farr(i, j, k, d) = x + y + z; + }); + } + } +} + +/** Mock test object to allow access to the data struct + */ +class TestActContainer : public amr_wind::actuator::ActuatorContainer +{ +public: + TestActContainer(amrex::AmrCore& mesh, const int nobj) + : amr_wind::actuator::ActuatorContainer(mesh, nobj) + {} + + // Accessor for the particle data holder object + amr_wind::actuator::ActuatorCloud& get_data_obj() { return point_data(); } +}; + +class ActuatorTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + + { + amrex::ParmParse pp("amr"); + amrex::Vector ncell{{32, 32, 64}}; + pp.add("max_level", 0); + pp.add("max_grid_size", 16); + pp.addarr("n_cell", ncell); + } + { + amrex::ParmParse pp("geometry"); + amrex::Vector problo{{0.0, 0.0, 0.0}}; + amrex::Vector probhi{{128.0, 128.0, 128.0}}; + + pp.addarr("prob_lo", problo); + pp.addarr("prob_hi", probhi); + } + } +}; + +} // namespace + +TEST_F(ActuatorTest, act_container) +{ + const int nprocs = amrex::ParallelDescriptor::NProcs(); + if (nprocs > 2) { + GTEST_SKIP(); + return; + } + + const int iproc = amrex::ParallelDescriptor::MyProc(); + initialize_mesh(); + auto& vel = sim().repo().declare_field("velocity", 3, 3); + init_field(vel); + + // Number of turbines in an MPI rank + const int num_turbines = 2; + // Number of points per turbine + const int num_nodes = 16; + + TestActContainer ac(mesh(), num_turbines); + auto& data = ac.get_data_obj(); + + for (int it = 0; it < num_turbines; ++it) { + data.num_pts[it] = num_nodes; + } + + ac.initialize_container(); + + { + const int lev = 0; + int idx = 0; + const amrex::Real dz = mesh().Geom(lev).CellSize(2); + const amrex::Real ypos = 32.0 * (iproc + 1); + auto& pvec = data.position; + for (int it = 0; it < num_turbines; ++it) { + const amrex::Real xpos = 32.0 * (it + 1); + for (int ni = 0; ni < num_nodes; ++ni) { + const amrex::Real zpos = (ni + 0.5) * dz; + + pvec[idx].x() = xpos; + pvec[idx].y() = ypos; + pvec[idx].z() = zpos; + ++idx; + } + } + ASSERT_EQ(idx, ac.num_actuator_points()); + } + + ac.update_positions(); + + // Check that the update positions scattered the particles to the MPI rank + // that contains the cell enclosing the particle location +#ifndef AMREX_USE_GPU + ASSERT_EQ( + ac.num_actuator_points() * nprocs, ac.NumberOfParticlesAtLevel(0)); +#endif + { + using ParIter = amr_wind::actuator::ActuatorContainer::ParIterType; + int counter = 0; + int total_particles = 0; + const int lev = 0; + for (ParIter pti(ac, lev); pti.isValid(); ++pti) { + ++counter; + total_particles += pti.numParticles(); + } + + if (iproc == 0) { + ASSERT_EQ(num_turbines * num_nodes * nprocs, total_particles); + } else { + ASSERT_EQ(total_particles, 0); + } + } + + ac.sample_velocities(vel); + + // Check to make sure that the velocity sampling gathered the particles back + // to their original MPI ranks +#ifndef AMREX_USE_GPU + ASSERT_EQ( + ac.num_actuator_points() * nprocs, ac.NumberOfParticlesAtLevel(0)); +#endif + { + using ParIter = amr_wind::actuator::ActuatorContainer::ParIterType; + int counter = 0; + int total_particles = 0; + const int lev = 0; + for (ParIter pti(ac, lev); pti.isValid(); ++pti) { + ++counter; + total_particles += pti.numParticles(); + } + + // All particles should have been recalled back to the same patch within + // each MPI rank + ASSERT_EQ(counter, 1); + // Total number of particles should be the same as what this MPI rank + // created + ASSERT_EQ(num_turbines * num_nodes, total_particles); + } + + // Check the interpolated velocity field + { + namespace vs = amr_wind::vs; + constexpr amrex::Real rtol = 1.0e-12; + amrex::Real rerr = 0.0; + const int npts = ac.num_actuator_points(); + const auto& pvec = data.position; + const auto& vvec = data.velocity; + for (int ip = 0; ip < npts; ++ip) { + const auto& pos = pvec[ip]; + const auto& pvel = vvec[ip]; + + const amrex::Real vval = pos.x() + pos.y() + pos.z(); + const vs::Vector vgold{vval, vval, vval}; + rerr += vs::mag_sqr(pvel - vgold); + } + EXPECT_NEAR(rerr, 0.0, rtol); + } +} + +} // namespace amr_wind_tests