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

#include <ERF_GeneralAD.H>

Inheritance diagram for GeneralAD:
Collaboration diagram for GeneralAD:

Public Member Functions

 GeneralAD ()
 
virtual ~GeneralAD ()=default
 
void advance (const amrex::Geometry &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
 
void compute_freestream_velocity (const amrex::MultiFab &cons_in, const amrex::MultiFab &U_old, const amrex::MultiFab &V_old, const amrex::MultiFab &mf_SMark)
 
void source_terms_cellcentered (const amrex::Geometry &geom, const amrex::MultiFab &cons_in, const amrex::MultiFab &mf_Smark, amrex::MultiFab &mf_vars_generalAD)
 
void update (const double &dt_advance, amrex::MultiFab &cons_in, amrex::MultiFab &U_old, amrex::MultiFab &V_old, amrex::MultiFab &W_old, const amrex::MultiFab &mf_vars)
 
void compute_power_output (const double &time)
 
- Public Member Functions inherited from NullWindFarm
 NullWindFarm ()
 
virtual ~NullWindFarm ()=default
 
virtual void set_turb_spec (const amrex::Real &rotor_rad, const amrex::Real &hub_height, const amrex::Real &thrust_coeff_standing, const amrex::Vector< amrex::Real > &wind_speed, const amrex::Vector< amrex::Real > &thrust_coeff, const amrex::Vector< amrex::Real > &power)
 
virtual void set_turb_loc (const amrex::Vector< amrex::Real > &xloc, const amrex::Vector< amrex::Real > &yloc)
 
virtual void set_turb_disk_angle (const amrex::Real &turb_disk_angle)
 
virtual void set_blade_spec (const amrex::Vector< amrex::Real > &bld_rad_loc, const amrex::Vector< amrex::Real > &bld_twist, const amrex::Vector< amrex::Real > &bld_chord)
 
virtual void set_blade_airfoil_spec (const amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_aoa, const amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_Cl, const amrex::Vector< amrex::Vector< amrex::Real >> &bld_airfoil_Cd)
 
virtual void set_turb_spec_extra (const amrex::Vector< amrex::Real > &velocity, const amrex::Vector< amrex::Real > &C_P, const amrex::Vector< amrex::Real > &C_T, const amrex::Vector< amrex::Real > &rotor_RPM, const amrex::Vector< amrex::Real > &blade_pitch)
 
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 Attributes

amrex::Vector< amrex::Realxloc
 
amrex::Vector< amrex::Realyloc
 
amrex::Real turb_disk_angle
 
amrex::Real hub_height
 
amrex::Real rotor_rad
 
amrex::Real thrust_coeff_standing
 
amrex::Real nominal_power
 
amrex::Vector< amrex::Realwind_speed
 
amrex::Vector< amrex::Realthrust_coeff
 
amrex::Vector< amrex::Realpower
 
amrex::Vector< amrex::Realfreestream_velocity
 
amrex::Vector< amrex::Realfreestream_phi
 
amrex::Vector< amrex::Realdisk_cell_count
 
amrex::Vector< amrex::Realbld_rad_loc
 
amrex::Vector< amrex::Realbld_twist
 
amrex::Vector< amrex::Realbld_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
 
amrex::Vector< amrex::Realvelocity
 
amrex::Vector< amrex::RealC_P
 
amrex::Vector< amrex::RealC_T
 
amrex::Vector< amrex::Realrotor_RPM
 
amrex::Vector< amrex::Realblade_pitch
 
- Protected Attributes inherited from NullWindFarm
amrex::Vector< amrex::Realm_xloc
 
amrex::Vector< amrex::Realm_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::Realm_wind_speed
 
amrex::Vector< amrex::Realm_thrust_coeff
 
amrex::Vector< amrex::Realm_power
 
amrex::Vector< amrex::Realm_bld_rad_loc
 
amrex::Vector< amrex::Realm_bld_twist
 
amrex::Vector< amrex::Realm_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::Realm_velocity
 
amrex::Vector< amrex::Realm_C_P
 
amrex::Vector< amrex::Realm_C_T
 
amrex::Vector< amrex::Realm_rotor_RPM
 
amrex::Vector< amrex::Realm_blade_pitch
 

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)
 

Constructor & Destructor Documentation

◆ GeneralAD()

GeneralAD::GeneralAD ( )
inline
12 {}

◆ ~GeneralAD()

