big update for ref temp

amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H

@@ -43,15 +43,15 @@ private:
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_ht;
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_vals;
 
-    //! Reference temperature (Kelvin)
-    amrex::Real m_ref_theta{300.0};
-
     //! Transport model
     const transport::TransportModel& m_transport;
 
     //! Thermal expansion coefficient
     std::unique_ptr<ScratchField> m_beta;
 
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
+
     int m_axis{2};
 
     bool m_const_profile{false};

amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp

@@ -8,14 +8,6 @@
 namespace amr_wind::pde::icns {
 
 /** Boussinesq buoyancy source term for ABL simulations
- *
- *  Reads in the following parameters from `ABLMeanBoussinesq` namespace:
- *
- *  - `reference_temperature` (Mandatory) temperature (`T0`) in Kelvin
- *  - `thermal_expansion_coeff` Optional, default = `1.0 / T0`
- *  - `gravity` acceleration due to gravity (m/s)
- *  - `read_temperature_profile`
- *  - `tprofile_filename`
  */
 ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim)
     : m_mesh(sim.mesh()), m_transport(sim.transport_model())
@@ -24,31 +16,67 @@ ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim)
     const auto& abl = sim.physics_manager().get<amr_wind::ABL>();
     abl.register_mean_boussinesq_term(this);
 
-    amrex::ParmParse pp_boussinesq_buoyancy("BoussinesqBuoyancy");
-    pp_boussinesq_buoyancy.get("reference_temperature", m_ref_theta);
-
     m_beta = m_transport.beta();
+    m_ref_theta = m_transport.ref_theta();
 
     // gravity in `incflo` namespace
     amrex::ParmParse pp_incflo("incflo");
     pp_incflo.queryarr("gravity", m_gravity);
 
-    bool read_temp_prof = false;
-    pp_boussinesq_buoyancy.query("read_temperature_profile", read_temp_prof);
-
-    if ((!pp_boussinesq_buoyancy.contains("read_temperature_profile") &&
-         pp_boussinesq_buoyancy.contains("tprofile_filename")) ||
-        read_temp_prof) {
+    // Backwards compatibility
+    amrex::ParmParse pp_boussinesq_buoyancy("BoussinesqBuoyancy");
+    amrex::ParmParse pp_abl("ABLMeanBoussinesq");
 
-        m_const_profile = true;
+    bool read_temp_prof = false;
+    if (pp_abl.contains("read_temperature_profile")) {
+        pp_abl.get("read_temperature_profile", read_temp_prof);
+        if (pp_boussinesq_buoyancy.contains("read_temperature_profile")) {
+            amrex::Print()
+                << "WARNING: BoussinesqBuoyancy.read_temperature_profile "
+                   "option has been deprecated in favor of "
+                   "ABLMeanBoussinesq.read_temperature_profile. Ignoring the"
+                   "BoussinesqBuoyancy option in favor of the "
+                   "ABLMeanBoussinesq "
+                   "option."
+                << std::endl;
+        }
+    } else if (pp_boussinesq_buoyancy.contains("read_temperature_profile")) {
+        amrex::Print()
+            << "WARNING: BoussinesqBuoyancy.read_temperature_profile option "
+               "has been deprecated in favor of "
+               "ABLMeanBoussinesq.read_temperature_profile. Please replace "
+               "this option."
+            << std::endl;
+        pp_boussinesq_buoyancy.get("read_temperature_profile", read_temp_prof);
+    }
 
-        std::string tprofile_filename;
+    std::string tprofile_filename;
+    if (pp_abl.contains("temperature_profile_filename")) {
+        pp_abl.get("temperature_profile_filename", tprofile_filename);
+        if (pp_boussinesq_buoyancy.contains("tprofile_filename")) {
+            amrex::Print() << "WARNING: BoussinesqBuoyancy.tprofile_filename "
+                              "option has been deprecated in favor of "
+                              "ABLMeanBoussinesq.temperature_profile_filename. "
+                              "Ignoring the"
+                              "BoussinesqBuoyancy option in favor of the "
+                              "ABLMeanBoussinesq "
+                              "option."
+                           << std::endl;
+        }
+    } else if (pp_boussinesq_buoyancy.contains("tprofile_filename")) {
+        amrex::Print()
+            << "WARNING: BoussinesqBuoyancy.tprofile_filename option "
+               "has been deprecated in favor of "
+               "ABLMeanBoussinesq.temperature_profile_filename. Please replace "
+               "this option."
+            << std::endl;
         pp_boussinesq_buoyancy.get("tprofile_filename", tprofile_filename);
+    }
 
+    if ((read_temp_prof) && (tprofile_filename.empty())) {
+        m_const_profile = true;
         read_temperature_profile(tprofile_filename);
-
     } else {
-
         mean_temperature_init(abl.abl_statistics().theta_profile());
     }
 }
