Some corrections to wind models (erf-model#1839)

* Correcting compilation error with TINY_PROFILE=TRUE * Changes to Make.ERF * Merging with dev * Getting ready for generalized actuator disk * Framework for reading in blade data * Removing print * Minor changes to error hooks and docs * Correcting errors in turbine positioning * Removing unwanted prints * Correcting error hook * Minor correction * Correcting CMake errors --------- Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]>
AMLattanzi · Oct 2, 2024 · 84ef97e · 84ef97e
1 parent c52ddce
commit 84ef97e
Show file tree

Hide file tree

Showing 16 changed files with 323 additions and 35 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -196,6 +196,7 @@ function(build_erf_lib erf_lib_name)
 	   ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
 	   ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp
 	   ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
+	   ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp
        ${SRC_DIR}/LandSurfaceModel/SLM/ERF_SLM.cpp
        ${SRC_DIR}/LandSurfaceModel/MM5/ERF_MM5.cpp
   )
@@ -251,6 +252,7 @@ endif()
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/WindFarmParametrization/Fitch>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/WindFarmParametrization/EWP>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/WindFarmParametrization/SimpleActuatorDisk>)
+  target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/WindFarmParametrization/GeneralActuatorDisk>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/LandSurfaceModel>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/LandSurfaceModel/Null>)
   target_include_directories(${erf_lib_name} PUBLIC $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/Source/LandSurfaceModel/SLM>)

diff --git a/Docs/sphinx_doc/theory/WindFarmModels.rst b/Docs/sphinx_doc/theory/WindFarmModels.rst
@@ -4,7 +4,7 @@ Wind farm models
 Introduction
 -------------
 
-ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) and Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) are supported.
+ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently the Fitch model (`Fitch et al. 2012`_), Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) and Simplified actuator disk model (See Section 3.2 in `Wind Energy Handbook 2nd edition`_) are supported.
 
 .. _Fitch model:
 

diff --git a/Exec/Make.ERF b/Exec/Make.ERF
@@ -152,12 +152,14 @@ ERF_WINDFARM_NULL_DIR   = $(ERF_WINDFARM_DIR)/Null
 ERF_WINDFARM_FITCH_DIR = $(ERF_WINDFARM_DIR)/Fitch
 ERF_WINDFARM_EWP_DIR   = $(ERF_WINDFARM_DIR)/EWP
 ERF_WINDFARM_SIMPLEAD_DIR   = $(ERF_WINDFARM_DIR)/SimpleActuatorDisk
+ERF_WINDFARM_GENERALAD_DIR   = $(ERF_WINDFARM_DIR)/GeneralActuatorDisk
 
 include $(ERF_WINDFARM_DIR)/Make.package
 include $(ERF_WINDFARM_NULL_DIR)/Make.package
 include $(ERF_WINDFARM_FITCH_DIR)/Make.package
 include $(ERF_WINDFARM_EWP_DIR)/Make.package
 include $(ERF_WINDFARM_SIMPLEAD_DIR)/Make.package
+include $(ERF_WINDFARM_GENERALAD_DIR)/Make.package
 
 VPATH_LOCATIONS   += $(ERF_WINDFARM_DIR)
 INCLUDE_LOCATIONS += $(ERF_WINDFARM_DIR)
@@ -173,6 +175,10 @@ INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR)
 
 VPATH_LOCATIONS   += $(ERF_WINDFARM_SIMPLEAD_DIR)
 INCLUDE_LOCATIONS += $(ERF_WINDFARM_SIMPLEAD_DIR)
+
+VPATH_LOCATIONS   += $(ERF_WINDFARM_GENERALAD_DIR)
+INCLUDE_LOCATIONS += $(ERF_WINDFARM_GENERALAD_DIR)
+
 endif
 
 ifeq ($(USE_HEFFTE),TRUE)

diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
@@ -38,7 +38,7 @@ enum struct MoistureType {
 };
 
 enum struct WindFarmType {
-    Fitch, EWP, SimpleAD, None
+    Fitch, EWP, SimpleAD, GeneralAD, None
 };
 
 enum struct WindFarmLocType{
@@ -359,9 +359,12 @@ struct SolverChoice {
         else if (windfarm_type_string == "SimpleActuatorDisk") {
             windfarm_type = WindFarmType::SimpleAD;
         }
+        else if (windfarm_type_string == "GeneralActuatorDisk") {
+            windfarm_type = WindFarmType::GeneralAD;
+        }
         else if (windfarm_type_string != "None") {
             amrex::Abort("Are you using windfarms? Dont know this windfarm_type. windfarm_type"
-                         " has to be Fitch or EWP or None.");
+                         " has to be Fitch, EWP, SimpleActuatorDisk, GeneralActuatorDisk or None.");
         }
 
         static std::string windfarm_loc_type_string = "None";
@@ -380,6 +383,8 @@ struct SolverChoice {
 
         pp.query("windfarm_loc_table",  windfarm_loc_table);
         pp.query("windfarm_spec_table", windfarm_spec_table);
+        pp.query("windfarm_blade_table", windfarm_blade_table);
+        pp.query("windfarm_airofil_tables", windfarm_airfoil_tables);
 
         // Sampling distance upstream of the turbine to find the
         // incoming free stream velocity as a factor of the diameter of the
@@ -424,7 +429,7 @@ struct SolverChoice {
                            " distance as a factor of the turbine diameter at which the incoming free stream"
                            " velocity will be computed at.");
         }
-        if (windfarm_type==WindFarmType::SimpleAD and turb_disk_angle < 0.0) {
+        if((windfarm_type==WindFarmType::SimpleAD or windfarm_type==WindFarmType::GeneralAD) and turb_disk_angle < 0.0) {
             amrex::Abort("To use simplified actuator disks, you need to provide a variable"
                           " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the"
                           " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction"
@@ -667,6 +672,7 @@ struct SolverChoice {
     int RhoQr_comp {-1};
 
     std::string windfarm_loc_table, windfarm_spec_table;
+    std::string windfarm_blade_table, windfarm_airfoil_tables;
     amrex::Real sampling_distance_by_D = -1.0;
     amrex::Real turb_disk_angle = -1.0;
     amrex::Real windfarm_x_shift = -1.0;

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -888,7 +888,7 @@ private:
                                                     "vorticity_x","vorticity_y","vorticity_z",
                                                     "magvel", "divU",
                                                     "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens",
-                                                    "eq_pot_temp", "num_turb", "SMark",
+                                                    "eq_pot_temp", "num_turb", "SMark0", "SMark1",
                                                     "dpdx", "dpdy", "pres_hse_x", "pres_hse_y",
                                                     "z_phys", "detJ" , "mapfac", "lat_m", "lon_m",
                                                     // Time averaged velocity

diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp
@@ -234,6 +234,9 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     }
     if (solverChoice.windfarm_type == WindFarmType::SimpleAD) {
         vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
+    }
+    if (solverChoice.windfarm_type == WindFarmType::GeneralAD) {
+        vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt
     }
         Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell
         SMark[lev].define(ba, dm, 2, ngrow_state); // Free stream velocity/source term

diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp
@@ -83,14 +83,38 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
                       derived_names[i] != "graup_accum") )
                 {
                     if ( (solverChoice.moisture_type != MoistureType::None)      ||
-                         (derived_names[i] != "qv" && derived_names[i] != "qc") ) {
+                         (derived_names[i] != "qv" && derived_names[i] != "qc" and
+                          derived_names[i] != "num_turb" and derived_names[i] != "SMark0" and
+                          derived_names[i] != "SMark1") ) {
                         tmp_plot_names.push_back(derived_names[i]);
                     }
                 } // moisture_type
             } // use_terrain?
         } // hasElement
     }
 
+#ifdef ERF_USE_WINDFARM
+    for (int i = 0; i < derived_names.size(); ++i) {
+        if ( containerHasElement(plot_var_names, derived_names[i]) ) {
+            if(solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP) {
+                if(derived_names[i] == "num_turb") {
+                    tmp_plot_names.push_back(derived_names[i]);
+                }
+            }
+            if(solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+                if(derived_names[i] == "num_turb" or derived_names[i] == "SMark0" or derived_names[i] == "Smark1") {
+                    tmp_plot_names.push_back(derived_names[i]);
+                }
+            }
+            if(solverChoice.windfarm_type == WindFarmType::GeneralAD) {
+                if(derived_names[i] == "num_turb" or derived_names[i] == "SMark1") {
+                    tmp_plot_names.push_back(derived_names[i]);
+                }
+            }
+        }
+    }
+#endif
+
 #ifdef ERF_USE_PARTICLES
     const auto& particles_namelist( particleData.getNamesUnalloc() );
     for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) {
@@ -444,7 +468,9 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         }
 
 #ifdef ERF_USE_WINDFARM