virtual GeneralAD::~GeneralAD ( )
virtualdefault

Member Function Documentation

◆ advance()

void GeneralAD::advance ( const amrex::Geometry &  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 
)
overridevirtual

Implements NullWindFarm.

19 {
20  AMREX_ALWAYS_ASSERT(W_old.nComp() > 0);
21  AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0);
22  AMREX_ALWAYS_ASSERT(mf_vars_generalAD.nComp() > 0);
23  AMREX_ALWAYS_ASSERT(time > -one);
24  compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark);
25  source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD);
26  update(dt_advance, cons_in, U_old, V_old, W_old, mf_vars_generalAD);
28 }
constexpr amrex::Real one
Definition: ERF_Constants.H:9
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
void update(const double &dt_advance, amrex::MultiFab &cons_in, amrex::MultiFab &U_old, amrex::MultiFab &V_old, amrex::MultiFab &W_old, const amrex::MultiFab &mf_vars)
Definition: ERF_AdvanceGeneralAD.cpp:60
void source_terms_cellcentered(const amrex::Geometry &geom, const amrex::MultiFab &cons_in, const amrex::MultiFab &mf_Smark, amrex::MultiFab &mf_vars_generalAD)
Definition: ERF_AdvanceGeneralAD.cpp:322
void compute_freestream_velocity(const amrex::MultiFab &cons_in, const amrex::MultiFab &U_old, const amrex::MultiFab &V_old, const amrex::MultiFab &mf_SMark)
Definition: ERF_AdvanceGeneralAD.cpp:96
void compute_power_output(const double &time)
Definition: ERF_AdvanceGeneralAD.cpp:31
Here is the call graph for this function:

◆ compute_freestream_velocity()

