Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Power output for all wind models #1991

Merged
merged 1 commit into from
Dec 3, 2024
Merged

Power output for all wind models #1991

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions Source/IO/ERF_Plotfile.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading