ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
WindFarm Class Reference

Manage a collection of wind turbine models. More...

#include <ERF_WindFarm.H>

Inheritance diagram for WindFarm:
Collaboration diagram for WindFarm:

Public Member Functions

 WindFarm ()
 
virtual ~WindFarm ()=default
 
 WindFarm (int nlev, const WindFarmType &a_windfarm_type)
 
void read_tables (std::string windfarm_loc_table, std::string windfarm_spec_table, bool x_y, bool lat_lon, const amrex::Real windfarm_x_shift=zero, const amrex::Real windfarm_y_shift=zero)
 Read turbine location and specification tables. More...
 
void init_windfarm_lat_lon (const std::string windfarm_loc_table, const amrex::Real windfarm_x_shift, const amrex::Real windfarm_y_shift)
 Initialize wind farm using latitude and longitude coordinates. More...
 
void init_windfarm_x_y (const std::string windfarm_loc_table)
 Initialize wind farm using Cartesian coordinates. More...
 
void read_windfarm_locations_table (const std::string windfarm_loc_table, bool x_y, bool lat_lon, const amrex::Real windfarm_x_shift=zero, const amrex::Real windfarm_y_shift=zero)
 Read the turbine location table. More...
 
void read_windfarm_spec_table (const std::string windfarm_spec_table)
 Read the turbine specification table. More...
 
void read_windfarm_blade_table (const std::string windfarm_blade_table)
 Read the turbine blade table. More...
 
void read_windfarm_airfoil_tables (const std::string windfarm_airfoil_tables, const std::string windfarm_blade_table)
 Read the airfoil tables. More...
 
void read_windfarm_spec_table_extra (const std::string windfarm_spec_table_extra)
 Read extra turbine specifications. More...
 
void fill_Nturb_multifab (const amrex::Geometry &geom, amrex::MultiFab &mf_Nturb, std::unique_ptr< amrex::MultiFab > &z_phys_nd)
 Fill the turbine count MultiFab. More...
 
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)
 Fill the SMark MultiFab. More...
 
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)
 Fill the SMark MultiFab for mesoscale models. More...
 
void write_turbine_locations_vtk ()
 Write turbine locations to a VTK file. More...
 
void write_actuator_disks_vtk (const amrex::Geometry &geom, const amrex::Real &sampling_distance_by_D)
 Write actuator disks to a VTK file. More...
 
void advance (const amrex::Geometry &a_geom, const double &dt_advance, amrex::MultiFab &cons_in, amrex::MultiFab &mf_vars_windfarm, amrex::MultiFab &U_old, amrex::MultiFab &V_old, amrex::MultiFab &W_old, const amrex::MultiFab &mf_Nturb, const amrex::MultiFab &mf_SMark, const double &time) override
 Advance the wind farm model in time. More...
 
void set_turb_spec (const amrex::Real &a_rotor_rad, const amrex::Real &a_hub_height, const amrex::Real &a_thrust_coeff_standing, const amrex::Vector< amrex::Real > &a_wind_speed, const amrex::Vector< amrex::Real > &a_thrust_coeff, const amrex::Vector< amrex::Real > &a_power) override
 Set general turbine specifications. More...
 
void set_turb_loc (const amrex::Vector< amrex::Real > &a_xloc, const amrex::Vector< amrex::Real > &a_yloc) override
 Set turbine locations. More...
 
void set_turb_disk_angle (const amrex::Real &a_turb_disk_angle) override
 Set the turbine disk angle. More...
 
void set_blade_spec (const amrex::Vector< amrex::Real > &a_bld_rad_loc, const amrex::Vector< amrex::Real > &a_bld_twist, const amrex::Vector< amrex::Real > &a_bld_chord) override
 Set blade specifications. More...
 
void set_blade_airfoil_spec (const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_aoa, const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_Cl, const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_Cd) override
 Set blade airfoil specifications. More...
 
void set_turb_spec_extra (const amrex::Vector< amrex::Real > &a_velocity, const amrex::Vector< amrex::Real > &a_C_P, const amrex::Vector< amrex::Real > &a_C_T, const amrex::Vector< amrex::Real > &a_rotor_RPM, const amrex::Vector< amrex::Real > &a_blade_pitch) override
 Set extra turbine specifications. More...
 
- Public Member Functions inherited from NullWindFarm
 NullWindFarm ()
 
virtual ~NullWindFarm ()=default
 
void get_turb_spec (amrex::Real &rotor_rad, amrex::Real &hub_height, amrex::Real &thrust_coeff_standing, amrex::Vector< amrex::Real > &wind_speed, amrex::Vector< amrex::Real > &thrust_coeff, amrex::Vector< amrex::Real > &power)
 
void get_turb_loc (amrex::Vector< amrex::Real > &xloc, amrex::Vector< amrex::Real > &yloc)
 
void get_turb_disk_angle (amrex::Real &turb_disk_angle)
 
void get_blade_spec (amrex::Vector< amrex::Real > &bld_rad_loc, amrex::Vector< amrex::Real > &bld_twist, amrex::Vector< amrex::Real > &bld_chord)
 
void get_blade_airfoil_spec (amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_aoa, amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_Cl, amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_Cd)
 
void get_turb_spec_extra (amrex::Vector< amrex::Real > &velocity, amrex::Vector< amrex::Real > &C_P, amrex::Vector< amrex::Real > &C_T, amrex::Vector< amrex::Real > &rotor_RPM, amrex::Vector< amrex::Real > &blade_pitch)
 

Protected Member Functions

template<class NewWindFarmModel >
void SetModel ()
 Create and set the specified windfarm model. More...
 

Protected Attributes

amrex::Vector< amrex::Real > xloc
 
amrex::Vector< amrex::Real > yloc
 
amrex::Vector< amrex::Real > zloc
 
amrex::Real my_turb_disk_angle
 
amrex::Real hub_height
 
amrex::Real rotor_rad
 
amrex::Real thrust_coeff_standing
 
amrex::Real nominal_power
 
amrex::Vector< amrex::Real > wind_speed
 
amrex::Vector< amrex::Real > thrust_coeff
 
amrex::Vector< amrex::Real > power
 
amrex::Vector< amrex::Real > bld_rad_loc
 
amrex::Vector< amrex::Real > bld_twist
 
amrex::Vector< amrex::Real > bld_chord
 
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_aoa
 
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cl
 
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cd
 
int n_bld_sections
 
amrex::Vector< amrex::Real > velocity
 
amrex::Vector< amrex::Real > C_P
 
amrex::Vector< amrex::Real > C_T
 
amrex::Vector< amrex::Real > rotor_RPM
 
amrex::Vector< amrex::Real > blade_pitch
 
- Protected Attributes inherited from NullWindFarm
amrex::Vector< amrex::Real > m_xloc
 
amrex::Vector< amrex::Real > m_yloc
 
amrex::Real m_turb_disk_angle
 
amrex::Real m_hub_height
 
amrex::Real m_rotor_rad
 
amrex::Real m_thrust_coeff_standing
 
amrex::Real m_nominal_power
 
amrex::Vector< amrex::Real > m_wind_speed
 
amrex::Vector< amrex::Real > m_thrust_coeff
 
amrex::Vector< amrex::Real > m_power
 
amrex::Vector< amrex::Real > m_bld_rad_loc
 
amrex::Vector< amrex::Real > m_bld_twist
 
amrex::Vector< amrex::Real > m_bld_chord
 
amrex::Vector< amrex::Vector< amrex::Real > > m_bld_airfoil_aoa
 
amrex::Vector< amrex::Vector< amrex::Real > > m_bld_airfoil_Cl
 
amrex::Vector< amrex::Vector< amrex::Real > > m_bld_airfoil_Cd
 
amrex::Vector< amrex::Real > m_velocity
 
amrex::Vector< amrex::Real > m_C_P
 
amrex::Vector< amrex::Real > m_C_T
 
amrex::Vector< amrex::Real > m_rotor_RPM
 
amrex::Vector< amrex::Real > m_blade_pitch
 

Private Attributes

amrex::Vector< std::unique_ptr< NullWindFarm > > m_windfarm_model
 

Additional Inherited Members