void GeneralAD::compute_freestream_velocity ( const amrex::MultiFab &  cons_in,
const amrex::MultiFab &  U_old,
const amrex::MultiFab &  V_old,
const amrex::MultiFab &  mf_SMark 
)
100 {
102  freestream_velocity.clear();
103  freestream_phi.clear();
104  disk_cell_count.clear();
105  freestream_velocity.resize(xloc.size(),zero);
106  freestream_phi.resize(xloc.size(),zero);
107  disk_cell_count.resize(xloc.size(),zero);
108 
109  Gpu::DeviceVector<Real> d_freestream_velocity(xloc.size());
110  Gpu::DeviceVector<Real> d_freestream_phi(yloc.size());
111  Gpu::DeviceVector<Real> d_disk_cell_count(yloc.size());
112  Gpu::copy(Gpu::hostToDevice, freestream_velocity.begin(), freestream_velocity.end(), d_freestream_velocity.begin());
113  Gpu::copy(Gpu::hostToDevice, freestream_phi.begin(), freestream_phi.end(), d_freestream_phi.begin());
114  Gpu::copy(Gpu::hostToDevice, disk_cell_count.begin(), disk_cell_count.end(), d_disk_cell_count.begin());
115 
116  Real* d_freestream_velocity_ptr = d_freestream_velocity.data();
117  Real* d_freestream_phi_ptr = d_freestream_phi.data();
118  Real* d_disk_cell_count_ptr = d_disk_cell_count.data();
119 
120 
121  for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
122 
123  auto SMark_array = mf_SMark.array(mfi);
124  auto u_vel = U_old.array(mfi);
125  auto v_vel = V_old.array(mfi);
126  // NOTE: this reduction is driven by the cell-centered SMark, so it must
127  // run over the cell-centered tilebox. Using nodaltilebox(0) here
128  // would include the plane at bigEnd(0)+1, which is a ghost cell of
129  // this box and a valid cell of its x-neighbor, and so would count
130  // that plane twice in the sums below.
131  Box tbx = mfi.tilebox();
132 
133  ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
134 
135  if(SMark_array(i,j,k,0) != -one) {
136  int turb_index = static_cast<int>(SMark_array(i,j,k,0));
137  Real phi = std::atan2(v_vel(i,j,k),u_vel(i,j,k)); // Wind direction w.r.t the x-direction
138  Gpu::Atomic::Add(&d_freestream_velocity_ptr[turb_index],std::pow(u_vel(i,j,k)*u_vel(i,j,k) + v_vel(i,j,k)*v_vel(i,j,k),myhalf));
139  Gpu::Atomic::Add(&d_disk_cell_count_ptr[turb_index],one);
140  Gpu::Atomic::Add(&d_freestream_phi_ptr[turb_index],phi);
141  }
142  });
143  }
144 
145  // Copy back to host
146  Gpu::copy(Gpu::deviceToHost, d_freestream_velocity.begin(), d_freestream_velocity.end(), freestream_velocity.begin());
147  Gpu::copy(Gpu::deviceToHost, d_freestream_phi.begin(), d_freestream_phi.end(), freestream_phi.begin());
148  Gpu::copy(Gpu::deviceToHost, d_disk_cell_count.begin(), d_disk_cell_count.end(), disk_cell_count.begin());
149 
150  // Reduce the data on every processor
151  amrex::ParallelAllReduce::Sum(freestream_velocity.data(),
152  static_cast<int>(freestream_velocity.size()),
153  amrex::ParallelContext::CommunicatorAll());
154 
155  amrex::ParallelAllReduce::Sum(freestream_phi.data(),
156  static_cast<int>(freestream_phi.size()),
157  amrex::ParallelContext::CommunicatorAll());
158 
159 
160  amrex::ParallelAllReduce::Sum(disk_cell_count.data(),
161  static_cast<int>(disk_cell_count.size()),
162  amrex::ParallelContext::CommunicatorAll());
163 
165 
166 
167  /*if (ParallelDescriptor::IOProcessor()){
168  for(int it=0; it<xloc.size(); it++){
169  std::cout << "turbine index, freestream velocity is " << it << " " << freestream_velocity[it] << " " <<
170  disk_cell_count[it] << " " <<
171  freestream_velocity[it]/(disk_cell_count[it] + 1e-10) << " " <<
172  freestream_phi[it]/(disk_cell_count[it] + 1e-10) << "\n";
173  }
174  }*/
175 }
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real myhalf
Definition: ERF_Constants.H:13
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);})
amrex::Real Real
Definition: ERF_ShocInterface.H:19
amrex::Vector< amrex::Real > freestream_velocity
Definition: ERF_GeneralAD.H:51
amrex::Vector< amrex::Real > xloc
Definition: ERF_GeneralAD.H:47
amrex::Vector< amrex::Real > disk_cell_count
Definition: ERF_GeneralAD.H:51
amrex::Vector< amrex::Real > freestream_phi
Definition: ERF_GeneralAD.H:51
amrex::Vector< amrex::Real > yloc
Definition: ERF_GeneralAD.H:47
void get_turb_loc(amrex::Vector< amrex::Real > &xloc, amrex::Vector< amrex::Real > &yloc)
Definition: ERF_NullWindFarm.H:96
Here is the call graph for this function:

◆ compute_power_output()

void GeneralAD::compute_power_output ( const double &  time)
32 {
36 
37  const int n_spec_table = static_cast<int>(wind_speed.size());
38  // Compute power based on the look-up table
39 
40  if (ParallelDescriptor::IOProcessor()){
41  static std::ofstream file("power_output_GeneralAD.txt", std::ios::app);
42  // Check if the file opened successfully
43  if (!file.is_open()) {
44  std::cerr << "Error opening file!" << std::endl;
45  Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp");
46  }
47  Real total_power = zero;
48  for(int it=0; it<static_cast<int>(xloc.size()); it++){
49  Real avg_vel = freestream_velocity[it]/(disk_cell_count[it] + Real(1e-10));
50  Real turb_power = interpolate_1d(wind_speed.data(), power.data(), avg_vel, n_spec_table);
51  total_power = total_power + turb_power;
52  }
53  file << time << " " << total_power << "\n";
54  file.flush();
55  }
56 }
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real interpolate_1d(const amrex::Real *alpha, const amrex::Real *beta, const amrex::Real alpha_interp, const int alpha_size)
Definition: ERF_Interpolation_1D.H:14
amrex::Vector< amrex::Real > wind_speed
Definition: ERF_GeneralAD.H:50
amrex::Vector< amrex::Real > power
Definition: ERF_GeneralAD.H:50
amrex::Vector< amrex::Real > thrust_coeff
Definition: ERF_GeneralAD.H:50
amrex::Real rotor_rad
Definition: ERF_GeneralAD.H:49
amrex::Real hub_height
Definition: ERF_GeneralAD.H:49
amrex::Real thrust_coeff_standing
Definition: ERF_GeneralAD.H:49
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)
Definition: ERF_NullWindFarm.H:84
Here is the call graph for this function:

