Reading in table for rpm and blade pitch

Source/DataStructs/ERF_DataStruct.H

@@ -385,6 +385,7 @@ struct SolverChoice {
         pp.query("windfarm_spec_table", windfarm_spec_table);
         pp.query("windfarm_blade_table", windfarm_blade_table);
         pp.query("windfarm_airfoil_tables", windfarm_airfoil_tables);
+        pp.query("windfarm_spec_table_extra", windfarm_spec_table_extra);
 
         // Sampling distance upstream of the turbine to find the
         // incoming free stream velocity as a factor of the diameter of the
@@ -671,7 +672,7 @@ struct SolverChoice {
     // if     SAM, then it will be set to RhoQ4
     int RhoQr_comp {-1};
 
-    std::string windfarm_loc_table, windfarm_spec_table;
+    std::string windfarm_loc_table, windfarm_spec_table, windfarm_spec_table_extra;
     std::string windfarm_blade_table, windfarm_airfoil_tables;
     amrex::Real sampling_distance_by_D = -1.0;
     amrex::Real turb_disk_angle = -1.0;

Source/Initialization/ERF_init_windfarm.cpp

@@ -46,6 +46,7 @@ ERF::init_windfarm (int lev)
         windfarm->read_windfarm_blade_table(solverChoice.windfarm_blade_table);
         windfarm->read_windfarm_airfoil_tables(solverChoice.windfarm_airfoil_tables,
                                                solverChoice.windfarm_blade_table);
+        windfarm->read_windfarm_spec_table_extra(solverChoice.windfarm_spec_table_extra);
     }
 }
 

Source/WindFarmParametrization/ERF_InitWindFarm.cpp

@@ -125,6 +125,14 @@ WindFarm::init_windfarm_lat_lon (const std::string windfarm_loc_table,
         xloc[it] = xloc[it] - xloc_min + windfarm_x_shift;
         yloc[it] = yloc[it] - yloc_min + windfarm_y_shift;
     }
+
+    FILE* file_xy_loc;
+    file_xy_loc = fopen("file_xy_loc_KingPlains.txt","w");
+
+    for(int it = 0;it<xloc.size(); it++){
+        fprintf(file_xy_loc,"%0.15g %0.15g %0.15g\n", xloc[it], yloc[it], 89.0);
+    }
+    fclose(file_xy_loc);
 }
 
 void
@@ -230,6 +238,45 @@ WindFarm::read_windfarm_blade_table(const std::string windfarm_blade_table)
     }
 }
 
+void
+WindFarm::read_windfarm_spec_table_extra(const std::string windfarm_spec_table_extra)
+{
+    // Open the file
+    std::ifstream file(windfarm_spec_table_extra);
+
+    // Check if file opened successfully
+    if (!file.is_open()) {
+        Abort("Error: You are using generalized wind farms option. This requires an input file erf.windfarm_spec_table_extra."
+              " Either this entry is missing in the inputs or the file specified -" + windfarm_spec_table_extra + " does"
+              " not exist. Exiting...");
+    } else {
+        printf("Reading in windfarm_spec_table_extra %s", windfarm_spec_table_extra.c_str());
+    }
+
+    // Ignore the first line (header)
+    std::string header;
+    std::getline(file, header);
+
+    // Variables to hold each row's values
+    double V, Cp, Ct, rpm, pitch, temp;
+
+    // Read the file row by row
+    while (file >> V) {
+        char comma;  // To ignore the commas
+        file >> comma >> Cp >> comma >> Ct >> comma >> temp >> comma >> temp >> comma
+             >> temp >> comma >> rpm >> comma >> pitch >> comma >> temp;
+
+        velocity.push_back(V);
+        C_P.push_back(Cp);
+        C_T.push_back(Ct);
+        rotor_RPM.push_back(rpm);
+        blade_pitch.push_back(pitch);
+    }
+
+    set_turb_spec_extra(velocity, C_P, C_T, rotor_RPM, blade_pitch);
+}
+
+
 void
 WindFarm::read_windfarm_airfoil_tables(const std::string windfarm_airfoil_tables,
                                        const std::string windfarm_blade_table)

Source/WindFarmParametrization/ERF_WindFarm.H

@@ -68,6 +68,8 @@ public:
     void read_windfarm_airfoil_tables(const std::string windfarm_airfoil_tables,
                                       const std::string windfarm_blade_table);
 
