FLLC nonuniform actuator point distribution (Exawind#1045)

* Added documentation for new options
* Modified FLLC unit test and added nonuniform test

amr-wind/wind_energy/actuator/FLLC.H

@@ -19,6 +19,7 @@ struct FLLCData
     amrex::Real epsilon;
     amrex::Real relaxation_factor{0.1};
     amrex::Real fllc_start_time{0.0};
+    RealList r;
     RealList dr;
     RealList optimal_epsilon;
 
@@ -33,6 +34,16 @@ struct FLLCData
     bool different_sizes;
     RealList span_distance_vel;
     RealList span_distance_force;
+
+    // non-uniform variables
+    bool nonuniform{true};      // non-uniform flag
+    amrex::Real eps_dr{1.};     // the ratio of epsilon to actuator width
+    RealList nonuniform_r;      // non-uniform radius
+    RealList nonuniform_dr;     // non-uniform spacing
+    VecList vel_rel;            // uniform relative velocity
+    VecList nonuniform_vel_rel; // non-uniform relative velocity
+    RealList nonuniform_optimal_epsilon; // non-uniform radius
+    VecList nonuniform_lift;             // non-uniform lift
 };
 
 /**

amr-wind/wind_energy/actuator/FLLC.cpp

@@ -1,10 +1,13 @@
 #include "amr-wind/wind_energy/actuator/FLLC.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 namespace amr_wind::actuator {
 
 void FLLCInit(
     FLLCData& data, const ComponentView& view, const amrex::Real eps_chord)
 {
+
+    namespace interp = ::amr_wind::interp;
     const int npts = static_cast<int>(view.pos.size());
     data.different_sizes = view.pos.size() != view.vel_pos.size();
     if (data.different_sizes) {
@@ -17,7 +20,9 @@ void FLLCInit(
             data.span_distance_vel[i] = vs::mag(view.vel_pos[i]);
         }
     }
+    data.r.resize(npts);
     data.dr.resize(npts);
+    data.vel_rel.resize(npts);
     data.optimal_epsilon.resize(npts);
     data.les_velocity.assign(npts, vs::Vector::zero());
     data.optimal_velocity.assign(npts, vs::Vector::zero());
@@ -27,7 +32,7 @@ void FLLCInit(
     data.force_point_velocity.assign(npts, vs::Vector::zero());
 
     for (int i = 0; i < npts; ++i) {
-
+        data.r[i] = vs::mag(view.pos[i] - view.pos[0]);
         data.optimal_epsilon[i] = view.chord[i] * eps_chord;
     }
 
@@ -38,6 +43,65 @@ void FLLCInit(
     data.dr[0] = vs::mag(view.pos[1] - view.pos[0]) / 2.;
     data.dr[npts - 1] = vs::mag(view.pos[npts - 1] - view.pos[npts - 2]) / 2.;
 
+    //
+    //  The following code is used to create a non-uniform distribution of
+    //  points. This should follow a spacing given by epsilon/dr and values
+    //  above 1 are recommended for converged solutions.
+    //
+    if (data.nonuniform) {
+        // Non-uniform radial distribution variables
+        amrex::Real dr;
+        amrex::Real dr1;
+        amrex::Real eps1;
+        amrex::Real dr2;
+        amrex::Real eps2;
+
+        // Initialize the first location
+        data.nonuniform_r.push_back(data.r[0]);
+        while (data.nonuniform_r.back() < data.r.back()) {
+
+            amrex::Real r_ =
+                data.nonuniform_r.back(); // The latest radial position
+
+            // Interpolate the value of epsilon to the current radial location
+            eps1 = interp::linear(data.r, data.optimal_epsilon, r_);
+            dr1 = eps1 / data.eps_dr; // spacing needed to match condition
+
+            // Interpolate the value of epsilon to the next radial location
+            eps2 = interp::linear(data.r, data.optimal_epsilon, r_ + dr1);
+            dr2 = eps2 / data.eps_dr;
+
+            // This will ensure that the spacing will always meet the
+            // requirement eps/dr
+            dr = std::min(dr1, dr2);
+            AMREX_ALWAYS_ASSERT(
+                dr > std::numeric_limits<amrex::Real>::epsilon());
+
+            // Append value to the array
+            // Ensure that the value is smaller than the tip
+            data.nonuniform_r.push_back(std::min(r_ + dr, data.r.back()));
+        }
+
+        int npts_n = static_cast<int>(data.nonuniform_r.size());
+        data.nonuniform_dr.resize(npts_n);
+        data.nonuniform_vel_rel.resize(npts_n);
+        data.nonuniform_optimal_epsilon.resize(npts_n);
+        data.nonuniform_lift.resize(npts_n);
+
+        // Central difference for non-uniform points
+        for (int i = 1; i < npts_n - 1; ++i) {
+            data.nonuniform_dr[i] =
+                std::abs(data.nonuniform_r[i + 1] - data.nonuniform_r[i - 1]) /
+                2.;
+        }
+        data.nonuniform_dr[0] =
+            std::abs(data.nonuniform_r[1] - data.nonuniform_r[0]) / 2.;
+        data.nonuniform_dr[npts_n - 1] =
+            std::abs(
+                data.nonuniform_r[npts_n - 1] - data.nonuniform_r[npts_n - 2]) /
+            2.;
+    }
+
     data.initialized = true;
 }
 
@@ -49,6 +113,8 @@ void FLLCParse(const utils::ActParser& pp, FLLCData& data)
     std::string typeString = "variable_chord";
     pp.query("fllc_type", typeString);
     data.correction_type = FLLCTypeMap.at(typeString);
+    pp.query("fllc_nonuniform", data.nonuniform);
+    pp.query("fllc_epsilon_dr_ratio", data.eps_dr);
 
     if (!pp.contains("epsilon") || !pp.contains("epsilon_chord")) {
         amrex::Abort(

amr-wind/wind_energy/actuator/FLLCOp.H

@@ -162,9 +162,11 @@ struct FLLCOp
         const int npts = static_cast<int>(fllc.correction_velocity.size());
         auto& du = fllc.correction_velocity;
 
+        auto& nonuniform_r = fllc.nonuniform_r;
         auto& u_les = fllc.les_velocity;
         auto& u_opt = fllc.optimal_velocity;
         auto& G = fllc.lift;
+        auto& nonuniform_G = fllc.nonuniform_lift;
         VecSlice u_force_pnt_slice =
             ::amr_wind::utils::slice(fllc.force_point_velocity, 0);
         vs::Vector* vel_ptr = data.vel.data();
@@ -191,12 +193,30 @@ struct FLLCOp
             const auto vmag =
                 std::max(vs::mag(vel), vs::DTraits<amrex::Real>::eps());
             const auto vmag2 = vmag * vmag;
-
+            fllc.vel_rel[ip] = vel;
             const auto fv = force & vel;
 
             G[ip] = (force - vel * fv / vmag2) / dr;
         }
 
+        // Get the non-uniform distribution of force
+        if (fllc.nonuniform) {
+            interp::linear(fllc.r, G, nonuniform_r, nonuniform_G);
+            interp::linear(
+                fllc.r, fllc.optimal_epsilon, nonuniform_r,
+                fllc.nonuniform_optimal_epsilon);
+            interp::linear(
+                fllc.r, fllc.vel_rel, nonuniform_r, fllc.nonuniform_vel_rel);
+        }
+
+        const auto optimal_epsilon = (fllc.nonuniform)
+                                         ? fllc.nonuniform_optimal_epsilon
+                                         : fllc.optimal_epsilon;
+        const auto dr = (fllc.nonuniform) ? fllc.nonuniform_dr : fllc.dr;
+        const auto vel_rel =
+            (fllc.nonuniform) ? fllc.nonuniform_vel_rel : fllc.vel_rel;
+        const auto G_new = (fllc.nonuniform) ? nonuniform_G : G;
+
         /**
          * Step 2
          * Compute the induced velocities
@@ -205,19 +225,24 @@ struct FLLCOp
 
             const auto coefficient = 1.0 / (2.0 * amr_wind::utils::pi());
 
-            for (int jp = 0; jp < npts; ++jp) {
+            for (int jp = 0; jp < dr.size(); ++jp) {
 
                 const auto eps_les = fllc.epsilon;
-                const auto eps_opt = fllc.optimal_epsilon[jp];
+                const auto eps_opt = optimal_epsilon[jp];
                 const auto eps_les2 = eps_les * eps_les;
                 const auto eps_opt2 = eps_opt * eps_opt;
-                const auto& vel = data.vel_rel[jp];
+                const auto& vel = vel_rel[jp];
                 const auto vmag =
                     std::max(vs::mag(vel), vs::DTraits<amrex::Real>::eps());
                 auto k_les = 0.;
                 auto k_opt = 0.;
 
-                if (ip == jp) {
+                const auto r_dis =
+                    (fllc.nonuniform)
+                        ? std::abs(fllc.r[ip] - fllc.nonuniform_r[jp])
+                        : vs::mag(data.vel_pos[ip] - data.vel_pos[jp]);
+
+                if (r_dis == 0.) {
 
                     k_les = .5 / (eps_les2);
                     k_opt = .5 / (eps_opt2);
@@ -226,8 +251,6 @@ struct FLLCOp
 
                 else {
 
-                    const auto r_dis =
-                        vs::mag(data.vel_pos[ip] - data.vel_pos[jp]);
                     auto exp_les = std::exp(-r_dis * r_dis / eps_les2);
                     auto exp_opt = std::exp(-r_dis * r_dis / eps_opt2);
 
@@ -237,10 +260,10 @@ struct FLLCOp
                             1. / (2. * r_dis * r_dis) * (exp_opt - 1.);
                 }
 
-                u_les[ip] = u_les[ip] +
-                            coefficient * G[jp] / vmag * k_les * fllc.dr[jp];
-                u_opt[ip] = u_opt[ip] +
-                            coefficient * G[jp] / vmag * k_opt * fllc.dr[jp];
+                u_les[ip] =
+                    u_les[ip] + coefficient * G_new[jp] / vmag * k_les * dr[jp];
+                u_opt[ip] =
+                    u_opt[ip] + coefficient * G_new[jp] / vmag * k_opt * dr[jp];
             }
 
             // Relaxation to compute the induced velocity

docs/sphinx/user/inputs_Actuator.rst

@@ -99,19 +99,36 @@ Example for ``FixedWingLine``::
 
 .. input_param:: Actuator.FixedWingLine.fllc_relaxation_factor
 
-  **type:** Double
+  **type:** Double, optional
 
   The relaxation factor to be applied to the updated velocity see:
   `Martinez-Tossas and Meneveau (2019) <https://doi.org/10.1017/jfm.2018.994>`_
   The default value is `0.1`.
 
 .. input_param:: Actuator.FixedWingLine.fllc_start_time
 
-  **type:** Double
+  **type:** Double, optional
 
   The time in the simulation from when to start using the correction.
   The default value is `0`.
 
+.. input_param:: Actuator.FixedWingLine.fllc_nonuniform
+
+  **type:** Bool
+
+  The flag to specify if the actuator points used to compute the correction should be
+  non-uniformly distributed. This helps in using less points for the fllc while
+  still maintaining the accuracy of the fllc.
+  The default value is `true`.
+
+.. input_param:: Actuator.FixedWingLine.fllc_epsilon_dr_ratio
+
+  **type:** Double, optional
+
+  The ratio of epsilon to actuator point spacing used to create a non-uniform distribution.
+  A value of `1` or greater is recommended.
+  The default value is `1`.
+
 .. input_param:: Actuator.FixedWingLine.pitch
 
    **type:** Real number, optional

unit_tests/wind_energy/actuator/test_FLLC.cpp

@@ -22,6 +22,8 @@ TEST(TestFLLCData, data_initializes_with_cviews)
     view.chord = ::amr_wind::utils::slice(dummy_real, 0, num_points);
     view.vel_rel = ::amr_wind::utils::slice(dummy_vec, 0, num_points);
 
+    data.nonuniform = false;
+
     FLLCInit(data, view, 1.0);
 
     ASSERT_EQ(num_points, data.les_velocity.size());
@@ -30,4 +32,38 @@ TEST(TestFLLCData, data_initializes_with_cviews)
     ASSERT_EQ(num_points, data.grad_lift.size());
 }
 
+TEST(TestFLLCData, data_initializes_with_cviews_nonuniform)
+{
+    const int num_points = 5;
+
+    VecList dummy_vec(num_points);
+    RealList dummy_real(num_points);
+    TensorList dummy_tensor(num_points);
+
+    FLLCData data;
+
+    ComponentView view;
+    view.pos = ::amr_wind::utils::slice(dummy_vec, 0, num_points);
+    view.force = ::amr_wind::utils::slice(dummy_vec, 0, num_points);
+    view.epsilon = ::amr_wind::utils::slice(dummy_vec, 0, num_points);
+    view.orientation = ::amr_wind::utils::slice(dummy_tensor, 0, num_points);
+    view.chord = ::amr_wind::utils::slice(dummy_real, 0, num_points);
+    view.vel_rel = ::amr_wind::utils::slice(dummy_vec, 0, num_points);
+
+    data.nonuniform = true;
+    data.eps_dr = 1.1;
+
+    for (int ip = 0; ip < num_points; ++ip) {
+        view.pos[ip] = {0, static_cast<amrex::Real>(ip), 0};
+        view.chord[ip] = 1.;
+    }
+
+    FLLCInit(data, view, 1.0);
+
+    const int npts_r = static_cast<int>(data.nonuniform_r.size());
+    const int npts_dr = static_cast<int>(data.nonuniform_dr.size());
+
+    ASSERT_EQ(npts_r, npts_dr);
+}
+
 } // namespace amr_wind::actuator