- Static Public Member Functions inherited from NullWindFarm
static AMREX_GPU_DEVICE bool find_if_marked (amrex::Real x1, amrex::Real x2, amrex::Real y1, amrex::Real y2, amrex::Real x0, amrex::Real y0, amrex::Real nx, amrex::Real ny, amrex::Real d_hub_height, amrex::Real d_rotor_rad, amrex::Real z)
 

Detailed Description

Manage a collection of wind turbine models.

Constructor & Destructor Documentation

◆ WindFarm() [1/2]

WindFarm::WindFarm ( )
inline
21 {}

◆ ~WindFarm()

virtual WindFarm::~WindFarm ( )
virtualdefault

◆ WindFarm() [2/2]

WindFarm::WindFarm ( int  nlev,
const WindFarmType &  a_windfarm_type 
)
inline
27  {
28  m_windfarm_model.resize(nlev);
29  if (a_windfarm_type == WindFarmType::Fitch) {
30  SetModel<Fitch>();
31  amrex::Print() << "Fitch windfarm model!\n";
32  }
33  else if (a_windfarm_type == WindFarmType::EWP) {
34  SetModel<EWP>();
35  amrex::Print() << "EWP windfarm model!\n";
36  }
37  else if (a_windfarm_type == WindFarmType::SimpleAD) {
38  SetModel<SimpleAD>();
39  amrex::Print() << "Simplified actuator disk windfarm model!\n";
40  }
41  else if (a_windfarm_type == WindFarmType::GeneralAD) {
42  SetModel<GeneralAD>();
43  amrex::Print() << "Generalized actuator disk windfarm model!\n";
44  }
45  else {
46  amrex::Abort("WindFarm: Dont know this windfarm_type!") ;
47  }
48  }
amrex::Vector< std::unique_ptr< NullWindFarm > > m_windfarm_model
Definition: ERF_WindFarm.H:299

Member Function Documentation

◆ advance()

void WindFarm::advance ( const amrex::Geometry &  a_geom,
const double &  dt_advance,
amrex::MultiFab &  cons_in,
amrex::MultiFab &  mf_vars_windfarm,
amrex::MultiFab &  U_old,
amrex::MultiFab &  V_old,
amrex::MultiFab &  W_old,
const amrex::MultiFab &  mf_Nturb,
const amrex::MultiFab &  mf_SMark,
const double &  time 
)
inlineoverridevirtual

Advance the wind farm model in time.

Parameters
[in]a_geomGeometry
[in]dt_advanceTimestep.
[in,out]cons_inConservative variables.
[in,out]mf_vars_windfarmWind farm variables MultiFab.
[in,out]U_oldOld x-velocity MultiFab.
[in,out]V_oldOld y-velocity MultiFab.
[in,out]W_oldOld z-velocity MultiFab.
[in]mf_NturbTurbine count MultiFab.
[in]mf_SMarkSMark MultiFab.
[in]timeCurrent simulation time.

Implements NullWindFarm.

192  {
193  m_windfarm_model[0]->advance(a_geom, dt_advance, cons_in, mf_vars_windfarm,
194  U_old, V_old, W_old, mf_Nturb, mf_SMark, time);
195  }

◆ fill_Nturb_multifab()

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

Fill the turbine count MultiFab.

Parameters
[in]geomGeometry
[out]mf_NturbMultiFab to fill with turbine counts.
[in,out]z_phys_ndPhysical height MultiFab.

Fill a MultiFab with the number of turbines per cell and determine base heights.

Parameters
[in]geomGeometry defining the domain.
[in,out]mf_NturbMultiFab to be filled with turbine counts per cell.
[in]z_phys_ndNode-centered physical height FAB.
447 {
448 
449  zloc.resize(xloc.size(),zero);
450  Vector<int> is_counted;
451  is_counted.resize(xloc.size(),0);
452 
453  amrex::Gpu::DeviceVector<Real> d_xloc(xloc.size());
454  amrex::Gpu::DeviceVector<Real> d_yloc(yloc.size());
455  amrex::Gpu::DeviceVector<Real> d_zloc(xloc.size());
456  amrex::Gpu::DeviceVector<int> d_is_counted(xloc.size());
457  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
458  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
459  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin());
460  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, is_counted.begin(), is_counted.end(), d_is_counted.begin());
461 
462  Real* d_xloc_ptr = d_xloc.data();
463  Real* d_yloc_ptr = d_yloc.data();
464  Real* d_zloc_ptr = d_zloc.data();
465  int* d_is_counted_ptr = d_is_counted.data();
466 
467  mf_Nturb.setVal(0);
468 
469  int i_lo = geom.Domain().smallEnd(0); int i_hi = geom.Domain().bigEnd(0);
470  int j_lo = geom.Domain().smallEnd(1); int j_hi = geom.Domain().bigEnd(1);
471  auto dx = geom.CellSizeArray();
472  if(dx[0]<= 1e-3 or dx[1]<=1e-3 or dx[2]<= 1e-3) {
473  Abort("The value of grid spacing for wind farm parametrization cannot be less than 1e-3 m. "
474  "It should be usually of order 1 m");
475  }
476  auto ProbLoArr = geom.ProbLoArray();
477  auto ProbHiArr = geom.ProbHiArray();
478  int num_turb = xloc.size();
479 
480  bool is_terrain = z_phys_nd ? true: false;
481 
482  // The ground is at the bottom of the domain, not at the bottom of each box --
483  // these differ for every box that doesn't touch the ground, which happens
484  // whenever the grids are decomposed in the vertical direction
485  const int k0 = geom.Domain().smallEnd(2);
486 
487  // Initialize wind farm
488  for ( MFIter mfi(mf_Nturb,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
489  const Box& bx = mfi.tilebox();
490  auto Nturb_array = mf_Nturb.array(mfi);
491  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
492  ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
493  int li = amrex::min(amrex::max(i, i_lo), i_hi);
494  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
495 
496  Real x1 = ProbLoArr[0] + li*dx[0];
497  Real x2 = ProbLoArr[0] + (li+1)*dx[0];
498  Real y1 = ProbLoArr[1] + lj*dx[1];
499  Real y2 = ProbLoArr[1] + (lj+1)*dx[1];
500 
501  for(int it=0; it<num_turb; it++){
502  if( d_xloc_ptr[it]+1e-3 > x1 and d_xloc_ptr[it]+1e-3 < x2 and
503  d_yloc_ptr[it]+1e-3 > y1 and d_yloc_ptr[it]+1e-3 < y2){
504  Nturb_array(i,j,k,0) = Nturb_array(i,j,k,0) + 1;
505  //
506  // Only the cell at the bottom of the domain counts this turbine and
507  // stores its surface elevation. d_is_counted is a per-rank array
508  // which is summed over all ranks below, so if we counted at every k
509  // then a column spread over more than one rank (which happens when
510  // the grids are decomposed in z) would be counted more than once.
511  //
512  if (is_terrain and k == k0) {
513  int expected = 0;
514  int desired = 1;
515  // Atomic Compare-And-Swap: Increment only if d_is_counted_ptr[it] was 0
516  if (Gpu::Atomic::CAS(&d_is_counted_ptr[it], expected, desired) == expected) {
517  // The current thread successfully set d_is_counted_ptr[it] from 0 to 1
518  Gpu::Atomic::Add(&d_zloc_ptr[it], z_nd_arr(i, j, k0));
519  }
520  }
521  }
522  }
523  });
524  }
525 
526  mf_Nturb.FillBoundary(geom.periodicity());
527 
528  Gpu::copy(Gpu::deviceToHost, d_zloc.begin(), d_zloc.end(), zloc.begin());
529  Gpu::copy(Gpu::deviceToHost, d_is_counted.begin(), d_is_counted.end(), is_counted.begin());
530 
531  amrex::ParallelAllReduce::Sum(zloc.data(),
532  zloc.size(),
533  amrex::ParallelContext::CommunicatorAll());
534 
535  amrex::ParallelAllReduce::Sum(is_counted.data(),
536  is_counted.size(),
537  amrex::ParallelContext::CommunicatorAll());
538 
539  for(int it=0;it<num_turb;it++) {
540  if(is_terrain and
541  xloc[it] > ProbLoArr[0] and
542  xloc[it] < ProbHiArr[0] and
543  yloc[it] > ProbLoArr[1] and
544  yloc[it] < ProbHiArr[1] ) {
545  if(is_counted[it] != 1) {
546  Abort("Wind turbine " + std::to_string(it) + "has been counted " + std::to_string(is_counted[it]) + " times" +
547  " It should have been counted only once. Aborting....");
548  }
549  }
550  }
551 
552  // Debugging
553  /*int my_rank = amrex::ParallelDescriptor::MyProc();
554 
555  for(int it=0;it<num_turb;it++) {
556  std::cout << "The value of zloc is " << my_rank << " " << zloc[it] << " " << is_counted[it] << "\n";
557  }*/
558 }
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
ParallelFor(fab_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:36
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Vector< amrex::Real > yloc
Definition: ERF_WindFarm.H:280
amrex::Vector< amrex::Real > zloc
Definition: ERF_WindFarm.H:280
amrex::Vector< amrex::Real > xloc
Definition: ERF_WindFarm.H:280
Here is the call graph for this function:

◆ fill_SMark_multifab()

void WindFarm::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 
)

Fill the SMark MultiFab.

Parameters
[in]geomGeometry
[out]mf_SMarkMultiFab to fill with SMark values.
[in]sampling_distance_by_DSampling distance normalized by rotor diameter.
[in]turb_disk_angleTurbine disk angle.
[in,out]z_phys_ccPhysical height MultiFab.

Fill the SMark MultiFab identifying turbine disks and sampling points.

Parameters
[in]geomGeometry defining the domain.
[in,out]mf_SMarkMultiFab to be filled with turbine markers.
[in]sampling_distance_by_DNormalized sampling distance.
[in]turb_disk_angleAngle of the turbine disk in degrees.
[in]z_phys_ccCell-centered physical height FAB.
635 {
636  // Copy turbine locations to device
637  Gpu::DeviceVector<Real> d_xloc(xloc.size());
638  Gpu::DeviceVector<Real> d_yloc(yloc.size());
639  Gpu::DeviceVector<Real> d_zloc(zloc.size());
640 
641  Gpu::copy(Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
642 
643  Gpu::copy(Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
644 
645  Gpu::copy(Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin());
646 
647  Real d_rotor_rad = rotor_rad;
648  Real d_hub_height = hub_height;
649  Real d_sampling_distance = sampling_distance_by_D * two * rotor_rad;
650 
651  Real* d_xloc_ptr = d_xloc.data();
652  Real* d_yloc_ptr = d_yloc.data();
653  Real* d_zloc_ptr = d_zloc.data();
654 
655  mf_SMark.setVal(-1.0);
656 
657  const int i_lo = geom.Domain().smallEnd(0);
658  const int i_hi = geom.Domain().bigEnd(0);
659  const int j_lo = geom.Domain().smallEnd(1);
660  const int j_hi = geom.Domain().bigEnd(1);
661  const int k_lo = geom.Domain().smallEnd(2);
662  const int k_hi = geom.Domain().bigEnd(2);
663 
664  auto dx = geom.CellSizeArray();
665  auto ProbLoArr = geom.ProbLoArray();
666 
667  const int num_turb = static_cast<int>(xloc.size());
668 
669  Real theta = turb_disk_angle * M_PI / Real(180.0) - myhalf * M_PI;
670 
673 
674  Real nx = -std::cos(theta);
675  Real ny = -std::sin(theta);
676 
677  // ------------------------------------------------------------
678  // Device-side overlap information
679  //
680  // overlap_info[0] = overlap flag
681  // overlap_info[1] = first turbine index
682  // overlap_info[2] = second turbine index
683  // ------------------------------------------------------------
684 
685  int h_overlap_info[3] = {0, -1, -1};
686 
687  Gpu::DeviceVector<int> d_overlap_info(3);
688 
689  Gpu::copy(Gpu::hostToDevice, h_overlap_info, h_overlap_info + 3, d_overlap_info.begin());
690 
691  int* d_overlap_info_ptr = d_overlap_info.data();
692 
693  // ------------------------------------------------------------
694  // Initialize wind farm
695  // ------------------------------------------------------------
696 
697  for (MFIter mfi(mf_SMark, TilingIfNotGPU()); mfi.isValid(); ++mfi)
698  {
699  const Box& gbx = mfi.growntilebox(1);
700 
701  auto SMark_array = mf_SMark.array(mfi);
702 
703  const Array4<const Real>& z_cc_arr = z_phys_cc->const_array(mfi);
704 
705  ParallelFor(
706  gbx,
707  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
708  {
709  int ii = amrex::min(amrex::max(i, i_lo), i_hi);
710  int jj = amrex::min(amrex::max(j, j_lo), j_hi);
711  int kk = amrex::min(amrex::max(k, k_lo), k_hi);
712 
713  // ------------------------------------------------
714  // Mesh-cell x/y extents
715  // ------------------------------------------------
716 
717  Real x1 = ProbLoArr[0] + ii * dx[0];
718  Real x2 = ProbLoArr[0] + (ii + 1) * dx[0];
719 
720  Real y1 = ProbLoArr[1] + jj * dx[1];
721  Real y2 = ProbLoArr[1] + (jj + 1) * dx[1];
722 
723  // ------------------------------------------------
724  // Cell-centered physical z
725  // ------------------------------------------------
726 
727  Real z = z_cc_arr(ii, jj, kk);
728 
729  int check_int = 0;
730  int first_turb = -1;
731 
732  for (int it = 0; it < num_turb; ++it)
733  {
734  Real z0 = d_zloc_ptr[it];
735 
736  // ------------------------------------------------
737  // Sampling point upstream/downstream of turbine
738  // ------------------------------------------------
739 
740  Real x0 = d_xloc_ptr[it] + d_sampling_distance * nx;
741 
742  Real y0 = d_yloc_ptr[it] + d_sampling_distance * ny;
743 
744  bool is_cell_marked =
746  x1, x2, y1, y2,
747  x0, y0,
748  nx, ny,
749  d_hub_height + z0,
750  d_rotor_rad,
751  z);
752 
753  if (is_cell_marked)
754  {
755  SMark_array(i, j, k, 0) = it;
756  }
757 
758  // ------------------------------------------------
759  // Actual turbine disk
760  // ------------------------------------------------
761 
762  x0 = d_xloc_ptr[it];
763  y0 = d_yloc_ptr[it];
764 
765  is_cell_marked =
767  x1, x2, y1, y2,
768  x0, y0,
769  nx, ny,
770  d_hub_height + z0,
771  d_rotor_rad,
772  z);
773 
774  if (is_cell_marked)
775  {
776  SMark_array(i, j, k, 1) = it;
777 
778  if (check_int == 0)
779  {
780  first_turb = it;
781  }
782 
783  ++check_int;
784 
785  // ------------------------------------------------
786  // More than one turbine disk marks this cell
787  // ------------------------------------------------
788 
789  if (check_int > 1)
790  {
791  // Set overlap flag.
792  amrex::Gpu::Atomic::Exch(&d_overlap_info_ptr[0], 1);
793 
794  // Save the turbine indices. Multiple threads
795  // may write these, so use atomics.
796  amrex::Gpu::Atomic::Exch(&d_overlap_info_ptr[1], first_turb);
797 
798  amrex::Gpu::Atomic::Exch(&d_overlap_info_ptr[2], it);
799  }
800  }
801  }
802  });
803  }
804 
805  // ------------------------------------------------------------
806  // Copy overlap information back to host.
807  // This also synchronizes the GPU work.
808  // ------------------------------------------------------------
809 
810  Gpu::copy(Gpu::deviceToHost, d_overlap_info.begin(), d_overlap_info.end(), h_overlap_info);
811 
812  // ------------------------------------------------------------
813  // Handle overlap on the HOST, not inside the GPU kernel.
814  // ------------------------------------------------------------
815 
816  if (h_overlap_info[0] != 0)
817  {
818  amrex::Error(
819  "Actuator disks with indices " + std::to_string(h_overlap_info[1]) +
820  " and " + std::to_string(h_overlap_info[2]) + " are overlapping. "
821  "Visualize actuator_disks.vtk and check the "
822  "windturbine locations input file. Exiting..");
823  }
824 }
#define M_PI
Definition: ERF_HurricaneDiagnostics.cpp:27
const int nx
Definition: ERF_InitCustomPertVels_CloudChamber.H:14
const int ny
Definition: ERF_InitCustomPertVels_CloudChamber.H:15
Real z0
Definition: ERF_InitCustomPertVels_ScalarAdvDiff.H:8
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:38
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:41
static AMREX_GPU_DEVICE bool find_if_marked(amrex::Real x1, amrex::Real x2, amrex::Real y1, amrex::Real y2, amrex::Real x0, amrex::Real y0, amrex::Real nx, amrex::Real ny, amrex::Real d_hub_height, amrex::Real d_rotor_rad, amrex::Real z)
Definition: ERF_NullWindFarm.H:140
amrex::Real my_turb_disk_angle
Definition: ERF_WindFarm.H:281
amrex::Real hub_height
Definition: ERF_WindFarm.H:282
amrex::Real rotor_rad
Definition: ERF_WindFarm.H:282
void set_turb_disk_angle(const amrex::Real &a_turb_disk_angle) override
Set the turbine disk angle.
Definition: ERF_WindFarm.H:230
@ theta
Definition: ERF_Kessler.H:26
Here is the call graph for this function:

◆ fill_SMark_multifab_mesoscale_models()

void WindFarm::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 
)

Fill the SMark MultiFab for mesoscale models.

Parameters
[in]geomGeometry
[out]mf_SMarkMultiFab to fill with SMark values.
[in]mf_NturbMultiFab containing turbine counts.
[in,out]z_phys_ccPhysical height MultiFab.

Fill the SMark MultiFab for mesoscale model parametrization.

Parameters
[in]geomGeometry defining the domain.
[in,out]mf_SMarkMultiFab to be filled with turbine markers.
[in]mf_NturbMultiFab containing turbine counts per cell.
[in]z_phys_ndNode-centered physical height FAB.
573 {
574  mf_SMark.setVal(-1.0);
575 
576  Real d_hub_height = hub_height;
577 
578  amrex::Gpu::DeviceVector<Real> d_xloc(xloc.size());
579  amrex::Gpu::DeviceVector<Real> d_yloc(yloc.size());
580  amrex::Gpu::DeviceVector<Real> d_zloc(xloc.size());
581  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
582  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
583  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zloc.begin(), zloc.end(), d_zloc.begin());
584 
585  int i_lo = geom.Domain().smallEnd(0); int i_hi = geom.Domain().bigEnd(0);
586  int j_lo = geom.Domain().smallEnd(1); int j_hi = geom.Domain().bigEnd(1);
587  int k_lo = geom.Domain().smallEnd(2); int k_hi = geom.Domain().bigEnd(2);
588 
589  //auto dx = geom.CellSizeArray();
590  //auto ProbLoArr = geom.ProbLoArray();
591 
592  // Initialize wind farm
593  for ( MFIter mfi(mf_SMark,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
594 
595  const Box& gbx = mfi.growntilebox(1);
596  auto SMark_array = mf_SMark.array(mfi);
597  auto Nturb_array = mf_Nturb.array(mfi);
598  const Array4<const Real>& z_nd_arr = z_phys_nd->const_array(mfi);
599  int k0 = std::max(gbx.smallEnd(2), geom.Domain().smallEnd(2));
600 
601  ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
602  if(Nturb_array(i,j,k,0) > 0) {
603  int li = amrex::min(amrex::max(i, i_lo), i_hi);
604  int lj = amrex::min(amrex::max(j, j_lo), j_hi);
605  int lk = amrex::min(amrex::max(k, k_lo), k_hi);
606 
607  Real z1 = z_nd_arr(li,lj,lk );
608  Real z2 = z_nd_arr(li,lj,lk+1);
609 
610  Real zturb = z_nd_arr(li,lj,k0) + d_hub_height;
611  if (zturb+1e-3 > z1 and zturb+1e-3 < z2) {
612  SMark_array(i,j,k,0) = one;
613  }
614  }
615  });
616  }
617 }
constexpr amrex::Real one
Definition: ERF_NumericalConstants.H:37
Here is the call graph for this function:

◆ init_windfarm_lat_lon()

void WindFarm::init_windfarm_lat_lon ( const std::string  windfarm_loc_table,
const amrex::Real  windfarm_x_shift,
const amrex::Real  windfarm_y_shift 
)

Initialize wind farm using latitude and longitude coordinates.

Parameters
[in]windfarm_loc_tablePath to the location table.
[in]windfarm_x_shiftX-coordinate shift.
[in]windfarm_y_shiftY-coordinate shift.

Initialize turbine locations from latitude and longitude coordinates.

Parameters
[in]windfarm_loc_tablePath to the turbine locations table.
[in]windfarm_x_shiftShift applied to the resulting x coordinates.
[in]windfarm_y_shiftShift applied to the resulting y coordinates.
74 {
75  // Re-reading the file must not append to a previous level's list.
76  xloc.clear();
77  yloc.clear();
78 
79  // Read turbine locations from windturbines.txt
80  std::ifstream file(windfarm_loc_table);
81  if (!file.is_open()) {
82  amrex::Error("Wind turbines location table not found. Either the inputs is missing the"
83  " erf.windfarm_loc_table entry or the file specified in the entry " + windfarm_loc_table + " is missing.");
84  }
85  // Vector of vectors to store the matrix
86  Vector<Real> lat, lon;
87  Real value1, value2, value3;
88 
89  while (file >> value1 >> value2 >> value3) {
90 
91  if(std::fabs(value1) > Real(90.0)) {
92  amrex::Error("The value of latitude for entry " + std::to_string(lat.size() + 1) +
93  " in " + windfarm_loc_table + " should be within -90 and 90");
94  }
95 
96  if(std::fabs(value2) > Real(180.0)) {
97  amrex::Error("The value of longitude for entry " + std::to_string(lat.size() + 1) +
98  " in " + windfarm_loc_table + " should be within -180 and 180");
99  }
100  lat.push_back(value1);
101  lon.push_back(value2);
102  }
103  file.close();
104 
105  Real rad_earth = Real(6371.0e3); // Radius of the earth
106  Real m_per_deg_lat = rad_earth*two*M_PI/(two*Real(180.0));
107 
108  // Find the coordinates of average of min and max of the farm
109  // Rotate about that point
110  ParmParse pp("erf");
111  std::string fname_usgs;
112  auto valid_fname_USGS = pp.queryAdd("terrain_file_name_USGS",fname_usgs);
113  Real lon_ref, lat_ref;
114 
115  if (valid_fname_USGS) {
116  std::ifstream file_usgs(fname_usgs);
117  if (!file_usgs.is_open()) {
118  amrex::Error("Could not open the file " + fname_usgs + " specified by the inputs file entry"
119  " erf.terrain_file_name_USGS. Either correct the entry or remove it, in which"
120  " case the reference point is taken from the turbine locations.");
121  }
122  file_usgs >> lon_ref >> lat_ref;
123  if (file_usgs.fail()) {
124  amrex::Error("Could not read the reference longitude and latitude from the first line of "
125  + fname_usgs + ", specified by the inputs file entry erf.terrain_file_name_USGS.");
126  }
127  file_usgs.close();
128  lon_ref = lon_ref*M_PI/Real(180.0);
129  lat_ref = lat_ref*M_PI/Real(180.0);
130  } else {
131  Real lat_min = *std::min_element(lat.begin(), lat.end());
132  Real lon_min = *std::min_element(lon.begin(), lon.end());
133 
134  lon_ref = lon_min*M_PI/Real(180.0);
135  lat_ref = lat_min*M_PI/Real(180.0);
136  }
137 
138 
139  for(int it=0;it<lat.size();it++){
140  lat[it] = lat[it]*M_PI/Real(180.0);
141  lon[it] = lon[it]*M_PI/Real(180.0);
142  Real delta_lat = (lat[it] - lat_ref);
143  Real delta_lon = (lon[it] - lon_ref);
144 
145  Real term1 = amrex::Math::powi<2>(sin(delta_lat/two));
146  Real term2 = cos(lat[it])*cos(lat_ref)*amrex::Math::powi<2>(sin(delta_lon/two));
147  Real dist = two*rad_earth*std::asin(std::sqrt(term1 + term2));
148  Real dy_turb = delta_lat * m_per_deg_lat * Real(180.0)/M_PI ;
149 
150  if(dist<dy_turb){
151  if(std::fabs(dist-dy_turb)<1e-8){
152  dist=dy_turb;
153  }
154  else{
155  Abort("The value of dist is less than dy_turb "+ std::to_string(dist) + " " + std::to_string(dy_turb));
156  }
157  }
158  Real tmp = amrex::Math::powi<2>(dist) - amrex::Math::powi<2>(dy_turb);
159 
160  if(std::fabs(tmp)<1e-8){
161  tmp = zero;
162  }
163  Real dx_turb = std::sqrt(tmp);
164 
165 
166  if(delta_lon >= zero) {
167  xloc.push_back(dx_turb);
168  }
169  else {
170  xloc.push_back(-dx_turb);
171  }
172  yloc.push_back(dy_turb);
173  }
174 
175  for(int it = 0;it<xloc.size(); it++){
176  xloc[it] = xloc[it] + windfarm_x_shift;
177  yloc[it] = yloc[it] + windfarm_y_shift;
178  }
179 }
ParmParse pp("prob")
@ tmp
Definition: ERF_AdvanceWSM6.cpp:116
Here is the call graph for this function:

