Actuator syntax changes (Exawind#635)

* Remove unused file

* Remove redundant point specification methods

* Default Jouwkowsky num blades to 3

* Output TSR

* Out put Joukowksy forces to netcdf

* Change default root correction constants

* Fix typo

* Fix unit test

* Add backwards compatibility checks

Also add error_stream checks outside the base op

* Format

amr-wind/wind_energy/actuator/disk/Joukowsky.H

@@ -13,14 +13,17 @@ struct JoukowskyData : public DiskBaseData
     RealList angular_velocity;
     RealList tip_correction;
     RealList root_correction;
+    RealList f_normal;
+    RealList f_theta;
     amrex::Real current_angular_velocity{0.0};
+    amrex::Real current_tip_speed_ratio{0.0};
     amrex::Real vortex_core_size;
     amrex::Real current_cp;
     // --- Sorenson 2020 equation 10 constants ----
-    amrex::Real root_correction_coefficient{1.256};
-    amrex::Real root_correction_exponent{2.0};
+    amrex::Real root_correction_coefficient{2.335};
+    amrex::Real root_correction_exponent{4.0};
     // --------------------------------------------
-    int num_blades;
+    int num_blades{3};
     vs::Vector disk_force{0.0, 0.0, 0.0};
     bool use_tip_correction{true};
     bool use_root_correction{true};

amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H

@@ -58,6 +58,8 @@ struct InitDataOp<Joukowsky, ActSrcDisk>
         grid.vel_pos.resize(meta.num_vel_pts);
         meta.tip_correction.resize(meta.num_vel_pts_r);
         meta.root_correction.resize(meta.num_vel_pts_r);
+        meta.f_normal.resize(meta.num_vel_pts_r);
+        meta.f_theta.resize(meta.num_force_pts);
 
         std::fill(meta.tip_correction.begin(), meta.tip_correction.end(), 1.0);
 
@@ -130,6 +132,9 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
 
         const amrex::Real lambda =
             ddata.radius() * ddata.current_angular_velocity / U_ref;
+
+        ddata.current_tip_speed_ratio = lambda;
+
         const amrex::Real lambda_2 = lambda * lambda;
 
         // step 0: compute tip correction based on current wind speed (eq 9)
@@ -191,6 +196,8 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
             const amrex::Real f_z =
                 q0 * g * F / x * (lambda * x + 0.5 * q0 * g * F / x);
 
+            ddata.f_normal[i] = f_z;
+
             for (int j = 0; j < ddata.num_vel_pts_t; j++, ip++) {
                 const auto point_current = grid.pos[ip];
 
@@ -204,6 +211,7 @@ struct ComputeForceOp<Joukowsky, ActSrcDisk>
                     -ddata.normal_vec & disk_velocity[ip];
 
                 const amrex::Real f_theta = u_disk_ij / U_ref * q0 * g * F / x;
+                ddata.f_theta[ip] = f_theta;
 
                 const amrex::Real dArea =
                     geometry_factor * (std::pow(i + 1.0, 2) - std::pow(i, 2));

amr-wind/wind_energy/actuator/disk/Joukowsky_ops.cpp

@@ -1,6 +1,7 @@
 #include "amr-wind/wind_energy/actuator/disk/Joukowsky_ops.H"
 #include "amr-wind/utilities/ncutils/nc_interface.H"
 #include "amr-wind/utilities/io_utils.H"
+#include "amr-wind/wind_energy/actuator/disk/disk_ops.H"
 
 namespace amr_wind {
 namespace actuator {
@@ -13,6 +14,7 @@ void check_for_parse_conflicts(const utils::ActParser& pp)
     base::collect_parse_dependencies(pp, "thrust_coeff", "angular_velocity", error_collector);
     base::collect_parse_dependencies(pp, "use_root_correction", "vortex_core_size", error_collector);
     // clang-format on
+    ops::base::check_error_stream(error_collector);
 }
 
 void optional_parameters(JoukowskyData& meta, const utils::ActParser& pp)
@@ -23,6 +25,7 @@ void optional_parameters(JoukowskyData& meta, const utils::ActParser& pp)
     pp.query("vortex_core_size", meta.vortex_core_size);
     pp.query("root_correction_exponent", meta.root_correction_exponent);
     pp.query("root_correction_coefficient", meta.root_correction_coefficient);
+    pp.query("num_blades", meta.num_blades);
 }
 
 void required_parameters(JoukowskyData& meta, const utils::ActParser& pp)
@@ -31,7 +34,6 @@ void required_parameters(JoukowskyData& meta, const utils::ActParser& pp)
     pp.get("num_points_t", meta.num_vel_pts_t);
     pp.get("num_points_r", meta.num_vel_pts_r);
     pp.getarr("angular_velocity", meta.angular_velocity);
-    pp.get("num_blades", meta.num_blades);
     meta.num_force_pts = meta.num_vel_pts_r * meta.num_vel_pts_t;
     meta.num_vel_pts = meta.num_force_pts * 2;
     meta.dr = 0.5 * meta.diameter / meta.num_vel_pts_r;
@@ -59,13 +61,15 @@ void prepare_netcdf_file(
     const std::string nt_name = "num_time_steps";
     const std::string np_name = "num_actuator_points";
     const std::string nv_name = "num_velocity_points";
+    const std::string nr_name = "num_radial_points";
 
     ncf.enter_def_mode();
     ncf.put_attr("title", "AMR-Wind ActuatorDisk actuator output");
     ncf.put_attr("version", ioutils::amr_wind_version());
     ncf.put_attr("created_on", ioutils::timestamp());
     ncf.def_dim(nt_name, NC_UNLIMITED);
     ncf.def_dim("ndim", AMREX_SPACEDIM);
+    ncf.def_dim(nr_name, data.num_vel_pts_r);
 
     auto grp = ncf.def_group(info.label);
     // clang-format off
@@ -84,11 +88,14 @@ void prepare_netcdf_file(
     grp.def_var("xyz_v", NC_DOUBLE, {nv_name, "ndim"});
     grp.def_var("vref", NC_DOUBLE, {nt_name, "ndim"});
     grp.def_var("vdisk", NC_DOUBLE, {nt_name, "ndim"});
+    grp.def_var("tsr", NC_DOUBLE, {nt_name});
     grp.def_var("ct", NC_DOUBLE, {nt_name});
     grp.def_var("cp", NC_DOUBLE, {nt_name});
     grp.def_var("density", NC_DOUBLE, {nt_name});
     grp.def_var("total_disk_force", NC_DOUBLE, {nt_name, "ndim"});
     grp.def_var("angular_velocity", NC_DOUBLE, {nt_name});
+    grp.def_var("f_normal", NC_DOUBLE, {nt_name, nr_name});
+    grp.def_var("f_theta", NC_DOUBLE, {nt_name, np_name});
     ncf.exit_def_mode();
 
     {
@@ -124,18 +131,23 @@ void write_netcdf(
     const std::string nt_name = "num_time_steps";
     // Index of next timestep
     const size_t nt = ncf.dim(nt_name).len();
+    const size_t nr = data.num_vel_pts_r;
+    const size_t np = data.num_force_pts;
     auto grp = ncf.group(info.label);
     grp.var("time").put(&time, {nt}, {1});
     grp.var("vref").put(
         &(data.reference_velocity[0]), {nt, 0}, {1, AMREX_SPACEDIM});
     grp.var("vdisk").put(
         &(data.mean_disk_velocity[0]), {nt, 0}, {1, AMREX_SPACEDIM});
+    grp.var("tsr").put(&data.current_tip_speed_ratio, {nt}, {1});
     grp.var("ct").put(&data.current_ct, {nt}, {1});
     grp.var("cp").put(&data.current_cp, {nt}, {1});
     grp.var("density").put(&data.density, {nt}, {1});
     grp.var("total_disk_force")
         .put(&data.disk_force[0], {nt}, {1, AMREX_SPACEDIM});
     grp.var("angular_velocity").put(&data.current_angular_velocity, {nt}, {1});
+    grp.var("f_normal").put(&(data.f_normal[0]), {nt, 0}, {1, nr});
+    grp.var("f_theta").put(&(data.f_theta[0]), {nt, 0}, {1, np});
 #else
     amrex::ignore_unused(name, data, info, time);
 #endif

amr-wind/wind_energy/actuator/disk/disk_ops.H

@@ -33,12 +33,19 @@ void collect_parse_conflicts(
     const std::string& p1,
     const std::string& p2,
     std::ostringstream& ss);
+void collect_parse_dependencies_one_way(
+    const utils::ActParser& pp,
+    const std::string& independent,
+    const std::string& dependent,
+    std::ostringstream& ss);
 void collect_parse_dependencies(
     const utils::ActParser& pp,
     const std::string& p1,
     const std::string& p2,
     std::ostringstream& ss);
 
+void check_error_stream(std::ostringstream& error_collector);
+
 void required_parameters(DiskBaseData& meta, const utils::ActParser& pp);
 
 void optional_parameters(DiskBaseData& meta, const utils::ActParser& pp);

amr-wind/wind_energy/actuator/disk/disk_ops.cpp

@@ -130,20 +130,27 @@ void collect_parse_conflicts(
         ss << "ActuatorDisk Conflict: " << p1 << " and " << p2 << std::endl;
     }
 }
+
+void collect_parse_dependencies_one_way(
+    const utils::ActParser& pp,
+    const std::string& independent,
+    const std::string& dependent,
+    std::ostringstream& ss)
+{
+    if (pp.contains(dependent) && !pp.contains(independent)) {
+        ss << "ActuatorDisk Dependency Missing: " << independent
+           << " required with " << dependent << std::endl;
+    }
+}
+
 void collect_parse_dependencies(
     const utils::ActParser& pp,
     const std::string& p1,
     const std::string& p2,
     std::ostringstream& ss)
 {
-    if (pp.contains(p1) && !pp.contains(p2)) {
-        ss << "ActuatorDisk Dependency Missing: " << p2 << " required with "
-           << p1 << std::endl;
-    }
-    if (!pp.contains(p1) && pp.contains(p2)) {
-        ss << "ActuatorDisk Dependency Missing: " << p1 << " required with "
-           << p2 << std::endl;
-    }
+    collect_parse_dependencies_one_way(pp, p1, p2, ss);
+    collect_parse_dependencies_one_way(pp, p2, p1, ss);
 }
 
 void required_parameters(DiskBaseData& meta, const utils::ActParser& pp)
@@ -175,7 +182,7 @@ void optional_parameters(DiskBaseData& meta, const utils::ActParser& pp)
     pp.query("disk_normal", meta.normal_vec);
     pp.query("density", meta.density);
     pp.query("diameters_to_sample", meta.diameters_to_sample);
-    pp.query("num_theta_force_points", meta.num_force_theta_pts);
+    pp.query("num_points_t", meta.num_force_theta_pts);
 
     // make sure we compute normal vec contribution from tilt before yaw
     // since we won't know a reference axis to rotate for tilt after
@@ -238,6 +245,14 @@ void optional_parameters(DiskBaseData& meta, const utils::ActParser& pp)
     }
 }
 
+void check_error_stream(std::ostringstream& error_collector)
+{
+    if (!error_collector.str().empty()) {
+        amrex::Abort(
+            "Errors found while parsing ActuatorDisk Inputs:\n" +
+            error_collector.str());
+    }
+}
 std::ostringstream check_for_parse_conflicts(const utils::ActParser& pp)
 {
     std::ostringstream error_collector;
@@ -273,11 +288,7 @@ std::ostringstream check_for_parse_conflicts(const utils::ActParser& pp)
         }
     }
 
-    if (!error_collector.str().empty()) {
-        amrex::Abort(
-            "Errors found while parsing ActuatorDisk Inputs:\n" +
-            error_collector.str());
-    }
+    check_error_stream(error_collector);
     return error_collector;
 }
 
@@ -307,8 +318,7 @@ void compute_and_normalize_coplanar_vector(DiskBaseData& meta)
 void final_checks(const DiskBaseData& meta)
 {
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-        meta.num_vel_pts > 0,
-        "num_vel_points_r and num_vel_points_t must both be >=1");
+        meta.num_vel_pts > 0, "num_points_r and num_points_t must both be >=1");
 }
 
 amrex::RealBox compute_bounding_box(const DiskBaseData& meta)