@@ -64,8 +92,8 @@ void ABLMeanBoussinesq::operator()(
 {
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
-    const amrex::Real T0 = m_ref_theta;
     const auto& beta = (*m_beta)(lev).const_array(mfi);
+    const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi);
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
 
@@ -79,6 +107,7 @@ void ABLMeanBoussinesq::operator()(
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::IntVect iv(i, j, k);
         const amrex::Real ht = problo[idir] + (iv[idir] + 0.5) * dx[idir];
+        const amrex::Real T0 = ref_theta(i, j, k);
         const amrex::Real temp =
             amr_wind::interp::linear(theights, theights_end, tvals, ht);
         const amrex::Real fac = beta(i, j, k) * (temp - T0);

amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H

@@ -40,6 +40,9 @@ private:
 
     //! Thermal expansion coefficient
     std::unique_ptr<ScratchField> m_beta;
+
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 };
 } // namespace amr_wind::pde::icns
 

amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp

@@ -14,6 +14,7 @@ BoussinesqBuoyancy::BoussinesqBuoyancy(const CFDSim& sim)
     , m_transport(sim.transport_model())
 {
     m_beta = m_transport.beta();
+    m_ref_theta = m_transport.ref_theta();
 
     // gravity in `incflo` namespace
     amrex::ParmParse pp_incflo("incflo");
@@ -29,16 +30,17 @@ void BoussinesqBuoyancy::operator()(
     const FieldState fstate,
     const amrex::Array4<amrex::Real>& src_term) const
 {
-    const amrex::Real T0 = m_transport.reference_temperature();
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
 
     const auto& temp =
         m_temperature.state(field_impl::phi_state(fstate))(lev).const_array(
             mfi);
     const auto& beta = (*m_beta)(lev).const_array(mfi);
+    const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi);
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const amrex::Real T = temp(i, j, k, 0);
+        const amrex::Real T0 = ref_theta(i, j, k);
         const amrex::Real fac = beta(i, j, k) * (T0 - T);
 
         src_term(i, j, k, 0) += gravity[0] * fac;

amr-wind/equation_systems/temperature/source_terms/DragTempForcing.H

@@ -4,6 +4,7 @@
 #include "amr-wind/equation_systems/temperature/TemperatureSource.H"
 #include "amr-wind/core/SimTime.H"
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/transport_models/TransportModel.H"
 
 namespace amr_wind::pde::temperature {
 
@@ -29,7 +30,12 @@ private:
     const Field& m_velocity;
     const Field& m_temperature;
     amrex::Real m_drag_coefficient{1.0};
-    amrex::Real m_reference_temperature{300.0};
+
+    //! Transport model
+    const transport::TransportModel& m_transport;
+
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 };
 
 } // namespace amr_wind::pde::temperature

amr-wind/equation_systems/temperature/source_terms/DragTempForcing.cpp

@@ -13,10 +13,11 @@ DragTempForcing::DragTempForcing(const CFDSim& sim)
     , m_mesh(sim.mesh())
     , m_velocity(sim.repo().get_field("velocity"))
     , m_temperature(sim.repo().get_field("temperature"))
+    , m_transport(sim.transport_model())
 {
     amrex::ParmParse pp("DragTempForcing");
     pp.query("drag_coefficient", m_drag_coefficient);
-    pp.query("reference_temperature", m_reference_temperature);
+    m_ref_theta = m_transport.ref_theta();
 }
 
 DragTempForcing::~DragTempForcing() = default;
