From 9160d419ee99be1f241849a77f2f08d7b661a80d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 18 Jun 2024 15:58:50 -0700 Subject: [PATCH 1/3] Making work on GPUs --- Source/Initialization/ERF_init_windfarm.cpp | 20 +++++++++---- .../Fitch/AdvanceFitch.cpp | 28 +++++++++++++------ Source/WindFarmParametrization/WindFarm.cpp | 6 ++-- 3 files changed, 37 insertions(+), 17 deletions(-) diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index fdaf02765..e0ff4135a 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/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 d_xloc(xloc.size()); + amrex::Gpu::DeviceVector 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 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 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; } } }); diff --git a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp index 1526a24f6..569145f17 100644 --- a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp +++ b/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) { @@ -134,7 +143,7 @@ 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 @@ -143,6 +152,7 @@ void fitch_source_terms_cellcentered (const Geometry& geom, 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_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); diff --git a/Source/WindFarmParametrization/WindFarm.cpp b/Source/WindFarmParametrization/WindFarm.cpp index 3e54d02c3..b487c91b6 100644 --- a/Source/WindFarmParametrization/WindFarm.cpp +++ b/Source/WindFarmParametrization/WindFarm.cpp @@ -1,6 +1,6 @@ #include #include -#include +//#include using namespace amrex; amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; @@ -63,7 +63,7 @@ void advance_windfarm (int lev, fitch_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, mf_vars_windfarm, mf_Nturb); } else if (solver_choice.windfarm_type == WindFarmType::EWP) { - ewp_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, - mf_vars_windfarm, mf_Nturb); + //ewp_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, + // mf_vars_windfarm, mf_Nturb); } } From 7810615f488298701fcdfe73cc20e78fdb462c69 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 20 Jun 2024 17:44:09 -0700 Subject: [PATCH 2/3] Making work on GPUs --- Source/IO/Plotfile.cpp | 1 - .../Fitch/AdvanceFitch.cpp | 2 +- Source/WindFarmParametrization/WindFarm.H | 2 ++ Source/WindFarmParametrization/WindFarm.cpp | 25 +++++++++++++------ 4 files changed, 20 insertions(+), 10 deletions(-) diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp index c9e0ff65e..b3315809c 100644 --- a/Source/IO/Plotfile.cpp +++ b/Source/IO/Plotfile.cpp @@ -427,7 +427,6 @@ ERF::WritePlotFile (int which, Vector 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 diff --git a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp index 569145f17..5590ac032 100644 --- a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp @@ -151,7 +151,7 @@ void fitch_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 = interpolate_1d(wind_speed.dataPtr(), thrust_coeff.dataPtr(), z, wind_speed.size()); + Real C_T = interpolate_1d(d_wind_speed.dataPtr(), d_thrust_coeff.dataPtr(), z, d_wind_speed.size()); Real C_TKE = 0.0; fitch_array(i,j,k,0) = Vabs; diff --git a/Source/WindFarmParametrization/WindFarm.H b/Source/WindFarmParametrization/WindFarm.H index a0376ae52..9f101175b 100644 --- a/Source/WindFarmParametrization/WindFarm.H +++ b/Source/WindFarmParametrization/WindFarm.H @@ -4,9 +4,11 @@ #include #include #include +#include extern amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; extern amrex::Vector wind_speed, thrust_coeff, power; +extern amrex::Gpu::DeviceVector d_wind_speed, d_thrust_coeff, d_power; void read_in_table(std::string windfarm_spec_table); diff --git a/Source/WindFarmParametrization/WindFarm.cpp b/Source/WindFarmParametrization/WindFarm.cpp index b487c91b6..854c9d44d 100644 --- a/Source/WindFarmParametrization/WindFarm.cpp +++ b/Source/WindFarmParametrization/WindFarm.cpp @@ -3,8 +3,9 @@ //#include using namespace amrex; -amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; -amrex::Vector wind_speed, thrust_coeff, power; +Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; +Vector wind_speed, thrust_coeff, power; +Gpu::DeviceVector 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> 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) { From 188821a2d617eeca3b4b2c858f75301789593acf Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Fri, 28 Jun 2024 15:29:56 -0700 Subject: [PATCH 3/3] Fitch and EWP works on GPUs now --- .../EWP/AdvanceEWP.cpp | 21 ++++++++++++++----- .../Fitch/AdvanceFitch.cpp | 6 +++++- Source/WindFarmParametrization/WindFarm.cpp | 6 +++--- 3 files changed, 24 insertions(+), 9 deletions(-) diff --git a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp index d92f3ded9..098868367 100644 --- a/Source/WindFarmParametrization/EWP/AdvanceEWP.cpp +++ b/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; diff --git a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp index 5590ac032..56a1a50a2 100644 --- a/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/AdvanceFitch.cpp @@ -130,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); @@ -151,7 +155,7 @@ void fitch_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 = interpolate_1d(d_wind_speed.dataPtr(), d_thrust_coeff.dataPtr(), z, d_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; diff --git a/Source/WindFarmParametrization/WindFarm.cpp b/Source/WindFarmParametrization/WindFarm.cpp index 854c9d44d..5468bc16e 100644 --- a/Source/WindFarmParametrization/WindFarm.cpp +++ b/Source/WindFarmParametrization/WindFarm.cpp @@ -1,6 +1,6 @@ #include #include -//#include +#include using namespace amrex; Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; @@ -72,7 +72,7 @@ void advance_windfarm (int lev, fitch_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, mf_vars_windfarm, mf_Nturb); } else if (solver_choice.windfarm_type == WindFarmType::EWP) { - //ewp_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, - // mf_vars_windfarm, mf_Nturb); + ewp_advance(lev, geom, dt_advance, cons_in, U_old, V_old, W_old, + mf_vars_windfarm, mf_Nturb); } }