343 Gpu::DeviceVector<Real> d_xloc(
xloc.size());
344 Gpu::DeviceVector<Real> d_yloc(
yloc.size());
345 Gpu::copy(Gpu::hostToDevice,
xloc.begin(),
xloc.end(), d_xloc.begin());
346 Gpu::copy(Gpu::hostToDevice,
yloc.begin(),
yloc.end(), d_yloc.begin());
348 auto dx = geom.CellSizeArray();
351 const amrex::Box& domain = geom.Domain();
352 auto ProbLoArr = geom.ProbLoArray();
353 int domlo_x = domain.smallEnd(0);
354 int domhi_x = domain.bigEnd(0) + 1;
355 int domlo_y = domain.smallEnd(1);
356 int domhi_y = domain.bigEnd(1) + 1;
357 int domlo_z = domain.smallEnd(2);
358 int domhi_z = domain.bigEnd(2) + 1;
361 mf_vars_generalAD.setVal(0.0);
363 long unsigned int nturbs =
static_cast<long unsigned int>(
xloc.size());
371 Gpu::DeviceVector<Real> d_freestream_velocity(nturbs);
372 Gpu::DeviceVector<Real> d_disk_cell_count(nturbs);
376 Real* d_xloc_ptr = d_xloc.data();
377 Real* d_yloc_ptr = d_yloc.data();
378 Real* d_freestream_velocity_ptr = d_freestream_velocity.data();
379 Real* d_disk_cell_count_ptr = d_disk_cell_count.data();
381 int n_bld_sections =
static_cast<int>(
bld_rad_loc.size());
383 Gpu::DeviceVector<Real> d_bld_rad_loc(n_bld_sections);
384 Gpu::DeviceVector<Real> d_bld_twist(n_bld_sections);
385 Gpu::DeviceVector<Real> d_bld_chord(n_bld_sections);
388 Gpu::copy(Gpu::hostToDevice,
bld_twist.begin(),
bld_twist.end(), d_bld_twist.begin());
389 Gpu::copy(Gpu::hostToDevice,
bld_chord.begin(),
bld_chord.end(), d_bld_chord.begin());
391 Real* bld_rad_loc_ptr = d_bld_rad_loc.data();
392 Real* bld_twist_ptr = d_bld_twist.data();
393 Real* bld_chord_ptr = d_bld_chord.data();
395 Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_aoa(n_bld_sections);
396 Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_Cl(n_bld_sections);
397 Vector<Gpu::DeviceVector<Real>> d_bld_airfoil_Cd(n_bld_sections);
399 Vector<int> h_n_pts_airfoil(n_bld_sections);
400 Gpu::DeviceVector<int> d_n_pts_airfoil(n_bld_sections);
402 for (
int i = 0; i < n_bld_sections; ++i) {
409 Gpu::copy(Gpu::hostToDevice,
411 d_bld_airfoil_aoa[i].begin());
413 Gpu::copy(Gpu::hostToDevice,
415 d_bld_airfoil_Cl[i].begin());
417 Gpu::copy(Gpu::hostToDevice,
419 d_bld_airfoil_Cd[i].begin());
422 Gpu::copy(Gpu::hostToDevice,
423 h_n_pts_airfoil.begin(), h_n_pts_airfoil.end(),
424 d_n_pts_airfoil.begin());
426 int* d_n_pts_airfoil_ptr = d_n_pts_airfoil.data();
428 Vector<Real*> hp_bld_airfoil_aoa, hp_bld_airfoil_Cl, hp_bld_airfoil_Cd;
429 for (
auto & v :d_bld_airfoil_aoa) {
430 hp_bld_airfoil_aoa.push_back(v.data());
432 for (
auto & v :d_bld_airfoil_Cl) {
433 hp_bld_airfoil_Cl.push_back(v.data());
435 for (
auto & v :d_bld_airfoil_Cd) {
436 hp_bld_airfoil_Cd.push_back(v.data());
439 Gpu::AsyncArray<Real*> aoa(hp_bld_airfoil_aoa.data(), n_bld_sections);
440 Gpu::AsyncArray<Real*> Cl(hp_bld_airfoil_Cl.data(), n_bld_sections);
441 Gpu::AsyncArray<Real*> Cd(hp_bld_airfoil_Cd.data(), n_bld_sections);
443 auto d_bld_airfoil_aoa_ptr = aoa.data();
444 auto d_bld_airfoil_Cl_ptr = Cl.data();
445 auto d_bld_airfoil_Cd_ptr = Cd.data();
447 int n_spec_extra =
static_cast<int>(
velocity.size());
449 Gpu::DeviceVector<Real> d_velocity(n_spec_extra);
450 Gpu::DeviceVector<Real> d_rotor_RPM(n_spec_extra);
451 Gpu::DeviceVector<Real> d_blade_pitch(n_spec_extra);
453 Gpu::copy(Gpu::hostToDevice,
velocity.begin(),
velocity.end(), d_velocity.begin());
454 Gpu::copy(Gpu::hostToDevice,
rotor_RPM.begin(),
rotor_RPM.end(), d_rotor_RPM.begin());
457 auto d_velocity_ptr = d_velocity.data();
458 auto d_rotor_RPM_ptr = d_rotor_RPM.data();
459 auto d_blade_pitch_ptr = d_blade_pitch.data();
461 for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
463 const Box& gbx = mfi.growntilebox(1);
464 auto SMark_array = mf_SMark.array(mfi);
465 auto generalAD_array = mf_vars_generalAD.array(mfi);
467 ParallelFor(gbx, [=] AMREX_GPU_DEVICE(
int i,
int j,
int k) noexcept {
468 int ii = amrex::min(amrex::max(i, domlo_x), domhi_x);
469 int jj = amrex::min(amrex::max(j, domlo_y), domhi_y);
470 int kk = amrex::min(amrex::max(k, domlo_z), domhi_z);
480 std::array<Real,2> Fn_and_Ft;
482 for(
long unsigned int it=0;it<nturbs;it++) {
483 Real avg_vel = d_freestream_velocity_ptr[it]/(d_disk_cell_count_ptr[it] +
Real(1e-10));
484 Real phi = d_turb_disk_angle;
487 if(SMark_array(ii,jj,kk,1) ==
static_cast<double>(it)) {
491 Real rad = std::pow( (
x-d_xloc_ptr[it])*(
x-d_xloc_ptr[it]) +
492 (
y-d_yloc_ptr[it])*(
y-d_yloc_ptr[it]) +
493 (
z-d_hub_height)*(
z-d_hub_height),
myhalf );
500 if(rad >=
two and rad <= d_rotor_rad) {
506 Real vec_proj = (
x-d_xloc_ptr[it])*(std::sin(phi)) +
507 (
y-d_yloc_ptr[it])*(-std::cos(phi));
510 Real zeta = std::atan2(
z-d_hub_height, vec_proj);
517 d_bld_airfoil_aoa_ptr[index],
518 d_bld_airfoil_Cl_ptr[index],
519 d_bld_airfoil_Cd_ptr[index],
520 d_n_pts_airfoil_ptr[index],
530 Real Fx = Fn*std::cos(phi) + Ft*std::sin(zeta)*std::sin(phi);
531 Real Fy = Fn*std::sin(phi) - Ft*std::sin(zeta)*std::cos(phi);
532 Real Fz = -Ft*std::cos(zeta);
547 amrex::Error(
"Actuator disks are overlapping. Visualize actuator_disks.vtk "
548 "and check the windturbine locations input file. Exiting..");
551 generalAD_array(i,j,k,0) = source_x;
552 generalAD_array(i,j,k,1) = source_y;
553 generalAD_array(i,j,k,2) = source_z;
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:180
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:213
const Real dx
Definition: ERF_InitCustomPert_ABL.H:44
constexpr amrex::Real three
Definition: ERF_NumericalConstants.H:32
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
constexpr amrex::Real PI
Definition: ERF_NumericalConstants.H:39
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