@@ -44,7 +45,7 @@ void DragTempForcing::operator()(
     const auto& geom = m_mesh.Geom(lev);
     const auto& dx = geom.CellSizeArray();
     const amrex::Real drag_coefficient = m_drag_coefficient / dx[2];
-    const amrex::Real reference_temperature = m_reference_temperature;
+    const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi);
     const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
     const amrex::Real cd_max = 10.0;
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -54,9 +55,9 @@ void DragTempForcing::operator()(
         const amrex::Real m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
         const amrex::Real Cd =
             std::min(drag_coefficient / (m + tiny), cd_max / dx[2]);
+        const amrex::Real T0 = ref_theta(i, j, k);
         src_term(i, j, k, 0) -=
-            (Cd * (temperature(i, j, k, 0) - reference_temperature) *
-             blank(i, j, k, 0));
+            (Cd * (temperature(i, j, k, 0) - T0) * blank(i, j, k, 0));
     });
 }
 

amr-wind/equation_systems/tke/source_terms/KransAxell.H

@@ -1,6 +1,7 @@
 #ifndef KRANSAXELL_H
 #define KRANSAXELL_H
 
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/equation_systems/tke/TKESource.H"
 
 namespace amr_wind::pde::tke {
@@ -35,7 +36,6 @@ private:
     Field& m_tke;
     amrex::Real m_Cmu{0.556};
     amrex::Real m_heat_flux{0.0};
-    amrex::Real m_ref_temp{300.0};
     amrex::Real m_z0{0.1};
     amrex::Real m_kappa{0.41};
     amrex::Real m_sponge_start{600};
@@ -45,6 +45,12 @@ private:
     const CFDSim& m_sim;
     const amrex::AmrCore& m_mesh;
     const Field& m_velocity;
+
+    //! Transport model
+    const transport::TransportModel& m_transport;
+
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 };
 
 } // namespace amr_wind::pde::tke

amr-wind/equation_systems/tke/source_terms/KransAxell.cpp

@@ -16,20 +16,21 @@ KransAxell::KransAxell(const CFDSim& sim)
     , m_sim(sim)
     , m_mesh(sim.mesh())
     , m_velocity(sim.repo().get_field("velocity"))
+    , m_transport(sim.transport_model())
 {
     AMREX_ALWAYS_ASSERT(sim.turbulence_model().model_name() == "KLAxell");
     auto coeffs = sim.turbulence_model().model_coeffs();
     amrex::ParmParse pp("ABL");
     pp.query("Cmu", m_Cmu);
     pp.query("kappa", m_kappa);
     pp.query("surface_roughness_z0", m_z0);
-    pp.query("reference_temperature", m_ref_temp);
     pp.query("surface_temp_flux", m_heat_flux);
     pp.query("meso_sponge_start", m_sponge_start);
     {
         amrex::ParmParse pp_incflow("incflo");
         pp_incflow.queryarr("gravity", m_gravity);
     }
+    m_ref_theta = m_transport.ref_theta();
 }
 
 KransAxell::~KransAxell() = default;
@@ -48,20 +49,21 @@ void KransAxell::operator()(
     const auto& buoy_prod_arr = (this->m_buoy_prod)(lev).array(mfi);
     const auto& dissip_arr = (this->m_dissip)(lev).array(mfi);
     const auto& tke_arr = m_tke(lev).array(mfi);
+    const auto& ref_theta_arr = (*m_ref_theta)(lev).const_array(mfi);
     const auto& geom = m_mesh.Geom(lev);
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& probhi = m_mesh.Geom(lev).ProbHiArray();
     const auto& dx = geom.CellSizeArray();
     const auto& dt = m_time.delta_t();
-    const amrex::Real ref_temp = m_ref_temp;
-    const amrex::Real heat_flux =
-        std::abs(m_gravity[2]) / ref_temp * m_heat_flux;
+    const amrex::Real heat_flux = m_heat_flux;
     const amrex::Real Cmu = m_Cmu;
     const amrex::Real sponge_start = m_sponge_start;
     const amrex::Real ref_tke = m_ref_tke;
     const auto tiny = std::numeric_limits<amrex::Real>::epsilon();
     const amrex::Real kappa = m_kappa;
     const amrex::Real z0 = m_z0;
+    const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
+        m_gravity[0], m_gravity[1], m_gravity[2]};
     amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         amrex::Real bcforcing = 0;
         const amrex::Real ux = vel(i, j, k, 0);
