From 84ef97e309c99bdf48656f811449b15e48a7c20d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 2 Oct 2024 07:08:53 -0700 Subject: [PATCH] Some corrections to wind models (#1839) * Correcting compilation error with TINY_PROFILE=TRUE * Changes to Make.ERF * Merging with dev * Getting ready for generalized actuator disk * Framework for reading in blade data * Removing print * Minor changes to error hooks and docs * Correcting errors in turbine positioning * Removing unwanted prints * Correcting error hook * Minor correction * Correcting CMake errors --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan --- CMake/BuildERFExe.cmake | 2 + Docs/sphinx_doc/theory/WindFarmModels.rst | 2 +- Exec/Make.ERF | 6 + Source/DataStructs/ERF_DataStruct.H | 12 +- Source/ERF.H | 2 +- Source/ERF_make_new_arrays.cpp | 3 + Source/IO/ERF_Plotfile.cpp | 48 +++++- Source/Initialization/ERF_init_windfarm.cpp | 9 +- .../ERF_InitWindFarm.cpp | 35 ++-- Source/WindFarmParametrization/ERF_WindFarm.H | 16 +- .../Fitch/ERF_AdvanceFitch.cpp | 2 +- .../ERF_AdvanceGeneralAD.cpp | 149 ++++++++++++++++++ .../GeneralActuatorDisk/ERF_GeneralAD.H | 46 ++++++ .../GeneralActuatorDisk/Make.package | 2 + .../Null/ERF_NullWindFarm.H | 21 +-- .../ERF_AdvanceSimpleAD.cpp | 3 +- 16 files changed, 323 insertions(+), 35 deletions(-) create mode 100644 Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp create mode 100644 Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H create mode 100644 Source/WindFarmParametrization/GeneralActuatorDisk/Make.package diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index e284a2115..dedb436bb 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -196,6 +196,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp + ${SRC_DIR}/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp ${SRC_DIR}/LandSurfaceModel/SLM/ERF_SLM.cpp ${SRC_DIR}/LandSurfaceModel/MM5/ERF_MM5.cpp ) @@ -251,6 +252,7 @@ endif() target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) + target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) diff --git a/Docs/sphinx_doc/theory/WindFarmModels.rst b/Docs/sphinx_doc/theory/WindFarmModels.rst index 952995547..ee041bed2 100644 --- a/Docs/sphinx_doc/theory/WindFarmModels.rst +++ b/Docs/sphinx_doc/theory/WindFarmModels.rst @@ -4,7 +4,7 @@ Wind farm models Introduction ------------- -ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) and Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) are supported. +ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently the Fitch model (`Fitch et al. 2012`_), Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) and Simplified actuator disk model (See Section 3.2 in `Wind Energy Handbook 2nd edition`_) are supported. .. _Fitch model: diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 5d129369a..6c8a8328e 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -152,12 +152,14 @@ ERF_WINDFARM_NULL_DIR = $(ERF_WINDFARM_DIR)/Null ERF_WINDFARM_FITCH_DIR = $(ERF_WINDFARM_DIR)/Fitch ERF_WINDFARM_EWP_DIR = $(ERF_WINDFARM_DIR)/EWP ERF_WINDFARM_SIMPLEAD_DIR = $(ERF_WINDFARM_DIR)/SimpleActuatorDisk +ERF_WINDFARM_GENERALAD_DIR = $(ERF_WINDFARM_DIR)/GeneralActuatorDisk include $(ERF_WINDFARM_DIR)/Make.package include $(ERF_WINDFARM_NULL_DIR)/Make.package include $(ERF_WINDFARM_FITCH_DIR)/Make.package include $(ERF_WINDFARM_EWP_DIR)/Make.package include $(ERF_WINDFARM_SIMPLEAD_DIR)/Make.package +include $(ERF_WINDFARM_GENERALAD_DIR)/Make.package VPATH_LOCATIONS += $(ERF_WINDFARM_DIR) INCLUDE_LOCATIONS += $(ERF_WINDFARM_DIR) @@ -173,6 +175,10 @@ INCLUDE_LOCATIONS += $(ERF_WINDFARM_EWP_DIR) VPATH_LOCATIONS += $(ERF_WINDFARM_SIMPLEAD_DIR) INCLUDE_LOCATIONS += $(ERF_WINDFARM_SIMPLEAD_DIR) + +VPATH_LOCATIONS += $(ERF_WINDFARM_GENERALAD_DIR) +INCLUDE_LOCATIONS += $(ERF_WINDFARM_GENERALAD_DIR) + endif ifeq ($(USE_HEFFTE),TRUE) diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 3c245b8bc..70f2c7ca6 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -38,7 +38,7 @@ enum struct MoistureType { }; enum struct WindFarmType { - Fitch, EWP, SimpleAD, None + Fitch, EWP, SimpleAD, GeneralAD, None }; enum struct WindFarmLocType{ @@ -359,9 +359,12 @@ struct SolverChoice { else if (windfarm_type_string == "SimpleActuatorDisk") { windfarm_type = WindFarmType::SimpleAD; } + else if (windfarm_type_string == "GeneralActuatorDisk") { + windfarm_type = WindFarmType::GeneralAD; + } else if (windfarm_type_string != "None") { amrex::Abort("Are you using windfarms? Dont know this windfarm_type. windfarm_type" - " has to be Fitch or EWP or None."); + " has to be Fitch, EWP, SimpleActuatorDisk, GeneralActuatorDisk or None."); } static std::string windfarm_loc_type_string = "None"; @@ -380,6 +383,8 @@ struct SolverChoice { pp.query("windfarm_loc_table", windfarm_loc_table); pp.query("windfarm_spec_table", windfarm_spec_table); + pp.query("windfarm_blade_table", windfarm_blade_table); + pp.query("windfarm_airofil_tables", windfarm_airfoil_tables); // Sampling distance upstream of the turbine to find the // incoming free stream velocity as a factor of the diameter of the @@ -424,7 +429,7 @@ struct SolverChoice { " distance as a factor of the turbine diameter at which the incoming free stream" " velocity will be computed at."); } - if (windfarm_type==WindFarmType::SimpleAD and turb_disk_angle < 0.0) { + if((windfarm_type==WindFarmType::SimpleAD or windfarm_type==WindFarmType::GeneralAD) and turb_disk_angle < 0.0) { amrex::Abort("To use simplified actuator disks, you need to provide a variable" " erf.turb_disk_angle_from_x in the inputs which is the angle of the face of the" " turbine disk from the x-axis. A turbine facing an oncoming flow in the x-direction" @@ -667,6 +672,7 @@ struct SolverChoice { int RhoQr_comp {-1}; std::string windfarm_loc_table, windfarm_spec_table; + std::string windfarm_blade_table, windfarm_airfoil_tables; amrex::Real sampling_distance_by_D = -1.0; amrex::Real turb_disk_angle = -1.0; amrex::Real windfarm_x_shift = -1.0; diff --git a/Source/ERF.H b/Source/ERF.H index 0d950437f..84a0e33d5 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -888,7 +888,7 @@ private: "vorticity_x","vorticity_y","vorticity_z", "magvel", "divU", "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", - "eq_pot_temp", "num_turb", "SMark", + "eq_pot_temp", "num_turb", "SMark0", "SMark1", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", "lat_m", "lon_m", // Time averaged velocity diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index fbb627bb5..6d485ebb0 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -234,6 +234,9 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, } if (solverChoice.windfarm_type == WindFarmType::SimpleAD) { vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt + } + if (solverChoice.windfarm_type == WindFarmType::GeneralAD) { + vars_windfarm[lev].define(ba, dm, 2, ngrow_state);// dudt, dvdt } Nturb[lev].define(ba, dm, 1, ngrow_state); // Number of turbines in a cell SMark[lev].define(ba, dm, 2, ngrow_state); // Free stream velocity/source term diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 7b9f76def..c87587be7 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -83,7 +83,9 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector derived_names[i] != "graup_accum") ) { if ( (solverChoice.moisture_type != MoistureType::None) || - (derived_names[i] != "qv" && derived_names[i] != "qc") ) { + (derived_names[i] != "qv" && derived_names[i] != "qc" and + derived_names[i] != "num_turb" and derived_names[i] != "SMark0" and + derived_names[i] != "SMark1") ) { tmp_plot_names.push_back(derived_names[i]); } } // moisture_type @@ -91,6 +93,28 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector } // hasElement } +#ifdef ERF_USE_WINDFARM + for (int i = 0; i < derived_names.size(); ++i) { + if ( containerHasElement(plot_var_names, derived_names[i]) ) { + if(solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP) { + if(derived_names[i] == "num_turb") { + tmp_plot_names.push_back(derived_names[i]); + } + } + if(solverChoice.windfarm_type == WindFarmType::SimpleAD) { + if(derived_names[i] == "num_turb" or derived_names[i] == "SMark0" or derived_names[i] == "Smark1") { + tmp_plot_names.push_back(derived_names[i]); + } + } + if(solverChoice.windfarm_type == WindFarmType::GeneralAD) { + if(derived_names[i] == "num_turb" or derived_names[i] == "SMark1") { + tmp_plot_names.push_back(derived_names[i]); + } + } + } + } +#endif + #ifdef ERF_USE_PARTICLES const auto& particles_namelist( particleData.getNamesUnalloc() ); for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) { @@ -444,7 +468,9 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } #ifdef ERF_USE_WINDFARM - if (containerHasElement(plot_var_names, "num_turb")) + if (containerHasElement(plot_var_names, "num_turb") and + (solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or + solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD)) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -461,7 +487,8 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } - if(containerHasElement(plot_var_names, "SMark") and solverChoice.windfarm_type == WindFarmType::SimpleAD) { + if(containerHasElement(plot_var_names, "SMark0") and + solverChoice.windfarm_type == WindFarmType::SimpleAD) { for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); @@ -473,6 +500,21 @@ ERF::WritePlotFile (int which, Vector plot_var_names) } mf_comp ++; } + + if(containerHasElement(plot_var_names, "SMark1") and + (solverChoice.windfarm_type == WindFarmType::SimpleAD or solverChoice.windfarm_type == WindFarmType::GeneralAD)) { + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& SMark_array = SMark[lev].const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = SMark_array(i,j,k,1); + }); + } + mf_comp ++; + } + #endif int klo = geom[lev].Domain().smallEnd(2); diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index ac239b6b0..bf7d3dcc2 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/Source/Initialization/ERF_init_windfarm.cpp @@ -29,16 +29,23 @@ ERF::init_windfarm (int lev) true, false); } + windfarm->fill_Nturb_multifab(geom[lev], Nturb[lev]); windfarm->write_turbine_locations_vtk(); - if(solverChoice.windfarm_type == WindFarmType::SimpleAD) { + if(solverChoice.windfarm_type == WindFarmType::SimpleAD or + solverChoice.windfarm_type == WindFarmType::GeneralAD) { windfarm->fill_SMark_multifab(geom[lev], SMark[lev], solverChoice.sampling_distance_by_D, solverChoice.turb_disk_angle); windfarm->write_actuator_disks_vtk(geom[lev]); } + + if(solverChoice.windfarm_type == WindFarmType::GeneralAD) { + //windfarm->read_windfarm_blade_table(solverChoice.windfarm_blade_table); + //windfarm->read_airfoil_tables + } } void diff --git a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp index 40d9a4009..32fae15ce 100644 --- a/Source/WindFarmParametrization/ERF_InitWindFarm.cpp +++ b/Source/WindFarmParametrization/ERF_InitWindFarm.cpp @@ -42,7 +42,7 @@ WindFarm::read_windfarm_locations_table(const std::string windfarm_loc_table, } else { amrex::Abort("Are you using windfarms? For windfarm simulations, the inputs need to have an" - " entry erf.windfarm_loc_table which should not be None. \n"); + " entry erf.windfarm_loc_type which should be either lat_lon or x_y. \n"); } set_turb_loc(xloc, yloc); @@ -138,6 +138,8 @@ WindFarm::init_windfarm_x_y (const std::string windfarm_loc_table) Real value1, value2; while (file >> value1 >> value2) { + value1 = value1 + 1e-3; + value2 = value2 + 1e-3; xloc.push_back(value1); yloc.push_back(value2); } @@ -193,6 +195,20 @@ WindFarm::read_windfarm_spec_table(const std::string windfarm_spec_table) } +void +WindFarm::read_windfarm_blade_table(const std::string windfarm_blade_table) +{ + std::ifstream filename(windfarm_blade_table); + if (!filename.is_open()) { + Error("You are using a generalized actuator disk model based on blade element theory. This needs info of blades." + " An entry erf.windfarm_blade_table is needed. Either the entry is missing or the file specified" + " in the entry - " + windfarm_blade_table + " is missing."); + } + else { + Print() << "Reading in wind farm blade table: " << windfarm_blade_table << "\n"; + } +} + void WindFarm::fill_Nturb_multifab(const Geometry& geom, MultiFab& mf_Nturb) @@ -369,22 +385,19 @@ WindFarm::write_actuator_disks_vtk(const Geometry& geom) } fprintf(file_actuator_disks_in_dom, "%s %ld %s\n", "POINTS", static_cast(num_turb_in_dom*npts), "float"); - Real nx = std::cos(my_turb_disk_angle); - Real ny = std::sin(my_turb_disk_angle); + Real nx = std::cos(my_turb_disk_angle+0.5*M_PI); + Real ny = std::sin(my_turb_disk_angle+0.5*M_PI); for(int it=0; it ProbLoArr[0] and xloc[it] < ProbHiArr[0] and yloc[it] > ProbLoArr[1] and yloc[it] < ProbHiArr[1]) { - Real xp = (x-xloc[it])*nx - (y-yloc[it])*ny; - Real yp = (x-xloc[it])*nx + (y-yloc[it])*ny; - fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", xloc[it]+xp, yloc[it]+yp, z); + fprintf(file_actuator_disks_in_dom, "%0.15g %0.15g %0.15g\n", x, y, z); } } } diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H index c2bdb33f6..ee7742de9 100644 --- a/Source/WindFarmParametrization/ERF_WindFarm.H +++ b/Source/WindFarmParametrization/ERF_WindFarm.H @@ -9,6 +9,7 @@ #include "ERF_Fitch.H" #include "ERF_EWP.H" #include "ERF_SimpleAD.H" +#include "ERF_GeneralAD.H" class WindFarm : public NullWindFarm { @@ -24,14 +25,21 @@ public: m_windfarm_model.resize(nlev); if (a_windfarm_type == WindFarmType::Fitch) { SetModel(); + amrex::Print() << "Fitch windfarm model!\n"; } else if (a_windfarm_type == WindFarmType::EWP) { SetModel(); + amrex::Print() << "EWP windfarm model!\n"; } else if (a_windfarm_type == WindFarmType::SimpleAD) { SetModel(); - amrex::Print() << "Simple actuator disk windfarm model!\n"; - } else { + amrex::Print() << "Simplified actuator disk windfarm model!\n"; + } + else if (a_windfarm_type == WindFarmType::GeneralAD) { + SetModel(); + amrex::Print() << "Generalized actuator disk windfarm model!\n"; + } + else { amrex::Abort("WindFarm: Dont know this windfarm_type!") ; } } @@ -55,6 +63,10 @@ public: void read_windfarm_spec_table(const std::string windfarm_spec_table); + void read_windfarm_blade_table(const std::string windfarm_blade_table); + + void read_windfarm_airofil_tables(const std::string windfarm_airfoils_tables); + void fill_Nturb_multifab(const amrex::Geometry& geom, amrex::MultiFab& mf_Nturb); diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp index 7770ea32d..5fbf55232 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp @@ -14,7 +14,7 @@ Real compute_A(const Real z, Real d = std::min(std::fabs(z - hub_height), rotor_rad); Real theta = std::acos(d/rotor_rad); - Real A_s = rotor_rad*rotor_rad*theta - d*std::pow(rotor_rad*rotor_rad - d*d, 0.5); + Real A_s = rotor_rad*rotor_rad*theta - d*std::pow(std::max(rotor_rad*rotor_rad - d*d,0.0), 0.5); Real A = PI*rotor_rad*rotor_rad/2.0 - A_s; return A; diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp new file mode 100644 index 000000000..96243d29b --- /dev/null +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp @@ -0,0 +1,149 @@ +#include +#include +#include + +using namespace amrex; + +void +GeneralAD::advance (const Geometry& geom, + const Real& dt_advance, + MultiFab& cons_in, + MultiFab& mf_vars_generalAD, + MultiFab& U_old, + MultiFab& V_old, + MultiFab& W_old, + const MultiFab& mf_Nturb, + const MultiFab& mf_SMark) +{ + AMREX_ALWAYS_ASSERT(W_old.nComp() > 0); + AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0); + AMREX_ALWAYS_ASSERT(mf_vars_generalAD.nComp() > 0); + source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD); + update(dt_advance, cons_in, U_old, V_old, mf_vars_generalAD); +} + +void +GeneralAD::update (const Real& dt_advance, + MultiFab& cons_in, + MultiFab& U_old, MultiFab& V_old, + const MultiFab& mf_vars_generalAD) +{ + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + Box tbx = mfi.nodaltilebox(0); + Box tby = mfi.nodaltilebox(1); + + auto generalAD_array = mf_vars_generalAD.array(mfi); + auto u_vel = U_old.array(mfi); + auto v_vel = V_old.array(mfi); + + ParallelFor(tbx, tby, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + u_vel(i,j,k) = u_vel(i,j,k) + (generalAD_array(i-1,j,k,0) + generalAD_array(i,j,k,0))/2.0*dt_advance; + }, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + v_vel(i,j,k) = v_vel(i,j,k) + (generalAD_array(i,j-1,k,1) + generalAD_array(i,j,k,1))/2.0*dt_advance; + }); + } +} + +void +GeneralAD::source_terms_cellcentered (const Geometry& geom, + const MultiFab& cons_in, + const MultiFab& mf_SMark, + MultiFab& mf_vars_generalAD) +{ + + get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, + wind_speed, thrust_coeff, power); + + Gpu::DeviceVector d_xloc(xloc.size()); + Gpu::DeviceVector d_yloc(yloc.size()); + Gpu::copy(Gpu::hostToDevice, xloc.begin(), xloc.end(), d_xloc.begin()); + Gpu::copy(Gpu::hostToDevice, yloc.begin(), yloc.end(), d_yloc.begin()); + + auto dx = geom.CellSizeArray(); + + // Domain valid box + const amrex::Box& domain = geom.Domain(); + int domlo_x = domain.smallEnd(0); + int domhi_x = domain.bigEnd(0) + 1; + int domlo_y = domain.smallEnd(1); + int domhi_y = domain.bigEnd(1) + 1; + int domlo_z = domain.smallEnd(2); + int domhi_z = domain.bigEnd(2) + 1; + + // The order of variables are - Vabs dVabsdt, dudt, dvdt, dTKEdt + mf_vars_generalAD.setVal(0.0); + + long unsigned int nturbs = xloc.size(); + + get_turb_disk_angle(turb_disk_angle); + Real nx = -std::cos(turb_disk_angle); + Real ny = -std::sin(turb_disk_angle); + Real d_turb_disk_angle = turb_disk_angle; + + Gpu::DeviceVector d_wind_speed(wind_speed.size()); + Gpu::DeviceVector d_thrust_coeff(thrust_coeff.size()); + + // Copy data from host vectors to device vectors + Gpu::copy(Gpu::hostToDevice, wind_speed.begin(), wind_speed.end(), d_wind_speed.begin()); + Gpu::copy(Gpu::hostToDevice, thrust_coeff.begin(), thrust_coeff.end(), d_thrust_coeff.begin()); + + const Real* wind_speed_d = d_wind_speed.dataPtr(); + const Real* thrust_coeff_d = d_thrust_coeff.dataPtr(); + const int n_spec_table = d_wind_speed.size(); + + for ( MFIter mfi(cons_in,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const Box& gbx = mfi.growntilebox(1); + auto SMark_array = mf_SMark.array(mfi); + auto generalAD_array = mf_vars_generalAD.array(mfi); + + ParallelFor(gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int ii = amrex::min(amrex::max(i, domlo_x), domhi_x); + int jj = amrex::min(amrex::max(j, domlo_y), domhi_y); + int kk = amrex::min(amrex::max(k, domlo_z), domhi_z); + + int check_int = 0; + + Real source_x = 0.0; + Real source_y = 0.0; + + for(long unsigned int it=0;it(it)) { + check_int++; + if(C_T <= 1) { + source_x = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); + source_y = -2.0*std::pow(Uinfty_dot_nhat, 2.0)*a*(1.0-a)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + } + else { + source_x = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::cos(phi); + source_y = -0.5*C_T*std::pow(Uinfty_dot_nhat, 2.0)*dx[1]*dx[2]*std::cos(d_turb_disk_angle)/(dx[0]*dx[1]*dx[2])*std::sin(phi); + } + } + } + if(check_int > 1){ + amrex::Error("Actuator disks are overlapping. Visualize actuator_disks.vtk " + "and check the windturbine locations input file. Exiting.."); + } + + generalAD_array(i,j,k,0) = source_x; + generalAD_array(i,j,k,1) = source_y; + }); + } +} diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H new file mode 100644 index 000000000..f3ab40227 --- /dev/null +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H @@ -0,0 +1,46 @@ +#ifndef ERF_GENERALAD_H +#define ERF_GENERALAD_H + +#include +#include +#include "ERF_NullWindFarm.H" + +class GeneralAD : public NullWindFarm { + +public: + + GeneralAD() {} + + virtual ~GeneralAD() = default; + + void advance (const amrex::Geometry& geom, + const amrex::Real& 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) override; + + 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 amrex::Real& dt_advance, + amrex::MultiFab& cons_in, + amrex::MultiFab& U_old, + amrex::MultiFab& V_old, + const amrex::MultiFab& mf_vars); + +protected: + amrex::Vector xloc, yloc; + amrex::Real turb_disk_angle; + amrex::Real hub_height, rotor_rad, thrust_coeff_standing, nominal_power; + amrex::Vector wind_speed, thrust_coeff, power; + amrex::Vector freestream_velocity, freestream_phi, disk_cell_count; +}; + +#endif + diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/Make.package b/Source/WindFarmParametrization/GeneralActuatorDisk/Make.package new file mode 100644 index 000000000..072a3d4f8 --- /dev/null +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += ERF_AdvanceGeneralAD.cpp +CEXE_headers += ERF_GeneralAD.H diff --git a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H index f2b0b23f5..934832b19 100644 --- a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H +++ b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H @@ -81,32 +81,33 @@ bool find_if_marked(amrex::Real x1, amrex::Real x2, amrex::Real y1, amrex::Real amrex::Real z) { - // Plane is (x-x0)*nx + (y-y0*ny = 0. And the planes to intersect are - // x = x1, x = x2, y = y1, y = y2 + // Plane is (x-x0)*nx + (y-y0)*ny = 0. And the planes to intersect are + // x = x1, x = x2, y = y1, y = y2. The intersections are + // (xval1,y1), (xval2,y2), (x1,yval1) and (x2,yval2) - amrex::Real xval1 = x0+1e-12 - (y1-y0)*ny/(nx+1e-10); - amrex::Real xval2 = x0+1e-12 - (y2-y0)*ny/(nx+1e-10); + amrex::Real xval1 = x0 - (y1-y0)*ny/(nx+1e-10); + amrex::Real xval2 = x0 - (y2-y0)*ny/(nx+1e-10); - amrex::Real yval1 = y0+1e-12 - (x1-x0)*nx/(ny+1e-10); - amrex::Real yval2 = y0+1e-12 - (x2-x0)*nx/(ny+1e-10); + amrex::Real yval1 = y0 - (x1-x0)*nx/(ny+1e-10); + amrex::Real yval2 = y0 - (x2-x0)*nx/(ny+1e-10); if( xval1 >=x1 and xval1 <=x2 ) { - if(std::pow((y1-y0)*(y1-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { + if(std::pow((xval1-x0)*(xval1-x0) + (y1-y0)*(y1-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { return true; } } if( xval2 >=x1 and xval2 <=x2 ) { - if(std::pow((y2-y0)*(y2-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { + if(std::pow((xval2-x0)*(xval2-x0) + (y2-y0)*(y2-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { return true; } } if (yval1 >=y1 and yval1 <=y2) { - if(std::pow((yval1-y0)*(yval1-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { + if(std::pow((x1-x0)*(x1-x0) + (yval1-y0)*(yval1-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { return true; } } if (yval2 >=y1 and yval2 <=y2) { - if(std::pow((yval2-y0)*(yval2-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { + if(std::pow((x2-x0)*(x2-x0) + (yval2-y0)*(yval2-y0) + (z-d_hub_height)*(z-d_hub_height),0.5) < d_rotor_rad ) { return true; } } diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp index 8aa569a03..38a10d65e 100644 --- a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp +++ b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -114,7 +114,6 @@ void SimpleAD::compute_freestream_velocity(const MultiFab& cons_in, amrex::ParallelContext::CommunicatorAll()); get_turb_loc(xloc, yloc); - std::cout << "xloc size is " << xloc.size() << "\n"; if (ParallelDescriptor::IOProcessor()){ for(int it=0; it