Thermal expansion coefficient and ref temp consolidation (#1365)

* Thermal expansion coeffcient consolidation

* fix shadow

* try this

* or this

* fix

* fix unit tests

* beta in turb models

* Ref temp return

* more beta

* format

* try to fix amd

* Move vof to const transport

* promote transport

* typo

* fix warns

* fix conditionals for multiphase

* Add warning

* fix unit test

* fix comment

* big update for ref temp

* tweak for multiphase

* ugh redo a bunch

* beta impl

* same for 2 phase

* fix naming

* documentation update

* input file update

* more input params

* unit tests

* format

* fix warn

amr-wind/CFDSim.H

@@ -20,6 +20,10 @@ class OversetManager;
 class ExtSolverMgr;
 class HelicsStorage;
 
+namespace transport {
+class TransportModel;
+}
+
 namespace turbulence {
 class TurbulenceModel;
 }
@@ -74,6 +78,12 @@ public:
         return m_physics_mgr.objects();
     }
 
+    transport::TransportModel& transport_model() { return *m_transport; }
+    const transport::TransportModel& transport_model() const
+    {
+        return *m_transport;
+    }
+
     turbulence::TurbulenceModel& turbulence_model() { return *m_turbulence; }
     const turbulence::TurbulenceModel& turbulence_model() const
     {
@@ -103,6 +113,9 @@ public:
 
     bool has_overset() const;
 
+    //! Instantiate the transport model based on user inputs
+    void create_transport_model();
+
     //! Instantiate the turbulence model based on user inputs
     void create_turbulence_model();
 
@@ -137,6 +150,8 @@ private:
 
     PhysicsMgr m_physics_mgr;
 
+    std::unique_ptr<transport::TransportModel> m_transport;
+
     std::unique_ptr<turbulence::TurbulenceModel> m_turbulence;
 
     std::unique_ptr<IOManager> m_io_mgr;

amr-wind/CFDSim.cpp

@@ -1,4 +1,5 @@
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/turbulence/TurbulenceModel.H"
 #include "amr-wind/utilities/IOManager.H"
 #include "amr-wind/utilities/PostProcessing.H"
@@ -22,6 +23,17 @@ CFDSim::CFDSim(amrex::AmrCore& mesh)
 
 CFDSim::~CFDSim() = default;
 
+void CFDSim::create_transport_model()
+{
+    std::string transport_model_name = "ConstTransport";
+    {
+        amrex::ParmParse pp("transport");
+        pp.query("model", transport_model_name);
+    }
+    m_transport =
+        transport::TransportModel::create(transport_model_name, *this);
+}
+
 void CFDSim::create_turbulence_model()
 {
     std::string transport_model_name = "ConstTransport";

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

@@ -2,6 +2,7 @@
 #define ABLMEANBOUSSINESQ_H
 
 #include "amr-wind/core/FieldRepo.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/equation_systems/icns/MomentumSource.H"
 #include "amr-wind/utilities/FieldPlaneAveraging.H"
 
@@ -42,11 +43,11 @@ 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
-    amrex::Real m_beta{0.0};
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 
     int m_axis{2};
 

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

@@ -8,50 +8,74 @@
 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())
-{
+ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim)
+    : m_mesh(sim.mesh()), m_transport(sim.transport_model())
 
+{
     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);
-
-    if (pp_boussinesq_buoyancy.contains("thermal_expansion_coeff")) {
-        pp_boussinesq_buoyancy.get("thermal_expansion_coeff", m_beta);
-    } else {
-        m_beta = 1.0 / m_ref_theta;
-    }
+    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());
     }
 }
@@ -60,15 +84,19 @@ ABLMeanBoussinesq::~ABLMeanBoussinesq() = default;
 
 void ABLMeanBoussinesq::operator()(
     const int lev,
-    const amrex::MFIter& /*mfi*/,
+    const amrex::MFIter& mfi,
     const amrex::Box& bx,
     const FieldState /*fstate*/,
     const amrex::Array4<amrex::Real>& src_term) const
 {
     const auto& problo = m_mesh.Geom(lev).ProbLoArray();
     const auto& dx = m_mesh.Geom(lev).CellSizeArray();
-    const amrex::Real T0 = m_ref_theta;
-    const amrex::Real beta = m_beta;
+
+    amrex::FArrayBox beta_fab(bx, 1, amrex::The_Async_Arena());
+    amrex::Array4<amrex::Real> const& beta_arr = beta_fab.array();
+    m_transport.beta_impl(lev, mfi, bx, beta_arr);
+
+    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]};
 
