Skip to content

Commit

Permalink
Power output for all wind models
Browse files Browse the repository at this point in the history
  • Loading branch information
Mahesh Natarajan committed Dec 2, 2024
1 parent fc0f5a9 commit ffc17e3
Show file tree
Hide file tree
Showing 11 changed files with 270 additions and 27 deletions.
5 changes: 3 additions & 2 deletions Source/IO/ERF_Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
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]);
}
}
Expand Down Expand Up @@ -510,7 +510,8 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> 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();
Expand Down
11 changes: 10 additions & 1 deletion Source/Initialization/ERF_init_windfarm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand Down
78 changes: 65 additions & 13 deletions Source/WindFarmParametrization/ERF_InitWindFarm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -352,7 +352,7 @@ WindFarm::read_windfarm_airfoil_tables (const std::string windfarm_airfoil_table

void
WindFarm::gatherKeyValuePairs(const std::vector<std::pair<int, double>>& localData,
std::vector<std::pair<int, double>>& globalData)
std::vector<std::pair<int, double>>& globalData)
{
int myRank = amrex::ParallelDescriptor::MyProc();
int nProcs = amrex::ParallelDescriptor::NProcs();
Expand Down Expand Up @@ -419,17 +419,13 @@ WindFarm::gatherKeyValuePairs(const std::vector<std::pair<int, double>>& 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<MultiFab>& z_phys_cc)
std::unique_ptr<MultiFab>& z_phys_nd)
{

zloc.resize(xloc.size(),-1.0);
Expand Down Expand Up @@ -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<const Real>& z_cc_arr = (z_phys_cc) ? z_phys_cc->const_array(mfi) : Array4<Real>{};
const Array4<const Real>& z_nd_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4<Real>{};
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);
Expand All @@ -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 {
Expand All @@ -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<std::pair<int, double>> turb_index_zloc;
for(int it=0;it<xloc.size();it++){
if(turb_index[it] != -1) {
Expand All @@ -513,9 +507,64 @@ WindFarm::fill_Nturb_multifab (const Geometry& geom,

// Each process now has the global array
for (const auto& kv : turb_index_zloc_glob) {
std::cout << "Rank " << amrex::ParallelDescriptor::MyProc() << "Global data" << kv.first << " " << kv.second << "\n";
//std::cout << "Rank " << amrex::ParallelDescriptor::MyProc() << "Global data" << kv.first << " " << kv.second << "\n";
zloc[kv.first] = kv.second;
}
}

void
WindFarm::fill_SMark_multifab_mesoscale_models (const Geometry& geom,
MultiFab& mf_SMark,
const MultiFab& mf_Nturb,
std::unique_ptr<MultiFab>& z_phys_nd)
{
mf_SMark.setVal(-1.0);

Real d_hub_height = hub_height;

amrex::Gpu::DeviceVector<Real> d_xloc(xloc.size());
amrex::Gpu::DeviceVector<Real> d_yloc(yloc.size());
amrex::Gpu::DeviceVector<Real> 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<const Real>& z_nd_arr = (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4<Real>{};
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;
}
}
});
}
}

Expand Down Expand Up @@ -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];
Expand Down
7 changes: 6 additions & 1 deletion Source/WindFarmParametrization/ERF_WindFarm.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,14 +72,19 @@ public:

void fill_Nturb_multifab (const amrex::Geometry& geom,
amrex::MultiFab& mf_Nturb,
std::unique_ptr<amrex::MultiFab>& z_phys_cc);
std::unique_ptr<amrex::MultiFab>& z_phys_nd);

void fill_SMark_multifab (const amrex::Geometry& geom,
amrex::MultiFab& mf_SMark,
const amrex::Real& sampling_distance_by_D,
const amrex::Real& turb_disk_angle,
std::unique_ptr<amrex::MultiFab>& 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<amrex::MultiFab>& z_phys_cc);

void write_turbine_locations_vtk ();

void write_actuator_disks_vtk (const amrex::Geometry& geom,
Expand Down
68 changes: 67 additions & 1 deletion Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real> d_wind_speed(wind_speed.size());
Gpu::DeviceVector<Real> 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<Real> 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,
Expand Down Expand Up @@ -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)*
Expand Down
11 changes: 10 additions & 1 deletion Source/WindFarmParametrization/EWP/ERF_EWP.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::Real> xloc, yloc;
amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power;
Expand Down
63 changes: 62 additions & 1 deletion Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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<Real> d_wind_speed(wind_speed.size());
Gpu::DeviceVector<Real> 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<Real> 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,
Expand Down
Loading

0 comments on commit ffc17e3

Please sign in to comment.