-        if (containerHasElement(plot_var_names, "num_turb"))
+        if (containerHasElement(plot_var_names, "num_turb") and
+            (solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or
+            solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD))
         {
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
@@ -461,7 +487,8 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             mf_comp ++;
         }
 
-        if(containerHasElement(plot_var_names, "SMark") and solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+        if(containerHasElement(plot_var_names, "SMark0") and
+           solverChoice.windfarm_type == WindFarmType::SimpleAD) {
              for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
             {
                 const Box& bx = mfi.tilebox();
@@ -473,6 +500,21 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             }
             mf_comp ++;
         }
+
+         if(containerHasElement(plot_var_names, "SMark1") and
+           (solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD)) {
+             for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
+            {
+                const Box& bx = mfi.tilebox();
+                const Array4<Real>& derdat  = mf[lev].array(mfi);
+                const Array4<Real const>& SMark_array = SMark[lev].const_array(mfi);
+                ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    derdat(i, j, k, mf_comp) = SMark_array(i,j,k,1);
+                });
+            }
+            mf_comp ++;
+        }
+
 #endif
 
         int klo = geom[lev].Domain().smallEnd(2);

diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp
@@ -29,16 +29,23 @@ ERF::init_windfarm (int lev)
                              true, false);
     }
 
+
     windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev]);
 
     windfarm->write_turbine_locations_vtk();
 
-    if(solverChoice.windfarm_type == WindFarmType::SimpleAD) {
+    if(solverChoice.windfarm_type == WindFarmType::SimpleAD or
+       solverChoice.windfarm_type == WindFarmType::GeneralAD) {
         windfarm->fill_SMark_multifab(geom[lev], SMark[lev],
                                       solverChoice.sampling_distance_by_D,
                                       solverChoice.turb_disk_angle);
         windfarm->write_actuator_disks_vtk(geom[lev]);
     }
+
+    if(solverChoice.windfarm_type == WindFarmType::GeneralAD) {
+        //windfarm->read_windfarm_blade_table(solverChoice.windfarm_blade_table);
+        //windfarm->read_airfoil_tables
+    }
 }
 
 void

diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp
@@ -42,7 +42,7 @@ WindFarm::read_windfarm_locations_table(const std::string windfarm_loc_table,
     }
     else {
         amrex::Abort("Are you using windfarms? For windfarm simulations, the inputs need to have an"
-                     " entry erf.windfarm_loc_table which should not be None. \n");
+                     " entry erf.windfarm_loc_type which should be either lat_lon or x_y. \n");
     }
 
     set_turb_loc(xloc, yloc);
@@ -138,6 +138,8 @@ WindFarm::init_windfarm_x_y (const std::string windfarm_loc_table)
     Real value1, value2;
 
     while (file >> value1 >> value2) {
+        value1 = value1 + 1e-3;
+        value2 = value2 + 1e-3;
         xloc.push_back(value1);
         yloc.push_back(value2);
     }
@@ -193,6 +195,20 @@ WindFarm::read_windfarm_spec_table(const std::string windfarm_spec_table)
 
 }
 
+void
+WindFarm::read_windfarm_blade_table(const std::string windfarm_blade_table)
+{
+    std::ifstream filename(windfarm_blade_table);
+    if (!filename.is_open()) {
+        Error("You are using a generalized actuator disk model based on blade element theory. This needs info of blades."
+                      " An entry erf.windfarm_blade_table is needed. Either the entry is missing or the file specified"
+                      " in the entry - " + windfarm_blade_table + " is missing.");
+    }
+    else {
+        Print() << "Reading in wind farm blade table: " << windfarm_blade_table << "\n";
+    }
+}
+
 void
 WindFarm::fill_Nturb_multifab(const Geometry& geom,
                               MultiFab& mf_Nturb)
