2D Pitching Fixed Wing

amr-wind/core/MultiParser.H

@@ -65,6 +65,51 @@ public:
         }
     }
 
+    /** Get a vector from the input file, but could be specified uniform instead
+     *
+     *  \param [in] name Keyword to search for
+     *  \param [out] value Vector value
+     */
+    void get_either(const std::string& name, vs::Vector& value) const
+    {
+        amrex::Vector<vs::Vector::value_type> val;
+        getarr(name, val);
+        AMREX_ALWAYS_ASSERT(val.size() == AMREX_SPACEDIM || val.size() == 1);
+        if (val.size() == 1) {
+            value.x() = val[0];
+            value.y() = val[0];
+            value.z() = val[0];
+        } else {
+            value.x() = val[0];
+            value.y() = val[1];
+            value.z() = val[2];
+        }
+    }
+
+    /** Query and return if a scalar (uniform) or vector exists in input file
+     *
+     *  \param [in] name Keyword to search for
+     *  \param [out] value Vector value
+     */
+    void query_either(const std::string& name, vs::Vector& value) const
+    {
+        amrex::Vector<vs::Vector::value_type> val;
+        queryarr(name, val);
+        if (!val.empty()) {
+            AMREX_ALWAYS_ASSERT(
+                val.size() == AMREX_SPACEDIM || val.size() == 1);
+            if (val.size() == 1) {
+                value.x() = val[0];
+                value.y() = val[0];
+                value.z() = val[0];
+            } else {
+                value.x() = val[0];
+                value.y() = val[1];
+                value.z() = val[2];
+            }
+        }
+    }
+
     //! Get the value for the keyword entry from either namespace
     template <typename T>
     void get(const std::string& name, T& value) const

amr-wind/core/vs/vector.H

@@ -111,6 +111,9 @@ struct VectorT
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VectorT<T> operator*=(const T val);
 
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VectorT<T>
+    operator*=(const VectorT<T>& valvec);
+
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VectorT<T> operator/=(const T val);
 
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T& operator[](size_type pos) &

amr-wind/core/vs/vectorI.H

@@ -24,6 +24,16 @@ VectorT<T>::operator*=(const T fac)
     return *this;
 }
 
+template <typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VectorT<T>
+VectorT<T>::operator*=(const VectorT<T>& vin)
+{
+    vv[0] *= vin.x();
+    vv[1] *= vin.y();
+    vv[2] *= vin.z();
+    return *this;
+}
+
 template <typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE VectorT<T>
 VectorT<T>::operator/=(const T fac)

amr-wind/wind_energy/actuator/ActSrcLineOp.H

