From deed504e3ae80fdda10661cf6999be8c95800ee6 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Thu, 10 Jul 2025 13:27:47 -0600 Subject: [PATCH 01/12] Initial code changes. --- amr-wind/wind_energy/ABLWallFunction.H | 16 ++- amr-wind/wind_energy/ABLWallFunction.cpp | 169 ++++++++++++++++++----- amr-wind/wind_energy/MOData.H | 4 + amr-wind/wind_energy/MOData.cpp | 6 + amr-wind/wind_energy/ShearStress.H | 5 + 5 files changed, 161 insertions(+), 39 deletions(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.H b/amr-wind/wind_energy/ABLWallFunction.H index 9d8e4b1e1c..d5ceb034a8 100644 --- a/amr-wind/wind_energy/ABLWallFunction.H +++ b/amr-wind/wind_energy/ABLWallFunction.H @@ -58,8 +58,20 @@ private: amrex::Vector m_surf_temp_time; amrex::Vector m_surf_temp_value; - bool m_tempflux{true}; - bool m_temp_table{false}; + //! Ability to read in a table of near-surface temperature verus time. + std::string m_near_surf_temp_timetable; + amrex::Vector m_near_surf_temp_time; + amrex::Vector m_near_surf_temp_value; + + bool m_temp_flux{true}; + bool m_surf_temp{false}; + bool m_nearsurf_temp{false}; + + bool m_surf_temp_use_table{false}; + bool m_surf_temp_use_rate{false}; + + bool m_nearsurf_temp_use_table{false}; + amrex::Real m_surf_temp_rate{0.0}; amrex::Real m_surf_temp_rate_tstart{0.0}; amrex::Real m_surf_temp_init{300.0}; diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 4438573751..28d9a5d734 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -60,19 +60,45 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) << std::endl; } + + + + // Check for type of surface heating and read in data. + // - Specified surface temperature flux mode: if (pp.contains("surface_temp_flux")) { + m_temp_flux = true; + m_surf_temp = false; + m_nearsurf_temp = false; + pp.query("surface_temp_flux", m_mo.surf_temp_flux); + amrex::Print() << "ABLWallFunction: Surface temperature flux mode is selected." << std::endl; - } else if (pp.contains("surface_temp_timetable")) { - pp.query("surface_temp_timetable", m_surf_temp_timetable); - m_tempflux = false; - m_temp_table = true; - amrex::Print() << "ABLWallFunction: Surface temperature time table " + + // - Specified near-surface temperature mode: + } else if (pp.contains("near_surface_temp_timetable")) { + m_temp_flux = false; + m_surf_temp = false; + m_nearsurf_temp = true; + + m_nearsurf_temp_use_table = true; + + if (pp.contains("near_surface_height")) { + pp.get("near_surface_height",m_mo.zNearSurf); + } else { + amrex::Abort( + "ABLWallFunction: near_surface_height must be specified when " + "specifying near-surface temperature."); + } + + pp.query("near_surface_temp_timetable", m_near_surf_temp_timetable); + + amrex::Print() << "ABLWallFunction: Near-surface temperature time table " "mode is selected." << std::endl; - if (!m_surf_temp_timetable.empty()) { + + if (!m_nearsurf_temp_timetable.empty()) { std::ifstream ifh(m_surf_temp_timetable, std::ios::in); if (!ifh.good()) { amrex::Abort( @@ -88,38 +114,83 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_surf_temp_value.push_back(data_value); } } - } else if (pp.contains("surface_temp_rate")) { - m_tempflux = false; - pp.get("surface_temp_rate", m_surf_temp_rate); - amrex::Print() - << "ABLWallFunction: Surface temperature rate mode is selected." - << std::endl; - if (pp.contains("surface_temp_init")) { - pp.get("surface_temp_init", m_surf_temp_init); - } else { + + // - Specified surface temperature mode: + } else if ((pp.contains("surface_temp_timetable")) || + (pp.contains("surface_temp_rate"))) { + m_temp_flux = false; + m_surf_temp = true; + m_nearsurf_temp = false; + + if (pp.contains("surface_temp_timetable") { + m_surf_temp_use_table = true; + m_surf_temp_use_rate = false; + + pp.query("surface_temp_timetable", m_surf_temp_timetable); + amrex::Print() << "ABLWallFunction: Surface temperature time table " + "mode is selected." + << std::endl; + if (!m_surf_temp_timetable.empty()) { + std::ifstream ifh(m_surf_temp_timetable, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find surface_temp_timetable file: " + + m_surf_temp_timetable); + } + amrex::Real data_time; + amrex::Real data_value; + ifh.ignore(std::numeric_limits::max(), '\n'); + while (ifh >> data_time) { + ifh >> data_value; + m_surf_temp_time.push_back(data_time); + m_surf_temp_value.push_back(data_value); + } + } + } else if (pp.contains("surface_temp_rate")) { + m_surf_temp_use_table = false; + m_surf_temp_use_rate = true; + + pp.get("surface_temp_rate", m_surf_temp_rate); + amrex::Print() - << "ABLWallFunction: Initial surface temperature not found for " - "ABL. Assuming to be equal to the reference temperature" + << "ABLWallFunction: Surface temperature rate mode is selected." << std::endl; - m_surf_temp_init = sim.transport_model().reference_temperature(); - } - if (pp.contains("surface_temp_rate_tstart")) { - pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); - } else { - amrex::Print() - << "ABLWallFunction: Surface temperature heating/cooling start " - "time (surface_temp_rate_tstart) not found for ABL. " - "Assuming zero." - << m_surf_temp_rate_tstart << std::endl; + + if (pp.contains("surface_temp_init")) { + pp.get("surface_temp_init", m_surf_temp_init); + } else { + amrex::Print() + << "ABLWallFunction: Initial surface temperature not found for " + "ABL. Assuming to be equal to the reference temperature" + << std::endl; + m_surf_temp_init = sim.transport_model().reference_temperature(); + } + + if (pp.contains("surface_temp_rate_tstart")) { + pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); + } else { + amrex::Print() + << "ABLWallFunction: Surface temperature heating/cooling start " + "time (surface_temp_rate_tstart) not found for ABL. " + "Assuming zero." + << m_surf_temp_rate_tstart << std::endl; + } } + // - If no surface heating mode is specified, default to no surface heating. } else { - amrex::Print() << "ABLWallFunction: Neither surface_temp_flux nor " - "surface_temp_rate specified for ABL physics. " - "Assuming Neutral Stratification" + amrex::Print() << "ABLWallFunction: None of 'surface_temp_flux,' " + "'surface_temp_rate', 'surface_temp_timetable,' " + "or 'near_surface_temp_timetable' are specified " + "for ABL physics. Assuming no heat flux at surface." << std::endl; } + + + + // Check to see if the simulation is non-periodic. If so, recommend + // using the wall model locally instead of planar averaging. if (pp.contains("inflow_outflow_mode")) { pp.query("inflow_outflow_mode", m_inflow_outflow); if (m_inflow_outflow) { @@ -127,14 +198,30 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) pp.query("wf_vmag", m_wf_vmag); pp.query("wf_theta", m_wf_theta); amrex::Print() << "ABLWallFunction: Inflow/Outflow mode is turned " - "on. Please make sure wall shear stress type is " - "set to local." + "on. It is recommended that the wall shear stress " + "type is set to local IF the flow is heterogeneous." << std::endl; } } - m_mo.alg_type = m_tempflux ? MOData::ThetaCalcType::HEAT_FLUX - : MOData::ThetaCalcType::SURFACE_TEMPERATURE; + + + + // Set the Monin-Obukhov algorith type. + if (m_temp_flux) { + m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; + } else if (m_surf_temp) { + m_mo.alg_type = MOData::ThetaCalcType::SURFACE_TEMPERATURE; + } else if (m_nearsurf_temp) { + m_mo.alg_type = MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE; + } else { + m_mo.alg_type = MODATA::ThetaCalcType::HEAT_FLUX; + } + + + + + // Set gravity strength in Monin-Obukhov. m_mo.gravity = utils::vec_mag(m_gravity.data()); } @@ -151,20 +238,28 @@ void ABLWallFunction::update_umean( { const auto& time = m_sim.time(); - if (!m_tempflux) { - if (!m_temp_table) { + // If not using a constant heat flux... + // - if specifying surface temperature + if (m_surf_temp) { + if (m_surf_temp_use_rate) { m_mo.surf_temp = m_surf_temp_init + m_surf_temp_rate * amrex::max( time.current_time() - m_surf_temp_rate_tstart, 0.0) / 3600.0; - } else { + } else if (m_surf_temp_use_table) { m_mo.surf_temp = amr_wind::interp::linear( m_surf_temp_time, m_surf_temp_value, time.current_time()); } amrex::Print() << "Current surface temperature: " << m_mo.surf_temp << std::endl; + // - if specifying near-surface temperature + } else if (m_nearsurf_temp) { + if (m_nearsurf_temp_use_table) { + m_mo.near_surf_temp = amr_wind::interp::linear( + m_nearsurf_temp_time, m_nearsurf_temp_value, time.current_time()); + } } if (m_inflow_outflow) { diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H index 101a1c627a..0aab220029 100644 --- a/amr-wind/wind_energy/MOData.H +++ b/amr-wind/wind_energy/MOData.H @@ -22,9 +22,12 @@ struct MOData enum class ThetaCalcType { HEAT_FLUX = 0, ///< Heat-flux specified SURFACE_TEMPERATURE ///< Surface temperature specified + NEAR_SURFACE_TEMPERATURE = 0, ///< Temperature specified above but + ///< near the surface }; amrex::Real zref{0.0}; ///< Reference height (m) + amrex::Real zNearSurf{2.0}; ///< Near-surface temperature driving height (m) amrex::Real z0{0.1}; ///< Aerodynamic roughness height (m) amrex::Real z0t{z0}; ///< Thermal roughness height (m) amrex::Real utau; ///< Friction velocity (m/s) @@ -40,6 +43,7 @@ struct MOData amrex::Real surf_temp_flux{0.0}; ///< Heat flux amrex::Real surf_temp; ///< Instantaneous surface temperature + amrex::Real near_surf_temp; ///< Instantaneous near-surface temperature (at zNearSurf). amrex::Real gamma_m{5.0}; amrex::Real gamma_h{5.0}; diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index 75241cbc1f..f9deb748e8 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -70,6 +70,12 @@ void MOData::update_fluxes(int max_iters) surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa / (std::log(zref / z0t) - psi_h); break; + + case ThetaCalcType::NEAR_SURFACE_TEMPERATURE: + amrex::Print() + << "To be implemented..." + << std::endl; + break; } if (std::abs(surf_temp_flux) > eps) { diff --git a/amr-wind/wind_energy/ShearStress.H b/amr-wind/wind_energy/ShearStress.H index 34d20f97cb..e921a6323b 100644 --- a/amr-wind/wind_energy/ShearStress.H +++ b/amr-wind/wind_energy/ShearStress.H @@ -216,6 +216,11 @@ struct ShearStressDonelan case amr_wind::MOData::ThetaCalcType::SURFACE_TEMPERATURE: flux = 0.0012 * wspd_mean * (theta_surface - theta_mean); break; + case amr_wind::MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE: + amrex::Print() + << "To be implemented..." + << std::endl; + break; } return flux; }; From 3a69d8028c7e91605d3e978288127169e5428e2e Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Thu, 10 Jul 2025 14:28:23 -0600 Subject: [PATCH 02/12] Added functionality to specify a near surface temperature as a way to drive surface heat flux. --- amr-wind/wind_energy/ABLWallFunction.H | 6 ++--- amr-wind/wind_energy/ABLWallFunction.cpp | 17 +++++++------ amr-wind/wind_energy/MOData.H | 9 ++++--- amr-wind/wind_energy/MOData.cpp | 32 +++++++++++++++--------- 4 files changed, 37 insertions(+), 27 deletions(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.H b/amr-wind/wind_energy/ABLWallFunction.H index d5ceb034a8..c897237177 100644 --- a/amr-wind/wind_energy/ABLWallFunction.H +++ b/amr-wind/wind_energy/ABLWallFunction.H @@ -59,9 +59,9 @@ private: amrex::Vector m_surf_temp_value; //! Ability to read in a table of near-surface temperature verus time. - std::string m_near_surf_temp_timetable; - amrex::Vector m_near_surf_temp_time; - amrex::Vector m_near_surf_temp_value; + std::string m_nearsurf_temp_timetable; + amrex::Vector m_nearsurf_temp_time; + amrex::Vector m_nearsurf_temp_value; bool m_temp_flux{true}; bool m_surf_temp{false}; diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 28d9a5d734..c3c8a97994 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -26,6 +26,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) amrex::ParmParse pp("ABL"); pp.query("kappa", m_mo.kappa); + pp.query("mo_alpha_h", m_mo.alpha_h); pp.query("mo_gamma_m", m_mo.gamma_m); pp.query("mo_gamma_h", m_mo.gamma_h); pp.query("mo_beta_m", m_mo.beta_m); @@ -92,26 +93,26 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) "specifying near-surface temperature."); } - pp.query("near_surface_temp_timetable", m_near_surf_temp_timetable); + pp.query("near_surface_temp_timetable", m_nearsurf_temp_timetable); amrex::Print() << "ABLWallFunction: Near-surface temperature time table " "mode is selected." << std::endl; if (!m_nearsurf_temp_timetable.empty()) { - std::ifstream ifh(m_surf_temp_timetable, std::ios::in); + std::ifstream ifh(m_nearsurf_temp_timetable, std::ios::in); if (!ifh.good()) { amrex::Abort( - "Cannot find surface_temp_timetable file: " + - m_surf_temp_timetable); + "Cannot find near_surface_temp_timetable file: " + + m_nearsurf_temp_timetable); } amrex::Real data_time; amrex::Real data_value; ifh.ignore(std::numeric_limits::max(), '\n'); while (ifh >> data_time) { ifh >> data_value; - m_surf_temp_time.push_back(data_time); - m_surf_temp_value.push_back(data_value); + m_nearsurf_temp_time.push_back(data_time); + m_nearsurf_temp_value.push_back(data_value); } } @@ -122,7 +123,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_surf_temp = true; m_nearsurf_temp = false; - if (pp.contains("surface_temp_timetable") { + if (pp.contains("surface_temp_timetable")) { m_surf_temp_use_table = true; m_surf_temp_use_rate = false; @@ -215,7 +216,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) } else if (m_nearsurf_temp) { m_mo.alg_type = MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE; } else { - m_mo.alg_type = MODATA::ThetaCalcType::HEAT_FLUX; + m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; } diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H index 0aab220029..1526360999 100644 --- a/amr-wind/wind_energy/MOData.H +++ b/amr-wind/wind_energy/MOData.H @@ -20,10 +20,10 @@ namespace amr_wind { struct MOData { enum class ThetaCalcType { - HEAT_FLUX = 0, ///< Heat-flux specified - SURFACE_TEMPERATURE ///< Surface temperature specified - NEAR_SURFACE_TEMPERATURE = 0, ///< Temperature specified above but - ///< near the surface + HEAT_FLUX = 0, ///< Heat-flux specified + SURFACE_TEMPERATURE = 1, ///< Surface temperature specified + NEAR_SURFACE_TEMPERATURE = 2, ///< Temperature specified above but + // near the surface }; amrex::Real zref{0.0}; ///< Reference height (m) @@ -45,6 +45,7 @@ struct MOData amrex::Real surf_temp; ///< Instantaneous surface temperature amrex::Real near_surf_temp; ///< Instantaneous near-surface temperature (at zNearSurf). + amrex::Real alpha_h{1.0}; amrex::Real gamma_m{5.0}; amrex::Real gamma_h{5.0}; amrex::Real beta_m{16.0}; diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index f9deb748e8..bf1024b6bc 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -52,8 +52,8 @@ void MOData::update_fluxes(int max_iters) amrex::Real utau_iter = 0.0; // Initialize variables - amrex::Real psi_m = 0.0; - amrex::Real psi_h = 0.0; + amrex::Real psi_m_zref = 0.0; + amrex::Real psi_h_zref = 0.0; utau = kappa * vmag_mean / (std::log(zref / z0)); int iter = 0; @@ -61,20 +61,28 @@ void MOData::update_fluxes(int max_iters) utau_iter = utau; switch (alg_type) { case ThetaCalcType::HEAT_FLUX: - surf_temp = surf_temp_flux * (std::log(zref / z0t) - psi_h) / + surf_temp = alpha_h * surf_temp_flux * (std::log(zref / z0t) - psi_h_zref) / (utau * kappa) + theta_mean; break; case ThetaCalcType::SURFACE_TEMPERATURE: surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa / - (std::log(zref / z0t) - psi_h); + (alpha_h * (std::log(zref / z0t) - psi_h_zref)); break; case ThetaCalcType::NEAR_SURFACE_TEMPERATURE: - amrex::Print() - << "To be implemented..." - << std::endl; + amrex::Real psi_h_zNearSurf = calc_psi_h(zNearSurf / obukhov_len); + amrex::Real a11 = 1.0; + amrex::Real a12 = -(alpha_h / (kappa * utau)) * (std::log(zref / z0t) - psi_h_zref); + amrex::Real a21 = 1.0; + amrex::Real a22 = -(alpha_h / (kappa * utau)) * (std::log(zNearSurf / z0t) - psi_h_zNearSurf); + + amrex::Real coeff = 1.0 / (a11*a21 - a12*a21); + + surf_temp = coeff * (a22*theta_mean - a12*near_surf_temp); + surf_temp_flux = coeff * (-a21*theta_mean + a11*near_surf_temp); + break; } @@ -88,9 +96,9 @@ void MOData::update_fluxes(int max_iters) obukhov_len = std::numeric_limits::max(); zeta = 0.0; } - psi_m = calc_psi_m(zeta); - psi_h = calc_psi_h(zeta); - utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m); + psi_m_zref = calc_psi_m(zeta); + psi_h_zref = calc_psi_h(zeta); + utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m_zref); ++iter; } while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters); @@ -98,8 +106,8 @@ void MOData::update_fluxes(int max_iters) amrex::Print() << "MOData::update_fluxes: Convergence criteria not met after " << max_iters << " iterations\nObuhov length = " << obukhov_len - << " zeta = " << zeta << "\npsi_m = " << psi_m - << " psi_h = " << psi_h << "\nutau = " << utau + << " zeta = (z_ref/L) = " << zeta << "\npsi_m(z_ref/L) = " << psi_m_zref + << " psi_h(z_ref/L) = " << psi_h_zref << "\nutau = " << utau << " Tsurf = " << surf_temp << " q = " << surf_temp_flux << std::endl; } From 85f40ac58d86c519a4d60b2793f228a2b0fd912e Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Mon, 28 Jul 2025 11:30:19 -0600 Subject: [PATCH 03/12] Small changes for writing to screen. --- amr-wind/wind_energy/ABLWallFunction.cpp | 6 ++++++ amr-wind/wind_energy/MOData.cpp | 9 ++++++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index c3c8a97994..a9744b3e39 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -34,6 +34,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) const char* z0_same = "surface_roughness_z0"; const char* z0_aero = "aerodynamic_roughness_length"; const char* z0_therm = "thermal_roughness_length"; + pp.query(z0_same, m_mo.z0); if (!pp.contains(z0_same)) { pp.query(z0_aero, m_mo.z0); @@ -114,6 +115,11 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_nearsurf_temp_time.push_back(data_time); m_nearsurf_temp_value.push_back(data_value); } + + for (int i =0; i < m_nearsurf_temp_time.size(); i++) + { + amrex::Print() << i << " " << m_nearsurf_temp_time[i] << " " << m_nearsurf_temp_value[i] << std::endl; + } } // - Specified surface temperature mode: diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index bf1024b6bc..64a94b6990 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -56,6 +56,8 @@ void MOData::update_fluxes(int max_iters) amrex::Real psi_h_zref = 0.0; utau = kappa * vmag_mean / (std::log(zref / z0)); + amrex::Print() << "Using z0 = " << z0 << ", z0t = " << z0t << std::endl; + int iter = 0; do { utau_iter = utau; @@ -78,7 +80,7 @@ void MOData::update_fluxes(int max_iters) amrex::Real a21 = 1.0; amrex::Real a22 = -(alpha_h / (kappa * utau)) * (std::log(zNearSurf / z0t) - psi_h_zNearSurf); - amrex::Real coeff = 1.0 / (a11*a21 - a12*a21); + amrex::Real coeff = 1.0 / (a11*a22 - a12*a21); surf_temp = coeff * (a22*theta_mean - a12*near_surf_temp); surf_temp_flux = coeff * (-a21*theta_mean + a11*near_surf_temp); @@ -100,6 +102,11 @@ void MOData::update_fluxes(int max_iters) psi_h_zref = calc_psi_h(zeta); utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m_zref); ++iter; + + amrex::Print() << "iteration: " << iter << ", T(" << zNearSurf << " m) = " << near_surf_temp + << " K, T(" << zref << " m) = " << theta_mean << "K , L = " << obukhov_len + << " m, uStar = " << utau << "m/s, Qs = " << surf_temp_flux + << " K-m/s, theta_0 = " << surf_temp << " K" << std::endl; } while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters); if (iter >= max_iters) { From c4e0aba0d93a2f94512907a1835503080e696569 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Fri, 1 Aug 2025 10:57:39 -0600 Subject: [PATCH 04/12] Small addition to code commenting. --- amr-wind/wind_energy/ABLWallFunction.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index a9744b3e39..64cc7df331 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -153,6 +153,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_surf_temp_value.push_back(data_value); } } + // - Specified surface temperature rate mode: } else if (pp.contains("surface_temp_rate")) { m_surf_temp_use_table = false; m_surf_temp_use_rate = true; From f2f3cc510f609a5a89fb6e715d93e57bd760b7bc Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Mon, 4 Aug 2025 12:49:54 -0600 Subject: [PATCH 05/12] Updates to be able to apply time-varying surface heat flux. --- amr-wind/wind_energy/ABLWallFunction.H | 9 +++- amr-wind/wind_energy/ABLWallFunction.cpp | 58 ++++++++++++++++++++---- 2 files changed, 57 insertions(+), 10 deletions(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.H b/amr-wind/wind_energy/ABLWallFunction.H index c897237177..79c5a0b163 100644 --- a/amr-wind/wind_energy/ABLWallFunction.H +++ b/amr-wind/wind_energy/ABLWallFunction.H @@ -53,6 +53,11 @@ private: amrex::Vector m_gravity{0.0, 0.0, -9.81}; + //! Ability to read in a table of surface temperature flux versus time. + std::string m_surf_temp_flux_timetable; + amrex::Vector m_surf_temp_flux_time; + amrex::Vector m_surf_temp_flux_value; + //! Ability to read in a table of surface temperature versus time. std::string m_surf_temp_timetable; amrex::Vector m_surf_temp_time; @@ -63,10 +68,12 @@ private: amrex::Vector m_nearsurf_temp_time; amrex::Vector m_nearsurf_temp_value; - bool m_temp_flux{true}; + bool m_surf_temp_flux{true}; bool m_surf_temp{false}; bool m_nearsurf_temp{false}; + bool m_surf_temp_flux_use_table{false}; + bool m_surf_temp_use_table{false}; bool m_surf_temp_use_rate{false}; diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 64cc7df331..769ba515c9 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -67,20 +67,54 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) // Check for type of surface heating and read in data. // - Specified surface temperature flux mode: - if (pp.contains("surface_temp_flux")) { - m_temp_flux = true; + if ((pp.contains("surface_temp_flux")) || + (pp.contains("surface_temp_flux_timetable"))) { + m_surf_temp_flux = true; m_surf_temp = false; m_nearsurf_temp = false; - pp.query("surface_temp_flux", m_mo.surf_temp_flux); + // - Fixed, time-invariant surface temperature flux mode: + if (pp.contains("surface_temp_flux")) { + m_surf_temp_flux_use_table = false; - amrex::Print() - << "ABLWallFunction: Surface temperature flux mode is selected." - << std::endl; + pp.query("surface_temp_flux", m_mo.surf_temp_flux); + + amrex::Print() + << "ABLWallFunction: Fixed, time-invariant surface temperature " + "flux mode is selected." + << std::endl; + } + // - surface temperature flux from time-table mode: + else if (pp.contains("surface_temp_flux_timetable")) { + m_surf_temp_flux_use_table = true; + + pp.query("surface_temp_flux_timetable", m_surf_temp_flux_timetable); + + amrex::Print() + << "ABLWallFunction: Surface temperature flux time table " + "mode is selected." + << std::endl; + if (!m_surf_temp_flux_timetable.empty()) { + std::ifstream ifh(m_surf_temp_flux_timetable, std::ios::in); + if (!ifh.good()) { + amrex::Abort( + "Cannot find surface_temp_flux_timetable file: " + + m_surf_temp_flux_timetable); + } + amrex::Real data_time; + amrex::Real data_value; + ifh.ignore(std::numeric_limits::max(), '\n'); + while (ifh >> data_time) { + ifh >> data_value; + m_surf_temp_flux_time.push_back(data_time); + m_surf_temp_flux_value.push_back(data_value); + } + } + } // - Specified near-surface temperature mode: } else if (pp.contains("near_surface_temp_timetable")) { - m_temp_flux = false; + m_surf_temp_flux = false; m_surf_temp = false; m_nearsurf_temp = true; @@ -125,7 +159,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) // - Specified surface temperature mode: } else if ((pp.contains("surface_temp_timetable")) || (pp.contains("surface_temp_rate"))) { - m_temp_flux = false; + m_surf_temp_flux = false; m_surf_temp = true; m_nearsurf_temp = false; @@ -216,7 +250,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) // Set the Monin-Obukhov algorith type. - if (m_temp_flux) { + if (m_surf_temp_flux) { m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; } else if (m_surf_temp) { m_mo.alg_type = MOData::ThetaCalcType::SURFACE_TEMPERATURE; @@ -247,6 +281,12 @@ void ABLWallFunction::update_umean( const auto& time = m_sim.time(); // If not using a constant heat flux... + if (m_surf_temp_flux && m_surf_temp_flux_use_table) { + m_mo.surf_temp_flux = amr_wind::interp::linear( + m_surf_temp_flux_time, m_surf_temp_flux_value, time.current_time()); + amrex::Print() << "Current surface temperature flux: " << m_mo.surf_temp_flux + << std::endl; + } // - if specifying surface temperature if (m_surf_temp) { if (m_surf_temp_use_rate) { From 7b9260b9389a49bf5f6d52552c29df2247350807 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Tue, 19 Aug 2025 12:02:55 -0600 Subject: [PATCH 06/12] Updates to ABLForcing, ABLStats, and ABL to implement free-atmosphere wind vector damper for use with ABLForcing. --- .../icns/source_terms/ABLForcing.H | 53 ++++++++ .../icns/source_terms/ABLForcing.cpp | 120 +++++++++++++++++- amr-wind/wind_energy/ABL.cpp | 6 + amr-wind/wind_energy/ABLStats.H | 5 + amr-wind/wind_energy/ABLStatsBase.H | 2 + 5 files changed, 182 insertions(+), 4 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H index 53a9f0fbe6..54f7cdb363 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H @@ -5,6 +5,7 @@ #include "amr-wind/core/SimTime.H" #include "amr-wind/utilities/trig_ops.H" #include "amr-wind/utilities/linear_interpolation.H" +#include "amr-wind/utilities/FieldPlaneAveraging.H" namespace amr_wind::pde::icns { @@ -29,12 +30,22 @@ public: const FieldState fstate, const amrex::Array4& src_term) const override; + void mean_velocity_init(const VelPlaneAveraging& /*vavg*/); + + void mean_velocity_update(const VelPlaneAveraging& /*vavg*/); + inline void set_target_velocities(amrex::Real ux, amrex::Real uy) { m_target_vel[0] = ux; m_target_vel[1] = uy; } + inline void set_boundary_layer_height(amrex::Real zi) + { + m_zi = zi; + amrex::Print() << "z_i = " << zi << std::endl; + } + inline void set_mean_velocities(amrex::Real ux, amrex::Real uy) { m_mean_vel[0] = ux; @@ -62,6 +73,15 @@ public: 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; + amrex::Print() << "U(zRef), V(zRef) = " << m_mean_vel[0] << " " << + m_mean_vel[1] << std::endl; + amrex::Print() << "dt = " << dt << std::endl; + amrex::Print() << "f_c = " << m_coriolis_factor << std::endl; + amrex::Print() << "Sx, Sy = " << m_abl_forcing[0] << " " << + m_abl_forcing[1] << std::endl; + amrex::Print() << "Ug, Vg = " << m_abl_forcing[1]/m_coriolis_factor << " " << + -m_abl_forcing[0]/m_coriolis_factor << std::endl; + if (m_write_force_timetable && amrex::ParallelDescriptor::IOProcessor() && (t_step % m_force_outfreq == 0) && @@ -83,9 +103,36 @@ private: const SimTime& m_time; const amrex::AmrCore& m_mesh; + amrex::Gpu::DeviceVector m_vel_ht; + amrex::Gpu::DeviceVector m_vel_vals; + + //! Axis over which averages are computed. + int m_axis{2}; + //! Activated when water is present in domain bool m_use_phase_ramp{false}; + //! Activate free-atmosphere damping + bool m_fa_damping{false}; + + //! Detect free-atmosphere base height (set to zi from ABLStats). + bool m_fa_detect_height{false}; + + //! Free-atmosphere base height. + amrex::Real m_fa_height{750.0}; + + //! Start time for free-atmosphere damping. + amrex::Real m_fa_time_start{0.0}; + + //! End time for free-atmosphere damping. + amrex::Real m_fa_time_end{1.0E9}; + + //! Time scale for free-atmosphere damping. + amrex::Real m_fa_tau{100.0}; + + //! Coriolis factor used for free-atmosphere damping. + amrex::Real m_coriolis_factor; + //! Number of cells in band to prevent forcing near liquid int m_n_band{2}; @@ -97,10 +144,13 @@ private: //! 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}; @@ -130,6 +180,9 @@ private: //! Local storage of interface location amrex::Real m_water_level; + //! boundary layer height. + amrex::Real m_zi; + //! VOF field, to avoid forcing on liquid above force-off height const Field* m_vof; }; diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp index 79ad1fbbd6..ab01fd4543 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp @@ -57,6 +57,33 @@ ABLForcing::ABLForcing(const CFDSim& sim) } } + // If free-atmosphere damping is used, read these inputs. + pp_abl.query("free_atmosphere_damping", m_fa_damping); + if (m_fa_damping) { + pp_abl.query("detect_free_atmosphere_height", m_fa_detect_height); + if (! m_fa_detect_height) { + pp_abl.get("free_atmosphere_height", m_fa_height); + } + pp_abl.query("free_atmosphere_damping_start_time", m_fa_time_start); + pp_abl.query("free_atmosphere_damping_end_time", m_fa_time_end); + pp_abl.query("free_atmosphere_damping_time_scale", m_fa_tau); + + amrex::ParmParse pp_coriolis("CoriolisForcing"); + amrex::Real rot_time_period = 86164.091; + pp_coriolis.query("rotational_time_period", rot_time_period); + + amrex::Real omega = utils::two_pi() / rot_time_period; + + amrex::Real latitude = 90.0; + pp_coriolis.query("latitude", latitude); + latitude = utils::radians(latitude); + amrex::Real sinphi = std::sin(latitude); + + m_coriolis_factor = 2.0 * omega * sinphi; + } + amrex::Print() << "free atmosphere damping = " << m_fa_damping << std::endl; + + for (int i = 0; i < AMREX_SPACEDIM; ++i) { m_mean_vel[i] = m_target_vel[i]; } @@ -80,6 +107,8 @@ ABLForcing::ABLForcing(const CFDSim& sim) // Point to something, will not be used m_vof = &sim.repo().get_field("velocity"); } + + mean_velocity_init(abl.abl_statistics().vel_profile_coarse()); } ABLForcing::~ABLForcing() = default; @@ -95,6 +124,20 @@ void ABLForcing::operator()( const amrex::Real dudt = m_abl_forcing[0]; const amrex::Real dvdt = m_abl_forcing[1]; + const bool fa_damping = m_fa_damping; + const bool fa_detect_height = m_fa_detect_height; + const amrex::Real fa_height = fa_detect_height ? m_zi : m_fa_height; + amrex::AllPrint() << "fa_height = " << fa_height << " fa_detect_height = " << fa_detect_height << " m_zi = " << m_zi << std::endl; + const amrex::Real fa_time_start = m_fa_time_start; + const amrex::Real fa_time_end = m_fa_time_end; + const amrex::Real fa_tau = m_fa_tau; + const amrex::Real fa_u = dvdt / m_coriolis_factor; + const amrex::Real fa_v = -dudt / m_coriolis_factor; + + 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 bool ph_ramp = m_use_phase_ramp; const int n_band = m_n_band; const amrex::Real wlev = m_water_level; @@ -103,11 +146,31 @@ void ABLForcing::operator()( const auto& problo = m_mesh.Geom(lev).ProbLoArray(); const auto& dx = m_mesh.Geom(lev).CellSizeArray(); + const int idir = m_axis; + const amrex::Real* heights = m_vel_ht.data(); + const amrex::Real* heights_end = m_vel_ht.end(); + const amrex::Real* vals = m_vel_vals.data(); + const auto& vof = (*m_vof)(lev).const_array(mfi); + + + //amrex::Print() << "fa_height = " << fa_height << std::endl; + //amrex::Print() << "Sx, Sy = " << dudt << ", " << dvdt << std::endl; + //amrex::Print() << "Ug, Vg = " << fa_u << ", " << fa_v << std::endl; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::IntVect iv(i, j, k); + const amrex::Real z = problo[idir] + (iv[idir] + 0.5) * dx[idir]; + + const amrex::Real umean = + amr_wind::interp::linear(heights, heights_end, vals, z, 3, 0); + const amrex::Real vmean = + amr_wind::interp::linear(heights, heights_end, vals, z, 3, 1); + + // This part deals with air and water phase and only masks ABL forcing + // on the water phase. amrex::Real fac = 1.0; if (ph_ramp) { - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; if (z - wlev < wrht0 + wrht1) { if (z - wlev < wrht0) { // Apply no forcing within first interval @@ -127,11 +190,60 @@ void ABLForcing::operator()( fac = 0.0; } } - src_term(i, j, k, 0) += fac * dudt; - src_term(i, j, k, 1) += fac * dvdt; - // No forcing in z-direction + // This part applies the free atmosphere damping, which nudges + // the solution in the free atmosphere toward geostrophic as + // computed by the wind speed controller. (i.e., the wind speed + // controller computes a forcing term that is really a driving + // pressure gradient, so there is a corresponding geostrophic + // wind. Nudge the free atmosphere toward that geostrophic wind.) + amrex::Real src_fa_damping_x = 0.0; + amrex::Real src_fa_damping_y = 0.0; + if (fa_damping && + ((nph_time >= fa_time_start) && + (nph_time <= fa_time_end)) && + (z >= fa_height)) { + src_fa_damping_x = (1.0/fa_tau) * (fa_u - umean); + src_fa_damping_y = (1.0/fa_tau) * (fa_v - vmean); + } + + amrex::AllPrint() << "z = " << z << ", fa_u = " << fa_u << ", umean = " << umean << ", src_fa_damping_x = " << src_fa_damping_x << ", dpdx = " << fac*dudt << std::endl; + amrex::AllPrint() << "z = " << z << ", fa_v = " << fa_v << ", vmean = " << vmean << ", src_fa_damping_y = " << src_fa_damping_y << ", dpdy = " << fac*dvdt << std::endl; + + // Sum up the source term + src_term(i, j, k, 0) += (fac * dudt) + src_fa_damping_x; + src_term(i, j, k, 1) += (fac * dvdt) + src_fa_damping_y; + src_term(i, j, k, 2) += 0.0; }); } +void ABLForcing::mean_velocity_init(const VelPlaneAveraging& vavg) +{ + m_axis = vavg.axis(); + + // The implementation depends the assumption that the ABL statistics class + // computes statistics at the cell-centeres only on level 0. If this + // assumption changes in future, the implementation will break... so put in + // a check here to catch this. + AMREX_ALWAYS_ASSERT( + m_mesh.Geom(0).Domain().length(m_axis) == + static_cast(vavg.line_centroids().size())); + + m_vel_ht.resize(vavg.line_centroids().size()); + m_vel_vals.resize(vavg.line_average().size()); + + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, vavg.line_centroids().begin(), + vavg.line_centroids().end(), m_vel_ht.begin()); + + mean_velocity_update(vavg); +} + +void ABLForcing::mean_velocity_update(const VelPlaneAveraging& vavg) +{ + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, vavg.line_average().begin(), + vavg.line_average().end(), m_vel_vals.begin()); +} + } // namespace amr_wind::pde::icns diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index 8ca419de57..eaa3dde427 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -173,6 +173,8 @@ void ABL::post_regrid_actions() { m_abl_anelastic->post_regrid_actions(); } */ void ABL::pre_advance_work() { + const auto& zi = m_stats->zi(); + amrex::Print() << "in ABL.H: zi = " << zi << std::endl; const auto& vel_pa = m_stats->vel_profile(); m_abl_wall_func.update_umean( m_stats->vel_profile(), m_stats->theta_profile_fine()); @@ -201,6 +203,10 @@ void ABL::pre_advance_work() #endif m_abl_forcing->set_mean_velocities(vx, vy); + m_abl_forcing->set_boundary_layer_height(zi); + + m_abl_forcing->mean_velocity_update( + m_stats->vel_profile_coarse()); } if (m_abl_mean_bous != nullptr) { diff --git a/amr-wind/wind_energy/ABLStats.H b/amr-wind/wind_energy/ABLStats.H index 1ba8269e3e..cdc14ffb7e 100644 --- a/amr-wind/wind_energy/ABLStats.H +++ b/amr-wind/wind_energy/ABLStats.H @@ -58,6 +58,11 @@ public: //! Compute height of capping inversion void compute_zi(); + amrex::Real zi() const override + { + return m_zi; + }; + //! Return vel plane averaging instance const VelPlaneAveragingFine& vel_profile() const override { diff --git a/amr-wind/wind_energy/ABLStatsBase.H b/amr-wind/wind_energy/ABLStatsBase.H index 7a10f5ab48..2fb08cd3ab 100644 --- a/amr-wind/wind_energy/ABLStatsBase.H +++ b/amr-wind/wind_energy/ABLStatsBase.H @@ -32,6 +32,8 @@ public: //! Flag indicating ABL simulation mode virtual ABLStatsMode abl_mode() const = 0; + virtual amrex::Real zi() const = 0; + //! Interpolating object for vertical velocity profile virtual const VelPlaneAveraging& vel_profile_coarse() const = 0; virtual const VelPlaneAveragingFine& vel_profile() const = 0; From d90e450b5ba46ab29cb5a819cc0e3c49cb6ba079 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Thu, 21 Aug 2025 15:27:18 -0600 Subject: [PATCH 07/12] Adding adequate comments to the function that computes z_i in ABLStats. --- amr-wind/wind_energy/ABLStats.cpp | 32 ++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index e855cdd78b..604611d1df 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -230,13 +230,29 @@ void ABLStats::post_advance_work() void ABLStats::compute_zi() { - + // This finds the location of the capping inversion (z_i) using the + // method of Sullivan et al. in "Structure of the entrainment zone + // capping the convective atmospheric boundary layer," JAS, Vol. 55, + // 1998. + // + // In this method, for every x,y location in the computational domain, + // the maximum d(\theta)/dz is found. You then have an x,y array + // of x,y local z_i. z_i is then averaged over the x,y array to give + // . + + // Compute the potential temperature gradient. auto gradT = (this->m_sim.repo()) .create_scratch_field(3, m_temperature.num_grow()[0]); fvm::gradient(*gradT, m_temperature); // Only compute zi using coarsest level const int lev = 0; + + // Using the AMReX ReduceToPlane function, find the maximum + // d(\theta)/dz over the mesh along lines in the normal direction, + // which typically would be z (if z is up). We end up with this + // 2D "planar" array of max d(\theta)/dz over the horizontal + // extent of the domain. const int dir = m_normal_dir; const auto& geom = (this->m_sim.repo()).mesh().Geom(lev); auto const& domain_box = geom.Domain(); @@ -260,11 +276,21 @@ void ABLStats::compute_zi() auto& pinned_tg_fab = device_tg_fab; #endif + // Because of the parallel decomposition, there may be multiple + // of these 2D "planar" arrays local to different processes handling + // boxes of mesh at different heights. So, we need to parallel + // reduce these by picking the max value at each x,y location to get + // the final resultant planar array of max temperature gradient values. + // This final array is collected only to the IO processor. amrex::ParallelReduce::Max( pinned_tg_fab.dataPtr(), static_cast(pinned_tg_fab.size()), amrex::ParallelDescriptor::IOProcessorNumber(), amrex::ParallelDescriptor::Communicator()); + // To get the planar average over the non-normal directions (typically + // x and y), we sum up all the max temperature gradient values. Again, + // this only happens on the IO processor. Once we have the sum, divide + // by the number of coarse grid points in x and y. if (amrex::ParallelDescriptor::IOProcessor()) { const auto dnval = m_dn; auto* p = pinned_tg_fab.dataPtr(); @@ -276,6 +302,10 @@ void ABLStats::compute_zi() 0.0); m_zi /= static_cast(pinned_tg_fab.size()); } + + // It is necessary for other processors to know , so broadcast + // it back out to all processors. + amrex::ParallelDescriptor::Bcast(&m_zi, 1, amrex::ParallelDescriptor::IOProcessorNumber()); } void ABLStats::process_output() From 0d31a60f9da9db0dbe857174823d9ab0ac71fbd4 Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Tue, 26 Aug 2025 10:12:59 -0600 Subject: [PATCH 08/12] Finalizing changes for free-atmosphere damping, and cleaned up some stuff about calculating surface heat flux from temperature. --- .../icns/source_terms/ABLForcing.H | 14 +++------- .../icns/source_terms/ABLForcing.cpp | 13 +++------ amr-wind/wind_energy/ABL.cpp | 1 - amr-wind/wind_energy/ABLWallFunction.cpp | 27 ++++++------------- amr-wind/wind_energy/MOData.cpp | 6 +++++ 5 files changed, 20 insertions(+), 41 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H index 54f7cdb363..287e46f6d1 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H @@ -43,7 +43,6 @@ public: inline void set_boundary_layer_height(amrex::Real zi) { m_zi = zi; - amrex::Print() << "z_i = " << zi << std::endl; } inline void set_mean_velocities(amrex::Real ux, amrex::Real uy) @@ -73,15 +72,6 @@ public: 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; - amrex::Print() << "U(zRef), V(zRef) = " << m_mean_vel[0] << " " << - m_mean_vel[1] << std::endl; - amrex::Print() << "dt = " << dt << std::endl; - amrex::Print() << "f_c = " << m_coriolis_factor << std::endl; - amrex::Print() << "Sx, Sy = " << m_abl_forcing[0] << " " << - m_abl_forcing[1] << std::endl; - amrex::Print() << "Ug, Vg = " << m_abl_forcing[1]/m_coriolis_factor << " " << - -m_abl_forcing[0]/m_coriolis_factor << std::endl; - if (m_write_force_timetable && amrex::ParallelDescriptor::IOProcessor() && (t_step % m_force_outfreq == 0) && @@ -91,7 +81,9 @@ public: outfile.open(m_force_timetable, std::ios::out | std::ios::app); outfile << std::setprecision(17) << nph_time << "\t" << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t" - << 0.0 << std::endl; + << 0.0 << "\t" << m_abl_forcing[1] / m_coriolis_factor + << "\t" << -m_abl_forcing[0] / m_coriolis_factor + << "\t" << m_zi << std::endl; } } diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp index ab01fd4543..6a3f345b93 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp @@ -47,13 +47,15 @@ ABLForcing::ABLForcing(const CFDSim& sim) m_write_force_timetable = pp_abl.contains("forcing_timetable_output_file"); if (m_write_force_timetable) { + amrex::Print() << "Here!!!!" << std::endl; 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()) { + amrex::Print() << m_force_timetable << std::endl; std::ofstream outfile; outfile.open(m_force_timetable, std::ios::out); - outfile << "time\tfx\tfy\tfz\n"; + outfile << "time\tfx\tfy\tfz\tUgx\tUgy\tzi\n"; } } @@ -81,7 +83,6 @@ ABLForcing::ABLForcing(const CFDSim& sim) m_coriolis_factor = 2.0 * omega * sinphi; } - amrex::Print() << "free atmosphere damping = " << m_fa_damping << std::endl; for (int i = 0; i < AMREX_SPACEDIM; ++i) { @@ -127,7 +128,6 @@ void ABLForcing::operator()( const bool fa_damping = m_fa_damping; const bool fa_detect_height = m_fa_detect_height; const amrex::Real fa_height = fa_detect_height ? m_zi : m_fa_height; - amrex::AllPrint() << "fa_height = " << fa_height << " fa_detect_height = " << fa_detect_height << " m_zi = " << m_zi << std::endl; const amrex::Real fa_time_start = m_fa_time_start; const amrex::Real fa_time_end = m_fa_time_end; const amrex::Real fa_tau = m_fa_tau; @@ -154,10 +154,6 @@ void ABLForcing::operator()( const auto& vof = (*m_vof)(lev).const_array(mfi); - //amrex::Print() << "fa_height = " << fa_height << std::endl; - //amrex::Print() << "Sx, Sy = " << dudt << ", " << dvdt << std::endl; - //amrex::Print() << "Ug, Vg = " << fa_u << ", " << fa_v << std::endl; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::IntVect iv(i, j, k); const amrex::Real z = problo[idir] + (iv[idir] + 0.5) * dx[idir]; @@ -207,9 +203,6 @@ void ABLForcing::operator()( src_fa_damping_y = (1.0/fa_tau) * (fa_v - vmean); } - amrex::AllPrint() << "z = " << z << ", fa_u = " << fa_u << ", umean = " << umean << ", src_fa_damping_x = " << src_fa_damping_x << ", dpdx = " << fac*dudt << std::endl; - amrex::AllPrint() << "z = " << z << ", fa_v = " << fa_v << ", vmean = " << vmean << ", src_fa_damping_y = " << src_fa_damping_y << ", dpdy = " << fac*dvdt << std::endl; - // Sum up the source term src_term(i, j, k, 0) += (fac * dudt) + src_fa_damping_x; src_term(i, j, k, 1) += (fac * dvdt) + src_fa_damping_y; diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index eaa3dde427..9aaae2e64b 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -174,7 +174,6 @@ void ABL::post_regrid_actions() { m_abl_anelastic->post_regrid_actions(); } void ABL::pre_advance_work() { const auto& zi = m_stats->zi(); - amrex::Print() << "in ABL.H: zi = " << zi << std::endl; const auto& vel_pa = m_stats->vel_profile(); m_abl_wall_func.update_umean( m_stats->vel_profile(), m_stats->theta_profile_fine()); diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 769ba515c9..f9dc62c880 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -101,13 +101,11 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) "Cannot find surface_temp_flux_timetable file: " + m_surf_temp_flux_timetable); } - amrex::Real data_time; - amrex::Real data_value; - ifh.ignore(std::numeric_limits::max(), '\n'); + amrex::Real data_time, data_value; while (ifh >> data_time) { - ifh >> data_value; m_surf_temp_flux_time.push_back(data_time); m_surf_temp_flux_value.push_back(data_value); + ifh.ignore(std::numeric_limits::max(), '\n'); } } } @@ -141,18 +139,11 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) "Cannot find near_surface_temp_timetable file: " + m_nearsurf_temp_timetable); } - amrex::Real data_time; - amrex::Real data_value; - ifh.ignore(std::numeric_limits::max(), '\n'); - while (ifh >> data_time) { - ifh >> data_value; + amrex::Real data_time, data_value; + while (ifh >> data_time >> data_value) { m_nearsurf_temp_time.push_back(data_time); m_nearsurf_temp_value.push_back(data_value); - } - - for (int i =0; i < m_nearsurf_temp_time.size(); i++) - { - amrex::Print() << i << " " << m_nearsurf_temp_time[i] << " " << m_nearsurf_temp_value[i] << std::endl; + ifh.ignore(std::numeric_limits::max(), '\n'); } } @@ -178,13 +169,11 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) "Cannot find surface_temp_timetable file: " + m_surf_temp_timetable); } - amrex::Real data_time; - amrex::Real data_value; - ifh.ignore(std::numeric_limits::max(), '\n'); - while (ifh >> data_time) { - ifh >> data_value; + amrex::Real data_time, data_value; + while (ifh >> data_time >> data_value) { m_surf_temp_time.push_back(data_time); m_surf_temp_value.push_back(data_value); + ifh.ignore(std::numeric_limits::max(), '\n'); } } // - Specified surface temperature rate mode: diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index 64a94b6990..c7b5718322 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -66,11 +66,17 @@ void MOData::update_fluxes(int max_iters) surf_temp = alpha_h * surf_temp_flux * (std::log(zref / z0t) - psi_h_zref) / (utau * kappa) + theta_mean; + near_surf_temp = alpha_h * surf_temp_flux * (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; break; case ThetaCalcType::SURFACE_TEMPERATURE: surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa / (alpha_h * (std::log(zref / z0t) - psi_h_zref)); + near_surf_temp = alpha_h * surf_temp_flux * (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; break; case ThetaCalcType::NEAR_SURFACE_TEMPERATURE: From 3cec074aa14639048747c762d655ff63bee30a1b Mon Sep 17 00:00:00 2001 From: Matt Churchfield Date: Tue, 26 Aug 2025 11:06:17 -0600 Subject: [PATCH 09/12] Removed a small print statement about what z0 and z0t are. --- amr-wind/wind_energy/MOData.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index c7b5718322..73a49758a0 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -56,8 +56,6 @@ void MOData::update_fluxes(int max_iters) amrex::Real psi_h_zref = 0.0; utau = kappa * vmag_mean / (std::log(zref / z0)); - amrex::Print() << "Using z0 = " << z0 << ", z0t = " << z0t << std::endl; - int iter = 0; do { utau_iter = utau; From 1edfbdb82aa2d340891517b561997806f9e1b61a Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 4 Sep 2025 13:33:50 -0600 Subject: [PATCH 10/12] formatting --- .../icns/source_terms/ABLForcing.H | 11 ++- .../icns/source_terms/ABLForcing.cpp | 19 +++-- amr-wind/wind_energy/ABL.cpp | 3 +- amr-wind/wind_energy/ABLStats.H | 5 +- amr-wind/wind_energy/ABLStats.cpp | 5 +- amr-wind/wind_energy/ABLWallFunction.cpp | 70 +++++++++---------- amr-wind/wind_energy/MOData.H | 19 ++--- amr-wind/wind_energy/MOData.cpp | 43 +++++++----- amr-wind/wind_energy/ShearStress.H | 4 +- 9 files changed, 88 insertions(+), 91 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H index 287e46f6d1..739983e64b 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.H @@ -40,10 +40,7 @@ public: m_target_vel[1] = uy; } - inline void set_boundary_layer_height(amrex::Real zi) - { - m_zi = zi; - } + inline void set_boundary_layer_height(amrex::Real zi) { m_zi = zi; } inline void set_mean_velocities(amrex::Real ux, amrex::Real uy) { @@ -81,9 +78,9 @@ public: outfile.open(m_force_timetable, std::ios::out | std::ios::app); outfile << std::setprecision(17) << nph_time << "\t" << m_abl_forcing[0] << "\t" << m_abl_forcing[1] << "\t" - << 0.0 << "\t" << m_abl_forcing[1] / m_coriolis_factor - << "\t" << -m_abl_forcing[0] / m_coriolis_factor - << "\t" << m_zi << std::endl; + << 0.0 << "\t" << m_abl_forcing[1] / m_coriolis_factor + << "\t" << -m_abl_forcing[0] / m_coriolis_factor << "\t" + << m_zi << std::endl; } } diff --git a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp index 6a3f345b93..bdc00fba77 100644 --- a/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/ABLForcing.cpp @@ -63,7 +63,7 @@ ABLForcing::ABLForcing(const CFDSim& sim) pp_abl.query("free_atmosphere_damping", m_fa_damping); if (m_fa_damping) { pp_abl.query("detect_free_atmosphere_height", m_fa_detect_height); - if (! m_fa_detect_height) { + if (!m_fa_detect_height) { pp_abl.get("free_atmosphere_height", m_fa_height); } pp_abl.query("free_atmosphere_damping_start_time", m_fa_time_start); @@ -75,7 +75,7 @@ ABLForcing::ABLForcing(const CFDSim& sim) pp_coriolis.query("rotational_time_period", rot_time_period); amrex::Real omega = utils::two_pi() / rot_time_period; - + amrex::Real latitude = 90.0; pp_coriolis.query("latitude", latitude); latitude = utils::radians(latitude); @@ -83,7 +83,6 @@ ABLForcing::ABLForcing(const CFDSim& sim) m_coriolis_factor = 2.0 * omega * sinphi; } - for (int i = 0; i < AMREX_SPACEDIM; ++i) { m_mean_vel[i] = m_target_vel[i]; @@ -129,9 +128,9 @@ void ABLForcing::operator()( const bool fa_detect_height = m_fa_detect_height; const amrex::Real fa_height = fa_detect_height ? m_zi : m_fa_height; const amrex::Real fa_time_start = m_fa_time_start; - const amrex::Real fa_time_end = m_fa_time_end; + const amrex::Real fa_time_end = m_fa_time_end; const amrex::Real fa_tau = m_fa_tau; - const amrex::Real fa_u = dvdt / m_coriolis_factor; + const amrex::Real fa_u = dvdt / m_coriolis_factor; const amrex::Real fa_v = -dudt / m_coriolis_factor; const auto& current_time = m_time.current_time(); @@ -153,7 +152,6 @@ void ABLForcing::operator()( const auto& vof = (*m_vof)(lev).const_array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { amrex::IntVect iv(i, j, k); const amrex::Real z = problo[idir] + (iv[idir] + 0.5) * dx[idir]; @@ -195,12 +193,11 @@ void ABLForcing::operator()( // wind. Nudge the free atmosphere toward that geostrophic wind.) amrex::Real src_fa_damping_x = 0.0; amrex::Real src_fa_damping_y = 0.0; - if (fa_damping && - ((nph_time >= fa_time_start) && - (nph_time <= fa_time_end)) && + if (fa_damping && + ((nph_time >= fa_time_start) && (nph_time <= fa_time_end)) && (z >= fa_height)) { - src_fa_damping_x = (1.0/fa_tau) * (fa_u - umean); - src_fa_damping_y = (1.0/fa_tau) * (fa_v - vmean); + src_fa_damping_x = (1.0 / fa_tau) * (fa_u - umean); + src_fa_damping_y = (1.0 / fa_tau) * (fa_v - vmean); } // Sum up the source term diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index 9aaae2e64b..7d1aa3b751 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -204,8 +204,7 @@ void ABL::pre_advance_work() m_abl_forcing->set_mean_velocities(vx, vy); m_abl_forcing->set_boundary_layer_height(zi); - m_abl_forcing->mean_velocity_update( - m_stats->vel_profile_coarse()); + m_abl_forcing->mean_velocity_update(m_stats->vel_profile_coarse()); } if (m_abl_mean_bous != nullptr) { diff --git a/amr-wind/wind_energy/ABLStats.H b/amr-wind/wind_energy/ABLStats.H index cdc14ffb7e..3bcc6490b5 100644 --- a/amr-wind/wind_energy/ABLStats.H +++ b/amr-wind/wind_energy/ABLStats.H @@ -58,10 +58,7 @@ public: //! Compute height of capping inversion void compute_zi(); - amrex::Real zi() const override - { - return m_zi; - }; + amrex::Real zi() const override { return m_zi; }; //! Return vel plane averaging instance const VelPlaneAveragingFine& vel_profile() const override diff --git a/amr-wind/wind_energy/ABLStats.cpp b/amr-wind/wind_energy/ABLStats.cpp index 604611d1df..b19b0fbde7 100644 --- a/amr-wind/wind_energy/ABLStats.cpp +++ b/amr-wind/wind_energy/ABLStats.cpp @@ -234,7 +234,7 @@ void ABLStats::compute_zi() // method of Sullivan et al. in "Structure of the entrainment zone // capping the convective atmospheric boundary layer," JAS, Vol. 55, // 1998. - // + // // In this method, for every x,y location in the computational domain, // the maximum d(\theta)/dz is found. You then have an x,y array // of x,y local z_i. z_i is then averaged over the x,y array to give @@ -305,7 +305,8 @@ void ABLStats::compute_zi() // It is necessary for other processors to know , so broadcast // it back out to all processors. - amrex::ParallelDescriptor::Bcast(&m_zi, 1, amrex::ParallelDescriptor::IOProcessorNumber()); + amrex::ParallelDescriptor::Bcast( + &m_zi, 1, amrex::ParallelDescriptor::IOProcessorNumber()); } void ABLStats::process_output() diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index f9dc62c880..2c33a3c6bc 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -62,9 +62,6 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) << std::endl; } - - - // Check for type of surface heating and read in data. // - Specified surface temperature flux mode: if ((pp.contains("surface_temp_flux")) || @@ -93,7 +90,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) amrex::Print() << "ABLWallFunction: Surface temperature flux time table " "mode is selected." - << std::endl; + << std::endl; if (!m_surf_temp_flux_timetable.empty()) { std::ifstream ifh(m_surf_temp_flux_timetable, std::ios::in); if (!ifh.good()) { @@ -105,12 +102,13 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) while (ifh >> data_time) { m_surf_temp_flux_time.push_back(data_time); m_surf_temp_flux_value.push_back(data_value); - ifh.ignore(std::numeric_limits::max(), '\n'); + ifh.ignore( + std::numeric_limits::max(), '\n'); } } } - // - Specified near-surface temperature mode: + // - Specified near-surface temperature mode: } else if (pp.contains("near_surface_temp_timetable")) { m_surf_temp_flux = false; m_surf_temp = false; @@ -119,7 +117,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_nearsurf_temp_use_table = true; if (pp.contains("near_surface_height")) { - pp.get("near_surface_height",m_mo.zNearSurf); + pp.get("near_surface_height", m_mo.zNearSurf); } else { amrex::Abort( "ABLWallFunction: near_surface_height must be specified when " @@ -128,9 +126,10 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) pp.query("near_surface_temp_timetable", m_nearsurf_temp_timetable); - amrex::Print() << "ABLWallFunction: Near-surface temperature time table " - "mode is selected." - << std::endl; + amrex::Print() + << "ABLWallFunction: Near-surface temperature time table " + "mode is selected." + << std::endl; if (!m_nearsurf_temp_timetable.empty()) { std::ifstream ifh(m_nearsurf_temp_timetable, std::ios::in); @@ -147,9 +146,10 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) } } - // - Specified surface temperature mode: - } else if ((pp.contains("surface_temp_timetable")) || - (pp.contains("surface_temp_rate"))) { + // - Specified surface temperature mode: + } else if ( + (pp.contains("surface_temp_timetable")) || + (pp.contains("surface_temp_rate"))) { m_surf_temp_flux = false; m_surf_temp = true; m_nearsurf_temp = false; @@ -173,10 +173,11 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) while (ifh >> data_time >> data_value) { m_surf_temp_time.push_back(data_time); m_surf_temp_value.push_back(data_value); - ifh.ignore(std::numeric_limits::max(), '\n'); + ifh.ignore( + std::numeric_limits::max(), '\n'); } } - // - Specified surface temperature rate mode: + // - Specified surface temperature rate mode: } else if (pp.contains("surface_temp_rate")) { m_surf_temp_use_table = false; m_surf_temp_use_rate = true; @@ -191,24 +192,28 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) pp.get("surface_temp_init", m_surf_temp_init); } else { amrex::Print() - << "ABLWallFunction: Initial surface temperature not found for " + << "ABLWallFunction: Initial surface temperature not found " + "for " "ABL. Assuming to be equal to the reference temperature" << std::endl; - m_surf_temp_init = sim.transport_model().reference_temperature(); + m_surf_temp_init = + sim.transport_model().reference_temperature(); } if (pp.contains("surface_temp_rate_tstart")) { pp.get("surface_temp_rate_tstart", m_surf_temp_rate_tstart); } else { amrex::Print() - << "ABLWallFunction: Surface temperature heating/cooling start " + << "ABLWallFunction: Surface temperature heating/cooling " + "start " "time (surface_temp_rate_tstart) not found for ABL. " "Assuming zero." << m_surf_temp_rate_tstart << std::endl; } } - // - If no surface heating mode is specified, default to no surface heating. + // - If no surface heating mode is specified, default to no surface + // heating. } else { amrex::Print() << "ABLWallFunction: None of 'surface_temp_flux,' " "'surface_temp_rate', 'surface_temp_timetable,' " @@ -217,9 +222,6 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) << std::endl; } - - - // Check to see if the simulation is non-periodic. If so, recommend // using the wall model locally instead of planar averaging. if (pp.contains("inflow_outflow_mode")) { @@ -228,16 +230,14 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) pp.query("wf_velocity", m_wf_vel); pp.query("wf_vmag", m_wf_vmag); pp.query("wf_theta", m_wf_theta); - amrex::Print() << "ABLWallFunction: Inflow/Outflow mode is turned " - "on. It is recommended that the wall shear stress " - "type is set to local IF the flow is heterogeneous." - << std::endl; + amrex::Print() + << "ABLWallFunction: Inflow/Outflow mode is turned " + "on. It is recommended that the wall shear stress " + "type is set to local IF the flow is heterogeneous." + << std::endl; } } - - - // Set the Monin-Obukhov algorith type. if (m_surf_temp_flux) { m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; @@ -249,9 +249,6 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; } - - - // Set gravity strength in Monin-Obukhov. m_mo.gravity = utils::vec_mag(m_gravity.data()); } @@ -273,8 +270,8 @@ void ABLWallFunction::update_umean( if (m_surf_temp_flux && m_surf_temp_flux_use_table) { m_mo.surf_temp_flux = amr_wind::interp::linear( m_surf_temp_flux_time, m_surf_temp_flux_value, time.current_time()); - amrex::Print() << "Current surface temperature flux: " << m_mo.surf_temp_flux - << std::endl; + amrex::Print() << "Current surface temperature flux: " + << m_mo.surf_temp_flux << std::endl; } // - if specifying surface temperature if (m_surf_temp) { @@ -291,11 +288,12 @@ void ABLWallFunction::update_umean( } amrex::Print() << "Current surface temperature: " << m_mo.surf_temp << std::endl; - // - if specifying near-surface temperature + // - if specifying near-surface temperature } else if (m_nearsurf_temp) { if (m_nearsurf_temp_use_table) { m_mo.near_surf_temp = amr_wind::interp::linear( - m_nearsurf_temp_time, m_nearsurf_temp_value, time.current_time()); + m_nearsurf_temp_time, m_nearsurf_temp_value, + time.current_time()); } } diff --git a/amr-wind/wind_energy/MOData.H b/amr-wind/wind_energy/MOData.H index 1526360999..b2a6342cd7 100644 --- a/amr-wind/wind_energy/MOData.H +++ b/amr-wind/wind_energy/MOData.H @@ -22,17 +22,17 @@ struct MOData enum class ThetaCalcType { HEAT_FLUX = 0, ///< Heat-flux specified SURFACE_TEMPERATURE = 1, ///< Surface temperature specified - NEAR_SURFACE_TEMPERATURE = 2, ///< Temperature specified above but + NEAR_SURFACE_TEMPERATURE = 2, ///< Temperature specified above but // near the surface }; - amrex::Real zref{0.0}; ///< Reference height (m) - amrex::Real zNearSurf{2.0}; ///< Near-surface temperature driving height (m) - amrex::Real z0{0.1}; ///< Aerodynamic roughness height (m) - amrex::Real z0t{z0}; ///< Thermal roughness height (m) - amrex::Real utau; ///< Friction velocity (m/s) - amrex::Real kappa{0.41}; ///< von Karman constant - amrex::Real gravity{9.81}; ///< Acceleration due to gravity (m/s^2) + amrex::Real zref{0.0}; ///< Reference height (m) + amrex::Real zNearSurf{2.0}; ///< Near-surface temperature driving height (m) + amrex::Real z0{0.1}; ///< Aerodynamic roughness height (m) + amrex::Real z0t{z0}; ///< Thermal roughness height (m) + amrex::Real utau; ///< Friction velocity (m/s) + amrex::Real kappa{0.41}; ///< von Karman constant + amrex::Real gravity{9.81}; ///< Acceleration due to gravity (m/s^2) amrex::Real obukhov_len{1.0e16}; ///< Non-dimensional Obukhov length amrex::RealArray vel_mean; ///< Mean velocity (at zref) @@ -43,7 +43,8 @@ struct MOData amrex::Real surf_temp_flux{0.0}; ///< Heat flux amrex::Real surf_temp; ///< Instantaneous surface temperature - amrex::Real near_surf_temp; ///< Instantaneous near-surface temperature (at zNearSurf). + amrex::Real near_surf_temp; ///< Instantaneous near-surface temperature (at + ///< zNearSurf). amrex::Real alpha_h{1.0}; amrex::Real gamma_m{5.0}; diff --git a/amr-wind/wind_energy/MOData.cpp b/amr-wind/wind_energy/MOData.cpp index 73a49758a0..bc31d46e3f 100644 --- a/amr-wind/wind_energy/MOData.cpp +++ b/amr-wind/wind_energy/MOData.cpp @@ -61,33 +61,38 @@ void MOData::update_fluxes(int max_iters) utau_iter = utau; switch (alg_type) { case ThetaCalcType::HEAT_FLUX: - surf_temp = alpha_h * surf_temp_flux * (std::log(zref / z0t) - psi_h_zref) / - (utau * kappa) + - theta_mean; - near_surf_temp = alpha_h * surf_temp_flux * (std::log(2.0 / z0t) - psi_h_zref) / + surf_temp = alpha_h * surf_temp_flux * + (std::log(zref / z0t) - psi_h_zref) / (utau * kappa) + theta_mean; + near_surf_temp = alpha_h * surf_temp_flux * + (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; break; case ThetaCalcType::SURFACE_TEMPERATURE: surf_temp_flux = -(theta_mean - surf_temp) * utau * kappa / (alpha_h * (std::log(zref / z0t) - psi_h_zref)); - near_surf_temp = alpha_h * surf_temp_flux * (std::log(2.0 / z0t) - psi_h_zref) / - (utau * kappa) + - theta_mean; + near_surf_temp = alpha_h * surf_temp_flux * + (std::log(2.0 / z0t) - psi_h_zref) / + (utau * kappa) + + theta_mean; break; case ThetaCalcType::NEAR_SURFACE_TEMPERATURE: amrex::Real psi_h_zNearSurf = calc_psi_h(zNearSurf / obukhov_len); amrex::Real a11 = 1.0; - amrex::Real a12 = -(alpha_h / (kappa * utau)) * (std::log(zref / z0t) - psi_h_zref); + amrex::Real a12 = -(alpha_h / (kappa * utau)) * + (std::log(zref / z0t) - psi_h_zref); amrex::Real a21 = 1.0; - amrex::Real a22 = -(alpha_h / (kappa * utau)) * (std::log(zNearSurf / z0t) - psi_h_zNearSurf); + amrex::Real a22 = -(alpha_h / (kappa * utau)) * + (std::log(zNearSurf / z0t) - psi_h_zNearSurf); - amrex::Real coeff = 1.0 / (a11*a22 - a12*a21); + amrex::Real coeff = 1.0 / (a11 * a22 - a12 * a21); - surf_temp = coeff * (a22*theta_mean - a12*near_surf_temp); - surf_temp_flux = coeff * (-a21*theta_mean + a11*near_surf_temp); + surf_temp = coeff * (a22 * theta_mean - a12 * near_surf_temp); + surf_temp_flux = coeff * (-a21 * theta_mean + a11 * near_surf_temp); break; } @@ -107,17 +112,21 @@ void MOData::update_fluxes(int max_iters) utau = kappa * vmag_mean / (std::log(zref / z0) - psi_m_zref); ++iter; - amrex::Print() << "iteration: " << iter << ", T(" << zNearSurf << " m) = " << near_surf_temp - << " K, T(" << zref << " m) = " << theta_mean << "K , L = " << obukhov_len - << " m, uStar = " << utau << "m/s, Qs = " << surf_temp_flux - << " K-m/s, theta_0 = " << surf_temp << " K" << std::endl; + amrex::Print() << "iteration: " << iter << ", T(" << zNearSurf + << " m) = " << near_surf_temp << " K, T(" << zref + << " m) = " << theta_mean << "K , L = " << obukhov_len + << " m, uStar = " << utau + << "m/s, Qs = " << surf_temp_flux + << " K-m/s, theta_0 = " << surf_temp << " K" + << std::endl; } while ((std::abs(utau_iter - utau) > 1e-5) && iter <= max_iters); if (iter >= max_iters) { amrex::Print() << "MOData::update_fluxes: Convergence criteria not met after " << max_iters << " iterations\nObuhov length = " << obukhov_len - << " zeta = (z_ref/L) = " << zeta << "\npsi_m(z_ref/L) = " << psi_m_zref + << " zeta = (z_ref/L) = " << zeta + << "\npsi_m(z_ref/L) = " << psi_m_zref << " psi_h(z_ref/L) = " << psi_h_zref << "\nutau = " << utau << " Tsurf = " << surf_temp << " q = " << surf_temp_flux << std::endl; diff --git a/amr-wind/wind_energy/ShearStress.H b/amr-wind/wind_energy/ShearStress.H index e921a6323b..76aff3e26c 100644 --- a/amr-wind/wind_energy/ShearStress.H +++ b/amr-wind/wind_energy/ShearStress.H @@ -217,9 +217,7 @@ struct ShearStressDonelan flux = 0.0012 * wspd_mean * (theta_surface - theta_mean); break; case amr_wind::MOData::ThetaCalcType::NEAR_SURFACE_TEMPERATURE: - amrex::Print() - << "To be implemented..." - << std::endl; + amrex::Print() << "To be implemented..." << std::endl; break; } return flux; From 932295a865e339bc10789822c230a3b7967253cc Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 4 Sep 2025 13:38:47 -0600 Subject: [PATCH 11/12] spelling --- amr-wind/wind_energy/ABLWallFunction.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 2c33a3c6bc..9369daef48 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -238,7 +238,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) } } - // Set the Monin-Obukhov algorith type. + // Set the Monin-Obukhov algorithm type. if (m_surf_temp_flux) { m_mo.alg_type = MOData::ThetaCalcType::HEAT_FLUX; } else if (m_surf_temp) { From 169785ac79c5abdcc3f35180741b421296cfcb15 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Thu, 4 Sep 2025 14:01:40 -0600 Subject: [PATCH 12/12] Update amr-wind/wind_energy/ABLWallFunction.cpp: read data_value --- amr-wind/wind_energy/ABLWallFunction.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/wind_energy/ABLWallFunction.cpp b/amr-wind/wind_energy/ABLWallFunction.cpp index 9369daef48..70a8736508 100644 --- a/amr-wind/wind_energy/ABLWallFunction.cpp +++ b/amr-wind/wind_energy/ABLWallFunction.cpp @@ -99,7 +99,7 @@ ABLWallFunction::ABLWallFunction(const CFDSim& sim) m_surf_temp_flux_timetable); } amrex::Real data_time, data_value; - while (ifh >> data_time) { + while (ifh >> data_time >> data_value) { m_surf_temp_flux_time.push_back(data_time); m_surf_temp_flux_value.push_back(data_value); ifh.ignore(