@@ -70,9 +72,10 @@ void KransAxell::operator()(
         if (k == 0) {
             const amrex::Real m = std::sqrt(ux * ux + uy * uy);
             const amrex::Real ustar = m * kappa / std::log(z / z0);
+            const amrex::Real T0 = ref_theta_arr(i, j, k);
+            const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux;
             const amrex::Real rans_b = std::pow(
-                std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3),
-                (2.0 / 3.0));
+                std::max(hf, 0.0) * kappa * z / std::pow(Cmu, 3), (2.0 / 3.0));
             bcforcing =
                 (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) / dt;
         }
@@ -105,8 +108,10 @@ void KransAxell::operator()(
                 const amrex::Real z = 0.5 * dx[2];
                 amrex::Real m = std::sqrt(ux * ux + uy * uy);
                 const amrex::Real ustar = m * kappa / std::log(z / z0);
+                const amrex::Real T0 = ref_theta_arr(i, j, k);
+                const amrex::Real hf = std::abs(gravity[2]) / T0 * heat_flux;
                 const amrex::Real rans_b = std::pow(
-                    std::max(heat_flux, 0.0) * kappa * z / std::pow(Cmu, 3),
+                    std::max(hf, 0.0) * kappa * z / std::pow(Cmu, 3),
                     (2.0 / 3.0));
                 terrainforcing =
                     (ustar * ustar / (Cmu * Cmu) + rans_b - tke_arr(i, j, k)) /

amr-wind/transport_models/ConstTransport.H

@@ -198,6 +198,20 @@ public:
         return m_reference_temperature;
     }
 
+    //! Return the reference temperature
+    inline std::unique_ptr<ScratchField> ref_theta() const override
+    {
+        if (m_reference_temperature < 0.0) {
+            amrex::Abort("reference temperature was not set");
+        }
+
+        auto ref_theta = m_repo.create_scratch_field(1, 1);
+        for (int lev = 0; lev < m_repo.num_active_levels(); ++lev) {
+            (*ref_theta)(lev).setVal(m_reference_temperature);
+        }
+        return ref_theta;
+    }
+
 private:
     //! Reference to the field repository (for creating scratch fields)
     FieldRepo& m_repo;

amr-wind/transport_models/TransportModel.H

@@ -48,6 +48,9 @@ public:
 
     //! Reference temperature
     virtual amrex::Real reference_temperature() const = 0;
+
+    //! Reference temperature
+    virtual std::unique_ptr<ScratchField> ref_theta() const = 0;
 };
 } // namespace amr_wind::transport
 

amr-wind/transport_models/TwoPhaseTransport.H

@@ -356,6 +356,14 @@ public:
         return m_reference_temperature;
     }
 
+    //! Return the reference temperature
+    inline std::unique_ptr<ScratchField> ref_theta() const override
+    {
+        amrex::Abort("ref temp not implemented");
+        auto ref_theta = m_repo.create_scratch_field(1, 1);
+        return ref_theta;
+    }
+
 private:
     //! Reference to the CFD sim
     const CFDSim& m_sim;

amr-wind/turbulence/LES/AMD.H

@@ -50,9 +50,6 @@ private:
     //! discretization)
     amrex::Real m_C{0.333333333333333};
 
-    //! Reference temperature (Kelvin)
-    amrex::Real m_ref_theta{300.0};
-
     //! Wall-normal direction axis
     int m_normal_dir{2};
 

amr-wind/turbulence/LES/OneEqKsgs.H

@@ -80,9 +80,6 @@ private:
     //! Gravity vector (m/s^2)
     amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 
-    //! Reference temperature (Kelvin)
-    amrex::Real m_ref_theta{300.0};
-
     //! Hybrid RANS-LES with Nalu-wind
     bool m_hybrid_rl{false};
     Field* m_sdr;