+    void read_windfarm_spec_table_extra(const std::string windfarm_spec_table_extra);
+
     void fill_Nturb_multifab(const amrex::Geometry& geom,
                              amrex::MultiFab& mf_Nturb);
 
@@ -128,6 +130,15 @@ public:
         m_windfarm_model[0]->set_blade_airfoil_spec(a_bld_airfoil_aoa, a_bld_airfoil_Cl, a_bld_airfoil_Cd);
     }
 
+    void set_turb_spec_extra (const amrex::Vector<amrex::Real>& a_velocity,
+                              const amrex::Vector<amrex::Real>& a_C_P,
+                              const amrex::Vector<amrex::Real>& a_C_T,
+                              const amrex::Vector<amrex::Real>& a_rotor_RPM,
+                              const amrex::Vector<amrex::Real>& a_blade_pitch) override
+    {
+        m_windfarm_model[0]->set_turb_spec_extra(a_velocity, a_C_P, a_C_T, a_rotor_RPM, a_blade_pitch);
+    }
+
 protected:
 
     amrex::Vector<amrex::Real> xloc, yloc;
@@ -137,6 +148,7 @@ protected:
     amrex::Vector<amrex::Real> bld_rad_loc, bld_twist, bld_chord;
     amrex::Vector<amrex::Vector<amrex::Real>> bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd;
     int n_bld_sections;
+    amrex::Vector<amrex::Real> velocity, C_P, C_T, rotor_RPM, blade_pitch;
 
    /*! \brief Create and set the specified windfarm model */
     template<class NewWindFarmModel>

Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp

@@ -181,9 +181,17 @@ compute_source_terms_Fn_Ft (const Real rad,
                             const Real* bld_airfoil_aoa,
                             const Real* bld_airfoil_Cl,
                             const Real* bld_airfoil_Cd,
-                            const int n_pts_airfoil)
+                            const int n_pts_airfoil,
+                            const Real* velocity,
+                            const Real* rotor_RPM,
+                            const Real* blade_pitch,
+                            const int n_spec_extra)
 {
-    Real Omega = 9.0/60.0*2.0*PI;
+
+    Real rpm = interpolate_1d(velocity, rotor_RPM, avg_vel, n_spec_extra);
+    Real pitch = interpolate_1d(velocity, blade_pitch, avg_vel, n_spec_extra);
+
+    Real Omega = rpm/60.0*2.0*PI;
     Real rho = 1.226;
 
     Real B = 3.0;
@@ -193,6 +201,8 @@ compute_source_terms_Fn_Ft (const Real rad,
     Real twist = interpolate_1d(bld_rad_loc, bld_twist, rad, n_bld_sections);
     Real c = interpolate_1d(bld_rad_loc, bld_chord, rad, n_bld_sections);
 
+    printf("The avg_vel, rpm, pitch are %0.15g %0.15g %0.15g\n", avg_vel, rpm, pitch);
+
     // Iteration procedure
 
     Real s = 0.5*c*B/(PI*rad);
@@ -289,10 +299,7 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
 
     get_blade_airfoil_spec(bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd);
 
-    for(int i=0;i<bld_rad_loc.size();i++) {
-        //std::cout << "I am here ...." << bld_rad_loc[i] << " " << bld_twist[i] << " " << bld_chord[i] << "\n";
-    }
-
+    get_turb_spec_extra(velocity, C_P, C_T, rotor_RPM, blade_pitch);
 
     Real d_hub_height = hub_height;
     Real d_rotor_rad = rotor_rad;
@@ -383,6 +390,20 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
     auto d_bld_airfoil_Cl_ptr  = Cl.data();
     auto d_bld_airfoil_Cd_ptr  = Cd.data();
 
+    int n_spec_extra = velocity.size();
+
+    Gpu::DeviceVector<Real> d_velocity(n_spec_extra);
+    Gpu::DeviceVector<Real> d_rotor_RPM(n_spec_extra);
+    Gpu::DeviceVector<Real> d_blade_pitch(n_spec_extra);
+
+    Gpu::copy(Gpu::hostToDevice, velocity.begin(), velocity.end(), d_velocity.begin());
+    Gpu::copy(Gpu::hostToDevice, rotor_RPM.begin(), rotor_RPM.end(), d_rotor_RPM.begin());
+    Gpu::copy(Gpu::hostToDevice, blade_pitch.begin(), blade_pitch.end(), d_blade_pitch.begin());
+
+    auto d_velocity_ptr = d_velocity.data();
+    auto d_rotor_RPM_ptr = d_rotor_RPM.data();
+    auto d_blade_pitch_ptr = d_blade_pitch.data();
+
     for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
 
         const Box& gbx      = mfi.growntilebox(1);
@@ -443,7 +464,11 @@ GeneralAD::source_terms_cellcentered (const Geometry& geom,
                                                                d_bld_airfoil_aoa_ptr[index],
                                                                d_bld_airfoil_Cl_ptr[index],
                                                                d_bld_airfoil_Cd_ptr[index],
-                                                               n_pts_airfoil);
+                                                               n_pts_airfoil,
+                                                               d_velocity_ptr,
+                                                               d_rotor_RPM_ptr,
+                                                               d_blade_pitch_ptr,
+                                                               n_spec_extra);
 
                         Real Fn = Fn_and_Ft[0];
                         Real Ft = Fn_and_Ft[1];

Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H

@@ -48,6 +48,7 @@ protected:
     amrex::Vector<amrex::Real> freestream_velocity, freestream_phi, disk_cell_count;
     amrex::Vector<amrex::Real> bld_rad_loc, bld_twist, bld_chord;
     amrex::Vector<amrex::Vector<amrex::Real>> bld_airfoil_aoa, bld_airfoil_Cl, bld_airfoil_Cd;
+    amrex::Vector<amrex::Real> velocity, C_P, C_T, rotor_RPM, blade_pitch;
 };
 
 #endif

