diff --git a/amr-wind/wind_energy/actuator/Actuator.H b/amr-wind/wind_energy/actuator/Actuator.H index 7d6c2bda36..1486c55083 100644 --- a/amr-wind/wind_energy/actuator/Actuator.H +++ b/amr-wind/wind_energy/actuator/Actuator.H @@ -75,6 +75,8 @@ private: std::vector> m_actuators; std::unique_ptr m_container; + + bool m_sample_nmhalf{false}; }; } // namespace actuator diff --git a/amr-wind/wind_energy/actuator/Actuator.cpp b/amr-wind/wind_energy/actuator/Actuator.cpp index d8cb7e5698..0ecc66ae2e 100644 --- a/amr-wind/wind_energy/actuator/Actuator.cpp +++ b/amr-wind/wind_energy/actuator/Actuator.cpp @@ -34,6 +34,12 @@ void Actuator::pre_init_actions() std::string type; pp.query("type", type); pp1.query("type", type); + if (type == "TurbineFastLine" || type == "TurbineFastDisk") { + // Only one kind of sampling can be chosen. If OpenFAST is involved + // in the Actuator type, the default behavior is to sample velocity + // at n-1/2 so that forcing is at n+1/2 + m_sample_nmhalf = true; + } AMREX_ALWAYS_ASSERT(!type.empty()); auto obj = ActuatorModel::create(type, m_sim, tname, i); @@ -44,6 +50,9 @@ void Actuator::pre_init_actions() obj->read_inputs(inp); m_actuators.emplace_back(std::move(obj)); } + + // Check if sampling should be modified aside from default behavior + pp.query("sample_vel_nmhalf", m_sample_nmhalf); } void Actuator::post_init_actions() @@ -179,9 +188,17 @@ void Actuator::update_positions() m_container->update_positions(); // Sample velocities at the new locations - const auto& vel = m_sim.repo().get_field("velocity"); - const auto& density = m_sim.repo().get_field("density"); - m_container->sample_fields(vel, density); + if (m_sample_nmhalf) { + const auto& umac = m_sim.repo().get_field("u_mac"); + const auto& vmac = m_sim.repo().get_field("v_mac"); + const auto& wmac = m_sim.repo().get_field("w_mac"); + const auto& density = m_sim.repo().get_field("density"); + m_container->sample_fields(umac, vmac, wmac, density); + } else { + const auto& vel = m_sim.repo().get_field("velocity"); + const auto& density = m_sim.repo().get_field("density"); + m_container->sample_fields(vel, density); + } } /** Provide updated velocities from container to actuator instances diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.H b/amr-wind/wind_energy/actuator/ActuatorContainer.H index 7cf57cf345..030862b2b2 100644 --- a/amr-wind/wind_energy/actuator/ActuatorContainer.H +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.H @@ -78,6 +78,12 @@ public: void sample_fields(const Field& vel, const Field& density); + void sample_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density); + int num_actuator_points() const { return static_cast(m_data.position.size()); @@ -87,6 +93,12 @@ public: void interpolate_fields(const Field& vel, const Field& density); + void interpolate_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density); + void populate_field_buffers(); void initialize_particles(const int total_pts); diff --git a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp index 8a448d493f..0b4db03e40 100644 --- a/amr-wind/wind_energy/actuator/ActuatorContainer.cpp +++ b/amr-wind/wind_energy/actuator/ActuatorContainer.cpp @@ -213,6 +213,29 @@ void ActuatorContainer::sample_fields(const Field& vel, const Field& density) m_is_scattered = false; } +void ActuatorContainer::sample_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density) +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::sample_velocities"); + AMREX_ALWAYS_ASSERT(m_container_initialized && m_is_scattered); + + // Sample velocity field + interpolate_fields(umac, vmac, wmac, density); + + // Recall particles to the MPI ranks that contains their corresponding + // turbines + // Redistribute(); + + // Populate the velocity buffer that all actuator instances can access + populate_field_buffers(); + + // Indicate that the particles have been restored to their original MPI rank + m_is_scattered = false; +} + /** Helper method for ActuatorContainer::sample_fields * * Loops over the particle tiles and copies velocity data from particles to the @@ -356,6 +379,128 @@ void ActuatorContainer::interpolate_fields( } } +void ActuatorContainer::interpolate_fields( + const Field& umac, + const Field& vmac, + const Field& wmac, + const Field& density) +{ + BL_PROFILE("amr-wind::actuator::ActuatorContainer::interpolate_velocities"); + 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 uarr = umac(lev).const_array(pti); + const auto varr = vmac(lev).const_array(pti); + const auto warr = wmac(lev).const_array(pti); + const auto darr = density(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]; + + const amrex::Real xf = (pp.pos(0) - plo[0]) * dxi[0]; + const amrex::Real yf = (pp.pos(1) - plo[1]) * dxi[1]; + const amrex::Real zf = (pp.pos(2) - plo[2]) * dxi[2]; + + // Index of the low corner + const int i = static_cast(std::floor(x)); + const int j = static_cast(std::floor(y)); + const int k = static_cast(std::floor(z)); + + const int i_f = static_cast(std::floor(xf)); + const int j_f = static_cast(std::floor(yf)); + const int k_f = static_cast(std::floor(zf)); + + // 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 amrex::Real wxf_hi = (xf - i_f); + const amrex::Real wyf_hi = (yf - j_f); + const amrex::Real wzf_hi = (zf - k_f); + + const amrex::Real wxf_lo = 1.0 - wxf_hi; + const amrex::Real wyf_lo = 1.0 - wyf_hi; + const amrex::Real wzf_lo = 1.0 - wzf_hi; + + const int iproc = pp.cpu(); + + // u + int ic = 0; + pp.rdata(ic) = + wxf_lo * wy_lo * wz_lo * uarr(i_f, j, k) + + wxf_lo * wy_lo * wz_hi * uarr(i_f, j, k + 1) + + wxf_lo * wy_hi * wz_lo * uarr(i_f, j + 1, k) + + wxf_lo * wy_hi * wz_hi * uarr(i_f, j + 1, k + 1) + + wxf_hi * wy_lo * wz_lo * uarr(i_f + 1, j, k) + + wxf_hi * wy_lo * wz_hi * uarr(i_f + 1, j, k + 1) + + wxf_hi * wy_hi * wz_lo * uarr(i_f + 1, j + 1, k) + + wxf_hi * wy_hi * wz_hi * uarr(i_f + 1, j + 1, k + 1); + + // Reset position vectors so that the particles return back + // to the MPI ranks with the turbines upon redistribution + pp.pos(ic) = dptr[iproc][ic]; + + // v + ic = 1; + pp.rdata(ic) = + wx_lo * wyf_lo * wz_lo * varr(i, j_f, k) + + wx_lo * wyf_lo * wz_hi * varr(i, j_f, k + 1) + + wx_lo * wyf_hi * wz_lo * varr(i, j_f + 1, k) + + wx_lo * wyf_hi * wz_hi * varr(i, j_f + 1, k + 1) + + wx_hi * wyf_lo * wz_lo * varr(i + 1, j_f, k) + + wx_hi * wyf_lo * wz_hi * varr(i + 1, j_f, k + 1) + + wx_hi * wyf_hi * wz_lo * varr(i + 1, j_f + 1, k) + + wx_hi * wyf_hi * wz_hi * varr(i + 1, j_f + 1, k + 1); + pp.pos(ic) = dptr[iproc][ic]; + + // w + ic = 2; + pp.rdata(ic) = + wx_lo * wy_lo * wzf_lo * warr(i, j, k_f) + + wx_lo * wy_lo * wzf_hi * warr(i, j, k_f + 1) + + wx_lo * wy_hi * wzf_lo * warr(i, j + 1, k_f) + + wx_lo * wy_hi * wzf_hi * warr(i, j + 1, k_f + 1) + + wx_hi * wy_lo * wzf_lo * warr(i + 1, j, k_f) + + wx_hi * wy_lo * wzf_hi * warr(i + 1, j, k_f + 1) + + wx_hi * wy_hi * wzf_lo * warr(i + 1, j + 1, k_f) + + wx_hi * wy_hi * wzf_hi * warr(i + 1, j + 1, k_f + 1); + pp.pos(ic) = dptr[iproc][ic]; + + // density + pp.rdata(AMREX_SPACEDIM) = + wx_lo * wy_lo * wz_lo * darr(i, j, k) + + wx_lo * wy_lo * wz_hi * darr(i, j, k + 1) + + wx_lo * wy_hi * wz_lo * darr(i, j + 1, k) + + wx_lo * wy_hi * wz_hi * darr(i, j + 1, k + 1) + + wx_hi * wy_lo * wz_lo * darr(i + 1, j, k) + + wx_hi * wy_lo * wz_hi * darr(i + 1, j, k + 1) + + wx_hi * wy_hi * wz_lo * darr(i + 1, j + 1, k) + + wx_hi * wy_hi * wz_hi * darr(i + 1, j + 1, k + 1); + }); + } + } +} + /** Determine position vector of a point within each MPI rank * * Loops over the boxArray and determines the first patch that belongs to the diff --git a/docs/sphinx/user/inputs_Actuator.rst b/docs/sphinx/user/inputs_Actuator.rst index 6c0d43d140..a14203639b 100644 --- a/docs/sphinx/user/inputs_Actuator.rst +++ b/docs/sphinx/user/inputs_Actuator.rst @@ -29,6 +29,21 @@ turbines as actuator disks and actuator line models. supported are: ``TurbineFastLine``, ``TurbineFastDisk``, and ``FixedWingLine``. +.. input_param:: Actuator.sample_vel_nmhalf + + **type:** Bool, optional + + This option sets the type of velocity sampling used to inform the actuator model. + Setting this variable to true (or `1`) makes the Actuator velocity sampling use + the face-centered velocity field at the time instant of `n-1/2`. For simulations + coupled to OpenFAST, this enables the actuator forces to be implemented at the + time instant of `n+1/2`, which fits best with the underlying numerical model of AMR-Wind + and has been shown to dramatically improve the accuracy of the Actuator Line Model. + This option defaults to true (`1`) when the Actuator type is ``TurbineFastLine`` or + ``TurbineFastDisk``, and it defaults to false (`0`) in all other cases. When the + option is false, the actuator routine samples the cell-centered velocity + at the time instant `n`. + FixedWingLine """"""""""""" diff --git a/unit_tests/wind_energy/actuator/test_act_utils.H b/unit_tests/wind_energy/actuator/test_act_utils.H index 9bfc058908..9a3e41b0a6 100644 --- a/unit_tests/wind_energy/actuator/test_act_utils.H +++ b/unit_tests/wind_energy/actuator/test_act_utils.H @@ -13,9 +13,25 @@ inline void init_field(amr_wind::Field& fld) const int nlevels = fld.repo().num_active_levels(); const int ncomp = fld.num_comp(); - amrex::Real offset = 0.0; + amrex::Real offsetx = 0.0; + amrex::Real offsety = 0.0; + amrex::Real offsetz = 0.0; if (fld.field_location() == amr_wind::FieldLoc::CELL) { - offset = 0.5; + offsetx = 0.5; + offsety = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::XFACE) { + offsety = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::YFACE) { + offsetx = 0.5; + offsetz = 0.5; + } + if (fld.field_location() == amr_wind::FieldLoc::ZFACE) { + offsetx = 0.5; + offsety = 0.5; } for (int lev = 0; lev < nlevels; ++lev) { @@ -27,9 +43,9 @@ inline void init_field(amr_wind::Field& fld) 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]; + const amrex::Real x = problo[0] + (i + offsetx) * dx[0]; + const amrex::Real y = problo[1] + (j + offsety) * dx[1]; + const amrex::Real z = problo[2] + (k + offsetz) * dx[2]; for (int d = 0; d < ncomp; d++) { farr(i, j, k, d) = x + y + z; diff --git a/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp index cad8ef04c8..c45b709ff9 100644 --- a/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp +++ b/unit_tests/wind_energy/actuator/test_actuator_sampling.cpp @@ -168,4 +168,129 @@ TEST_F(ActuatorTest, act_container) } } +TEST_F(ActuatorTest, act_container_macvels) +{ + const int nprocs = amrex::ParallelDescriptor::NProcs(); + if (nprocs > 2) { + GTEST_SKIP(); + } + + const int iproc = amrex::ParallelDescriptor::MyProc(); + initialize_mesh(); + auto& umac = + sim().repo().declare_field("u_mac", 1, 3, 1, amr_wind::FieldLoc::XFACE); + auto& vmac = + sim().repo().declare_field("v_mac", 1, 3, 1, amr_wind::FieldLoc::YFACE); + auto& wmac = + sim().repo().declare_field("w_mac", 1, 3, 1, amr_wind::FieldLoc::ZFACE); + auto& density = sim().repo().declare_field("density", 1, 3); + init_field(umac); + init_field(vmac); + init_field(wmac); + density.setVal(1.0); + + // 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 total_particles = 0; + const int lev = 0; + for (ParIter pti(ac, lev); pti.isValid(); ++pti) { + total_particles += pti.numParticles(); + } + + if (iproc == 0) { + ASSERT_EQ(num_turbines * num_nodes * nprocs, total_particles); + } else { + ASSERT_EQ(total_particles, 0); + } + } + + ac.sample_fields(umac, vmac, wmac, density); + ac.Redistribute(); + + // 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