From ffc17e3efe738fcc766de7fc8a2ce42535363a31 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Mon, 2 Dec 2024 11:55:03 -0800 Subject: [PATCH] Power output for all wind models --- Source/IO/ERF_Plotfile.cpp | 5 +- Source/Initialization/ERF_init_windfarm.cpp | 11 ++- .../ERF_InitWindFarm.cpp | 78 +++++++++++++++---- Source/WindFarmParametrization/ERF_WindFarm.H | 7 +- .../EWP/ERF_AdvanceEWP.cpp | 68 +++++++++++++++- Source/WindFarmParametrization/EWP/ERF_EWP.H | 11 ++- .../Fitch/ERF_AdvanceFitch.cpp | 63 ++++++++++++++- .../WindFarmParametrization/Fitch/ERF_Fitch.H | 8 ++ .../ERF_AdvanceGeneralAD.cpp | 42 ++++++++-- .../GeneralActuatorDisk/ERF_GeneralAD.H | 2 + .../ERF_AdvanceSimpleAD.cpp | 2 +- 11 files changed, 270 insertions(+), 27 deletions(-) diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index ea7ef0287..01a79c1cd 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -97,7 +97,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector 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") { + if(derived_names[i] == "num_turb" or derived_names[i] == "SMark0") { tmp_plot_names.push_back(derived_names[i]); } } @@ -510,7 +510,8 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p } if( containerHasElement(plot_var_names, "SMark0") and - (solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD) ) { + (solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or + 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(); diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index bdaae7bdf..ab46ed327 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/Source/Initialization/ERF_init_windfarm.cpp @@ -29,10 +29,19 @@ ERF::init_windfarm (int lev) true, false); } - windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev], z_phys_cc[lev]); + windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev], z_phys_nd[lev]); windfarm->write_turbine_locations_vtk(); + + if(solverChoice.windfarm_type == WindFarmType::Fitch or + solverChoice.windfarm_type == WindFarmType::EWP) { + windfarm->fill_SMark_multifab_mesoscale_models(geom[lev], + SMark[lev], + Nturb[lev], + z_phys_nd[lev]); + } + if(solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD) { windfarm->fill_SMark_multifab(geom[lev], SMark[lev], diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp index 482a39379..e4dea885e 100644 --- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp +++ b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp @@ -352,7 +352,7 @@ WindFarm::read_windfarm_airfoil_tables (const std::string windfarm_airfoil_table void WindFarm::gatherKeyValuePairs(const std::vector>& localData, - std::vector>& globalData) + std::vector>& globalData) { int myRank = amrex::ParallelDescriptor::MyProc(); int nProcs = amrex::ParallelDescriptor::NProcs(); @@ -419,17 +419,13 @@ WindFarm::gatherKeyValuePairs(const std::vector>& localDa amrex::ParallelDescriptor::Bcast(&globalData[i].first, 1, 0); amrex::ParallelDescriptor::Bcast(&globalData[i].second, 1, 0); } - - for (const auto& kv : globalData) { - std::cout << "Rank " << myRank << "Key: " << kv.first << ", Value: " << kv.second << std::endl; - } } void WindFarm::fill_Nturb_multifab (const Geometry& geom, MultiFab& mf_Nturb, - std::unique_ptr& z_phys_cc) + std::unique_ptr& z_phys_nd) { zloc.resize(xloc.size(),-1.0); @@ -461,13 +457,13 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom, auto ProbLoArr = geom.ProbLoArray(); int num_turb = xloc.size(); - bool is_terrain = z_phys_cc ? true: false; + bool is_terrain = z_phys_nd ? true: false; // Initialize wind farm for ( MFIter mfi(mf_Nturb,TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); auto Nturb_array = mf_Nturb.array(mfi); - const Array4& z_cc_arr = (z_phys_cc) ? z_phys_cc->const_array(mfi) : Array4{}; + const Array4& z_nd_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4{}; int k0 = bx.smallEnd()[2]; ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { int li = amrex::min(amrex::max(i, i_lo), i_hi); @@ -483,7 +479,7 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom, d_yloc_ptr[it]+1e-3 > y1 and d_yloc_ptr[it]+1e-3 < y2){ Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1; if(is_terrain) { - d_zloc_ptr[it] = z_cc_arr(i,j,k0); + d_zloc_ptr[it] = z_nd_arr(i,j,k0); d_turb_index_ptr[it] = it; } else { @@ -498,8 +494,6 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom, Gpu::copy(Gpu::deviceToHost, d_zloc.begin(), d_zloc.end(), zloc.begin()); Gpu::copy(Gpu::deviceToHost, d_turb_index.begin(), d_turb_index.end(), turb_index.begin()); - if(is_terrain) { - std::vector> turb_index_zloc; for(int it=0;it& z_phys_nd) +{ + mf_SMark.setVal(-1.0); + + Real d_hub_height = hub_height; + + amrex::Gpu::DeviceVector d_xloc(xloc.size()); + amrex::Gpu::DeviceVector d_yloc(yloc.size()); + amrex::Gpu::DeviceVector d_zloc(xloc.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()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin()); + + int i_lo = geom.Domain().smallEnd(0); int i_hi = geom.Domain().bigEnd(0); + int j_lo = geom.Domain().smallEnd(1); int j_hi = geom.Domain().bigEnd(1); + int k_lo = geom.Domain().smallEnd(2); int k_hi = geom.Domain().bigEnd(2); + + auto dx = geom.CellSizeArray(); + auto ProbLoArr = geom.ProbLoArray(); + + // Initialize wind farm + for ( MFIter mfi(mf_SMark,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.tilebox(); + auto SMark_array = mf_SMark.array(mfi); + auto Nturb_array = mf_Nturb.array(mfi); + const Array4& z_nd_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4{}; + int k0 = bx.smallEnd()[2]; + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if(Nturb_array(i,j,k,0) > 0) { + int li = amrex::min(amrex::max(i, i_lo), i_hi); + int lj = amrex::min(amrex::max(j, j_lo), j_hi); + int lk = amrex::min(amrex::max(k, k_lo), k_hi); + + Real z1 = (z_nd_arr) ? z_nd_arr(li,lj,lk) : ProbLoArr[2] + lk * dx[2]; + Real z2 = (z_nd_arr) ? z_nd_arr(li,lj,lk+1) : ProbLoArr[2] + (lk+1) * dx[2]; + + Real zturb; + if(z_nd_arr) { + zturb = z_nd_arr(li,lj,k0) + d_hub_height; + } else { + zturb = d_hub_height; + } + if(zturb+1e-3 > z1 and zturb+1e-3 < z2) { + SMark_array(i,j,k,0) = 1.0; + } + } + }); } } @@ -571,12 +620,15 @@ WindFarm::fill_SMark_multifab (const Geometry& geom, int jj = amrex::min(amrex::max(j, j_lo), j_hi); int kk = amrex::min(amrex::max(k, k_lo), k_hi); + // The x and y extents of the current mesh cell + Real x1 = ProbLoArr[0] + ii*dx[0]; Real x2 = ProbLoArr[0] + (ii+1)*dx[0]; Real y1 = ProbLoArr[1] + jj*dx[1]; Real y2 = ProbLoArr[1] + (jj+1)*dx[1]; - //Real z = ProbLoArr[2] + (kk+0.5) * dx[2]; + // The mesh cell centered z value + Real z = (z_cc_arr) ? z_cc_arr(ii,jj,kk) : ProbLoArr[2] + (kk+0.5) * dx[2]; int turb_indices_overlap[2]; diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H index 68265a6db..71edd8acf 100644 --- a/Source/WindFarmParametrization/ERF_WindFarm.H +++ b/Source/WindFarmParametrization/ERF_WindFarm.H @@ -72,7 +72,7 @@ public: void fill_Nturb_multifab (const amrex::Geometry& geom, amrex::MultiFab& mf_Nturb, - std::unique_ptr& z_phys_cc); + std::unique_ptr& z_phys_nd); void fill_SMark_multifab (const amrex::Geometry& geom, amrex::MultiFab& mf_SMark, @@ -80,6 +80,11 @@ public: const amrex::Real& turb_disk_angle, std::unique_ptr& z_phys_cc); + void fill_SMark_multifab_mesoscale_models (const amrex::Geometry& geom, + amrex::MultiFab& mf_SMark, + const amrex::MultiFab& mf_Nturb, + std::unique_ptr& z_phys_cc); + void write_turbine_locations_vtk (); void write_actuator_disks_vtk (const amrex::Geometry& geom, diff --git a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp index 8ba23af86..8bcff6486 100644 --- a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp +++ b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp @@ -21,9 +21,75 @@ EWP::advance (const Geometry& geom, AMREX_ALWAYS_ASSERT(time > -1.0); source_terms_cellcentered(geom, cons_in, mf_vars_ewp, U_old, V_old, W_old, mf_Nturb); update(dt_advance, cons_in, U_old, V_old, mf_vars_ewp); + compute_power_output(cons_in, U_old, V_old, W_old, mf_SMark, mf_Nturb, time); } +void +EWP::compute_power_output (const MultiFab& cons_in, + const MultiFab& U_old, + const MultiFab& V_old, + const MultiFab& W_old, + const MultiFab& mf_SMark, + const MultiFab& mf_Nturb, + const Real& time) +{ + get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, + wind_speed, thrust_coeff, power); + + const int n_spec_table = wind_speed.size(); + + Gpu::DeviceVector d_wind_speed(wind_speed.size()); + Gpu::DeviceVector d_power(wind_speed.size()); + Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); + Gpu::copy(Gpu::hostToDevice, power.begin(), power.end(), d_power.begin()); + + Gpu::DeviceScalar d_total_power(0.0); + Real* d_total_power_ptr = d_total_power.dataPtr(); + + const Real* d_wind_speed_ptr = d_wind_speed.dataPtr(); + const Real* d_power_ptr = d_power.dataPtr(); + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + auto SMark_array = mf_SMark.array(mfi); + auto Nturb_array = mf_Nturb.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + auto w_vel = W_old.array(mfi); + Box tbx = mfi.nodaltilebox(0); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + if(SMark_array(i,j,k,0) == 1.0) { + Real avg_vel = 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,k)*w_vel(i,j,k),0.5); + Real turb_power = interpolate_1d(d_wind_speed_ptr, d_power_ptr, avg_vel, n_spec_table); + turb_power = turb_power*Nturb_array(i,j,k,0); + Gpu::Atomic::Add(d_total_power_ptr,turb_power); + } + }); + } + + Real h_total_power = 0.0; + Gpu::copy(Gpu::deviceToHost, d_total_power.dataPtr(), d_total_power.dataPtr()+1, &h_total_power); + + amrex::ParallelAllReduce::Sum(&h_total_power, 1, amrex::ParallelContext::CommunicatorAll()); + + if (ParallelDescriptor::IOProcessor()){ + static std::ofstream file("power_output_EWP.txt", std::ios::app); + // Check if the file opened successfully + if (!file.is_open()) { + std::cerr << "Error opening file!" << std::endl; + Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp"); + } + file << time << " " << h_total_power << "\n"; + file.flush(); + } +} + void EWP::update (const Real& dt_advance, MultiFab& cons_in, @@ -121,7 +187,7 @@ EWP::source_terms_cellcentered (const Geometry& geom, 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 K_turb = 6.0; Real L_wake = std::pow(dx[0]*dx[1],0.5)/2.0; Real sigma_e = Vabs/(3.0*K_turb*L_wake)* diff --git a/Source/WindFarmParametrization/EWP/ERF_EWP.H b/Source/WindFarmParametrization/EWP/ERF_EWP.H index 119590a27..f930704ed 100644 --- a/Source/WindFarmParametrization/EWP/ERF_EWP.H +++ b/Source/WindFarmParametrization/EWP/ERF_EWP.H @@ -35,9 +35,18 @@ public: void update (const amrex::Real& dt_advance, amrex::MultiFab& cons_in, - amrex::MultiFab& U_old, amrex::MultiFab& V_old, + amrex::MultiFab& U_old, + amrex::MultiFab& V_old, const amrex::MultiFab& mf_vars_ewp); + void compute_power_output (const amrex::MultiFab& cons_in, + const amrex::MultiFab& U_old, + const amrex::MultiFab& V_old, + const amrex::MultiFab& W_old, + const amrex::MultiFab& mf_SMark, + const amrex::MultiFab& mf_Nturb, + const amrex::Real& time); + protected: amrex::Vector xloc, yloc; amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp index 6f1036339..d96740fda 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp @@ -61,9 +61,9 @@ Fitch::advance (const Geometry& geom, AMREX_ALWAYS_ASSERT(time > -1.0); source_terms_cellcentered(geom, cons_in, mf_vars_fitch, U_old, V_old, W_old, mf_Nturb); update(dt_advance, cons_in, U_old, V_old, mf_vars_fitch); + compute_power_output(cons_in, U_old, V_old, mf_SMark, mf_Nturb, time); } - void Fitch::update (const Real& dt_advance, MultiFab& cons_in, @@ -98,6 +98,67 @@ Fitch::update (const Real& dt_advance, } } +void +Fitch::compute_power_output (const MultiFab& cons_in, + const MultiFab& U_old, + const MultiFab& V_old, + const MultiFab& mf_SMark, + const MultiFab& mf_Nturb, + const Real& time) +{ + get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, + wind_speed, thrust_coeff, power); + + const int n_spec_table = wind_speed.size(); + + Gpu::DeviceVector d_wind_speed(wind_speed.size()); + Gpu::DeviceVector d_power(wind_speed.size()); + Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); + Gpu::copy(Gpu::hostToDevice, power.begin(), power.end(), d_power.begin()); + + Gpu::DeviceScalar d_total_power(0.0); + Real* d_total_power_ptr = d_total_power.dataPtr(); + + const Real* d_wind_speed_ptr = d_wind_speed.dataPtr(); + const Real* d_power_ptr = d_power.dataPtr(); + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + auto SMark_array = mf_SMark.array(mfi); + auto Nturb_array = mf_Nturb.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + Box tbx = mfi.nodaltilebox(0); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + if(SMark_array(i,j,k,0) == 1.0) { + Real avg_vel = std::pow(u_vel(i,j,k)*u_vel(i,j,k) + v_vel(i,j,k)*v_vel(i,j,k),0.5); + Real turb_power = interpolate_1d(d_wind_speed_ptr, d_power_ptr, avg_vel, n_spec_table); + turb_power = turb_power*Nturb_array(i,j,k,0); + Gpu::Atomic::Add(d_total_power_ptr,turb_power); + } + }); + } + + Real h_total_power = 0.0; + Gpu::copy(Gpu::deviceToHost, d_total_power.dataPtr(), d_total_power.dataPtr()+1, &h_total_power); + + amrex::ParallelAllReduce::Sum(&h_total_power, 1, amrex::ParallelContext::CommunicatorAll()); + + if (ParallelDescriptor::IOProcessor()){ + static std::ofstream file("power_output_Fitch.txt", std::ios::app); + // Check if the file opened successfully + if (!file.is_open()) { + std::cerr << "Error opening file!" << std::endl; + Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp"); + } + file << time << " " << h_total_power << "\n"; + file.flush(); + } +} + void Fitch::source_terms_cellcentered (const Geometry& geom, const MultiFab& cons_in, diff --git a/Source/WindFarmParametrization/Fitch/ERF_Fitch.H b/Source/WindFarmParametrization/Fitch/ERF_Fitch.H index baa557e22..3476cfc6b 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_Fitch.H +++ b/Source/WindFarmParametrization/Fitch/ERF_Fitch.H @@ -38,7 +38,15 @@ public: amrex::MultiFab& U_old, amrex::MultiFab& V_old, const amrex::MultiFab& mf_vars_fitch); + void compute_power_output (const amrex::MultiFab& cons_in, + const amrex::MultiFab& U_old, + const amrex::MultiFab& V_old, + const amrex::MultiFab& mf_SMark, + const amrex::MultiFab& mf_Nturb, + const amrex::Real& time); + protected: + amrex::Vector hub_height_velocity; amrex::Vector xloc, yloc; amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; amrex::Vector wind_speed, thrust_coeff, power; diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp index 2f5873f63..acbdddfe9 100644 --- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp @@ -24,8 +24,38 @@ GeneralAD::advance (const Geometry& geom, compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark); source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD); update(dt_advance, cons_in, U_old, V_old, W_old, mf_vars_generalAD); + compute_power_output(time); } +void +GeneralAD::compute_power_output (const Real& time) +{ + get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, + wind_speed, thrust_coeff, power); + + const int n_spec_table = wind_speed.size(); + // Compute power based on the look-up table + + if (ParallelDescriptor::IOProcessor()){ + static std::ofstream file("power_output_GeneralAD.txt", std::ios::app); + // Check if the file opened successfully + if (!file.is_open()) { + std::cerr << "Error opening file!" << std::endl; + Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp"); + } + Real total_power = 0.0; + for(int it=0; it xloc, yloc; amrex::Real turb_disk_angle; diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp index 116de0b48..9938cbce7 100644 --- a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp +++ b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -35,7 +35,7 @@ SimpleAD::compute_power_output (const Real& time) // Compute power based on the look-up table if (ParallelDescriptor::IOProcessor()){ - static std::ofstream file("power_output.txt", std::ios::app); + static std::ofstream file("power_output_SimpleAD.txt", std::ios::app); // Check if the file opened successfully if (!file.is_open()) { std::cerr << "Error opening file!" << std::endl;