Source/WindFarmParametrization/Null/ERF_NullWindFarm.H

@@ -67,6 +67,19 @@ public:
         m_bld_airfoil_Cd = bld_airfoil_Cd;
     }
 
+    virtual void set_turb_spec_extra(const amrex::Vector<amrex::Real>& velocity,
+                              const amrex::Vector<amrex::Real>& C_P,
+                              const amrex::Vector<amrex::Real>& C_T,
+                              const amrex::Vector<amrex::Real>& rotor_RPM,
+                              const amrex::Vector<amrex::Real>& blade_pitch)
+    {
+        m_velocity = velocity;
+        m_C_P = C_P;
+        m_C_T = C_T;
+        m_rotor_RPM = rotor_RPM;
+        m_blade_pitch = blade_pitch;
+    }
+
     void get_turb_spec (amrex::Real& rotor_rad, amrex::Real& hub_height,
                         amrex::Real& thrust_coeff_standing, amrex::Vector<amrex::Real>& wind_speed,
                         amrex::Vector<amrex::Real>& thrust_coeff, amrex::Vector<amrex::Real>& power)
@@ -109,6 +122,19 @@ public:
         bld_airfoil_Cd = m_bld_airfoil_Cd;
     }
 
+    void get_turb_spec_extra(amrex::Vector<amrex::Real>& velocity,
+                             amrex::Vector<amrex::Real>& C_P,
+                                amrex::Vector<amrex::Real>& C_T,
+                             amrex::Vector<amrex::Real>& rotor_RPM,
+                             amrex::Vector<amrex::Real>& blade_pitch)
+    {
+        velocity = m_velocity;
+        C_P = m_C_P;
+        C_T = m_C_T;
+        rotor_RPM = m_rotor_RPM;
+        blade_pitch = m_blade_pitch;
+    }
+
 static AMREX_GPU_DEVICE
 bool find_if_marked(amrex::Real x1, amrex::Real x2, amrex::Real y1, amrex::Real y2,
                     amrex::Real x0, amrex::Real y0, amrex::Real nx, amrex::Real ny,
@@ -158,6 +184,7 @@ protected:
     amrex::Vector<amrex::Real> m_wind_speed, m_thrust_coeff, m_power;
     amrex::Vector<amrex::Real> m_bld_rad_loc, m_bld_twist, m_bld_chord;
     amrex::Vector<amrex::Vector<amrex::Real>> m_bld_airfoil_aoa, m_bld_airfoil_Cl, m_bld_airfoil_Cd;
+    amrex::Vector<amrex::Real> m_velocity, m_C_P, m_C_T, m_rotor_RPM, m_blade_pitch;
 };