@@ -369,22 +385,19 @@ WindFarm::write_actuator_disks_vtk(const Geometry& geom)
         }
         fprintf(file_actuator_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast<long int>(num_turb_in_dom*npts), "float");
 
-        Real nx = std::cos(my_turb_disk_angle);
-        Real ny = std::sin(my_turb_disk_angle);
+        Real nx = std::cos(my_turb_disk_angle+0.5*M_PI);
+        Real ny = std::sin(my_turb_disk_angle+0.5*M_PI);
 
         for(int it=0; it<xloc.size(); it++){
-            Real x;
-            x = xloc[it]+1e-12;
             for(int pt=0;pt<100;pt++){
-                Real y, z;
+                Real x, y, z;
                 Real theta = 2.0*M_PI/npts*pt;
-                y = yloc[it]+rotor_rad*cos(theta);
-                z = hub_height+rotor_rad*sin(theta);
+                x = xloc[it] + rotor_rad*cos(theta)*nx;
+                y = yloc[it] + rotor_rad*cos(theta)*ny;
+                z = hub_height + rotor_rad*sin(theta);
                 fprintf(file_actuator_disks_all, "%0.15g %0.15g %0.15g\n", x, y, z);
                 if(xloc[it] > ProbLoArr[0] and xloc[it] < ProbHiArr[0] and yloc[it] > ProbLoArr[1] and yloc[it] < ProbHiArr[1]) {
-                    Real xp = (x-xloc[it])*nx - (y-yloc[it])*ny;
-                    Real yp = (x-xloc[it])*nx + (y-yloc[it])*ny;
-                    fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", xloc[it]+xp, yloc[it]+yp, z);
+                    fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", x, y, z);
                 }
             }
         }

diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H
@@ -9,6 +9,7 @@
 #include "ERF_Fitch.H"
 #include "ERF_EWP.H"
 #include "ERF_SimpleAD.H"
+#include "ERF_GeneralAD.H"
 
 class WindFarm : public NullWindFarm {
 
@@ -24,14 +25,21 @@ public:
         m_windfarm_model.resize(nlev);
         if (a_windfarm_type == WindFarmType::Fitch) {
             SetModel<Fitch>();
+            amrex::Print() << "Fitch windfarm model!\n";
         }
         else if (a_windfarm_type == WindFarmType::EWP) {
             SetModel<EWP>();
+            amrex::Print() << "EWP windfarm model!\n";
         }
         else if (a_windfarm_type == WindFarmType::SimpleAD) {
             SetModel<SimpleAD>();
-            amrex::Print() << "Simple actuator disk windfarm model!\n";
-        } else {
+            amrex::Print() << "Simplified actuator disk windfarm model!\n";
+        }
+        else if (a_windfarm_type == WindFarmType::GeneralAD) {
+            SetModel<GeneralAD>();
+            amrex::Print() << "Generalized actuator disk windfarm model!\n";
+        }
+         else {
             amrex::Abort("WindFarm: Dont know this windfarm_type!") ;
         }
     }
@@ -55,6 +63,10 @@ public:
 
     void read_windfarm_spec_table(const std::string windfarm_spec_table);
 
+    void read_windfarm_blade_table(const std::string windfarm_blade_table);
+
+    void read_windfarm_airofil_tables(const std::string windfarm_airfoils_tables);
+
     void fill_Nturb_multifab(const amrex::Geometry& geom,
                              amrex::MultiFab& mf_Nturb);
 

diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
@@ -14,7 +14,7 @@ Real compute_A(const Real z,
 
     Real d  = std::min(std::fabs(z - hub_height), rotor_rad);
     Real theta = std::acos(d/rotor_rad);
-    Real A_s = rotor_rad*rotor_rad*theta - d*std::pow(rotor_rad*rotor_rad - d*d, 0.5);
+    Real A_s = rotor_rad*rotor_rad*theta - d*std::pow(std::max(rotor_rad*rotor_rad - d*d,0.0), 0.5);
     Real A = PI*rotor_rad*rotor_rad/2.0 - A_s;
 
     return A;