Temperature support for non-neutral on terrain

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

@@ -49,7 +49,7 @@ private:
     int m_sponge_north{1};
     bool m_is_laminar{false};
     std::string m_wall_het_model{"none"};
-    amrex::Real m_mol_length{10000};
+    amrex::Real m_mol_length{-1e30};
     amrex::Real m_z0{0.1};
     amrex::Real m_kappa{0.41};
     amrex::Real m_gamma_m{5.0};

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

@@ -4,7 +4,6 @@
 #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 {
 
@@ -25,17 +24,54 @@ public:
         const amrex::Array4<amrex::Real>& src_term) const override;
 
 private:
+    const SimTime& m_time;
     const CFDSim& m_sim;
     const amrex::AmrCore& m_mesh;
     const Field& m_velocity;
     const Field& m_temperature;
+    amrex::Real m_internal_temperature{300};
     amrex::Real m_drag_coefficient{1.0};
+    std::string m_wall_het_model{"none"};
+    amrex::Real m_mol_length{-1e30};
+    amrex::Real m_z0{0.1};
+    amrex::Real m_kappa{0.41};
+    amrex::Real m_gamma_m{5.0};
+    amrex::Real m_beta_m{16.0};
 
-    //! Transport model
-    const transport::TransportModel& m_transport;
+    static amrex::Real stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_m,
+        const amrex::Real beta_m)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_m = 0;
+        if (zeta > 0) {
+            psi_m = -gamma_m * zeta;
+        } else {
+            amrex::Real x = std::sqrt(std::sqrt(1 - beta_m * zeta));
+            psi_m = 2.0 * std::log(0.5 * (1.0 + x)) + log(0.5 * (1 + x * x)) -
+                    2.0 * std::atan(x) + 2.0 * std::atan(1.0);
+        }
+        return psi_m;
+    }
 
-    //! Reference temperature
-    std::unique_ptr<ScratchField> m_ref_theta;
+    static amrex::Real thermal_stability(
+        const amrex::Real z,
+        const amrex::Real mol_length,
+        const amrex::Real gamma_h,
+        const amrex::Real beta_h)
+    {
+        const amrex::Real zeta = z / mol_length;
+        amrex::Real psi_m = 0;
+        if (zeta > 0) {
+            psi_h = -gamma_h * zeta;
+        } else {
+            amrex::Real x = std::sqrt(1 - beta_h * zeta);
+            psi_h = 2.0 * std::log(0.5 * (1 + x));
+        }
+        return psi_h;
+    }
 };
 
 } // namespace amr_wind::pde::temperature

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

@@ -9,20 +9,23 @@
 namespace amr_wind::pde::temperature {
 
 DragTempForcing::DragTempForcing(const CFDSim& sim)
-    : m_sim(sim)
+    : m_time(sim.time())
+    , m_sim(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);
-    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.");
-    }
+    pp.query("internal_temperature", m_internal_temperature);
+    amrex::ParmParse pp_abl("ABL");
+    pp_abl.query("wall_het_model", m_wall_het_model);
+    pp_abl.query("mol_length", m_mol_length);
+    pp_abl.query("surface_roughness_z0", m_z0);
+    pp_abl.query("kappa", m_kappa);
+    pp_abl.query("mo_gamma_m", m_gamma_m);
+    pp_abl.query("mo_beta_m", m_beta_m);
 }
 
 DragTempForcing::~DragTempForcing() = default;
@@ -47,22 +50,55 @@ void DragTempForcing::operator()(
     auto* const m_terrain_blank =
         &this->m_sim.repo().get_int_field("terrain_blank");
     const auto& blank = (*m_terrain_blank)(lev).const_array(mfi);
+    auto* const m_terrain_drag =
+        &this->m_sim.repo().get_int_field("terrain_drag");
+    const auto& drag = (*m_terrain_drag)(lev).const_array(mfi);
     const auto& geom = m_mesh.Geom(lev);
     const auto& dx = geom.CellSizeArray();
     const amrex::Real drag_coefficient = m_drag_coefficient / dx[2];
-    const auto& ref_theta = (*m_ref_theta)(lev).const_array(mfi);
+    const amrex::Real internal_temperature = m_internal_temperature;
+    const amrex::Real gravity_mod = 9.81;
+    const amrex::Real psi_m = 0.0;
+    const amrex::Real psi_h_neighbour = 0.0;
+    const amrex::Real psi_h_cell = 0.0;
+    const amrex::Real m_kappa = kappa;
+    const amrex::Real z0 = m_z0;
+    const amrex::Real mol_length = m_mol_length;
+    const auto& dt = m_time.delta_t();
+    if (m_wall_het_model == "mol") {
+        psi_m = stability(1.5 * dx[2], m_mol_length, m_gamma_m, m_beta_m);
+        psi_h_neighbour =
+            thermal_stability(1.5 * dx[2], m_mol_length, m_gamma_h, m_beta_h);
+        psi_h_cell =
+            thermal_stability(0.5 * dx[2], m_mol_length, m_gamma_h, m_beta_h);
+    }
     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 {
         const amrex::Real ux1 = vel(i, j, k, 0);
         const amrex::Real uy1 = vel(i, j, k, 1);
         const amrex::Real uz1 = vel(i, j, k, 2);
+        const amrex::Real theta = temperature(i, j, k, 0);
+        const amrex::Real theta2 = temperature(i, j, k + 1, 0);
+        const amrex::Real wspd = std::sqrt(ux1 * ux1 + uy1 * uy1);
+        const amrex::Real ustar =
+            wspd * kappa / (std::log(1.5 * dx[2] / z0) - psi_m);
+        //! We do not know the actual temperature so use cell above
+        const amrex::Real thetastar =
+            theta * ustar * ustar / (kappa * gravity_mod * mol_length);
+        const amrex::Real surf_temp =
+            theta2 -
+            thetastar / kappa * (std::log(1.5 * dx[2] / z0) - psi_h_neighbour);
+        const amrex::Real tTarget =
+            surf_temp +
+            thetastar / kappa * (std::log(0.5 * dx[2] / z0) - psi_h_cell);
+        const amrex::Real bc_forcing_t = -(tTarget - theta) / dt;
         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) - T0) * blank(i, j, k, 0));
+            (Cd * (theta - internal_temperature) * blank(i, j, k, 0) +
+             bc_forcing_t * drag(i, j, k));
     });
 }
 

amr-wind/wind_energy/ABLWallFunction.cpp

@@ -431,7 +431,7 @@ void ABLTempWallFunc::wall_model(
                                 -theta2 * ustar * ustar /
                                 (kappa * gravity_mod * mol_length);
                             const amrex::Real surf_temp_flux =
-                                ustar * thetastar;
+                                -ustar * thetastar;
                             const amrex::Real surf_temp =
                                 surf_temp_flux * (std::log(z / z0) - psi_h) /
                                     (ustar * kappa) +