Wind models on GPU

Source/IO/Plotfile.cpp

@@ -427,7 +427,6 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
 #ifdef ERF_USE_WINDFARM
         if (containerHasElement(plot_var_names, "num_turb"))
         {
-        std::cout << "Plotting num_turb" << "\n";
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif

Source/Initialization/ERF_init_windfarm.cpp

@@ -13,6 +13,7 @@ using namespace amrex;
  *
  * @param lev Integer specifying the current level
  */
+
 void
 ERF::init_windfarm (int lev)
 {
@@ -44,10 +45,20 @@ ERF::init_windfarm (int lev)
         fclose(file_turbloc_vtk);
     }
 
+    amrex::Gpu::DeviceVector<Real> d_xloc(xloc.size());
+    amrex::Gpu::DeviceVector<Real> d_yloc(yloc.size());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
+    amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
+
+    Real* d_xloc_ptr     = d_xloc.data();
+    Real* d_yloc_ptr     = d_yloc.data();
+
     Nturb[lev].setVal(0);
 
     int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0);
     int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1);
+    auto dx = geom[lev].CellSizeArray();
+    int num_turb = xloc.size();
 
      // Initialize wind farm
     for ( MFIter mfi(Nturb[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
@@ -57,14 +68,13 @@ ERF::init_windfarm (int lev)
             int li = amrex::min(amrex::max(i, i_lo), i_hi);
             int lj = amrex::min(amrex::max(j, j_lo), j_hi);
 
-            auto dx = geom[lev].CellSizeArray();
             Real x1 = li*dx[0], x2 = (li+1)*dx[0];
             Real y1 = lj*dx[1], y2 = (lj+1)*dx[1];
 
-            for(int it=0; it<xloc.size(); it++){
-                if( xloc[it]+1e-12 > x1 and xloc[it]+1e-12 < x2 and
-                    yloc[it]+1e-12 > y1 and yloc[it]+1e-12 < y2){
-                    Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1;
+            for(int it=0; it<num_turb; it++){
+                if( d_xloc_ptr[it]+1e-12 > x1 and d_xloc_ptr[it]+1e-12 < x2 and
+                    d_yloc_ptr[it]+1e-12 > y1 and d_yloc_ptr[it]+1e-12 < y2){
+                       Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1;
                 }
             }
         });

Source/WindFarmParametrization/EWP/AdvanceEWP.cpp

@@ -5,8 +5,6 @@
 
 using namespace amrex;
 
-Real C_TKE_ewp = 0.0;
-Real K_turb = 1.0;
 
 void ewp_advance (int lev,
                   const Geometry& geom,
@@ -64,6 +62,9 @@ void ewp_source_terms_cellcentered (const Geometry& geom,
   auto ProbLoArr = geom.ProbLoArray();
   Real sigma_0 = 1.7*rotor_rad;
 
+  Real d_rotor_rad = rotor_rad;
+  Real d_hub_height = hub_height;
+
   // Domain valid box
   const amrex::Box& domain = geom.Domain();
   int domlo_x = domain.smallEnd(0);
@@ -89,6 +90,10 @@ void ewp_source_terms_cellcentered (const Geometry& geom,
 
         amrex::IntVect lo = bx.smallEnd();
 
+        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();
+
         ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
             int ii = amrex::min(amrex::max(i, domlo_x), domhi_x);
             int jj = amrex::min(amrex::max(j, domlo_y), domhi_y);
@@ -105,13 +110,19 @@ void ewp_source_terms_cellcentered (const Geometry& geom,
                                  v_vel(i,j,k)*v_vel(i,j,k) +
                                  w_vel(i,j,kk)*w_vel(i,j,kk), 0.5);
 
-            Real C_T_ewp = interpolate_1d(wind_speed.dataPtr(), thrust_coeff.dataPtr(), z, wind_speed.size());
+            Real C_T = interpolate_1d(wind_speed_d, thrust_coeff_d, Vabs, n_spec_table);
+
+            Real C_TKE = 0.0;
+            Real K_turb = 1.0;
 
             Real L_wake = std::pow(dx[0]*dx[1],0.5)/2.0;
-            Real sigma_e = Vabs/(3.0*K_turb*L_wake)*(std::pow(2.0*K_turb*L_wake/Vabs + std::pow(sigma_0,2),3.0/2.0) - std::pow(sigma_0,3));
+            Real sigma_e = Vabs/(3.0*K_turb*L_wake)*
+                           (std::pow(2.0*K_turb*L_wake/Vabs + std::pow(sigma_0,2),3.0/2.0) - std::pow(sigma_0,3));
 
             Real phi     = std::atan2(v_vel(i,j,k),u_vel(i,j,k)); // Wind direction w.r.t the x-dreiction
-            Real fac = -std::pow(M_PI/8.0,0.5)*C_T_ewp*std::pow(rotor_rad,2)*std::pow(Vabs,2)/(dx[0]*dx[1]*sigma_e)*std::exp(-0.5*std::pow((z - hub_height)/sigma_e,2));
+            Real fac = -std::pow(M_PI/8.0,0.5)*C_T*std::pow(d_rotor_rad,2)*
+                        std::pow(Vabs,2)/(dx[0]*dx[1]*sigma_e)*
+                        std::exp(-0.5*std::pow((z - d_hub_height)/sigma_e,2));
             ewp_array(i,j,k,0) = fac*std::cos(phi)*Nturb_array(i,j,k);
             ewp_array(i,j,k,1) = fac*std::sin(phi)*Nturb_array(i,j,k);
             ewp_array(i,j,k,2) = 0.0;

Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp

@@ -5,9 +5,11 @@
 
 using namespace amrex;
 
-Real C_TKE = 0.0;
-
-Real compute_A(Real z)
+AMREX_FORCE_INLINE
+AMREX_GPU_DEVICE
+Real compute_A(const Real z,
+               const Real hub_height,
+               const Real rotor_rad)
 {
 
     Real d  = std::min(std::fabs(z - hub_height), rotor_rad);
@@ -18,11 +20,16 @@ Real compute_A(Real z)
     return A;
 }
 
-Real compute_Aijk(Real z_k, Real z_kp1)
+AMREX_FORCE_INLINE
+AMREX_GPU_DEVICE
+Real compute_Aijk(const Real z_k,
+                  const Real z_kp1,
+                  const Real hub_height,
+                  const Real rotor_rad)
 {
 
-    Real A_k   = compute_A(z_k);
-    Real A_kp1 = compute_A(z_kp1);
+    Real A_k   = compute_A(z_k, hub_height, rotor_rad);
+    Real A_kp1 = compute_A(z_kp1, hub_height, rotor_rad);
 
     Real check = (z_k - hub_height)*(z_kp1 - hub_height);
     Real A_ijk;
@@ -102,11 +109,13 @@ void fitch_source_terms_cellcentered (const Geometry& geom,
   int domhi_z = domain.bigEnd(2) + 1;
 
 
-  Real sum = 0.0;
-  Real *sum_area = ∑
+  //Real sum = 0.0;
+  //Real *sum_area = ∑
 
   // The order of variables are - Vabs dVabsdt, dudt, dvdt, dTKEdt
   mf_vars_fitch.setVal(0.0);
+  Real d_hub_height = hub_height;
+  Real d_rotor_rad = rotor_rad;
 
   for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
 
@@ -121,6 +130,10 @@ void fitch_source_terms_cellcentered (const Geometry& geom,
 
         amrex::IntVect lo = bx.smallEnd();
 
+        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();
+
         ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
             int ii = amrex::min(amrex::max(i, domlo_x), domhi_x);
             int jj = amrex::min(amrex::max(j, domlo_y), domhi_y);
@@ -134,15 +147,16 @@ void fitch_source_terms_cellcentered (const Geometry& geom,
             Real z_k   = kk*dx[2];
             Real z_kp1 = (kk+1)*dx[2];
 
-            Real A_ijk = compute_Aijk(z_k, z_kp1);
+            Real A_ijk = compute_Aijk(z_k, z_kp1, d_hub_height, d_rotor_rad);
 
             // Compute Fitch source terms
 
             Real Vabs = std::pow(u_vel(i,j,k)*u_vel(i,j,k) +
                                  v_vel(i,j,k)*v_vel(i,j,k) +
                                  w_vel(i,j,kk)*w_vel(i,j,kk), 0.5);
 
-            Real C_T = interpolate_1d(wind_speed.dataPtr(), thrust_coeff.dataPtr(), z, wind_speed.size());
+            Real C_T = interpolate_1d(wind_speed_d, thrust_coeff_d, Vabs, n_spec_table);
+            Real C_TKE = 0.0;
 
             fitch_array(i,j,k,0) = Vabs;
             fitch_array(i,j,k,1) =  -0.5*Nturb_array(i,j,k)/(dx[0]*dx[1])*C_T*Vabs*Vabs*A_ijk/(z_kp1 - z_k);

Source/WindFarmParametrization/WindFarm.H

@@ -4,9 +4,11 @@
 #include <DataStruct.H>
 #include <AMReX_Geometry.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_Gpu.H>
 
 extern amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power;
 extern amrex::Vector<amrex::Real> wind_speed, thrust_coeff, power;
+extern amrex::Gpu::DeviceVector<amrex::Real> d_wind_speed, d_thrust_coeff, d_power;
 
 void read_in_table(std::string windfarm_spec_table);
 

Source/WindFarmParametrization/WindFarm.cpp

@@ -3,8 +3,9 @@
 #include <EWP.H>
 using namespace amrex;
 
-amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power;
-amrex::Vector<amrex::Real> wind_speed, thrust_coeff, power;
+Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power;
+Vector<Real> wind_speed, thrust_coeff, power;
+Gpu::DeviceVector<Real> d_wind_speed, d_thrust_coeff, d_power;
 
 void read_in_table(std::string windfarm_spec_table)
 {
@@ -18,11 +19,11 @@ void read_in_table(std::string windfarm_spec_table)
      // Read turbine data from wind-turbine-1.tbl
     std::ifstream file_turb_table(windfarm_spec_table);
     if (!file_turb_table.is_open()) {
-        amrex::Error("Wind farm specifications table not found. Either the inputs is missing the "
+        Error("Wind farm specifications table not found. Either the inputs is missing the "
                       "erf.windfarm_spec_table entry or the file specified in the entry - " + windfarm_spec_table + " is missing.");
     }
     else {
-        amrex::Print() << "Reading in wind farm specifications table: " << windfarm_spec_table << "\n";
+        Print() << "Reading in wind farm specifications table: " << windfarm_spec_table << "\n";
     }
 
     int nlines;
@@ -31,23 +32,31 @@ void read_in_table(std::string windfarm_spec_table)
     thrust_coeff.resize(nlines);
     power.resize(nlines);
 
+    d_wind_speed.resize(nlines);
+    d_thrust_coeff.resize(nlines);
+    d_power.resize(nlines);
+
     Real rotor_dia;
     file_turb_table >> hub_height >> rotor_dia >> thrust_coeff_standing >> nominal_power;
     rotor_rad = rotor_dia*0.5;
     if(rotor_rad > hub_height) {
-        amrex::Abort("The blade length is more than the hub height. Check the second line in wind-turbine-1.tbl. Aborting.....");
+        Abort("The blade length is more than the hub height. Check the second line in wind-turbine-1.tbl. Aborting.....");
     }
     if(thrust_coeff_standing > 1.0) {
-        amrex::Abort("The standing thrust coefficient is greater than 1. Check the second line in wind-turbine-1.tbl. Aborting.....");
+        Abort("The standing thrust coefficient is greater than 1. Check the second line in wind-turbine-1.tbl. Aborting.....");
     }
 
     for(int iline=0;iline<nlines;iline++){
         file_turb_table >> wind_speed[iline] >> thrust_coeff[iline] >> power[iline];
         if(thrust_coeff[iline] > 1.0) {
-            amrex::Abort("The thrust coefficient is greater than 1. Check wind-turbine-1.tbl. Aborting.....");
+            Abort("The thrust coefficient is greater than 1. Check wind-turbine-1.tbl. Aborting.....");
         }
     }
     file_turb_table.close();
+
+    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, power.begin(), power.end(), d_power.begin());
 }
 
 
@@ -56,7 +65,7 @@ void advance_windfarm (int lev,
                        const Real& dt_advance,
                        MultiFab& cons_in,
                        MultiFab& U_old, MultiFab& V_old, MultiFab& W_old,
-                       MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb,
+                       MultiFab& mf_vars_windfarm, const MultiFab& mf_Nturb,
                        SolverChoice& solver_choice)
 {
     if (solver_choice.windfarm_type == WindFarmType::Fitch) {