◆ init_windfarm_x_y()

void WindFarm::init_windfarm_x_y ( const std::string  windfarm_loc_table)

Initialize wind farm using Cartesian coordinates.

Parameters
[in]windfarm_loc_tablePath to the location table.

Initialize turbine locations from Cartesian x and y coordinates.

Parameters
[in]windfarm_loc_tablePath to the turbine locations table.
188 {
189  // Read turbine locations from windturbines.txt
190  std::ifstream file(windfarm_loc_table);
191  if (!file.is_open()) {
192  amrex::Error("Wind turbines location table not found. Either the inputs is missing the"
193  " erf.windfarm_loc_table entry or the file specified in the entry " + windfarm_loc_table + " is missing.");
194  }
195  // Vector of vectors to store the matrix
196  Real value1, value2;
197 
198  // Re-reading the file must not append to a previous level's list.
199  xloc.clear();
200  yloc.clear();
201 
202  while (file >> value1 >> value2) {
203  value1 = value1 + 1e-3;
204  value2 = value2 + 1e-3;
205  xloc.push_back(value1);
206  yloc.push_back(value2);
207  }
208  file.close();
209 }

◆ read_tables()

void WindFarm::read_tables ( std::string  windfarm_loc_table,
std::string  windfarm_spec_table,
bool  x_y,
bool  lat_lon,
const amrex::Real  windfarm_x_shift = zero,
const amrex::Real  windfarm_y_shift = zero 
)

Read turbine location and specification tables.

Parameters
[in]windfarm_loc_tablePath to the location table.
[in]windfarm_spec_tablePath to the specification table.
[in]x_yWhether to use Cartesian coordinates.
[in]lat_lonWhether to use latitude and longitude coordinates.
[in]windfarm_x_shiftX-coordinate shift.
[in]windfarm_y_shiftY-coordinate shift.

Read in the turbine locations in latitude-longitude from windturbines.txt and convert it into x and y coordinates in metres

Parameters
levInteger specifying the current level
24 {
25  amrex::Print() << "Reading wind turbine locations table" << "\n";
26  read_windfarm_locations_table(windfarm_loc_table,
27  x_y, lat_lon,
28  windfarm_x_shift, windfarm_y_shift);
29 
30  amrex::Print() << "Reading wind turbine specifications table" << "\n";
31  read_windfarm_spec_table(windfarm_spec_table);
32 }
void read_windfarm_spec_table(const std::string windfarm_spec_table)
Read the turbine specification table.
Definition: ERF_WindFarm.cpp:218
void read_windfarm_locations_table(const std::string windfarm_loc_table, bool x_y, bool lat_lon, const amrex::Real windfarm_x_shift=zero, const amrex::Real windfarm_y_shift=zero)
Read the turbine location table.
Definition: ERF_WindFarm.cpp:44

◆ read_windfarm_airfoil_tables()

void WindFarm::read_windfarm_airfoil_tables ( const std::string  windfarm_airfoil_tables,
const std::string  windfarm_blade_table 
)

Read the airfoil tables.

Parameters
[in]windfarm_airfoil_tablesPath to the airfoil tables.
[in]windfarm_blade_tablePath to the blade table.

Read airfoil data from a directory of tables.

Parameters
[in]windfarm_airfoil_tablesPath to the directory containing airfoil tables.
[in]windfarm_blade_tablePath to the turbine blade properties table.
356 {
357  DIR* dir;
358  struct dirent* entry;
359  std::vector<std::string> files;
360 
361  // Check if directory exists
362  if ((dir = opendir(windfarm_airfoil_tables.c_str())) == nullptr) {
363  Abort("You are using a generalized actuator disk model based on blade element theory. This needs info of airfoil"
364  " cross sections over the span of the blade. There needs to be an entry erf.airfoil_tables which is the directory that"
365  " contains the angle of attack, Cl, Cd data for each airfoil cross-section. Either the entry is missing or the directory specified"
366  " in the entry - " + windfarm_airfoil_tables + " is missing. Exiting...");
367  }
368 
369  // Loop through directory entries and collect filenames
370  while ((entry = readdir(dir)) != nullptr) {
371  // Skip special directory entries "." and ".."
372  if (std::string(entry->d_name) == "." || std::string(entry->d_name) == "..") {
373  continue;
374  }
375  files.emplace_back(windfarm_airfoil_tables + "/" + entry->d_name); // Add file path to vector
376  }
377 
378  // Close the directory
379  closedir(dir);
380 
381  if (files.empty()) {
382  Abort("It seems the directory containing the info of airfoil cross sections of the blades - " + windfarm_airfoil_tables +
383  " is empty. Exiting...");
384  }
385 
386  if(files.size() != static_cast<long double>(n_bld_sections)) {
387  printf("There are %d airfoil sections in the last column of %s. But the number"
388  " of files in %s is only %ld.\n", n_bld_sections, windfarm_blade_table.c_str(),
389  windfarm_airfoil_tables.c_str(), files.size());
390  Abort("The number of blade sections from " + windfarm_blade_table + " should match the number of"
391  " files in " + windfarm_airfoil_tables + ". Exiting...");
392  }
393 
394  // Sort filenames in lexicographical (alphabetical) order
395  std::sort(files.begin(), files.end());
396 
397  // Process each file
398  int count = 0;
402  for (const auto& filePath : files) {
403  std::ifstream filename(filePath.c_str());
404 
405  // Every file in this directory is one blade section, and the sections are
406  // filled in sorted order, so one that cannot be read cannot be skipped over:
407  // doing so would shift every later section into the wrong slot.
408  if (!filename.is_open()) {
409  Abort("Could not open the airfoil section file " + filePath + " found in the directory " +
410  windfarm_airfoil_tables + ". Exiting...");
411  }
412 
413  std::cout << "Reading file: " << filePath << std::endl;
414 
415  std::string line;
416  for (int i = 0; i < 54; ++i) {
417  if (std::getline(filename, line)) { // Read one line into the array
418  }
419  }
420 
421  Real var1, var2, var3, temp;
422 
423  while(filename >> var1 >> var2 >> var3 >> temp) {
424  bld_airfoil_aoa[count].push_back(var1);
425  bld_airfoil_Cl[count].push_back(var2);
426  bld_airfoil_Cd[count].push_back(var3);
427  //int idx = bld_airfoil_aoa.size()-1;
428  //printf("Values are = %0.15g %0.15g %0.15g\n", bld_airfoil_aoa[idx], bld_airfoil_Cl[idx], bld_airfoil_Cd[idx]);
429  }
430  count++;
431  }
432 
434 }
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cl
Definition: ERF_WindFarm.H:285
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_aoa
Definition: ERF_WindFarm.H:285
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cd
Definition: ERF_WindFarm.H:285
int n_bld_sections
Definition: ERF_WindFarm.H:286
void set_blade_airfoil_spec(const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_aoa, const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_Cl, const amrex::Vector< amrex::Vector< amrex::Real >> &a_bld_airfoil_Cd) override
Set blade airfoil specifications.
Definition: ERF_WindFarm.H:254

◆ read_windfarm_blade_table()

