From bdddf133e41a9b7b4c8ce28f1ea1bebec47678f5 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Tue, 30 Apr 2024 13:21:33 -0600 Subject: [PATCH] Time varying abl bodyforce (#1025) * changes to write abl forcing to text file, add to existing reg test * reading in body force from file * setting up reg test * formatting * using start time and frequency arguments * write to file with Geostrophic Forcing * formatting * reverting last change, because Geostrophic Forcing is reversible * time-varying target velocity for geostrophic wind * reg test for time-varying geostrophic wind * updating to use underscores, but keeping backwards-compatible * removing unused variable * proliferate m_coriolis_factor (fix issue from last commit) * input file documentation * unit test start, better error feedback * test to write forces file * body forcing unit test + better precision * geostrophic forcing unit test * Correct source term implementations for temporal discretization * make body force test self-contained * add headers to time table files in reg tests, note in documentation * unused variables * avoiding angle interpolation issues * Unit tests for new angle interpolation feature --- .../icns/source_terms/ABLForcing.H | 41 +- .../icns/source_terms/ABLForcing.cpp | 16 +- .../icns/source_terms/BodyForce.H | 11 +- .../icns/source_terms/BodyForce.cpp | 67 ++- .../icns/source_terms/GeostrophicForcing.H | 16 + .../icns/source_terms/GeostrophicForcing.cpp | 55 ++- .../temperature/source_terms/BodyForce.H | 2 +- .../temperature/source_terms/BodyForce.cpp | 10 +- amr-wind/utilities/linear_interpolation.H | 65 +++ amr-wind/wind_energy/ABLFieldInit.cpp | 4 +- docs/sphinx/user/inputs_Momentum_Sources.rst | 89 +++- test/CMakeLists.txt | 2 + .../abl_godunov_forcetimetable.inp | 91 ++++ .../abl_godunov_geostrophic_timetable.inp | 77 ++++ .../time_table.txt | 4 + .../abl_godunov_timetable.inp | 11 + .../abl_godunov_timetable/time_table.txt | 1 + .../utilities/test_linear_interpolation.cpp | 36 ++ unit_tests/wind_energy/CMakeLists.txt | 1 + unit_tests/wind_energy/test_abl_src.cpp | 12 +- .../wind_energy/test_abl_src_timetable.cpp | 394 ++++++++++++++++++ 21 files changed, 969 insertions(+), 36 deletions(-) create mode 100644 test/test_files/abl_godunov_forcetimetable/abl_godunov_forcetimetable.inp create mode 100644 test/test_files/abl_godunov_geostrophic_timetable/abl_godunov_geostrophic_timetable.inp create mode 100644 test/test_files/abl_godunov_geostrophic_timetable/time_table.txt create mode 100644 unit_tests/wind_energy/test_abl_src_timetable.cpp diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H index e24892434d..9df4f08b77 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H @@ -41,20 +41,38 @@ public: m_mean_vel[1] = uy; const auto& current_time = m_time.current_time(); + const auto& new_time = m_time.new_time(); + const auto& nph_time = 0.5 * (current_time + new_time); const auto& dt = m_time.deltaT(); + const auto& t_step = m_time.time_index(); if (!m_vel_timetable.empty()) { - const amrex::Real current_spd = ::amr_wind::interp::linear( - m_time_table, m_speed_table, current_time); - const amrex::Real current_dir = ::amr_wind::interp::linear( - m_time_table, m_direction_table, current_time); - - m_target_vel[0] = current_spd * std::cos(current_dir); - m_target_vel[1] = current_spd * std::sin(current_dir); + // Forces should be applied at n+1/2. Because ABL forcing is a + // delta, the difference between the target velocity (at n+1) and + // the current velocity (at n) puts the force term at n+1/2 + const amrex::Real new_spd = ::amr_wind::interp::linear( + m_time_table, m_speed_table, new_time); + const amrex::Real new_dir = ::amr_wind::interp::linear_angle( + m_time_table, m_direction_table, new_time, 2.0 * M_PI); + + m_target_vel[0] = new_spd * std::cos(new_dir); + m_target_vel[1] = new_spd * std::sin(new_dir); } m_abl_forcing[0] = (m_target_vel[0] - m_mean_vel[0]) / dt; m_abl_forcing[1] = (m_target_vel[1] - m_mean_vel[1]) / dt; + + if (m_write_force_timetable && + amrex::ParallelDescriptor::IOProcessor() && + (t_step % m_force_outfreq == 0) && + (current_time >= m_force_outstart)) { + std::ofstream outfile; + // Forces are recorded at n+1/2 + outfile.open(m_force_timetable, std::ios::out | std::ios::app); + outfile << std::fixed << std::setprecision(15) << nph_time << "\t" + << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t" + << 0.0 << std::endl; + } } amrex::RealArray abl_forcing() const { return m_abl_forcing; } @@ -77,6 +95,15 @@ private: //! File name for velocity forcing time table std::string m_vel_timetable; + //! Bool for writing forcing time table + bool m_write_force_timetable{false}; + //! File name for forcing time table output + std::string m_force_timetable; + //! Output frequency for forces + int m_force_outfreq{1}; + //! Output start time for force + amrex::Real m_force_outstart{0.0}; + //! Velocity forcing time table amrex::Vector m_time_table; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp index 3da67364ab..79ad1fbbd6 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp @@ -26,7 +26,9 @@ ABLForcing::ABLForcing(const CFDSim& sim) if (!m_vel_timetable.empty()) { std::ifstream ifh(m_vel_timetable, std::ios::in); if (!ifh.good()) { - amrex::Abort("Cannot find input file: " + m_vel_timetable); + amrex::Abort( + "Cannot find ABLForcing velocity_timetable file: " + + m_vel_timetable); } amrex::Real data_time; amrex::Real data_speed; @@ -43,6 +45,18 @@ ABLForcing::ABLForcing(const CFDSim& sim) pp_incflo.getarr("velocity", m_target_vel); } + m_write_force_timetable = pp_abl.contains("forcing_timetable_output_file"); + if (m_write_force_timetable) { + pp_abl.get("forcing_timetable_output_file", m_force_timetable); + pp_abl.query("forcing_timetable_frequency", m_force_outfreq); + pp_abl.query("forcing_timetable_start_time", m_force_outstart); + if (amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream outfile; + outfile.open(m_force_timetable, std::ios::out); + outfile << "time\tfx\tfy\tfz\n"; + } + } + for (int i = 0; i < AMREX_SPACEDIM; ++i) { m_mean_vel[i] = m_target_vel[i]; } diff --git a/amr-wind/equation_systems/icns/source_terms/BodyForce.H b/amr-wind/equation_systems/icns/source_terms/BodyForce.H index 30592044a5..eae6960871 100644 --- a/amr-wind/equation_systems/icns/source_terms/BodyForce.H +++ b/amr-wind/equation_systems/icns/source_terms/BodyForce.H @@ -33,6 +33,7 @@ public: const amrex::Array4& src_term) const override; void read_bforce_profile(const std::string& filename); + void read_bforce_timetable(const std::string& filename); private: //! Time @@ -44,7 +45,9 @@ private: amrex::Vector m_body_force{{0.0, 0.0, 0.0}}; //! Body Force Type - std::string m_type{"constant"}; + std::string m_type{"uniform_constant"}; + //! Uniform time table file + std::string m_utt_file; //! Angular frequency used in the oscillatory forcing amrex::Real m_omega{0.0}; @@ -54,6 +57,12 @@ private: amrex::Gpu::DeviceVector m_prof_x; amrex::Gpu::DeviceVector m_prof_y; amrex::Gpu::DeviceVector m_ht; + + //! Vectors for storing uniform_timetable inputs + amrex::Vector m_time_table; + amrex::Vector m_fx_table; + amrex::Vector m_fy_table; + amrex::Vector m_fz_table; }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp b/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp index 871ca255b4..66e160d20d 100644 --- a/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp +++ b/amr-wind/equation_systems/icns/source_terms/BodyForce.cpp @@ -1,6 +1,7 @@ #include "amr-wind/equation_systems/icns/source_terms/BodyForce.H" #include "amr-wind/CFDSim.H" #include "amr-wind/utilities/trig_ops.H" +#include "amr-wind/utilities/linear_interpolation.H" #include "AMReX_ParmParse.H" #include "AMReX_Gpu.H" @@ -16,18 +17,38 @@ namespace amr_wind::pde::icns { BodyForce::BodyForce(const CFDSim& sim) : m_time(sim.time()), m_mesh(sim.mesh()) { - // Read the geostrophic wind speed vector (in m/s) + // Body Force arguments amrex::ParmParse pp("BodyForce"); pp.query("type", m_type); m_type = amrex::toLower(m_type); - - if (m_type == "height-varying") { - pp.get("bodyforce-file", m_bforce_file); + bool no_type_specified = !pp.contains("type"); + bool file_specified = pp.contains("uniform_timetable_file"); + + // Prepare type of body force distribution + if (m_type == "height_varying" || m_type == "height-varying") { + // Constant in time, varies with z + // Using underscores is preferred, remains backwards compatible + pp.query("bodyforce-file", m_bforce_file); + if (m_bforce_file.empty()) { + pp.get("bodyforce_file", m_bforce_file); + } read_bforce_profile(m_bforce_file); + } else if ( + m_type == "uniform_timetable" || + (no_type_specified && file_specified)) { + // Still used if type not specified but file is + // Uniform in space, varies with time + pp.get("uniform_timetable_file", m_utt_file); + read_bforce_timetable(m_utt_file); } else { pp.getarr("magnitude", m_body_force); if (m_type == "oscillatory") { pp.get("angular_frequency", m_omega); + } else if (m_type != "uniform_constant") { + amrex::Abort( + "BodyForce type not supported. Please choose uniform_constant " + "(default), height_varying, oscillatory, or " + "uniform_timetable.\n"); } } } @@ -68,6 +89,28 @@ void BodyForce::read_bforce_profile(const std::string& filename) m_ht.begin()); } +void BodyForce::read_bforce_timetable(const std::string& filename) +{ + std::ifstream ifh(filename, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find BodyForce uniform_timetable_file: " + filename + "\n"); + } + amrex::Real data_time; + amrex::Real data_fx; + amrex::Real data_fy; + amrex::Real data_fz; + // Skip first line (header) + ifh.ignore(std::numeric_limits::max(), '\n'); + while (ifh >> data_time) { + ifh >> data_fx >> data_fy >> data_fz; + m_time_table.push_back(data_time); + m_fx_table.push_back(data_fx); + m_fy_table.push_back(data_fy); + m_fz_table.push_back(data_fz); + } +} + void BodyForce::operator()( const int lev, const amrex::MFIter& /*mfi*/, @@ -75,7 +118,7 @@ void BodyForce::operator()( const FieldState /*fstate*/, const amrex::Array4& src_term) const { - const auto& time = m_time.current_time(); + const auto& nph_time = 0.5 * (m_time.current_time() + m_time.new_time()); const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& dx = m_mesh.Geom(lev).CellSizeArray(); const int lp1 = lev + 1; @@ -85,7 +128,7 @@ void BodyForce::operator()( const amrex::Real* force_x = m_prof_x.data(); const amrex::Real* force_y = m_prof_y.data(); - if (m_type == "height-varying") { + if (m_type == "height_varying" || m_type == "height-varying") { amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -113,8 +156,18 @@ void BodyForce::operator()( amrex::GpuArray forcing{ {m_body_force[0], m_body_force[1], m_body_force[2]}}; + if (!m_utt_file.empty()) { + // Populate forcing from file if supplied + forcing[0] = + amr_wind::interp::linear(m_time_table, m_fx_table, nph_time); + forcing[1] = + amr_wind::interp::linear(m_time_table, m_fy_table, nph_time); + forcing[2] = + amr_wind::interp::linear(m_time_table, m_fz_table, nph_time); + } + amrex::Real coeff = - (m_type == "oscillatory") ? std::cos(m_omega * time) : 1.0; + (m_type == "oscillatory") ? std::cos(m_omega * nph_time) : 1.0; amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { src_term(i, j, k, 0) += coeff * forcing[0]; diff --git a/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H b/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H index badb846f79..59f5546f8c 100644 --- a/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H @@ -26,8 +26,24 @@ public: const amrex::Array4& src_term) const override; private: + const SimTime& m_time; const amrex::AmrCore& m_mesh; + //! File name for target velocity time table + std::string m_vel_timetable; + + //! Velocity forcing time table + amrex::Vector m_time_table; + + //! Velocity forcing speed table + amrex::Vector m_speed_table; + + //! Velocity forcing direction table + amrex::Vector m_direction_table; + + //! Coriolis factor + amrex::Real m_coriolis_factor; + //! Activated when water is present in domain bool m_use_phase_ramp{false}; diff --git a/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.cpp b/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.cpp index bb7003d4a1..6ccff2847f 100644 --- a/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.cpp @@ -4,6 +4,7 @@ #include "amr-wind/core/vs/vstraits.H" #include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/equation_systems/vof/volume_fractions.H" +#include "amr-wind/utilities/linear_interpolation.H" #include "AMReX_ParmParse.H" #include "AMReX_Gpu.H" @@ -22,10 +23,9 @@ namespace amr_wind::pde::icns { * GeostrophicForcing namespace * */ -GeostrophicForcing::GeostrophicForcing(const CFDSim& sim) : m_mesh(sim.mesh()) +GeostrophicForcing::GeostrophicForcing(const CFDSim& sim) + : m_time(sim.time()), m_mesh(sim.mesh()) { - amrex::Real coriolis_factor = 0.0; - // Read the rotational time period (in seconds) amrex::ParmParse ppc("CoriolisForcing"); // Read the rotational time period (in seconds) -- This is 23hrs and 56 @@ -39,18 +39,40 @@ GeostrophicForcing::GeostrophicForcing(const CFDSim& sim) : m_mesh(sim.mesh()) latitude = utils::radians(latitude); amrex::Real sinphi = std::sin(latitude); - coriolis_factor = 2.0 * omega * sinphi; + m_coriolis_factor = 2.0 * omega * sinphi; ppc.query("is_horizontal", m_is_horizontal); amrex::Print() << "Geostrophic forcing: Coriolis factor = " - << coriolis_factor << std::endl; + << m_coriolis_factor << std::endl; // Read the geostrophic wind speed vector (in m/s) amrex::ParmParse ppg("GeostrophicForcing"); - ppg.getarr("geostrophic_wind", m_target_vel); + ppg.query("geostrophic_wind_timetable", m_vel_timetable); + if (!m_vel_timetable.empty()) { + std::ifstream ifh(m_vel_timetable, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find GeostrophicForcing geostrophic_wind_timetable " + "file: " + + m_vel_timetable); + } + amrex::Real data_time; + amrex::Real data_speed; + amrex::Real data_deg; + ifh.ignore(std::numeric_limits::max(), '\n'); + while (ifh >> data_time) { + ifh >> data_speed >> data_deg; + amrex::Real data_rad = utils::radians(data_deg); + m_time_table.push_back(data_time); + m_speed_table.push_back(data_speed); + m_direction_table.push_back(data_rad); + } + } else { + ppg.getarr("geostrophic_wind", m_target_vel); + } m_g_forcing = { - -coriolis_factor * m_target_vel[1], coriolis_factor * m_target_vel[0], - 0.0}; + -m_coriolis_factor * m_target_vel[1], + m_coriolis_factor * m_target_vel[0], 0.0}; // Set up relaxation toward 0 forcing near the air-water interface if (sim.repo().field_exists("vof")) { @@ -83,6 +105,8 @@ void GeostrophicForcing::operator()( const amrex::Array4& src_term) const { amrex::Real hfac = (m_is_horizontal) ? 0. : 1.; + // Forces applied at n+1/2 + const auto& nph_time = 0.5 * (m_time.current_time() + m_time.new_time()); const bool ph_ramp = m_use_phase_ramp; const int n_band = m_n_band; @@ -94,6 +118,21 @@ void GeostrophicForcing::operator()( amrex::GpuArray forcing{ {m_g_forcing[0], m_g_forcing[1], m_g_forcing[2]}}; + // Calculate forcing values if target velocity is a function of time + if (!m_vel_timetable.empty()) { + const amrex::Real nph_spd = + amr_wind::interp::linear(m_time_table, m_speed_table, nph_time); + const amrex::Real nph_dir = amr_wind::interp::linear_angle( + m_time_table, m_direction_table, nph_time, 2.0 * M_PI); + + const amrex::Real target_u = nph_spd * std::cos(nph_dir); + const amrex::Real target_v = nph_spd * std::sin(nph_dir); + + forcing[0] = -m_coriolis_factor * target_v; + forcing[1] = m_coriolis_factor * target_u; + forcing[2] = 0.0; + } + const auto& vof = (*m_vof)(lev).const_array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::Real wfac = 1.0; diff --git a/amr-wind/equation_systems/temperature/source_terms/BodyForce.H b/amr-wind/equation_systems/temperature/source_terms/BodyForce.H index 403de38d46..0d8f07f114 100644 --- a/amr-wind/equation_systems/temperature/source_terms/BodyForce.H +++ b/amr-wind/equation_systems/temperature/source_terms/BodyForce.H @@ -33,7 +33,7 @@ private: const SimTime& m_time; const amrex::AmrCore& m_mesh; - std::string m_type{"height-varying"}; + std::string m_type{"height_varying"}; std::string m_bforce_file; size_t m_bforce_profile_nhts; diff --git a/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp b/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp index bdcbdaa9c0..ea198c5d76 100644 --- a/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp +++ b/amr-wind/equation_systems/temperature/source_terms/BodyForce.cpp @@ -23,8 +23,12 @@ BodyForce::BodyForce(const CFDSim& sim) : m_time(sim.time()), m_mesh(sim.mesh()) pp.query("type", m_type); m_type = amrex::toLower(m_type); - if (m_type == "height-varying") { - pp.get("bodyforce-file", m_bforce_file); + // Updated to underscores (more common) but still backwards-compatible + if (m_type == "height_varying" || m_type == "height-varying") { + pp.query("bodyforce-file", m_bforce_file); + if (m_bforce_file.empty()) { + pp.get("bodyforce_file", m_bforce_file); + } read_bforce_profile(m_bforce_file); } } @@ -74,7 +78,7 @@ void BodyForce::operator()( const amrex::Real* force_ht = m_ht.data(); const amrex::Real* force_theta = m_prof_theta.data(); - if (m_type == "height-varying") { + if (m_type == "height_varying" || m_type == "height-varying") { amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { diff --git a/amr-wind/utilities/linear_interpolation.H b/amr-wind/utilities/linear_interpolation.H index efe50af72d..40c5994a83 100644 --- a/amr-wind/utilities/linear_interpolation.H +++ b/amr-wind/utilities/linear_interpolation.H @@ -153,6 +153,71 @@ inline void linear(const C1& xinp, const C2& yinp, const C1& xout, C2& yout) } } +// Angles in yinp should be defined between [-upper_bound, +upper_bound] +// Use cases would be degrees (upper_bound =360) and radians (upper_bound =2pi) +template +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + typename std::iterator_traits::value_type + linear_angle( + const C1 xbegin, + const C1 xend, + const C2 yinp, + const typename std::iterator_traits::value_type& xout, + const typename std::iterator_traits::value_type& upper_bound) +{ + using DType1 = typename std::iterator_traits::value_type; + const auto idx = bisection_search(xbegin, xend, xout); + + if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) { + return yinp[idx.idx]; + } + static constexpr DType1 eps = 1.0e-8; + const int j = idx.idx; + const auto denom = (xbegin[j + 1] - xbegin[j]); + /* https://math.stackexchange.com/questions/2144234/interpolating-between-2-angles + */ + const auto R = (denom > eps) ? ((xout - xbegin[j]) / denom) : 0.0; + const auto ub = upper_bound; + const auto hb = 0.5 * upper_bound; + const auto ohb = 1.5 * upper_bound; + return yinp[j] + + (std::fmod((std::fmod(yinp[j + 1] - yinp[j], ub) + ohb), ub) - hb) * + R; +} + +template +inline typename C2::value_type linear_angle( + const C1& xinp, + const C2& yinp, + const typename C1::value_type& xout, + const typename C1::value_type& upper_bound) +{ + return linear_angle( + xinp.data(), (xinp.data() + xinp.size()), yinp.data(), xout, + upper_bound); +} + +template +inline void linear_angle( + const C1& xinp, + const C2& yinp, + const C1& xout, + C2& yout, + const typename C1::value_type& upper_bound) +{ + AMREX_ASSERT( + static_cast(xinp.size()) == + static_cast(yinp.size())); + AMREX_ASSERT( + static_cast(xout.size()) == + static_cast(yout.size())); + + int npts = xout.size(); + for (int i = 0; i < npts; ++i) { + yout[i] = linear_angle(xinp, yinp, xout[i], upper_bound); + } +} + } // namespace amr_wind::interp #endif /* LINEAR_INTERPOLATION_H */ diff --git a/amr-wind/wind_energy/ABLFieldInit.cpp b/amr-wind/wind_energy/ABLFieldInit.cpp index 2f67b97d29..3063a946d7 100644 --- a/amr-wind/wind_energy/ABLFieldInit.cpp +++ b/amr-wind/wind_energy/ABLFieldInit.cpp @@ -69,7 +69,9 @@ void ABLFieldInit::initialize_from_inputfile() if (!m_vel_timetable.empty()) { std::ifstream ifh(m_vel_timetable, std::ios::in); if (!ifh.good()) { - amrex::Abort("Cannot find input file: " + m_vel_timetable); + amrex::Abort( + "Cannot find ABLForcing velocity_timetable file: " + + m_vel_timetable); } amrex::Real m_vel_time; amrex::Real m_vel_ang; diff --git a/docs/sphinx/user/inputs_Momentum_Sources.rst b/docs/sphinx/user/inputs_Momentum_Sources.rst index 041aaef56b..15e20386ea 100644 --- a/docs/sphinx/user/inputs_Momentum_Sources.rst +++ b/docs/sphinx/user/inputs_Momentum_Sources.rst @@ -9,7 +9,10 @@ Section: Momentum Sources equations. These strings can be entered in any order with a space between each. Please consult `AMR-Wind developer documentation `_ for a - comprehensive list of all momentum source terms available. + comprehensive list of all momentum source terms available. Note that the + following input arguments specific to each source term will only be active + if the corresponding source term (the root name) is listed in + :input_param:`ICNS.source_terms`. .. input_param:: BoussinesqBuoyancy.reference_temperature @@ -60,6 +63,19 @@ Section: Momentum Sources The user has to choose between GeostrophicForcing and ABLForcing. CoriolisForcing input must be present when using GeostrophicForcing. These checks are not enforced for now. + +.. input_param:: GeostrophicForcing.geostrophic_wind_timetable + + **type:** String, optional + + Input file name for table that lists time in seconds, wind speed + in meters per second, and horizontal wind direction in degrees of the Geostrophic + forcing velocity. Each line in the file should be a sequence of + three floats specifying the inputs in that order (e.g., 0.0 8.0 -5.0). If this + argument is present, the :input_param:`GeostrophicForcing.geostrophic_wind` + will be ignored. Note that the code expects there to be a single-line header + at the beginning of the geostrophic wind timetable file; if no header exists, + the first line of data will be ignored. .. input_param:: ABLForcing.abl_forcing_height @@ -73,6 +89,73 @@ Section: Momentum Sources **type:** String, optional Input file name for table that lists time in seconds, wind speed - in meters per second, and wind direction in degrees of the ABL + in meters per second, and horizontal wind direction in degrees of the ABL forcing velocity. Each line in the file should be a sequence of - three floats specifying the inputs in that order (e.g., 0.0 8.0 -5.0). \ No newline at end of file + three floats specifying the inputs in that order (e.g., 0.0 8.0 -5.0). If this + argument is present, the :input_param:`incflo.velocity` argument + will be ignored. Note that the code expects there to be a single-line header + at the beginning of the velocity timetable file; if no header exists, the first + line of data will be ignored. + +.. input_param:: ABLForcing.forcing_timetable_output_file + + **type:** String, optional + + Output file name for writing the ABL forcing vector to a text file over the course + of a simulation. This output is primarily intended for replicating the ABL forcing + from a precursor simulation in a subsequent inflow-outflow simulation by providing + an input file for Body Forcing. The output file will contain the time and three vector + components of the force. + +.. input_param:: ABLForcing.forcing_timetable_frequency + + **type:** Int, optional + + The interval of timesteps for writing to the forcing timetable output file. The default + is 1, i.e., writing every step, which is also the default of the boundary plane output feature. + +.. input_param:: ABLForcing.forcing_timetable_start_time + + **type:** Real, optional + + The start time for writing to the forcing timetable output file. The default is 0. + +.. input_param:: BodyForce.type + + **type:** String, optional + + The type of body force being used. The default is uniform_constant, which applies a single constant + force vector over the entire domain. Other available types are height_varying, oscillatory, and + uniform_constant. + +.. input_param:: BodyForce.magnitude + + **type:** List of 3 reals, conditionally mandatory + + The force vector to be applied as a body force. This argument is mandatory for uniform_constant + (default) and oscillatory body force types. + +.. input_param:: BodyForce.angular_frequency + + **type:** Real, conditionally mandatory + + The angular frequency to be used for applying sinusoidal time variation to the body force. This + argument is mandatory for the oscillatory body force type and is only active for the oscillatory type. + +.. input_param:: BodyForce.bodyforce_file + + **type:** String, conditionally mandatory + + The text file for specifying the body force vector as a function of height. This text file must contain + heights (z coordinate values), force components in x, and force components in y. This argument is mandatory for + the height_varying body force type and is only active for the height_varying type. + +.. input_param:: BodyForce.uniform_timetable_file + + **type:** String, conditionally mandatory + + The text file for specifying the body force vector as a function of time. This text file must contain + times, force components in x, force components in y, and force components in z. This argument is mandatory for + the uniform_timetable body force type and is only active for the uniform_timetable type. Note that the code + expects there to be a single-line header at the beginning of the uniform timetable file; if no header exists, + the first line of data will be ignored. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8b92fffd86..84744fdcf1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -174,6 +174,7 @@ add_test_re(abl_godunov_static_refinement_rr4) add_test_re(abl_godunov_scalar_velocity_solve) add_test_re(abl_godunov_segregated_velocity_solve) add_test_re(abl_godunov_timetable) +add_test_re(abl_godunov_geostrophic_timetable) add_test_re(abl_ksgsm84_godunov) add_test_re(abl_mol_cn) add_test_re(abl_mol_explicit) @@ -290,6 +291,7 @@ add_test_red(abl_bndry_input_native abl_bndry_output_native) add_test_red(abl_godunov_restart abl_godunov) add_test_red(abl_bndry_input_amr_native abl_bndry_output_native) add_test_red(abl_bndry_input_amr_native_mlbc abl_bndry_output_amr_native) +add_test_red(abl_godunov_forcetimetable abl_godunov_timetable) if(AMR_WIND_ENABLE_NETCDF) add_test_red(abl_bndry_input abl_bndry_output) diff --git a/test/test_files/abl_godunov_forcetimetable/abl_godunov_forcetimetable.inp b/test/test_files/abl_godunov_forcetimetable/abl_godunov_forcetimetable.inp new file mode 100644 index 0000000000..4e0383008e --- /dev/null +++ b/test/test_files/abl_godunov_forcetimetable/abl_godunov_forcetimetable.inp @@ -0,0 +1,91 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 5.0 # Max (simulated) time to evolve +time.max_step = 10 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.4 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +io.restart_file = "../abl_godunov_timetable/chk00005" +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = 5 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density + +incflo.use_godunov = 1 +incflo.diffusion_type = 2 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + + +incflo.physics = ABL +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing BodyForce +BoussinesqBuoyancy.reference_temperature = 300.0 +ABL.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 + +# ABL forcing is not turned on, but ABL Field Init needs velocity values +# This is required in the absence of incflo.velocity +ABLForcing.velocity_timetable = "../abl_godunov_timetable/time_table.txt" +# Note: the use of forcing timetable does not depend on velocity timetable; +# however, if a precursor uses velocity timetable, force timetable is +# recommended for subsequent simulations +BodyForce.type = uniform_timetable +BodyForce.uniform_timetable_file = "../abl_godunov_timetable/abl_forces.txt" + +ABL.temperature_heights = 650.0 750.0 1000.0 +ABL.temperature_values = 300.0 308.0 308.75 + +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.15 + +ABL.bndry_file = "../abl_godunov_timetable/bndry_files" +ABL.bndry_io_mode = 1 +ABL.bndry_var_names = velocity temperature +ABL.bndry_output_format = native + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1) + +# Boundary conditions +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.temperature = 0.0 +xhi.type = "pressure_outflow" +ylo.type = "mass_inflow" +ylo.density = 1.0 +ylo.temperature = 0.0 +yhi.type = "pressure_outflow" +zlo.type = "wall_model" +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level diff --git a/test/test_files/abl_godunov_geostrophic_timetable/abl_godunov_geostrophic_timetable.inp b/test/test_files/abl_godunov_geostrophic_timetable/abl_godunov_geostrophic_timetable.inp new file mode 100644 index 0000000000..8060bd806c --- /dev/null +++ b/test/test_files/abl_godunov_geostrophic_timetable/abl_godunov_geostrophic_timetable.inp @@ -0,0 +1,77 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 5.0 # Max (simulated) time to evolve +time.max_step = 10 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.5 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = 5 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density + +incflo.use_godunov = 1 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + + +incflo.physics = ABL +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing +BoussinesqBuoyancy.reference_temperature = 300.0 +ABL.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 +GeostrophicForcing.geostrophic_wind_timetable = "time_table.txt" +incflo.velocity = 8.0 0.0 0.0 + +ABL.temperature_heights = 650.0 750.0 1000.0 +ABL.temperature_values = 300.0 308.0 308.75 + +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.15 + +ABL.bndry_file = "bndry_files" +ABL.bndry_io_mode = 0 +ABL.bndry_planes = ylo xlo +ABL.bndry_output_start_time = 2.0 +ABL.bndry_var_names = velocity temperature +ABL.bndry_output_format = native + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) + +# Boundary conditions +zlo.type = "wall_model" + +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level diff --git a/test/test_files/abl_godunov_geostrophic_timetable/time_table.txt b/test/test_files/abl_godunov_geostrophic_timetable/time_table.txt new file mode 100644 index 0000000000..784835f3fb --- /dev/null +++ b/test/test_files/abl_godunov_geostrophic_timetable/time_table.txt @@ -0,0 +1,4 @@ +Time WindSpeed HorzAngle +0 8 0 +2 8 10 +5 8 20 \ No newline at end of file diff --git a/test/test_files/abl_godunov_timetable/abl_godunov_timetable.inp b/test/test_files/abl_godunov_timetable/abl_godunov_timetable.inp index 5702886f89..c618d7344e 100644 --- a/test/test_files/abl_godunov_timetable/abl_godunov_timetable.inp +++ b/test/test_files/abl_godunov_timetable/abl_godunov_timetable.inp @@ -39,6 +39,10 @@ CoriolisForcing.latitude = 41.3 ABLForcing.abl_forcing_height = 90 ABLForcing.velocity_timetable = "time_table.txt" +# Note: the use of forcing timetable does not depend on velocity timetable; +# however, if a precursor uses velocity timetable, force timetable is +# recommended for subsequent simulations +ABLForcing.forcing_timetable_output_file = "abl_forces.txt" ABL.temperature_heights = 650.0 750.0 1000.0 ABL.temperature_values = 300.0 308.0 308.75 @@ -46,6 +50,13 @@ ABL.temperature_values = 300.0 308.0 308.75 ABL.kappa = .41 ABL.surface_roughness_z0 = 0.15 +ABL.bndry_file = "bndry_files" +ABL.bndry_io_mode = 0 +ABL.bndry_planes = ylo xlo +ABL.bndry_output_start_time = 2.0 +ABL.bndry_var_names = velocity temperature +ABL.bndry_output_format = native + #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # ADAPTIVE MESH REFINEMENT # #.......................................# diff --git a/test/test_files/abl_godunov_timetable/time_table.txt b/test/test_files/abl_godunov_timetable/time_table.txt index 0178bbd48d..784835f3fb 100644 --- a/test/test_files/abl_godunov_timetable/time_table.txt +++ b/test/test_files/abl_godunov_timetable/time_table.txt @@ -1,3 +1,4 @@ +Time WindSpeed HorzAngle 0 8 0 2 8 10 5 8 20 \ No newline at end of file diff --git a/unit_tests/utilities/test_linear_interpolation.cpp b/unit_tests/utilities/test_linear_interpolation.cpp index fcc89708c3..a18cde5eb1 100644 --- a/unit_tests/utilities/test_linear_interpolation.cpp +++ b/unit_tests/utilities/test_linear_interpolation.cpp @@ -1,6 +1,7 @@ #include "aw_test_utils/AmrexTest.H" #include "amr-wind/utilities/linear_interpolation.H" +#include "amr-wind/utilities/trig_ops.H" #include "AMReX_Random.H" #include @@ -147,4 +148,39 @@ TEST(LinearInterpolation, lin_monotonic) } } +TEST(LinearInterpolation, lin_interp_angle) +{ + namespace interp = amr_wind::interp; + + const int vecsize = 8; + std::vector xvec(vecsize); + std::iota(xvec.begin(), xvec.end(), 0.0); + // Set up with angles that are awkward to interpolate between + std::vector yvec_deg{0.0, 30.0, 330.0, 300.0, + -15.0, -180.0, 120.0, -150.0}; + // Create duplicate list in radians + std::vector yvec_rad(vecsize); + for (size_t i = 0; i < yvec_rad.size(); ++i) { + yvec_rad[i] = amr_wind::utils::radians(yvec_deg[i]); + } + + // Set up output vectors, interp locations, and golds + std::vector xtest(vecsize), ytest_deg(vecsize), + ytest_rad(vecsize); + std::iota(xtest.begin(), xtest.end(), 0.5); + xtest[7] = 4.0; + std::vector ygold_deg{15., 0., 315., 322.5, + -97.5, -210, 165., 360. - 15.}; + + // Upper bound is 360 because y is in degrees + interp::linear_angle(xvec, yvec_deg, xtest, ytest_deg, 360.); + // Upper bound is 2pi because y is in radians + interp::linear_angle(xvec, yvec_rad, xtest, ytest_rad, 2. * M_PI); + for (size_t i = 0; i < xtest.size(); ++i) { + EXPECT_NEAR(ytest_deg[i], ygold_deg[i], 1.0e-12); + EXPECT_NEAR( + ytest_rad[i], amr_wind::utils::radians(ygold_deg[i]), 1.0e-12); + } +} + } // namespace amr_wind_tests diff --git a/unit_tests/wind_energy/CMakeLists.txt b/unit_tests/wind_energy/CMakeLists.txt index 7cbc153c5f..27b67c2541 100644 --- a/unit_tests/wind_energy/CMakeLists.txt +++ b/unit_tests/wind_energy/CMakeLists.txt @@ -7,6 +7,7 @@ target_sources(${amr_wind_unit_test_exe_name} PRIVATE test_abl_src.cpp test_abl_stats.cpp test_abl_bc.cpp + test_abl_src_timetable.cpp ) if (AMR_WIND_ENABLE_NETCDF) diff --git a/unit_tests/wind_energy/test_abl_src.cpp b/unit_tests/wind_energy/test_abl_src.cpp index d2684e9900..5c8f096989 100644 --- a/unit_tests/wind_energy/test_abl_src.cpp +++ b/unit_tests/wind_energy/test_abl_src.cpp @@ -151,7 +151,8 @@ TEST_F(ABLMeshTest, body_force) // Mimic source term at later timesteps { auto& time = sim().time(); - time.current_time() = 0.1; + time.new_timestep(); + time.current_time() = 0.05; src_term.setVal(0.0); run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { @@ -162,9 +163,12 @@ TEST_F(ABLMeshTest, body_force) }); const amrex::Array golds{ - {static_cast(1.0 * std::cos(0.1)), - static_cast(2.0 * std::cos(0.1)), - static_cast(3.0 * std::cos(0.1))}}; + {static_cast( + 1.0 * std::cos(0.5 * (0.05 + time.new_time()))), + static_cast( + 2.0 * std::cos(0.5 * (0.05 + time.new_time()))), + static_cast( + 3.0 * std::cos(0.5 * (0.05 + time.new_time())))}}; const auto valx2 = utils::field_max(src_term, 0); const auto valy2 = utils::field_max(src_term, 1); const auto valz2 = utils::field_max(src_term, 2); diff --git a/unit_tests/wind_energy/test_abl_src_timetable.cpp b/unit_tests/wind_energy/test_abl_src_timetable.cpp new file mode 100644 index 0000000000..70fefb7ef6 --- /dev/null +++ b/unit_tests/wind_energy/test_abl_src_timetable.cpp @@ -0,0 +1,394 @@ +#include "aw_test_utils/MeshTest.H" +#include "aw_test_utils/iter_tools.H" +#include "aw_test_utils/test_utils.H" + +#include "amr-wind/wind_energy/ABL.H" +#include "amr-wind/equation_systems/icns/icns.H" +#include "amr-wind/equation_systems/icns/icns_ops.H" +#include "amr-wind/equation_systems/icns/MomentumSource.H" +#include "amr-wind/equation_systems/icns/source_terms/BodyForce.H" +#include "amr-wind/equation_systems/icns/source_terms/ABLForcing.H" +#include "amr-wind/equation_systems/icns/source_terms/GeostrophicForcing.H" +#include "amr-wind/equation_systems/icns/source_terms/CoriolisForcing.H" + +namespace { +void write_target_velocity_file(const std::string& fname) +{ + std::ofstream os(fname); + // Write header **??** + os << "Time\tWindSpeed\tHorzAngle\n"; + // Write time table + os << "0.0\t8.0\t0.0\n"; + os << "0.1\t8.0\t360.0\n"; + os << "0.3\t8.0\t5.0\n"; + os << "0.5\t8.0\t12.0\n"; +} +void write_body_force_file(const std::string& fname) +{ + std::ofstream os(fname); + // Write header **??** + os << "time\tfx\tfy\tfz\n"; + // Write time table + os << "0.000000000000000\t0.000000000000000\t0.000000000000000\t0." + "000000000000000\n"; + os << "0.050000000000000\t0.000000000000000\t0.000000000000000\t0." + "000000000000000\n"; + os << "0.150000000000000\t-0.076142273451376\t3.489550989226880\t0." + "000000000000000\n"; + os << "0.250000000000000\t-0.304424152660356\t6.972459419812656\t0." + "000000000000000\n"; + os << "0.350000000000000\t-0.878730931046654\t11.824752890368849\t0." + "000000000000000\n"; +} +} // namespace + +namespace amr_wind_tests { + +using ICNSFields = + amr_wind::pde::FieldRegOp; + +// Tests in this file involve ABLForcing, BodyForce, and GeostrophicWind +class ABLSrcTimeTableTest : public MeshTest +{ +protected: + void populate_parameters() override + { + MeshTest::populate_parameters(); + + // Make computational domain like ABL mesh + { + amrex::ParmParse pp("amr"); + amrex::Vector ncell{{8, 8, 64}}; + pp.addarr("n_cell", ncell); + } + + { + amrex::ParmParse pp("geometry"); + amrex::Vector probhi{{120.0, 120.0, 1000.0}}; + pp.addarr("prob_hi", probhi); + } + + // Parameters for initializing ABL case + { + amrex::ParmParse pp("ABL"); + amrex::Vector theights{{0.0, 650.0, 750.0, 1000.0}}; + amrex::Vector tvalues{{300.0, 300.0, 308.0, 308.75}}; + pp.addarr("temperature_heights", theights); + pp.addarr("temperature_values", tvalues); + pp.add("perturb_ref_height", 50.0); + pp.add("reference_temperature", 300.0); + pp.add("kappa", 0.41); + pp.add("surface_roughness_z0", 0.1); + } + + { + // Physics + amrex::ParmParse pp("incflo"); + amrex::Vector phystr{"ABL"}; + pp.addarr("physics", phystr); + pp.add("density", 1.0); + } + + // Timestep size + { + amrex::ParmParse pp("time"); + pp.add("fixed_dt", 0.1); + } + } + std::string tvel_fname = "target_velocities.txt"; + std::string forces_fname = "abl_forces.txt"; + amrex::Real dt{0.1}; + int nsteps = 5; +}; + +TEST_F(ABLSrcTimeTableTest, abl) +{ + constexpr amrex::Real tol = 1.0e-12; + + // Write target wind file + write_target_velocity_file(tvel_fname); + + // Set up simulation parameters and mesh + populate_parameters(); + // ABL Forcing + { + amrex::ParmParse pp("ABLForcing"); + pp.add("abl_forcing_height", 90.0); + pp.add("velocity_timetable", tvel_fname); + } + initialize_mesh(); + + // Set up PDEs and physics objects + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + + // Get icns source term, which will be tested + auto& src_term = pde_mgr.icns().fields().src_term; + src_term.setVal(0.0); + // Source term object for ABLForcing + amr_wind::pde::icns::ABLForcing abl_forcing(sim()); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + abl_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + + // Initial velocity should be the same as target velocity, so force is 0 + amrex::Vector target_force = {0.0, 0.0, 0.0}; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + // Change mean velocity and recalculate + const amrex::Vector init_vel{8.0, 0.0, 0.0}; + const amrex::Vector new_vel{5.0, 3.0, 0.0}; + // This function actually calculates the forcing terms, applied later + // This function also is the one that writes forces to file + abl_forcing.set_mean_velocities(new_vel[0], new_vel[1]); + // Doing this through ABL physics is complicated and would require + // calc_averages in ABLStats, then pre_advance_work in ABL physics + src_term.setVal(0.0); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + abl_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + // Velocity at hub height is (5, 3) and target is (8, 0) + target_force[0] = (init_vel[0] - new_vel[0]) / dt; + target_force[1] = (init_vel[1] - new_vel[1]) / dt; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + // Advance time (twice to get to interpolation interval) + sim().time().new_timestep(); + sim().time().new_timestep(); + + // Go back to original mean values on mesh + abl_forcing.set_mean_velocities(init_vel[0], init_vel[1]); + // Recalculate forcing and check + src_term.setVal(0.0); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + abl_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + // Velocity at hub height is 8 at 0deg and target is 8 at 2.5deg + target_force[0] = (8.0 * std::cos(M_PI / 180.0 * 2.5) - init_vel[0]) / dt; + target_force[1] = (8.0 * std::sin(M_PI / 180.0 * 2.5) - init_vel[1]) / dt; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + // Delete target wind file + const char* fname = tvel_fname.c_str(); + { + std::ifstream f(fname); + if (f.good()) { + remove(fname); + } + // Check that file is removed + std::ifstream ff(fname); + EXPECT_FALSE(ff.good()); + } +} + +TEST_F(ABLSrcTimeTableTest, bodyforce) +{ + + constexpr amrex::Real tol = 1.0e-12; + + // Write bodyforce file + write_body_force_file(forces_fname); + + // Set up simulation parameters and mesh + populate_parameters(); + // Body Forcing + { + amrex::ParmParse pp("BodyForce"); + pp.add("type", (std::string) "uniform_timetable"); + pp.add("uniform_timetable_file", forces_fname); + } + // Velocity + { + amrex::ParmParse pp("incflo"); + pp.addarr("velocity", amrex::Vector{8.0, 0.0, 0.0}); + } + initialize_mesh(); + + // Set up PDEs and physics objects + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + + // Get icns source term, which will be tested + auto& src_term = pde_mgr.icns().fields().src_term; + src_term.setVal(0.0); + // Source term object for BodyForce + amr_wind::pde::icns::BodyForce body_forcing(sim()); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + body_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + + // Force is 0 initially + amrex::Vector target_force = {0.0, 0.0, 0.0}; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + amrex::Vector angles{0.0, 2.5, 5.0, 8.5}; + + // Loop through timesteps to check force outputs + for (int n = 0; n < nsteps - 1; ++n) { + // Advance time + sim().time().new_timestep(); + // Recalculate forcing and check + src_term.setVal(0.0); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + body_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + // Forces correspond to ABL Forcing from other test + const amrex::Vector init_vel{8.0, 0.0, 0.0}; + target_force[0] = + (8.0 * std::cos(M_PI / 180.0 * angles[n]) - init_vel[0]) / dt; + target_force[1] = + (8.0 * std::sin(M_PI / 180.0 * angles[n]) - init_vel[1]) / dt; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + } + + // Delete body force file + const char* fname = forces_fname.c_str(); + { + std::ifstream f(fname); + if (f.good()) { + remove(fname); + } + // Check that file is removed + std::ifstream ff(fname); + EXPECT_FALSE(ff.good()); + } +} + +TEST_F(ABLSrcTimeTableTest, geostrophic) +{ + constexpr amrex::Real tol = 1.0e-12; + // Default Coriolis parameters + const amrex::Real cf = + 2.0 * (2.0 * M_PI / 86164.091) * std::sin(M_PI / 180.0 * 90.0); + + // Write target wind file + write_target_velocity_file(tvel_fname); + + // Set up simulation parameters and mesh + populate_parameters(); + // Geostrophic Forcing + { + amrex::ParmParse pp("GeostrophicForcing"); + pp.add("geostrophic_wind_timetable", tvel_fname); + } + // Coriolis Forcing (for parameters) + { + amrex::ParmParse pp("CoriolisForcing"); + pp.add("latitude", 90.); + } + // Velocity + { + amrex::ParmParse pp("incflo"); + pp.addarr("velocity", amrex::Vector{8.0, 0.0, 0.0}); + } + initialize_mesh(); + + // Set up PDEs and physics objects + auto& pde_mgr = sim().pde_manager(); + pde_mgr.register_icns(); + sim().init_physics(); + + // Get icns source term, which will be tested + auto& src_term = pde_mgr.icns().fields().src_term; + src_term.setVal(0.0); + // Source term object for geostrophic forcing + amr_wind::pde::icns::GeostrophicForcing gstr_forcing(sim()); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + gstr_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + + // Initial velocity should be the same as target velocity + const amrex::Vector init_vel{8.0, 0.0, 0.0}; + amrex::Vector target_force = { + -cf * init_vel[1], cf * init_vel[0], 0.0}; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + // Advance time (twice to get to interpolation interval) + sim().time().new_timestep(); + sim().time().new_timestep(); + + // Recalculate forcing and check + src_term.setVal(0.0); + run_algorithm(src_term, [&](const int lev, const amrex::MFIter& mfi) { + const auto& bx = mfi.tilebox(); + const auto& src_arr = src_term(lev).array(mfi); + + gstr_forcing(lev, mfi, bx, amr_wind::FieldState::New, src_arr); + }); + // New target velocity is 8 at 1.25deg (nph velocity) + const amrex::Vector targ_vel{ + 8.0 * std::cos(M_PI / 180.0 * 1.25), + 8.0 * std::sin(M_PI / 180.0 * 1.25), 0.0}; + target_force[0] = -cf * targ_vel[1]; + target_force[1] = cf * targ_vel[0]; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + const auto min_val = utils::field_min(src_term, i); + const auto max_val = utils::field_max(src_term, i); + EXPECT_NEAR(min_val, target_force[i], tol); + EXPECT_NEAR(min_val, max_val, tol); + } + + // Delete target wind file + const char* fname = tvel_fname.c_str(); + { + std::ifstream f(fname); + if (f.good()) { + remove(fname); + } + // Check that file is removed + std::ifstream ff(fname); + EXPECT_FALSE(ff.good()); + } +} + +} // namespace amr_wind_tests