amr-wind/wind_energy/actuator/disk/disk_spreading.H

@@ -2,7 +2,6 @@
 #define DISK_SPREADING_H_
 
 #include "amr-wind/wind_energy/actuator/actuator_utils.H"
-#include "amr-wind/wind_energy/actuator/disk/spreading_kernels.h"
 #include "amr-wind/wind_energy/actuator/disk/UniformCt.H"
 #include "amr-wind/core/FieldRepo.H"
 

amr-wind/wind_energy/actuator/disk/uniform_ct_ops.cpp

@@ -1,29 +1,57 @@
 #include "amr-wind/wind_energy/actuator/disk/uniform_ct_ops.H"
+#include "amr-wind/wind_energy/actuator/ActParser.H"
+#include <ostream>
 
 namespace amr_wind {
 namespace actuator {
 namespace ops {
 namespace uniformct {
+
+void check_for_removed_syntax(
+    const utils::ActParser& pp, std::ostringstream& stream)
+{
+    const std::string error_start = "ERROR ActuatorDisk:: ";
+    if (pp.contains("num_vel_points_r")) {
+        stream << error_start
+               << "'num_vel_points_r' has been replaced with 'num_points_r'\n";
+    }
+    if (pp.contains("num_vel_points_t")) {
+        stream << error_start
+               << "'num_vel_points_t' has been replaced with 'num_points_t'\n";
+    }
+    if (pp.contains("num_force_points")) {
+        stream << error_start
+               << "'num_force_points' has been repalced with 'num_points_r'\n";
+    }
+    if (pp.contains("num_theta_force_points")) {
+        stream << error_start
+               << "'num_theta_force_points' has been replaced with "
+                  "'num_points_t'\n";
+    }
+}
+
 void check_for_parse_conflicts(const utils::ActParser& pp)
 {
     auto error_collector = ops::base::check_for_parse_conflicts(pp);
-    ops::base::collect_parse_dependencies(
-        pp, "num_theta_force_points", "spreading", error_collector);
+    check_for_removed_syntax(pp, error_collector);
+    ops::base::collect_parse_dependencies_one_way(
+        pp, "num_points_t", "spreading", error_collector);
+    ops::base::check_error_stream(error_collector);
 }
 void optional_parameters(DiskBaseData& meta, const utils::ActParser& pp)
 {
     ops::base::optional_parameters(meta, pp);
     // params for the sampling disk
-    pp.query("num_vel_points_r", meta.num_vel_pts_r);
-    pp.query("num_vel_points_t", meta.num_vel_pts_t);
+    pp.query("num_points_t", meta.num_vel_pts_t);
     // 2x 1 for sampling up stream and one for sampling at the disk
     meta.num_vel_pts = meta.num_vel_pts_r * meta.num_vel_pts_t * 2;
     pp.query("spreading_type", meta.spreading_type);
 }
 void required_parameters(DiskBaseData& meta, const utils::ActParser& pp)
 {
     ops::base::required_parameters(meta, pp);
-    pp.get("num_force_points", meta.num_force_pts);
+    pp.get("num_points_r", meta.num_force_pts);
+    meta.num_vel_pts_r = meta.num_force_pts;
     meta.dr = 0.5 * meta.diameter / meta.num_force_pts;
 }
 void parse_and_gather_params(const utils::ActParser& pp, UniformCtData& data)

test/test_files/joukowsky_disk/joukowsky_disk.i

@@ -28,7 +28,6 @@ Actuator.type = JoukowskyDisk
 
 Actuator.JoukowskyDisk.rotor_diameter = 126.0
 Actuator.JoukowskyDisk.disk_center = 0.0 0.0 0.0
-Actuator.JoukowskyDisk.num_blades = 3
 Actuator.JoukowskyDisk.yaw = 270.0 # degrees (yaw is relative to north which defaults to {0,1,0})
 Actuator.JoukowskyDisk.sample_yaw = 270.0 # set velocity sampling to be in the normal flow direction
 Actuator.JoukowskyDisk.thrust_coeff = 0.0 0.7 1.2

test/test_files/uniform_ct_disk/uniform_ct_disk.i

@@ -31,14 +31,13 @@ Actuator.UniformCtDisk.base_position = 0.0 0.0 0.0
 Actuator.UniformCtDisk.hub_height = 0.0
 Actuator.UniformCtDisk.yaw = 315.0 # degrees (yaw is relative to north which defaults to {0,1,0})
 Actuator.UniformCtDisk.sample_yaw = 270.0 # set velocity sampling to be in the normal flow direction
-Actuator.UniformCtDisk.num_force_points = 5
 Actuator.UniformCtDisk.thrust_coeff = 0.0 0.7 1.2
 Actuator.UniformCtDisk.wind_speed = 0.0 10.0 12.0
 Actuator.UniformCtDisk.epsilon = 10.0
 Actuator.UniformCtDisk.density = 1.225
 Actuator.UniformCtDisk.diameters_to_sample = 1.0
-Actuator.UniformCtDisk.num_vel_points_r = 3
-Actuator.UniformCtDisk.num_vel_points_t = 3
+Actuator.UniformCtDisk.num_points_r = 5
+Actuator.UniformCtDisk.num_points_t = 3
 
 ICNS.source_terms = ActuatorForcing
 

test/test_files/uniform_ct_disk_gaussian/uniform_ct_disk_gaussian.i

@@ -32,15 +32,13 @@ Actuator.UniformCtDisk.base_position = 0.0 0.0 0.0
 Actuator.UniformCtDisk.hub_height = 0.0
 Actuator.UniformCtDisk.yaw = 315.0 # degrees (yaw is relative to north which defaults to {0,1,0})
 Actuator.UniformCtDisk.sample_yaw = 270.0 # set velocity sampling to be in the normal flow direction
-Actuator.UniformCtDisk.num_force_points = 5
 Actuator.UniformCtDisk.thrust_coeff = 0.0 0.7 1.2
 Actuator.UniformCtDisk.wind_speed = 0.0 10.0 12.0
 Actuator.UniformCtDisk.epsilon = 10.0
 Actuator.UniformCtDisk.density = 1.225
 Actuator.UniformCtDisk.diameters_to_sample = 1.0
-Actuator.UniformCtDisk.num_vel_points_r = 3
-Actuator.UniformCtDisk.num_vel_points_t = 3
-Actuator.UniformCtDisk.num_theta_force_points = 20
+Actuator.UniformCtDisk.num_points_r = 5
+Actuator.UniformCtDisk.num_points_t = 20
 
 ICNS.source_terms = ActuatorForcing
 

unit_tests/wind_energy/actuator/test_disk_uniform_ct.cpp

@@ -15,7 +15,7 @@ class UniformCtTest : public AmrexTest
         AmrexTest::SetUp();
         {
             amrex::ParmParse pu("Actuator.UniformCtDisk");
-            pu.add("num_force_points", 3);
+            pu.add("num_points_r", 3);
             pu.add("epsilon", 1);
             pu.add("rotor_diameter", 1);
             std::vector<double> ct{1};