Some changes to the simplified actuator disk model (erf-model#1830)

* Generalize simplified actuator disk model * Removing unwatned files --------- Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Aaron M. Lattanzi <[email protected]>
AMLattanzi · Sep 27, 2024 · a4b7da7 · a4b7da7
1 parent 1c1c737
commit a4b7da7
Show file tree

Hide file tree

Showing 7 changed files with 41 additions and 304 deletions.
diff --git a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py b/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qc.py
diff --git a/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py b/Exec/SuperCell_3D/PostProcessing/ParaViewExportSTLMacro_qr.py
diff --git a/Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py b/Exec/SuperCell_3D/PostProcessing/SuperCell3D_Blender2p80.py
diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
@@ -402,17 +402,17 @@ struct SolverChoice {
         // by the diameter of the turbine
         pp.query("sampling_distance_by_D", sampling_distance_by_D);
         if(windfarm_type==WindFarmType::SimpleAD and sampling_distance_by_D < 0.0) {
-             amrex::Abort("To use simplified actuator disks, you need to provide a variable "
-                           "erf.sampling_distance_by_D in the inputs which specifies the upstream "
-                           "distance as a factor of the turbine diameter at which the incoming free stream "
-                           "velocity will be computed at");
+             amrex::Abort("To use simplified actuator disks, you need to provide a variable"
+                           " erf.sampling_distance_by_D in the inputs which specifies the upstream"
+                           " distance as a factor of the turbine diameter at which the incoming free stream"
+                           " velocity will be computed at.");
         }
         pp.query("turb_disk_angle_from_x", turb_disk_angle);
         if(windfarm_type==WindFarmType::SimpleAD 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 "
-                          " will have turb_disk_angle value of 90 deg ");
+            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"
+                          " will have turb_disk_angle value of 90 deg.");
         }
 
         pp.query("windfarm_x_shift",windfarm_x_shift);

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -411,7 +411,9 @@ public:
                        const amrex::Real& dt_advance,
                        amrex::MultiFab& cons_in,
                        amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old,
-                       amrex::MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb);
+                       amrex::MultiFab& mf_vars_windfarm,
+                       const amrex::MultiFab& mf_Nturb,
+                       const amrex::MultiFab& mf_SMark);
 #endif
 
 #ifdef ERF_USE_EB

diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
@@ -124,12 +124,12 @@ Fitch::source_terms_cellcentered (const Geometry& geom,
   mf_vars_fitch.setVal(0.0);
   Real d_hub_height = hub_height;
   Real d_rotor_rad = rotor_rad;
-     Gpu::DeviceVector<Real> d_wind_speed(wind_speed.size());
-    Gpu::DeviceVector<Real> d_thrust_coeff(thrust_coeff.size());
+  Gpu::DeviceVector<Real> d_wind_speed(wind_speed.size());
+  Gpu::DeviceVector<Real> d_thrust_coeff(thrust_coeff.size());
 
   // Copy data from host vectors to device vectors
-    Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin());
-    Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin());
+  Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin());
+  Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin());
 
   for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
 

diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp
@@ -1,5 +1,6 @@
 #include <ERF_SimpleAD.H>
 #include <ERF_IndexDefines.H>
+#include <ERF_Interpolation_1D.H>
 
 using namespace amrex;
 
@@ -173,6 +174,17 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
     Real ny = -std::sin(turb_disk_angle);
     Real d_turb_disk_angle = turb_disk_angle;
 
+    Gpu::DeviceVector<Real> d_wind_speed(wind_speed.size());
+    Gpu::DeviceVector<Real> d_thrust_coeff(thrust_coeff.size());
+
+    // Copy data from host vectors to device vectors
+    Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin());
+    Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin());
+
+    const Real* wind_speed_d     = d_wind_speed.dataPtr();
+    const Real* thrust_coeff_d   = d_thrust_coeff.dataPtr();
+    const int n_spec_table = d_wind_speed.size();
+
     for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
 
         const Box& gbx      = mfi.growntilebox(1);
@@ -193,11 +205,22 @@ SimpleAD::source_terms_cellcentered (const Geometry& geom,
                 Real avg_vel  = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);
                 Real phi      = d_freestream_phi_ptr[it]/(d_disk_cell_count_ptr[it] + 1e-10);
 
+                Real C_T = interpolate_1d(wind_speed_d, thrust_coeff_d, avg_vel, n_spec_table);
+                Real a;
+                if(C_T <= 1) {
+                    a = 0.5 - 0.5*std::pow(1.0-C_T,0.5);
+                }
                 Real Uinfty_dot_nhat = avg_vel*(std::cos(phi)*nx + std::sin(phi)*ny);
                 if(SMark_array(ii,jj,kk,1) == static_cast<double>(it)) {
-                        check_int++;
-                        source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*0.25*(1.0-0.25)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi);
-                        source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*0.25*(1.0-0.25)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi);
+                    check_int++;
+                    if(C_T <= 1) {
+                        source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi);
+                           source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi);
+                    }
+                    else {
+                        source_x = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi);
+                        source_y = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi);
+                    }
                 }
             }
             if(check_int > 1){