Actuators: sample mac velocity instead of cell-centered velocity (#868)

* new sampling is default for TurbineFastLine and TurbineFastDisk
* big improvement shown for ALM
* new input parameter available for Actuator, added to documentation
* unit test confirms implementation of face velocity sampling

amr-wind/wind_energy/actuator/Actuator.H

@@ -75,6 +75,8 @@ private:
     std::vector<std::unique_ptr<ActuatorModel>> m_actuators;
 
     std::unique_ptr<ActuatorContainer> m_container;
+
+    bool m_sample_nmhalf{false};
 };
 
 } // namespace actuator

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

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<int>(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);

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<int>(std::floor(x));
+                const int j = static_cast<int>(std::floor(y));
+                const int k = static_cast<int>(std::floor(z));
+
+                const int i_f = static_cast<int>(std::floor(xf));
+                const int j_f = static_cast<int>(std::floor(yf));
+                const int k_f = static_cast<int>(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

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
 """""""""""""
 

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;