@@ -102,6 +102,8 @@ void ActSrcOp<ActTrait, ActSrcLine>::operator()(
     const auto* eps = m_epsilon.data();
     const auto* tmat = m_orientation.data();
 
+    const auto dcoord_flags = m_data.grid().dcoord_flags;
+
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const vs::Vector cc{
             problo[0] + (i + 0.5) * dx[0],
@@ -115,7 +117,10 @@ void ActSrcOp<ActTrait, ActSrcLine>::operator()(
             constexpr amrex::Real wt = 0.5;
             const auto pos_ip = wt * pos[ip] + (1.0 - wt) * opos[ip];
             const auto dist = cc - pos_ip;
-            const auto dist_local = tmat[ip] & dist;
+            // Convert to local (chord, span, thickness) coords
+            const auto dist_local_3D = tmat[ip] & dist;
+            // In local coords, zero disabled distances (e.g., for 2D)
+            const auto dist_local = dist_local_3D * dcoord_flags;
             const auto gauss_fac = utils::gaussian3d(dist_local, eps[ip]);
             const auto& pforce = force[ip];
 

amr-wind/wind_energy/actuator/actuator_types.H

@@ -93,6 +93,9 @@ struct ActGrid
     //! Density at the sampled locations
     RealList density;
 
+    //! Switch for turning off distance components in Gaussian calculation
+    vs::Vector dcoord_flags{1.0, 1.0, 1.0};
+
     /** Helper method to resize the data arrays defined on the grid
      *
      *  \params Number of nodes that contain forcing data

amr-wind/wind_energy/actuator/turbine/turbine_utils.cpp

@@ -21,10 +21,10 @@ void read_inputs(
             "of 'epsilon' or 'epsilon_chord'");
     }
 
-    pp.query("epsilon", tdata.eps_inp);
-    pp.query("epsilon_chord", tdata.eps_chord);
-    pp.query("epsilon_min", tdata.eps_min);
-    pp.query("epsilon_tower", tdata.eps_tower);
+    pp.query_either("epsilon", tdata.eps_inp);
+    pp.query_either("epsilon_chord", tdata.eps_chord);
+    pp.query_either("epsilon_min", tdata.eps_min);
+    pp.query_either("epsilon_tower", tdata.eps_tower);
 
     pp.get("base_position", tinfo.base_pos);
     pp.get("rotor_diameter", tinfo.rotor_diameter);

amr-wind/wind_energy/actuator/wing/ActuatorWing.H

@@ -37,6 +37,27 @@ struct WingBaseData
     //! Pitch angle for the wing
     amrex::Real pitch{0.0};
 
+    //! Switch for turning on 2D Gaussian (off in spanwise)
+    bool gauss_2D{false};
+
+    //! Switch for spanwise normalization in 2D setting
+    bool normalize_2D_spanwise{true};
+
+    //! Switch for turning off force components of actuator
+    amrex::Vector<int> force_coord_flags{1, 1, 1};
+
+    //! User-specified velocity magnitude, only active when nonnegative
+    amrex::Real prescribed_uinf{-1.0};
+
+    //! File name for pitch actuation table
+    std::string pitch_timetable_file;
+
+    //! Pitch actuation: time table
+    amrex::Vector<amrex::Real> time_table;
+
+    //! Pitch actuation: pitch angle table
+    amrex::Vector<amrex::Real> pitch_table;
+
     //! Relative velocity at the actuator node
     VecList vel_rel;
 

amr-wind/wind_energy/actuator/wing/fixed_wing_ops.H

@@ -8,24 +8,6 @@
 
 namespace amr_wind::actuator::ops {
 
-ComponentView make_component_view(FixedWing::DataType& data)
-{
-    auto& grid = data.grid();
-    auto& meta = data.meta();
-    ComponentView view;
-    view.pos = ::amr_wind::utils::slice(grid.pos, 0, meta.num_pts);
-    view.vel_pos = ::amr_wind::utils::slice(grid.vel_pos, 0, meta.num_pts);
-    view.force = ::amr_wind::utils::slice(grid.force, 0, meta.num_pts);
-    view.epsilon = ::amr_wind::utils::slice(grid.epsilon, 0, meta.num_pts);
-    view.orientation =
-        ::amr_wind::utils::slice(grid.orientation, 0, meta.num_pts);
-    view.chord = ::amr_wind::utils::slice(meta.chord, 0, meta.num_pts);
-    view.vel_rel = ::amr_wind::utils::slice(meta.vel_rel, 0, meta.num_pts);
-    view.vel = ::amr_wind::utils::slice(grid.vel, 0, meta.num_pts);
-
-    return view;
-}
-
 template <>
 struct ReadInputsOp<FixedWing, ActSrcLine>
 {
@@ -36,16 +18,56 @@ struct ReadInputsOp<FixedWing, ActSrcLine>
         pp.get("num_points", wdata.num_pts);
         pp.get("start", wdata.start);
         pp.get("end", wdata.end);
-        pp.query("epsilon", wdata.eps_inp);
-        pp.query("epsilon_chord", wdata.epsilon_chord);
-        pp.get("pitch", wdata.pitch);
+        pp.query_either("epsilon", wdata.eps_inp);
+        pp.query_either("epsilon_chord", wdata.epsilon_chord);
+        if (!pp.contains("pitch_timetable")) {
+            // Constant pitch is required without timetable
+            pp.get("pitch", wdata.pitch);
+        }
         pp.get("airfoil_table", wdata.airfoil_file);
         pp.query("airfoil_type", wdata.airfoil_type);
         pp.queryarr("span_locs", wdata.span_locs);
         pp.queryarr("chord", wdata.chord_inp);
         bool use_fllc = false;
         pp.query("fllc", use_fllc);
 
+        // If spanwise components of Gaussian should be ignored
+        pp.query("disable_spanwise_gaussian", wdata.gauss_2D);
+        // If spanwise components of Gaussian and forces should be normalized
+        // (on by default, not used unless gauss_2D = true)
+        pp.query("normalize_spanwise", wdata.normalize_2D_spanwise);
+        // If certain components of force should be ignored
+        // (check that correct number of components are specified)
+        amrex::Vector<int> force_coord_flags_query_inp;
+        pp.queryarr("active_force_dirs", force_coord_flags_query_inp);
+        if (!force_coord_flags_query_inp.empty()) {
+            AMREX_ALWAYS_ASSERT(
+                force_coord_flags_query_inp.size() == AMREX_SPACEDIM);
+            wdata.force_coord_flags = force_coord_flags_query_inp;
+        }
+
+        // If velocity magnitude for force calculation should be prescribed by
+        // user, not measured from the flow
+        pp.query("prescribed_uinf", wdata.prescribed_uinf);
+
+        // Initialize tables for pitch actuation (in degrees)
+        pp.query("pitch_timetable", wdata.pitch_timetable_file);
+        if (!wdata.pitch_timetable_file.empty()) {
+            std::ifstream ifh(wdata.pitch_timetable_file, std::ios::in);
+            if (!ifh.good()) {
+                amrex::Abort(
+                    "Cannot find input file: " + wdata.pitch_timetable_file);
+            }
+            amrex::Real data_time;
+            amrex::Real data_pitch_deg;
+            ifh.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+            while (ifh >> data_time) {
+                ifh >> data_pitch_deg;
+                wdata.time_table.push_back(data_time);
+                wdata.pitch_table.push_back(data_pitch_deg);
+            }
+        }
+
         pp.query("motion_type", wdata.motion_type);
         {
             // do nothing for default, "none"
@@ -85,17 +107,53 @@ struct ReadInputsOp<FixedWing, ActSrcLine>
             amrex::max(max_eps, max_epsc) * max_chord * 3.0;
         const auto& p1 = wdata.start;
         const auto& p2 = wdata.end;
+        // Bounding box quantities used for all motion types
+        amrex::Real minpx = amrex::min(p1.x(), p2.x());
+        amrex::Real maxpx = amrex::max(p1.x(), p2.x());
+        amrex::Real minpy = amrex::min(p1.y(), p2.y());
+        amrex::Real maxpy = amrex::max(p1.y(), p2.y());
+        amrex::Real minpz = amrex::min(p1.z(), p2.z());
+        amrex::Real maxpz = amrex::max(p1.z(), p2.z());
+        // Bounding box limits for 2D actuator
+        if (wdata.gauss_2D) {
+            // Assume largest component of span needs entire domain, but the
+            // others do not need the bounding box to grow
+            const auto wspan = wdata.end - wdata.start;
+            // Flag for checking periodicity in spanwise direction
+            bool warn_per{false};
+            std::string sdir;
+            auto per = data.sim().mesh().Geom()[0].periodicity();
+            if (std::abs(wspan.x()) >=
+                std::max(std::abs(wspan.y()), std::abs(wspan.z()))) {
+                minpx = data.sim().mesh().Geom(0).ProbLoArray()[0];
+                maxpx = data.sim().mesh().Geom(0).ProbHiArray()[0];
+                warn_per = !per.isPeriodic(0);
+                sdir = "x";
+            } else if (std::abs(wspan.y()) >= std::abs(wspan.z())) {
+                minpy = data.sim().mesh().Geom(0).ProbLoArray()[1];
+                maxpy = data.sim().mesh().Geom(0).ProbHiArray()[1];
+                warn_per = !per.isPeriodic(1);
+                sdir = "y";
+            } else {
+                minpz = data.sim().mesh().Geom(0).ProbLoArray()[2];
+                maxpz = data.sim().mesh().Geom(0).ProbHiArray()[2];
+                warn_per = !per.isPeriodic(2);
+                sdir = "z";
+            }
+            if (warn_per) {
+                amrex::Print()
+                    << "\nWARNING: fixed_wing_ops: Detected spanwise direction "
+                    << sdir
+                    << " is not periodic, though 2D Gaussian is being "
+                       "used.\n\n";
+            }
+        }
         // Set up bounding box depending on motion type
         if (amrex::toLower(wdata.motion_type) == "none") {
-            // clang-format off
             info.bound_box = amrex::RealBox(
-                amrex::min(p1.x(), p2.x()) - search_radius,
-                amrex::min(p1.y(), p2.y()) - search_radius,
-                amrex::min(p1.z(), p2.z()) - search_radius,
-                amrex::max(p1.x(), p2.x()) + search_radius,
-                amrex::max(p1.y(), p2.y()) + search_radius,
-                amrex::max(p1.z(), p2.z()) + search_radius);
-            // clang-format on
+                minpx - search_radius, minpy - search_radius,
+                minpz - search_radius, maxpx + search_radius,
+                maxpy - search_radius, maxpz + search_radius);
         } else if (amrex::toLower(wdata.motion_type) == "linear") {
             // Extend bounding box in case of velocity
             constexpr amrex::Real tiny = 1e-8;
@@ -119,22 +177,20 @@ struct ReadInputsOp<FixedWing, ActSrcLine>
             const amrex::Real wm_ext =
                 data.sim().mesh().Geom(0).ProbLoArray()[2];
             info.bound_box = amrex::RealBox(
-                um_ ? um_ext : amrex::min(p1.x(), p2.x()) - search_radius,
-                vm_ ? vm_ext : amrex::min(p1.y(), p2.y()) - search_radius,
-                wm_ ? wm_ext : amrex::min(p1.z(), p2.z()) - search_radius,
-                up_ ? up_ext : amrex::max(p1.x(), p2.x()) + search_radius,
-                vp_ ? vp_ext : amrex::max(p1.y(), p2.y()) + search_radius,
-                wp_ ? wp_ext : amrex::max(p1.z(), p2.z()) + search_radius);
+                um_ ? um_ext : minpx - search_radius,
+                vm_ ? vm_ext : minpy - search_radius,
+                wm_ ? wm_ext : minpz - search_radius,
+                up_ ? up_ext : maxpx + search_radius,
+                vp_ ? vp_ext : maxpy + search_radius,
+                wp_ ? wp_ext : maxpz + search_radius);
         } else if (amrex::toLower(wdata.motion_type) == "sine") {
-            // clang-format off
             info.bound_box = amrex::RealBox(
-                amrex::min(p1.x(), p2.x()) - wdata.s_vector.x() - search_radius,
-                amrex::min(p1.y(), p2.y()) - wdata.s_vector.y() - search_radius,
-                amrex::min(p1.z(), p2.z()) - wdata.s_vector.z() - search_radius,
-                amrex::max(p1.x(), p2.x()) + wdata.s_vector.x() + search_radius,
-                amrex::max(p1.y(), p2.y()) + wdata.s_vector.y() + search_radius,
-                amrex::max(p1.z(), p2.z()) + wdata.s_vector.z() + search_radius);
-            // clang-format on
+                minpx - wdata.s_vector.x() - search_radius,
+                minpy - wdata.s_vector.y() - search_radius,
+                minpz - wdata.s_vector.z() - search_radius,
+                maxpx + wdata.s_vector.x() + search_radius,
+                maxpy + wdata.s_vector.y() + search_radius,
+                maxpz + wdata.s_vector.z() + search_radius);
         }
     }
 };
@@ -159,17 +215,22 @@ struct InitDataOp<FixedWing, ActSrcLine>
             meta.chord.resize(npts);
             ::amr_wind::interp::linear_monotonic(
                 meta.span_locs, meta.chord_inp, wx, meta.chord);
-            // clang-format off
             meta.epsilon_chord = {
                 meta.epsilon_chord.x(), meta.epsilon_chord.z(),
                 meta.epsilon_chord.y()};
-            // clang-format on
+            meta.eps_inp = {
+                meta.eps_inp.x(), meta.eps_inp.z(), meta.eps_inp.y()};
             for (int i = 0; i < npts; ++i) {
                 for (int n = 0; n < AMREX_SPACEDIM; ++n) {
                     const auto eps = meta.epsilon_chord[n] * meta.chord[i];
                     grid.epsilon[i][n] = amrex::max(meta.eps_inp[n], eps);
                 }
             }
+            // Copy Gaussian flags to grid struct if 2D
+            if (meta.gauss_2D) {
+                // Local coords are chord, span, thickness internally
+                grid.dcoord_flags = vs::Vector(1.0, 0.0, 1.0);
+            }
         }
 
         meta.aflookup =

amr-wind/wind_energy/actuator/wing/flat_plate_ops.H

@@ -25,7 +25,7 @@ struct ReadInputsOp<FlatPlate, ActSrcLine>
         pp.get("num_points", wdata.num_pts);
         pp.get("start", wdata.start);
         pp.get("end", wdata.end);
-        pp.get("epsilon", wdata.eps_inp);
+        pp.get_either("epsilon", wdata.eps_inp);
         pp.get("pitch", wdata.pitch);
 
         pp.query("chord", wdata.chord_inp);