void WindFarm::read_windfarm_blade_table ( const std::string  windfarm_blade_table)

Read the turbine blade table.

Parameters
[in]windfarm_blade_tablePath to the blade table.

Read turbine blade geometry and properties from a table.

Parameters
[in]windfarm_blade_tablePath to the turbine blade properties table.
272 {
273  std::ifstream filename(windfarm_blade_table);
274  std::string line;
275  Real temp, var1, var2, var3;
276  if (!filename.is_open()) {
277  Error("You are using a generalized actuator disk model based on blade element theory. This needs info of blades."
278  " An entry erf.windfarm_blade_table is needed. Either the entry is missing or the file specified"
279  " in the entry - " + windfarm_blade_table + " is missing.");
280  }
281  else {
282  Print() << "Reading in wind farm blade table: " << windfarm_blade_table << "\n";
283 
284  // First 6 lines are comments
285 
286  for (int i = 0; i < 6; ++i) {
287  if (std::getline(filename, line)) { // Read one line into the array
288  }
289  }
290 
291  while(filename >> var1 >> temp >> temp >> temp >> var2 >> var3 >> temp) {
292  bld_rad_loc.push_back(var1);
293  bld_twist.push_back(var2);
294  bld_chord.push_back(var3);
295  //int idx = bld_rad_loc.size()-1;
296  //printf("Values are = %0.15g %0.15g %0.15g\n", bld_rad_loc[idx], bld_twist[idx], bld_chord[idx]);
297  }
299  n_bld_sections = bld_rad_loc.size();
300  }
301 }
void set_blade_spec(const amrex::Vector< amrex::Real > &a_bld_rad_loc, const amrex::Vector< amrex::Real > &a_bld_twist, const amrex::Vector< amrex::Real > &a_bld_chord) override
Set blade specifications.
Definition: ERF_WindFarm.H:241
amrex::Vector< amrex::Real > bld_twist
Definition: ERF_WindFarm.H:284
amrex::Vector< amrex::Real > bld_chord
Definition: ERF_WindFarm.H:284
amrex::Vector< amrex::Real > bld_rad_loc
Definition: ERF_WindFarm.H:284

◆ read_windfarm_locations_table()

void WindFarm::read_windfarm_locations_table ( const std::string  windfarm_loc_table,
bool  x_y,
bool  lat_lon,
const amrex::Real  windfarm_x_shift = zero,
const amrex::Real  windfarm_y_shift = zero 
)

Read the turbine location table.

Parameters
[in]windfarm_loc_tablePath to the location table.
[in]x_yWhether to use Cartesian coordinates.
[in]lat_lonWhether to use latitude and longitude coordinates.
[in]windfarm_x_shiftX-coordinate shift.
[in]windfarm_y_shiftY-coordinate shift.

Read turbine locations from a table and convert them to Cartesian coordinates.

Parameters
[in]windfarm_loc_tablePath to the turbine locations table.
[in]x_yBoolean indicating if locations are provided as x and y coordinates.
[in]lat_lonBoolean indicating if locations are provided as latitude and longitude.
[in]windfarm_x_shiftShift applied to the resulting x coordinates.
[in]windfarm_y_shiftShift applied to the resulting y coordinates.
48 {
49  if(x_y) {
50  init_windfarm_x_y(windfarm_loc_table);
51  }
52  else if(lat_lon) {
53  init_windfarm_lat_lon(windfarm_loc_table, windfarm_x_shift, windfarm_y_shift);
54  }
55  else {
56  amrex::Abort("Are you using windfarms? For windfarm simulations, the inputs need to have an"
57  " entry erf.windfarm_loc_type which should be either lat_lon or x_y. \n");
58  }
59 
61 }
void set_turb_loc(const amrex::Vector< amrex::Real > &a_xloc, const amrex::Vector< amrex::Real > &a_yloc) override
Set turbine locations.
Definition: ERF_WindFarm.H:220
void init_windfarm_x_y(const std::string windfarm_loc_table)
Initialize wind farm using Cartesian coordinates.
Definition: ERF_WindFarm.cpp:187
void init_windfarm_lat_lon(const std::string windfarm_loc_table, const amrex::Real windfarm_x_shift, const amrex::Real windfarm_y_shift)
Initialize wind farm using latitude and longitude coordinates.
Definition: ERF_WindFarm.cpp:71

◆ read_windfarm_spec_table()

void WindFarm::read_windfarm_spec_table ( const std::string  windfarm_spec_table)

Read the turbine specification table.

Parameters
[in]windfarm_spec_tablePath to the specification table.

Read turbine performance specifications from a table.

Parameters
[in]windfarm_spec_tablePath to the turbine specifications table.
219 {
220  //The first line is the number of pairs entries for the power curve and thrust coefficient.
221  //The second line gives first the height in meters of the turbine hub, second, the diameter in
222  //meters of the rotor, third the standing thrust coefficient, and fourth the nominal power of
223  //the turbine in MW.
224  //The remaining lines contain the three values of: wind speed, thrust coefficient, and power production in kW.
225 
226  // Read turbine data from wind-turbine-1.tbl
227  std::ifstream file_turb_table(windfarm_spec_table);
228  if (!file_turb_table.is_open()) {
229  Error("Wind farm specifications table not found. Either the inputs is missing the "
230  "erf.windfarm_spec_table entry or the file specified in the entry - " + windfarm_spec_table + " is missing.");
231  }
232  else {
233  Print() << "Reading in wind farm specifications table: " << windfarm_spec_table << "\n";
234  }
235 
236  int nlines;
237  file_turb_table >> nlines;
238  wind_speed.resize(nlines);
239  thrust_coeff.resize(nlines);
240  power.resize(nlines);
241 
242  Real rotor_dia;
243  file_turb_table >> hub_height >> rotor_dia >> thrust_coeff_standing >> nominal_power;
244  rotor_rad = rotor_dia*myhalf;
245  if(rotor_rad > hub_height) {
246  Abort("The blade length is more than the hub height. Check the second line in wind-turbine-1.tbl. Aborting.....");
247  }
248  if(thrust_coeff_standing > one) {
249  Abort("The standing thrust coefficient is greater than one. Check the second line in wind-turbine-1.tbl. Aborting.....");
250  }
251 
252  for(int iline=0;iline<nlines;iline++){
253  file_turb_table >> wind_speed[iline] >> thrust_coeff[iline] >> power[iline];
254  if(thrust_coeff[iline] > one) {
255  Abort("The thrust coefficient is greater than one. Check wind-turbine-1.tbl. Aborting.....");
256  }
257  }
258  file_turb_table.close();
259 
262 
263 }
amrex::Vector< amrex::Real > wind_speed
Definition: ERF_WindFarm.H:283
amrex::Vector< amrex::Real > thrust_coeff
Definition: ERF_WindFarm.H:283
amrex::Real nominal_power
Definition: ERF_WindFarm.H:282
void set_turb_spec(const amrex::Real &a_rotor_rad, const amrex::Real &a_hub_height, const amrex::Real &a_thrust_coeff_standing, const amrex::Vector< amrex::Real > &a_wind_speed, const amrex::Vector< amrex::Real > &a_thrust_coeff, const amrex::Vector< amrex::Real > &a_power) override
Set general turbine specifications.
Definition: ERF_WindFarm.H:206
amrex::Real thrust_coeff_standing
Definition: ERF_WindFarm.H:282
amrex::Vector< amrex::Real > power
Definition: ERF_WindFarm.H:283

◆ read_windfarm_spec_table_extra()

void WindFarm::read_windfarm_spec_table_extra ( const std::string  windfarm_spec_table_extra)

Read extra turbine specifications.

Parameters
[in]windfarm_spec_table_extraPath to the extra specification table.

Read extended turbine performance specifications from a table.

