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());
347 auto dx = geom.CellSizeArray();
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;
360 mf_vars_generalAD.setVal(0.0);
362 long unsigned int nturbs =
static_cast<long unsigned int>(
xloc.size());
370 Gpu::DeviceVector<Real> d_freestream_velocity(nturbs);
371 Gpu::DeviceVector<Real> d_disk_cell_count(nturbs);
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();
380 int n_bld_sections =
static_cast<int>(
bld_rad_loc.size());
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);
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());
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();
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);
398 Vector<int> h_n_pts_airfoil(n_bld_sections);
399 Gpu::DeviceVector<int> d_n_pts_airfoil(n_bld_sections);
401 for (
int i = 0; i < n_bld_sections; ++i) {
408 Gpu::copy(Gpu::hostToDevice,
410 d_bld_airfoil_aoa[i].begin());
412 Gpu::copy(Gpu::hostToDevice,
414 d_bld_airfoil_Cl[i].begin());
416 Gpu::copy(Gpu::hostToDevice,
418 d_bld_airfoil_Cd[i].begin());
421 Gpu::copy(Gpu::hostToDevice,
422 h_n_pts_airfoil.begin(), h_n_pts_airfoil.end(),
423 d_n_pts_airfoil.begin());
425 int* d_n_pts_airfoil_ptr = d_n_pts_airfoil.data();
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());
431 for (
auto & v :d_bld_airfoil_Cl) {
432 hp_bld_airfoil_Cl.push_back(v.data());
434 for (
auto & v :d_bld_airfoil_Cd) {
435 hp_bld_airfoil_Cd.push_back(v.data());
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);
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();
446 int n_spec_extra =
static_cast<int>(
velocity.size());
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);
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());
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();
460 for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
462 const Box& gbx = mfi.growntilebox(1);
463 auto SMark_array = mf_SMark.array(mfi);
464 auto generalAD_array = mf_vars_generalAD.array(mfi);
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);
479 std::array<Real,2> Fn_and_Ft;
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;
486 if(SMark_array(ii,jj,kk,1) ==
static_cast<double>(it)) {
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 );
499 if(rad >=
two and rad <= d_rotor_rad) {
505 Real vec_proj = (
x-d_xloc_ptr[it])*(std::sin(phi)) +
506 (
y-d_yloc_ptr[it])*(-std::cos(phi));
509 Real zeta = std::atan2(
z-d_hub_height, vec_proj);
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],
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);
546 amrex::Error(
"Actuator disks are overlapping. Visualize actuator_disks.vtk "
547 "and check the windturbine locations input file. Exiting..");
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;
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