From 1d7b18b18477a02cca29311d71246cb20a520bfe Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 29 Oct 2024 14:22:21 -0700 Subject: [PATCH] Computing power output for wind farms (#1915) * Adding power output for simple actuator disk * Adding time to power output --------- Co-authored-by: Mahesh Natarajan Co-authored-by: Mahesh Natarajan --- Exec/GeneralActuatorDisk/ERF_prob.cpp | 164 ++++++++++++++++++ Source/ERF.H | 3 +- Source/Initialization/ERF_init_windfarm.cpp | 5 +- Source/TimeIntegration/ERF_Advance.cpp | 2 +- Source/WindFarmParametrization/ERF_WindFarm.H | 5 +- .../EWP/ERF_AdvanceEWP.cpp | 4 +- Source/WindFarmParametrization/EWP/ERF_EWP.H | 3 +- .../Fitch/ERF_AdvanceFitch.cpp | 4 +- .../WindFarmParametrization/Fitch/ERF_Fitch.H | 3 +- .../ERF_AdvanceGeneralAD.cpp | 4 +- .../GeneralActuatorDisk/ERF_GeneralAD.H | 3 +- .../Null/ERF_NullWindFarm.H | 3 +- .../ERF_AdvanceSimpleAD.cpp | 33 +++- .../SimpleActuatorDisk/ERF_SimpleAD.H | 5 +- 14 files changed, 226 insertions(+), 15 deletions(-) create mode 100644 Exec/GeneralActuatorDisk/ERF_prob.cpp diff --git a/Exec/GeneralActuatorDisk/ERF_prob.cpp b/Exec/GeneralActuatorDisk/ERF_prob.cpp new file mode 100644 index 000000000..3eeff25fa --- /dev/null +++ b/Exec/GeneralActuatorDisk/ERF_prob.cpp @@ -0,0 +1,164 @@ +#include "ERF_prob.H" +#include "AMReX_Random.H" + +using namespace amrex; + +std::unique_ptr +amrex_probinit(const amrex_real* problo, const amrex_real* probhi) +{ + return std::make_unique(problo, probhi); +} + +Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi) +{ + // Parse params + ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("T_0", parms.T_0); + pp.query("A_0", parms.A_0); + pp.query("QKE_0", parms.QKE_0); + + pp.query("U_0", parms.U_0); + pp.query("V_0", parms.V_0); + pp.query("W_0", parms.W_0); + pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag); + pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag); + pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag); + pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag); + + pp.query("pert_deltaU", parms.pert_deltaU); + pp.query("pert_deltaV", parms.pert_deltaV); + pp.query("pert_periods_U", parms.pert_periods_U); + pp.query("pert_periods_V", parms.pert_periods_V); + pp.query("pert_ref_height", parms.pert_ref_height); + parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]); + parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]); + parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height; + parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height; + + init_base_parms(parms.rho_0, parms.T_0); +} + +void +Problem::init_custom_pert( + const amrex::Box& bx, + const amrex::Box& xbx, + const amrex::Box& ybx, + const amrex::Box& zbx, + amrex::Array4 const& /*state*/, + amrex::Array4 const& state_pert, + amrex::Array4 const& x_vel_pert, + amrex::Array4 const& y_vel_pert, + amrex::Array4 const& z_vel_pert, + amrex::Array4 const& /*r_hse*/, + amrex::Array4 const& /*p_hse*/, + amrex::Array4 const& /*z_nd*/, + amrex::Array4 const& /*z_cc*/, + amrex::GeometryData const& geomdata, + amrex::Array4 const& /*mf_m*/, + amrex::Array4 const& /*mf_u*/, + amrex::Array4 const& /*mf_v*/, + const SolverChoice& sc) +{ + const bool use_moisture = (sc.moisture_type != MoistureType::None); + + ParallelForRNG(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = 0.5 * (prob_lo[0] + prob_hi[0]); + const Real yc = 0.5 * (prob_lo[1] + prob_hi[1]); + const Real zc = 0.5 * (prob_lo[2] + prob_hi[2]); + + const Real r = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + + // Add temperature perturbations + if ((z <= parms_d.pert_ref_height) && (parms_d.T_0_Pert_Mag != 0.0)) { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms_d.T_0_Pert_Mag; + } + + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-10.*r*r); + + // Set an initial value for QKE + //state_pert(i, j, k, RhoQKE_comp) = parms_d.QKE_0; + + if (use_moisture) { + state_pert(i, j, k, RhoQ1_comp) = 0.0; + state_pert(i, j, k, RhoQ2_comp) = 0.0; + } + }); + + // Set the x-velocity + ParallelForRNG(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the x-velocity + x_vel_pert(i, j, k) = parms_d.U_0; + if ((z <= parms_d.pert_ref_height) && (parms_d.U_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real x_vel_prime = (rand_double*2.0 - 1.0)*parms_d.U_0_Pert_Mag; + x_vel_pert(i, j, k) += x_vel_prime; + } + if (parms_d.pert_deltaU != 0.0) + { + const amrex::Real yl = y - prob_lo[1]; + const amrex::Real zl = z / parms_d.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + x_vel_pert(i, j, k) += parms_d.ufac * damp * z * std::cos(parms_d.aval * yl); + } + }); + + // Set the y-velocity + ParallelForRNG(ybx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the y-velocity + y_vel_pert(i, j, k) = parms_d.V_0; + if ((z <= parms_d.pert_ref_height) && (parms_d.V_0_Pert_Mag != 0.0)) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real y_vel_prime = (rand_double*2.0 - 1.0)*parms_d.V_0_Pert_Mag; + y_vel_pert(i, j, k) += y_vel_prime; + } + if (parms_d.pert_deltaV != 0.0) + { + const amrex::Real xl = x - prob_lo[0]; + const amrex::Real zl = z / parms_d.pert_ref_height; + const amrex::Real damp = std::exp(-0.5 * zl * zl); + y_vel_pert(i, j, k) += parms_d.vfac * damp * z * std::cos(parms_d.bval * xl); + } + }); + + // Set the z-velocity + ParallelForRNG(zbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept { + const int dom_lo_z = geomdata.Domain().smallEnd()[2]; + const int dom_hi_z = geomdata.Domain().bigEnd()[2]; + + // Set the z-velocity + if (k == dom_lo_z || k == dom_hi_z+1) + { + z_vel_pert(i, j, k) = 0.0; + } + else if (parms_d.W_0_Pert_Mag != 0.0) + { + Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0 + Real z_vel_prime = (rand_double*2.0 - 1.0)*parms_d.W_0_Pert_Mag; + z_vel_pert(i, j, k) = parms_d.W_0 + z_vel_prime; + } + }); +} diff --git a/Source/ERF.H b/Source/ERF.H index cafc46489..3c3dda751 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -435,7 +435,8 @@ public: amrex::MultiFab& U_old, amrex::MultiFab& V_old, amrex::MultiFab& W_old, amrex::MultiFab& mf_vars_windfarm, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark); + const amrex::MultiFab& mf_SMark, + const amrex::Real& time); #endif #ifdef ERF_USE_EB diff --git a/Source/Initialization/ERF_init_windfarm.cpp b/Source/Initialization/ERF_init_windfarm.cpp index bc41acf57..b87ededd4 100644 --- a/Source/Initialization/ERF_init_windfarm.cpp +++ b/Source/Initialization/ERF_init_windfarm.cpp @@ -59,8 +59,9 @@ ERF::advance_windfarm (const Geometry& a_geom, MultiFab& W_old, MultiFab& mf_vars_windfarm, const MultiFab& mf_Nturb, - const MultiFab& mf_SMark) + const MultiFab& mf_SMark, + const Real& time) { windfarm->advance(a_geom, dt_advance, cons_in, mf_vars_windfarm, - U_old, V_old, W_old, mf_Nturb, mf_SMark); + U_old, V_old, W_old, mf_Nturb, mf_SMark, time); } diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index db7cbf1a0..87fc9ee79 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -110,7 +110,7 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/) #if defined(ERF_USE_WINDFARM) if (solverChoice.windfarm_type != WindFarmType::None) { advance_windfarm(Geom(lev), dt_lev, S_old, - U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], SMark[lev]); + U_old, V_old, W_old, vars_windfarm[lev], Nturb[lev], SMark[lev], time); } #endif diff --git a/Source/WindFarmParametrization/ERF_WindFarm.H b/Source/WindFarmParametrization/ERF_WindFarm.H index 100a81f0a..945483bbe 100644 --- a/Source/WindFarmParametrization/ERF_WindFarm.H +++ b/Source/WindFarmParametrization/ERF_WindFarm.H @@ -90,10 +90,11 @@ public: amrex::MultiFab& V_old, amrex::MultiFab& W_old, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark) override + const amrex::MultiFab& mf_SMark, + const amrex::Real& time) override { m_windfarm_model[0]->advance(a_geom, dt_advance, cons_in, mf_vars_windfarm, - U_old, V_old, W_old, mf_Nturb, mf_SMark); + U_old, V_old, W_old, mf_Nturb, mf_SMark, time); } void set_turb_spec(const amrex::Real& a_rotor_rad, const amrex::Real& a_hub_height, diff --git a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp index 0522183f5..8ba23af86 100644 --- a/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp +++ b/Source/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp @@ -14,9 +14,11 @@ EWP::advance (const Geometry& geom, MultiFab& V_old, MultiFab& W_old, const MultiFab& mf_Nturb, - const MultiFab& mf_SMark) + const MultiFab& mf_SMark, + const Real& time) { AMREX_ALWAYS_ASSERT(mf_SMark.nComp() > 0); + AMREX_ALWAYS_ASSERT(time > -1.0); source_terms_cellcentered(geom, cons_in, mf_vars_ewp, U_old, V_old, W_old, mf_Nturb); update(dt_advance, cons_in, U_old, V_old, mf_vars_ewp); } diff --git a/Source/WindFarmParametrization/EWP/ERF_EWP.H b/Source/WindFarmParametrization/EWP/ERF_EWP.H index c62dd33e4..fbf10efca 100644 --- a/Source/WindFarmParametrization/EWP/ERF_EWP.H +++ b/Source/WindFarmParametrization/EWP/ERF_EWP.H @@ -22,7 +22,8 @@ public: amrex::MultiFab& V_old, amrex::MultiFab& W_old, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark) override; + const amrex::MultiFab& mf_SMark, + const amrex::Real& time) override; void source_terms_cellcentered (const amrex::Geometry& geom, const amrex::MultiFab& cons_in, diff --git a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp index 7947c847a..68cce5558 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp +++ b/Source/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp @@ -53,10 +53,12 @@ Fitch::advance (const Geometry& geom, MultiFab& V_old, MultiFab& W_old, const MultiFab& mf_Nturb, - const MultiFab& mf_SMark) + const MultiFab& mf_SMark, + const Real& time) { AMREX_ALWAYS_ASSERT(W_old.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_SMark.nComp() > 0); + AMREX_ALWAYS_ASSERT(time > -1.0); source_terms_cellcentered(geom, cons_in, mf_vars_fitch, U_old, V_old, W_old, mf_Nturb); update(dt_advance, cons_in, U_old, V_old, mf_vars_fitch); } diff --git a/Source/WindFarmParametrization/Fitch/ERF_Fitch.H b/Source/WindFarmParametrization/Fitch/ERF_Fitch.H index ca450ecb9..906a55835 100644 --- a/Source/WindFarmParametrization/Fitch/ERF_Fitch.H +++ b/Source/WindFarmParametrization/Fitch/ERF_Fitch.H @@ -22,7 +22,8 @@ public: amrex::MultiFab& V_old, amrex::MultiFab& W_old, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark) override; + const amrex::MultiFab& mf_SMark, + const amrex::Real& time) override; void source_terms_cellcentered (const amrex::Geometry& geom, const amrex::MultiFab& cons_in, diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp index aaa64ceb2..3c4da9eb1 100644 --- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_AdvanceGeneralAD.cpp @@ -14,11 +14,13 @@ GeneralAD::advance (const Geometry& geom, MultiFab& V_old, MultiFab& W_old, const MultiFab& mf_Nturb, - const MultiFab& mf_SMark) + const MultiFab& mf_SMark, + const Real& time) { AMREX_ALWAYS_ASSERT(W_old.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_vars_generalAD.nComp() > 0); + AMREX_ALWAYS_ASSERT(time > -1.0); compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark); source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_generalAD); update(dt_advance, cons_in, U_old, V_old, W_old, mf_vars_generalAD); diff --git a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H index 593b16a48..9092ea312 100644 --- a/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H +++ b/Source/WindFarmParametrization/GeneralActuatorDisk/ERF_GeneralAD.H @@ -21,7 +21,8 @@ public: amrex::MultiFab& V_old, amrex::MultiFab& W_old, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark) override; + const amrex::MultiFab& mf_SMark, + const amrex::Real& time) override; void compute_freestream_velocity (const amrex::MultiFab& cons_in, const amrex::MultiFab& U_old, diff --git a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H index 2daea0aa5..f0714eef3 100644 --- a/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H +++ b/Source/WindFarmParametrization/Null/ERF_NullWindFarm.H @@ -22,7 +22,8 @@ public: amrex::MultiFab& V_old, amrex::MultiFab& W_old, const amrex::MultiFab& mf_Nturb, - const amrex::MultiFab& mf_SMark) = 0; + const amrex::MultiFab& mf_SMark, + const amrex::Real& time) = 0; virtual void set_turb_spec(const amrex::Real& rotor_rad, const amrex::Real& hub_height, const amrex::Real& thrust_coeff_standing, const amrex::Vector& wind_speed, diff --git a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp index 7709f6355..038022c7b 100644 --- a/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp +++ b/Source/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -13,13 +13,44 @@ SimpleAD::advance (const Geometry& geom, MultiFab& V_old, MultiFab& W_old, const MultiFab& mf_Nturb, - const MultiFab& mf_SMark) + const MultiFab& mf_SMark, + const Real& time) { AMREX_ALWAYS_ASSERT(W_old.nComp() > 0); AMREX_ALWAYS_ASSERT(mf_Nturb.nComp() > 0); compute_freestream_velocity(cons_in, U_old, V_old, mf_SMark); source_terms_cellcentered(geom, cons_in, mf_SMark, mf_vars_simpleAD); update(dt_advance, cons_in, U_old, V_old, mf_vars_simpleAD); + compute_power_output(time); +} + +void +SimpleAD::compute_power_output(const Real& time) +{ + get_turb_loc(xloc, yloc); + get_turb_spec(rotor_rad, hub_height, thrust_coeff_standing, + wind_speed, thrust_coeff, power); + + const int n_spec_table = wind_speed.size(); + // Compute power based on the look-up table + + if (ParallelDescriptor::IOProcessor()){ + static std::ofstream file("power_output.txt", std::ios::app); + // Check if the file opened successfully + if (!file.is_open()) { + std::cerr << "Error opening file!" << std::endl; + Abort("Could not open file to write power output in ERF_AdvanceSimpleAD.cpp"); + } + Real total_power = 0.0; + for(int it=0; it xloc, yloc; amrex::Real turb_disk_angle;