◆ source_terms_cellcentered()

void GeneralAD::source_terms_cellcentered ( const amrex::Geometry &  geom,
const amrex::MultiFab &  cons_in,
const amrex::MultiFab &  mf_Smark,
amrex::MultiFab &  mf_vars_generalAD 
)
326 {
327 
329 
332 
334 
336 
338 
339  Real d_hub_height = hub_height;
340  Real d_rotor_rad = rotor_rad;
341 
342  Gpu::DeviceVector<Real> d_xloc(xloc.size());
343  Gpu::DeviceVector<Real> d_yloc(yloc.size());
344  Gpu::copy(Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin());
345  Gpu::copy(Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin());
346 
347  auto dx = geom.CellSizeArray();
348 
349  // Domain valid box
350  const amrex::Box& domain = geom.Domain();
351  auto ProbLoArr = geom.ProbLoArray();
352  int domlo_x = domain.smallEnd(0);
353  int domhi_x = domain.bigEnd(0) + 1;
354  int domlo_y = domain.smallEnd(1);
355  int domhi_y = domain.bigEnd(1) + 1;
356  int domlo_z = domain.smallEnd(2);
357  int domhi_z = domain.bigEnd(2) + 1;
358 
359  // The order of variables are - Vabs dVabsdt, dudt, dvdt, dTKEdt
360  mf_vars_generalAD.setVal(0.0);
361 
362  long unsigned int nturbs = static_cast<long unsigned int>(xloc.size());
363 
364  // This is the angle phi in Fig. 10 in Mirocha et. al. 2014
365  // set_turb_disk angle in ERF_InitWindFarm.cpp sets this phi as
366  // the turb_disk_angle
368  Real d_turb_disk_angle = turb_disk_angle;
369 
370  Gpu::DeviceVector<Real> d_freestream_velocity(nturbs);
371  Gpu::DeviceVector<Real> d_disk_cell_count(nturbs);
372  Gpu::copy(Gpu::hostToDevice, freestream_velocity.begin(), freestream_velocity.end(), d_freestream_velocity.begin());
373  Gpu::copy(Gpu::hostToDevice, disk_cell_count.begin(), disk_cell_count.end(), d_disk_cell_count.begin());
374 
375  Real* d_xloc_ptr = d_xloc.data();
376  Real* d_yloc_ptr = d_yloc.data();
377  Real* d_freestream_velocity_ptr = d_freestream_velocity.data();
378  Real* d_disk_cell_count_ptr = d_disk_cell_count.data();
379 
380  int n_bld_sections = static_cast<int>(bld_rad_loc.size());
381 
382  Gpu::DeviceVector<Real> d_bld_rad_loc(n_bld_sections);
383  Gpu::DeviceVector<Real> d_bld_twist(n_bld_sections);
384  Gpu::DeviceVector<Real> d_bld_chord(n_bld_sections);
385 
386  Gpu::copy(Gpu::hostToDevice, bld_rad_loc.begin(), bld_rad_loc.end(), d_bld_rad_loc.begin());
387  Gpu::copy(Gpu::hostToDevice, bld_twist.begin(), bld_twist.end(), d_bld_twist.begin());
388  Gpu::copy(Gpu::hostToDevice, bld_chord.begin(), bld_chord.end(), d_bld_chord.begin());
389 
390  Real* bld_rad_loc_ptr = d_bld_rad_loc.data();
391  Real* bld_twist_ptr = d_bld_twist.data();
392  Real* bld_chord_ptr = d_bld_chord.data();
393 
394  Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_aoa(n_bld_sections);
395  Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_Cl(n_bld_sections);
396  Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_Cd(n_bld_sections);
397 
398  Vector<int> h_n_pts_airfoil(n_bld_sections);
399  Gpu::DeviceVector<int> d_n_pts_airfoil(n_bld_sections);
400 
401  for (int i = 0; i < n_bld_sections; ++i) {
402  h_n_pts_airfoil[i] = static_cast<int>(bld_airfoil_aoa[i].size());
403 
404  d_bld_airfoil_aoa[i].resize(bld_airfoil_aoa[i].size());
405  d_bld_airfoil_Cl[i].resize(bld_airfoil_Cl[i].size());
406  d_bld_airfoil_Cd[i].resize(bld_airfoil_Cd[i].size());
407 
408  Gpu::copy(Gpu::hostToDevice,
409  bld_airfoil_aoa[i].begin(), bld_airfoil_aoa[i].end(),
410  d_bld_airfoil_aoa[i].begin());
411 
412  Gpu::copy(Gpu::hostToDevice,
413  bld_airfoil_Cl[i].begin(), bld_airfoil_Cl[i].end(),
414  d_bld_airfoil_Cl[i].begin());
415 
416  Gpu::copy(Gpu::hostToDevice,
417  bld_airfoil_Cd[i].begin(), bld_airfoil_Cd[i].end(),
418  d_bld_airfoil_Cd[i].begin());
419  }
420 
421  Gpu::copy(Gpu::hostToDevice,
422  h_n_pts_airfoil.begin(), h_n_pts_airfoil.end(),
423  d_n_pts_airfoil.begin());
424 
425  int* d_n_pts_airfoil_ptr = d_n_pts_airfoil.data();
426 
427  Vector<Real*> hp_bld_airfoil_aoa, hp_bld_airfoil_Cl, hp_bld_airfoil_Cd;
428  for (auto & v :d_bld_airfoil_aoa) {
429  hp_bld_airfoil_aoa.push_back(v.data());
430  }
431  for (auto & v :d_bld_airfoil_Cl) {
432  hp_bld_airfoil_Cl.push_back(v.data());
433  }
434  for (auto & v :d_bld_airfoil_Cd) {
435  hp_bld_airfoil_Cd.push_back(v.data());
436  }
437 
438  Gpu::AsyncArray<Real*> aoa(hp_bld_airfoil_aoa.data(), n_bld_sections);
439  Gpu::AsyncArray<Real*> Cl(hp_bld_airfoil_Cl.data(), n_bld_sections);
440  Gpu::AsyncArray<Real*> Cd(hp_bld_airfoil_Cd.data(), n_bld_sections);
441 
442  auto d_bld_airfoil_aoa_ptr = aoa.data();
443  auto d_bld_airfoil_Cl_ptr = Cl.data();
444  auto d_bld_airfoil_Cd_ptr = Cd.data();
445 
446  int n_spec_extra = static_cast<int>(velocity.size());
447 
448  Gpu::DeviceVector<Real> d_velocity(n_spec_extra);
449  Gpu::DeviceVector<Real> d_rotor_RPM(n_spec_extra);
450  Gpu::DeviceVector<Real> d_blade_pitch(n_spec_extra);
451 
452  Gpu::copy(Gpu::hostToDevice, velocity.begin(), velocity.end(), d_velocity.begin());
453  Gpu::copy(Gpu::hostToDevice, rotor_RPM.begin(), rotor_RPM.end(), d_rotor_RPM.begin());
454  Gpu::copy(Gpu::hostToDevice, blade_pitch.begin(), blade_pitch.end(), d_blade_pitch.begin());
455 
456  auto d_velocity_ptr = d_velocity.data();
457  auto d_rotor_RPM_ptr = d_rotor_RPM.data();
458  auto d_blade_pitch_ptr = d_blade_pitch.data();
459 
460  for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
461 
462  const Box& gbx = mfi.growntilebox(1);
463  auto SMark_array = mf_SMark.array(mfi);
464  auto generalAD_array = mf_vars_generalAD.array(mfi);
465 
466  ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
467  int ii = amrex::min(amrex::max(i, domlo_x), domhi_x);
468  int jj = amrex::min(amrex::max(j, domlo_y), domhi_y);
469  int kk = amrex::min(amrex::max(k, domlo_z), domhi_z);
470 
471  Real x = ProbLoArr[0] + (ii+myhalf)*dx[0];
472  Real y = ProbLoArr[1] + (jj+myhalf)*dx[1];
473  Real z = ProbLoArr[2] + (kk+myhalf)*dx[2];
474  // ?? Density needed here
475 
476  int check_int = 0;
477 
478  Real source_x = zero, source_y = zero, source_z = zero;
479  std::array<Real,2> Fn_and_Ft;
480 
481  for(long unsigned int it=0;it<nturbs;it++) {
482  Real avg_vel = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] + Real(1e-10));
483  Real phi = d_turb_disk_angle;
484 
485  // This if check makes sure it is a point on the actuator disk
486  if(SMark_array(ii,jj,kk,1) == static_cast<double>(it)) {
487  check_int++;
488 
489  // Find radial distance of the point and the zeta angle
490  Real rad = std::pow( (x-d_xloc_ptr[it])*(x-d_xloc_ptr[it]) +
491  (y-d_yloc_ptr[it])*(y-d_yloc_ptr[it]) +
492  (z-d_hub_height)*(z-d_hub_height), myhalf );
493 
494  int index = find_rad_loc_index(rad, bld_rad_loc_ptr, n_bld_sections);
495 
496  // This if check makes sure it is a point with radial distance
497  // between the hub radius and the rotor radius.
498  // ?? hub radius needed here
499  if(rad >= two and rad <= d_rotor_rad) {
500  //AMREX_ASSERT( (z-d_hub_height) <= rad );
501  // Consider the vector that joines the point and the turbine center.
502  // Dot it on to the vector that joins the turbine center and along
503  // the plane of the disk. See fig. 10 in Mirocha et. al. Real(2014.)
504 
505  Real vec_proj = (x-d_xloc_ptr[it])*(std::sin(phi)) +
506  (y-d_yloc_ptr[it])*(-std::cos(phi));
507 
508 
509  Real zeta = std::atan2(z-d_hub_height, vec_proj);
510  //printf("zeta val is %0.15g\n", zeta*Real(180.0)/PI);
511  Fn_and_Ft = compute_source_terms_Fn_Ft(rad, avg_vel,
512  bld_rad_loc_ptr,
513  bld_twist_ptr,
514  bld_chord_ptr,
515  n_bld_sections,
516  d_bld_airfoil_aoa_ptr[index],
517  d_bld_airfoil_Cl_ptr[index],
518  d_bld_airfoil_Cd_ptr[index],
519  d_n_pts_airfoil_ptr[index],
520  d_velocity_ptr,
521  d_rotor_RPM_ptr,
522  d_blade_pitch_ptr,
523  n_spec_extra);
524 
525  Real Fn = three*Fn_and_Ft[0];
526  Real Ft = three*Fn_and_Ft[1];
527  // Compute the source terms - pass in radial distance, free stream velocity
528 
529  Real Fx = Fn*std::cos(phi) + Ft*std::sin(zeta)*std::sin(phi);
530  Real Fy = Fn*std::sin(phi) - Ft*std::sin(zeta)*std::cos(phi);
531  Real Fz = -Ft*std::cos(zeta);
532 
533  //Real dn = (
534 
535  source_x = -Fx/(two*PI*rad*dx[0])*one/std::pow(two*PI,myhalf);
536  source_y = -Fy/(two*PI*rad*dx[0])*one/std::pow(two*PI,myhalf);
537  source_z = -Fz/(two*PI*rad*dx[0])*one/std::pow(two*PI,myhalf);
538 
539 
540  //printf("Val source_x, is %0.15g, %0.15g, %0.15g %0.15g %0.15g %0.15g\n", rad, Fn, Ft, source_x, source_y, source_z);
541  }
542  }
543  }
544 
545  if(check_int > 1){
546  amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk "
547  "and check the windturbine locations input file. Exiting..");
548  }
549 
550  generalAD_array(i,j,k,0) = source_x;
551  generalAD_array(i,j,k,1) = source_y;
552  generalAD_array(i,j,k,2) = source_z;
553  });
554  }
555 }
AMREX_FORCE_INLINE AMREX_GPU_DEVICE int find_rad_loc_index(const Real rad, const Real *bld_rad_loc, const int n_bld_sections)
Definition: ERF_AdvanceGeneralAD.cpp:179
AMREX_FORCE_INLINE AMREX_GPU_DEVICE std::array< Real, 2 > compute_source_terms_Fn_Ft(const Real rad, const Real avg_vel, const Real *bld_rad_loc, const Real *bld_twist, const Real *bld_chord, int n_bld_sections, const Real *bld_airfoil_aoa, const Real *bld_airfoil_Cl, const Real *bld_airfoil_Cd, const int n_pts_airfoil, const Real *velocity, const Real *rotor_RPM, const Real *blade_pitch, const int n_spec_extra)
Definition: ERF_AdvanceGeneralAD.cpp:212
constexpr amrex::Real three
Definition: ERF_Constants.H:11
constexpr amrex::Real two
Definition: ERF_Constants.H:10
constexpr amrex::Real PI
Definition: ERF_Constants.H:42
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cl
Definition: ERF_GeneralAD.H:53
amrex::Vector< amrex::Real > velocity
Definition: ERF_GeneralAD.H:54
amrex::Vector< amrex::Real > blade_pitch
Definition: ERF_GeneralAD.H:54
amrex::Vector< amrex::Real > C_P
Definition: ERF_GeneralAD.H:54
amrex::Vector< amrex::Real > bld_twist
Definition: ERF_GeneralAD.H:52
amrex::Vector< amrex::Real > bld_chord
Definition: ERF_GeneralAD.H:52
amrex::Vector< amrex::Real > bld_rad_loc
Definition: ERF_GeneralAD.H:52
amrex::Real turb_disk_angle
Definition: ERF_GeneralAD.H:48
amrex::Vector< amrex::Real > C_T
Definition: ERF_GeneralAD.H:54
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_aoa
Definition: ERF_GeneralAD.H:53
amrex::Vector< amrex::Real > rotor_RPM
Definition: ERF_GeneralAD.H:54
amrex::Vector< amrex::Vector< amrex::Real > > bld_airfoil_Cd
Definition: ERF_GeneralAD.H:53
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)
Definition: ERF_NullWindFarm.H:126
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)
Definition: ERF_NullWindFarm.H:117
void get_turb_disk_angle(amrex::Real &turb_disk_angle)
Definition: ERF_NullWindFarm.H:103
void get_blade_spec(amrex::Vector< amrex::Real > &bld_rad_loc, amrex::Vector< amrex::Real > &bld_twist, amrex::Vector< amrex::Real > &bld_chord)
Definition: ERF_NullWindFarm.H:108
Here is the call graph for this function:

◆ update()

void GeneralAD::update ( const double &  dt_advance,
amrex::MultiFab &  cons_in,
amrex::MultiFab &  U_old,
amrex::MultiFab &  V_old,
amrex::MultiFab &  W_old,
const amrex::MultiFab &  mf_vars 
)
66 {
67 
68  for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
69 
70  Box tbx = mfi.nodaltilebox(0);
71  Box tby = mfi.nodaltilebox(1);
72  Box tbz = mfi.nodaltilebox(2);
73 
74  auto generalAD_array = mf_vars_generalAD.array(mfi);
75  auto u_vel = U_old.array(mfi);
76  auto v_vel = V_old.array(mfi);
77  auto w_vel = W_old.array(mfi);
78 
79  ParallelFor(tbx, tby, tbz,
80  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
81  {
82  u_vel(i,j,k) = u_vel(i,j,k) + (generalAD_array(i-1,j,k,0) + generalAD_array(i,j,k,0))/two*dt_advance;
83  },
84  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
85  {
86  v_vel(i,j,k) = v_vel(i,j,k) + (generalAD_array(i,j-1,k,1) + generalAD_array(i,j,k,1))/two*dt_advance;
87  },
88  [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
89  {
90  w_vel(i,j,k) = w_vel(i,j,k) + (generalAD_array(i,j,k-1,2) + generalAD_array(i,j,k,2))/two*dt_advance;
91  });
92 
93  }
94 }
Here is the call graph for this function:

Member Data Documentation

◆ blade_pitch

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

◆ bld_airfoil_aoa

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

◆ bld_airfoil_Cd

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

◆ bld_airfoil_Cl

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

◆ bld_chord

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

◆ bld_rad_loc

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

◆ bld_twist

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

◆ C_P

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

◆ C_T

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

◆ disk_cell_count

amrex::Vector<amrex::Real> GeneralAD::disk_cell_count
protected

◆ freestream_phi

amrex::Vector<amrex::Real> GeneralAD::freestream_phi
protected

◆ freestream_velocity

amrex::Vector<amrex::Real> GeneralAD::freestream_velocity
protected

◆ hub_height

amrex::Real GeneralAD::hub_height
protected

◆ nominal_power

amrex::Real GeneralAD::nominal_power
protected

◆ power

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

◆ rotor_rad

amrex::Real GeneralAD::rotor_rad
protected

◆ rotor_RPM

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

◆ thrust_coeff

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

◆ thrust_coeff_standing

amrex::Real GeneralAD::thrust_coeff_standing
protected

◆ turb_disk_angle

amrex::Real GeneralAD::turb_disk_angle
protected

◆ velocity

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

◆ wind_speed

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

◆ xloc

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

◆ yloc

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

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