Parameters
[in]windfarm_spec_table_extraPath to the extended turbine specifications table.
310 {
311  // Open the file
312  std::ifstream file(windfarm_spec_table_extra);
313 
314  // Check if file opened successfully
315  if (!file.is_open()) {
316  Abort("Error: You are using generalized wind farms option. This requires an input file erf.windfarm_spec_table_extra."
317  " Either this entry is missing in the inputs or the file specified -" + windfarm_spec_table_extra + " does"
318  " not exist. Exiting...");
319  } else {
320  printf("Reading in windfarm_spec_table_extra %s", windfarm_spec_table_extra.c_str());
321  }
322 
323  // Ignore the first line (header)
324  std::string header;
325  std::getline(file, header);
326 
327  // Variables to hold each row's values
328  double V, Cp, Ct, rpm, pitch, temp;
329 
330  // Read the file row by row
331  while (file >> V) {
332  char comma; // To ignore the commas
333  file >> comma >> Cp >> comma >> Ct >> comma >> temp >> comma >> temp >> comma
334  >> temp >> comma >> rpm >> comma >> pitch >> comma >> temp;
335 
336  velocity.push_back(V);
337  C_P.push_back(Cp);
338  C_T.push_back(Ct);
339  rotor_RPM.push_back(rpm);
340  blade_pitch.push_back(pitch);
341  }
342 
344 }
amrex::Vector< amrex::Real > blade_pitch
Definition: ERF_WindFarm.H:287
void set_turb_spec_extra(const amrex::Vector< amrex::Real > &a_velocity, const amrex::Vector< amrex::Real > &a_C_P, const amrex::Vector< amrex::Real > &a_C_T, const amrex::Vector< amrex::Real > &a_rotor_RPM, const amrex::Vector< amrex::Real > &a_blade_pitch) override
Set extra turbine specifications.
Definition: ERF_WindFarm.H:269
amrex::Vector< amrex::Real > velocity
Definition: ERF_WindFarm.H:287
amrex::Vector< amrex::Real > C_P
Definition: ERF_WindFarm.H:287
amrex::Vector< amrex::Real > C_T
Definition: ERF_WindFarm.H:287
amrex::Vector< amrex::Real > rotor_RPM
Definition: ERF_WindFarm.H:287
@ V
Definition: ERF_IndexDefines.H:130

◆ set_blade_airfoil_spec()

void WindFarm::set_blade_airfoil_spec ( const amrex::Vector< amrex::Vector< amrex::Real >> &  a_bld_airfoil_aoa,
const amrex::Vector< amrex::Vector< amrex::Real >> &  a_bld_airfoil_Cl,
const amrex::Vector< amrex::Vector< amrex::Real >> &  a_bld_airfoil_Cd 
)
inlineoverridevirtual

Set blade airfoil specifications.

Parameters
[in]a_bld_airfoil_aoaAngle of attack.
[in]a_bld_airfoil_ClLift coefficient.
[in]a_bld_airfoil_CdDrag coefficient.

Reimplemented from NullWindFarm.

257  {
258  m_windfarm_model[0]->set_blade_airfoil_spec(a_bld_airfoil_aoa, a_bld_airfoil_Cl, a_bld_airfoil_Cd);
259  }

◆ set_blade_spec()

void WindFarm::set_blade_spec ( const amrex::Vector< amrex::Real > &  a_bld_rad_loc,
const amrex::Vector< amrex::Real > &  a_bld_twist,
const amrex::Vector< amrex::Real > &  a_bld_chord 
)
inlineoverridevirtual

Set blade specifications.

Parameters
[in]a_bld_rad_locBlade radial locations.
[in]a_bld_twistBlade twist.
[in]a_bld_chordBlade chord.

Reimplemented from NullWindFarm.

244  {
245  m_windfarm_model[0]->set_blade_spec(a_bld_rad_loc, a_bld_twist, a_bld_chord);
246  }

◆ set_turb_disk_angle()

void WindFarm::set_turb_disk_angle ( const amrex::Real &  a_turb_disk_angle)
inlineoverridevirtual

Set the turbine disk angle.

Parameters
[in]a_turb_disk_angleDisk angle.

Reimplemented from NullWindFarm.

231  {
232  m_windfarm_model[0]->set_turb_disk_angle(a_turb_disk_angle);
233  }

◆ set_turb_loc()

void WindFarm::set_turb_loc ( const amrex::Vector< amrex::Real > &  a_xloc,
const amrex::Vector< amrex::Real > &  a_yloc 
)
inlineoverridevirtual

Set turbine locations.

Parameters
[in]a_xlocX locations.
[in]a_ylocY locations.

Reimplemented from NullWindFarm.

222  {
223  m_windfarm_model[0]->set_turb_loc(a_xloc, a_yloc);
224  }

◆ set_turb_spec()

void WindFarm::set_turb_spec ( const amrex::Real &  a_rotor_rad,
const amrex::Real &  a_hub_height,
const amrex::Real &  a_thrust_coeff_standing,
const amrex::Vector< amrex::Real > &  a_wind_speed,
const amrex::Vector< amrex::Real > &  a_thrust_coeff,
const amrex::Vector< amrex::Real > &  a_power 
)
inlineoverridevirtual

Set general turbine specifications.

Parameters
[in]a_rotor_radRotor radius.
[in]a_hub_heightHub height.
[in]a_thrust_coeff_standingStanding thrust coefficient.
[in]a_wind_speedWind speed vector.
[in]a_thrust_coeffThrust coefficient vector.
[in]a_powerPower vector.

Reimplemented from NullWindFarm.

210  {
211  m_windfarm_model[0]->set_turb_spec(a_rotor_rad, a_hub_height, a_thrust_coeff_standing,
212  a_wind_speed, a_thrust_coeff, a_power);
213  }

◆ set_turb_spec_extra()

void WindFarm::set_turb_spec_extra ( const amrex::Vector< amrex::Real > &  a_velocity,
const amrex::Vector< amrex::Real > &  a_C_P,
const amrex::Vector< amrex::Real > &  a_C_T,
const amrex::Vector< amrex::Real > &  a_rotor_RPM,
const amrex::Vector< amrex::Real > &  a_blade_pitch 
)
inlineoverridevirtual

Set extra turbine specifications.

Parameters
[in]a_velocityVelocity.
[in]a_C_PPower coefficient.
[in]a_C_TThrust coefficient.
[in]a_rotor_RPMRotor RPM.
[in]a_blade_pitchBlade pitch.

Reimplemented from NullWindFarm.

274  {
275  m_windfarm_model[0]->set_turb_spec_extra(a_velocity, a_C_P, a_C_T, a_rotor_RPM, a_blade_pitch);
276  }

◆ SetModel()

template<class NewWindFarmModel >
void WindFarm::SetModel ( )
inlineprotected

Create and set the specified windfarm model.

292  {
293  for (int lev(0); lev<m_windfarm_model.size(); ++lev) {
294  m_windfarm_model[lev] = std::make_unique<NewWindFarmModel>();
295  }
296  }

◆ write_actuator_disks_vtk()

void WindFarm::write_actuator_disks_vtk ( const amrex::Geometry &  geom,
const amrex::Real &  sampling_distance_by_D 
)

Write actuator disks to a VTK file.

Parameters
[in]geomGeometry
[in]sampling_distance_by_DSampling distance normalized by rotor diameter.

Write the actuator disks and sampling disks to VTK files.