@@ -82,9 +110,10 @@ 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 * (temp - T0);
+        const amrex::Real fac = beta_arr(i, j, k) * (temp - T0);
         src_term(i, j, k, 0) += gravity[0] * fac;
         src_term(i, j, k, 1) += gravity[1] * fac;
         src_term(i, j, k, 2) += gravity[2] * fac;

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

@@ -2,6 +2,7 @@
 #define BOUSSINESQBUOYANCY_H
 
 #include "amr-wind/core/FieldRepo.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/equation_systems/icns/MomentumSource.H"
 
 namespace amr_wind::pde::icns {
@@ -31,18 +32,14 @@ public:
 
 private:
     const Field& m_temperature;
-    const Field* m_vof;
 
     amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 
-    //! Reference temperature (Kelvin)
-    amrex::Real m_ref_theta{300.0};
+    //! Transport model
+    const transport::TransportModel& m_transport;
 
-    //! Thermal expansion coefficient
-    amrex::Real m_beta{0.0};
-
-    //! Check for VOF
-    bool m_is_vof{false};
+    //! Reference temperature
+    std::unique_ptr<ScratchField> m_ref_theta;
 };
 } // namespace amr_wind::pde::icns
 

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

@@ -8,32 +8,12 @@ namespace amr_wind::pde::icns {
 
 /** Boussinesq buoyancy source term for ABL simulations
  *
- *  Reads in the following parameters from `BoussinesqBuoyancy` namespace:
- *
- *  - `reference_temperature` (Mandatory) temperature (`T0`) in Kelvin
- *  - `thermal_expansion_coeff` Optional, default = `1.0 / T0`
- *  - `gravity` acceleration due to gravity (m/s)
  */
 BoussinesqBuoyancy::BoussinesqBuoyancy(const CFDSim& sim)
     : m_temperature(sim.repo().get_field("temperature"))
+    , m_transport(sim.transport_model())
 {
-    amrex::ParmParse pp_boussinesq_buoyancy(
-        amr_wind::pde::icns::BoussinesqBuoyancy::identifier());
-    pp_boussinesq_buoyancy.get("reference_temperature", m_ref_theta);
-
-    if (pp_boussinesq_buoyancy.contains("thermal_expansion_coeff")) {
-        pp_boussinesq_buoyancy.get("thermal_expansion_coeff", m_beta);
-    } else {
-        m_beta = 1.0 / m_ref_theta;
-    }
-
-    m_is_vof = sim.repo().field_exists("vof");
-    if (m_is_vof) {
-        m_vof = &sim.repo().get_field("vof");
-    } else {
-        // Point to something, will not be used
-        m_vof = &m_temperature;
-    }
+    m_ref_theta = m_transport.ref_theta();
 
     // gravity in `incflo` namespace
     amrex::ParmParse pp_incflo("incflo");
@@ -49,25 +29,20 @@ void BoussinesqBuoyancy::operator()(
     const FieldState fstate,
     const amrex::Array4<amrex::Real>& src_term) const
 {
-    const amrex::Real T0 = m_ref_theta;
-    const amrex::Real beta = m_beta;
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
 
-    const bool ivf = m_is_vof;
-    const auto& vof_arr = (*m_vof)(lev).const_array(mfi);
-    constexpr amrex::Real tol = 1e-12;
-
     const auto& temp =
         m_temperature.state(field_impl::phi_state(fstate))(lev).const_array(
             mfi);
+    amrex::FArrayBox beta_fab(bx, 1, amrex::The_Async_Arena());
+    amrex::Array4<amrex::Real> const& beta_arr = beta_fab.array();
+    m_transport.beta_impl(lev, mfi, bx, beta_arr);
+    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 fac_air = beta * (T0 - T);
-        // If vof exists, ignore Boussinesq term in cells with liquid
-        // If no vof, assume single phase and use the term for air everywhere
-        const amrex::Real fac =
-            ivf ? (vof_arr(i, j, k) > tol ? 0.0 : fac_air) : fac_air;
+        const amrex::Real T0 = ref_theta(i, j, k);
+        const amrex::Real fac = beta_arr(i, j, k) * (T0 - T);
 
         src_term(i, j, k, 0) += gravity[0] * fac;
         src_term(i, j, k, 1) += gravity[1] * 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,16 @@ 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();
+    if (pp.contains("reference_temperature")) {
+        amrex::Abort(
+            "DragTempForcing.reference_temperature has been deprecated. Please "
+            "replace with transport.reference_temperature.");
+    }
 }
 
 DragTempForcing::~DragTempForcing() = default;
@@ -44,7 +50,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 +60,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));
     });
 }