Parameters
[in]geomGeometry defining the domain.
[in]sampling_distance_by_DNormalized sampling distance.
857 {
858 
859  Real sampling_distance = sampling_distance_by_D*two*rotor_rad;
860 
861  if (ParallelDescriptor::IOProcessor()){
862  FILE *file_actuator_disks_all, *file_actuator_disks_in_dom, *file_averaging_disks_in_dom;
863  file_actuator_disks_all = fopen("actuator_disks_all.vtk","w");
864  fprintf(file_actuator_disks_all, "%s\n","# vtk DataFile Version 3.0");
865  fprintf(file_actuator_disks_all, "%s\n","Actuator Disks");
866  fprintf(file_actuator_disks_all, "%s\n","ASCII");
867  fprintf(file_actuator_disks_all, "%s\n","DATASET POLYDATA");
868 
869  file_actuator_disks_in_dom = fopen("actuator_disks_in_dom.vtk","w");
870  fprintf(file_actuator_disks_in_dom, "%s\n","# vtk DataFile Version 3.0");
871  fprintf(file_actuator_disks_in_dom, "%s\n","Actuator Disks");
872  fprintf(file_actuator_disks_in_dom, "%s\n","ASCII");
873  fprintf(file_actuator_disks_in_dom, "%s\n","DATASET POLYDATA");
874 
875  file_averaging_disks_in_dom = fopen("averaging_disks_in_dom.vtk","w");
876  fprintf(file_averaging_disks_in_dom, "%s\n","# vtk DataFile Version 3.0");
877  fprintf(file_averaging_disks_in_dom, "%s\n","Actuator Disks");
878  fprintf(file_averaging_disks_in_dom, "%s\n","ASCII");
879  fprintf(file_averaging_disks_in_dom, "%s\n","DATASET POLYDATA");
880 
881 
882  int npts = 100;
883  fprintf(file_actuator_disks_all, "%s %ld %s\n", "POINTS", xloc.size()*npts, "float");
884  auto ProbLoArr = geom.ProbLoArray();
885  auto ProbHiArr = geom.ProbHiArray();
886  int num_turb_in_dom = 0;
887 
888  // Find the number of turbines inside the specified computational domain
889 
890  for(int it=0; it<xloc.size(); it++){
891  Real x = xloc[it];
892  Real y = yloc[it];
893  if(x > ProbLoArr[0] and x < ProbHiArr[0] and y > ProbLoArr[1] and y < ProbHiArr[1]) {
894  num_turb_in_dom++;
895  }
896  }
897  fprintf(file_actuator_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast<long int>(num_turb_in_dom*npts), "float");
898  fprintf(file_averaging_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast<long int>(num_turb_in_dom*npts), "float");
899 
900  Real nx = std::cos(my_turb_disk_angle+myhalf*M_PI);
901  Real ny = std::sin(my_turb_disk_angle+myhalf*M_PI);
902 
903  Real nx1 = -std::cos(my_turb_disk_angle);
904  Real ny1 = -std::sin(my_turb_disk_angle);
905 
906  for(int it=0; it<xloc.size(); it++){
907  for(int pt=0;pt<100;pt++){
908  Real x, y, z, xavg, yavg;
909  Real theta = two*M_PI/npts*pt;
910  x = xloc[it] + rotor_rad*cos(theta)*nx;
911  y = yloc[it] + rotor_rad*cos(theta)*ny;
912  z = hub_height + zloc[it] + rotor_rad*sin(theta);
913 
914  xavg = xloc[it] + sampling_distance*nx1 + rotor_rad*cos(theta)*nx;
915  yavg = yloc[it] + sampling_distance*ny1 + rotor_rad*cos(theta)*ny;
916 
917  fprintf(file_actuator_disks_all, "%0.15g %0.15g %0.15g\n", x, y, z);
918  if(xloc[it] > ProbLoArr[0] and xloc[it] < ProbHiArr[0] and yloc[it] > ProbLoArr[1] and yloc[it] < ProbHiArr[1]) {
919  fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", x, y, z);
920  fprintf(file_averaging_disks_in_dom, "%0.15g %0.15g %0.15g\n", xavg, yavg, z);
921  }
922  }
923  }
924  fprintf(file_actuator_disks_all, "%s %ld %ld\n", "LINES", xloc.size()*(npts-1), static_cast<long int>(xloc.size()*(npts-1)*3));
925  fprintf(file_actuator_disks_in_dom, "%s %ld %ld\n", "LINES", static_cast<long int>(num_turb_in_dom*(npts-1)), static_cast<long int>(num_turb_in_dom*(npts-1)*3));
926  fprintf(file_averaging_disks_in_dom, "%s %ld %ld\n", "LINES", static_cast<long int>(num_turb_in_dom*(npts-1)), static_cast<long int>(num_turb_in_dom*(npts-1)*3));
927  for(int it=0; it<xloc.size(); it++){
928  for(int pt=0;pt<99;pt++){
929  fprintf(file_actuator_disks_all, "%ld %ld %ld\n",
930  static_cast<long int>(2),
931  static_cast<long int>(it*npts+pt),
932  static_cast<long int>(it*npts+pt+1));
933  }
934  }
935  for(int it=0; it<num_turb_in_dom; it++){
936  for(int pt=0;pt<99;pt++){
937  fprintf(file_actuator_disks_in_dom, "%ld %ld %ld\n",
938  static_cast<long int>(2),
939  static_cast<long int>(it*npts+pt),
940  static_cast<long int>(it*npts+pt+1));
941  }
942  }
943 
944  for(int it=0; it<num_turb_in_dom; it++){
945  for(int pt=0;pt<99;pt++){
946  fprintf(file_averaging_disks_in_dom, "%ld %ld %ld\n",
947  static_cast<long int>(2),
948  static_cast<long int>(it*npts+pt),
949  static_cast<long int>(it*npts+pt+1));
950  }
951  }
952 
953  fclose(file_actuator_disks_all);
954  fclose(file_actuator_disks_in_dom);
955  fclose(file_averaging_disks_in_dom);
956  }
957 }

◆ write_turbine_locations_vtk()

void WindFarm::write_turbine_locations_vtk ( )

Write turbine locations to a VTK file.

Write the turbine hub locations to a VTK file.

831 {
832  if (ParallelDescriptor::IOProcessor()){
833  FILE* file_turbloc_vtk;
834  file_turbloc_vtk = fopen("turbine_locations.vtk","w");
835  fprintf(file_turbloc_vtk, "%s\n","# vtk DataFile Version 3.0");
836  fprintf(file_turbloc_vtk, "%s\n","Wind turbine locations");
837  fprintf(file_turbloc_vtk, "%s\n","ASCII");
838  fprintf(file_turbloc_vtk, "%s\n","DATASET POLYDATA");
839  fprintf(file_turbloc_vtk, "%s %ld %s\n", "POINTS", xloc.size(), "float");
840  for(int it=0; it<xloc.size(); it++){
841  fprintf(file_turbloc_vtk, "%0.15g %0.15g %0.15g\n", xloc[it], yloc[it], hub_height);
842  }
843  fclose(file_turbloc_vtk);
844  }
845 }

Member Data Documentation

◆ blade_pitch

amrex::Vector<amrex::Real> WindFarm::blade_pitch
protected

◆ bld_airfoil_aoa

amrex::Vector<amrex::Vector<amrex::Real> > WindFarm::bld_airfoil_aoa
protected

◆ bld_airfoil_Cd

amrex::Vector<amrex::Vector<amrex::Real> > WindFarm::bld_airfoil_Cd
protected

◆ bld_airfoil_Cl

amrex::Vector<amrex::Vector<amrex::Real> > WindFarm::bld_airfoil_Cl
protected

◆ bld_chord

amrex::Vector<amrex::Real> WindFarm::bld_chord
protected

◆ bld_rad_loc

amrex::Vector<amrex::Real> WindFarm::bld_rad_loc
protected

◆ bld_twist

amrex::Vector<amrex::Real> WindFarm::bld_twist
protected

◆ C_P

amrex::Vector<amrex::Real> WindFarm::C_P
protected

◆ C_T

amrex::Vector<amrex::Real> WindFarm::C_T
protected

◆ hub_height

amrex::Real WindFarm::hub_height
protected

◆ m_windfarm_model

amrex::Vector<std::unique_ptr<NullWindFarm> > WindFarm::m_windfarm_model
private

◆ my_turb_disk_angle

amrex::Real WindFarm::my_turb_disk_angle
protected

◆ n_bld_sections

int WindFarm::n_bld_sections
protected

◆ nominal_power

amrex::Real WindFarm::nominal_power
protected

◆ power

amrex::Vector<amrex::Real> WindFarm::power
protected

◆ rotor_rad

amrex::Real WindFarm::rotor_rad
protected

◆ rotor_RPM

amrex::Vector<amrex::Real> WindFarm::rotor_RPM
protected

◆ thrust_coeff

amrex::Vector<amrex::Real> WindFarm::thrust_coeff
protected

◆ thrust_coeff_standing

amrex::Real WindFarm::thrust_coeff_standing
protected

◆ velocity

amrex::Vector<amrex::Real> WindFarm::velocity
protected

◆ wind_speed

amrex::Vector<amrex::Real> WindFarm::wind_speed
protected

◆ xloc

amrex::Vector<amrex::Real> WindFarm::xloc
protected

◆ yloc

amrex::Vector<amrex::Real> WindFarm::yloc
protected

◆ zloc

amrex::Vector<amrex::Real> WindFarm::zloc
protected

The documentation for this class was generated from the following files: