Compatible with mesoscale zone

amr-wind/equation_systems/icns/source_terms/CMakeLists.txt

@@ -16,4 +16,5 @@ target_sources(${amr_wind_lib_name} PRIVATE
   NonLinearSGSTerm.cpp
   DragForcing.cpp
   ForestForcing.cpp
+  WindSpongeForcing.cpp
   )

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

@@ -118,28 +118,28 @@ void DragForcing::operator()(
         yi_end = sponge_north * std::max(yi_end, 0.0);
         ystart_damping = sponge_strength * yi_start * yi_start;
         yend_damping = sponge_strength * yi_end * yi_end;
-
+        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 auto idx =
             interp::bisection_search(device_vel_ht, device_vel_ht + vsize, z);
         const amrex::Real spongeVelX =
             (vsize > 0) ? interp::linear_impl(
                               device_vel_ht, device_vel_vals, z, idx, 3, 0)
-                        : 0.0;
+                        : ux1;
         const amrex::Real spongeVelY =
             (vsize > 0) ? interp::linear_impl(
                               device_vel_ht, device_vel_vals, z, idx, 3, 1)
-                        : 0.0;
+                        : uy1;
         const amrex::Real spongeVelZ =
             (vsize > 0) ? interp::linear_impl(
                               device_vel_ht, device_vel_vals, z, idx, 3, 2)
-                        : 0.0;
+                        : uz1;
         amrex::Real Dxz = 0.0;
         amrex::Real Dyz = 0.0;
         amrex::Real bc_forcing_x = 0;
         amrex::Real bc_forcing_y = 0;
-        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 m = std::sqrt(ux1 * ux1 + uy1 * uy1 + uz1 * uz1);
         if (drag(i, j, k) == 1 && (!is_laminar)) {
             const amrex::Real ux2 = vel(i, j, k + 1, 0);

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

@@ -0,0 +1,45 @@
+#ifndef WINDSPONGEFORCING_H
+#define WINDSPONGEFORCING_H
+
+#include "amr-wind/equation_systems/icns/MomentumSource.H"
+#include "amr-wind/core/SimTime.H"
+#include "amr-wind/CFDSim.H"
+
+namespace amr_wind::pde::icns {
+
+class WindSpongeForcing : public MomentumSource::Register<WindSpongeForcing>
+{
+public:
+    static std::string identifier() { return "WindSpongeForcing"; }
+
+    explicit WindSpongeForcing(const CFDSim& sim);
+
+    ~WindSpongeForcing() override;
+
+    void operator()(
+        const int lev,
+        const amrex::MFIter& mfi,
+        const amrex::Box& bx,
+        const FieldState /*fstate*/,
+        const amrex::Array4<amrex::Real>& src_term) const override;
+
+private:
+    const SimTime& m_time;
+    const amrex::AmrCore& m_mesh;
+    std::string m_1d_rans;
+    const Field& m_velocity;
+    amrex::Vector<amrex::Real> m_wind_heights;
+    amrex::Vector<amrex::Real> m_u_values;
+    amrex::Vector<amrex::Real> m_v_values;
+    amrex::Vector<amrex::Real> m_w_values;
+    amrex::Gpu::DeviceVector<amrex::Real> m_wind_heights_d;
+    amrex::Gpu::DeviceVector<amrex::Real> m_u_values_d;
+    amrex::Gpu::DeviceVector<amrex::Real> m_v_values_d;
+    amrex::Gpu::DeviceVector<amrex::Real> m_w_values_d;
+    amrex::Real m_meso_start{600};
+    amrex::Real m_meso_timescale{30};
+    const CFDSim& m_sim;
+};
+
+} // namespace amr_wind::pde::icns
+#endif

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

@@ -0,0 +1,153 @@
+
+#include "amr-wind/equation_systems/icns/source_terms/WindSpongeForcing.H"
+#include "amr-wind/utilities/IOManager.H"
+#include "amr-wind/utilities/linear_interpolation.H"
+
+#include "AMReX_ParmParse.H"
+#include "AMReX_Gpu.H"
+#include "AMReX_Random.H"
+
+namespace amr_wind::pde::icns {
+
+WindSpongeForcing::WindSpongeForcing(const CFDSim& sim)
+    : m_time(sim.time())
+    , m_mesh(sim.mesh())
+    , m_velocity(sim.repo().get_field("velocity"))
+    , m_sim(sim)
+{
+    amrex::ParmParse pp_abl("ABL");
+    pp_abl.query("rans_1dprofile_file", m_1d_rans);
+    if (!m_1d_rans.empty()) {
+        std::ifstream ransfile(m_1d_rans, std::ios::in);
+        if (!ransfile.good()) {
+            amrex::Abort("Cannot find 1-D RANS profile file " + m_1d_rans);
+        }
+        amrex::Real value1, value2, value3, value4, value5;
+        while (ransfile >> value1 >> value2 >> value3 >> value4 >> value5) {
+            m_wind_heights.push_back(value1);
+            m_u_values.push_back(value2);
+            m_v_values.push_back(value3);
+            m_w_values.push_back(value4);
+        }
+    } else {
+        amrex::Abort("Cannot find 1-D RANS profile file " + m_1d_rans);
+    }
+    pp_abl.query("meso_sponge_start", m_meso_start);
+    pp_abl.query("meso_timescale", m_meso_timescale);
+    int num_wind_values = static_cast<int>(m_wind_heights.size());
+    m_wind_heights_d.resize(num_wind_values);
+    m_u_values_d.resize(num_wind_values);
+    m_v_values_d.resize(num_wind_values);
+    m_w_values_d.resize(num_wind_values);
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_wind_heights.begin(), m_wind_heights.end(),
+        m_wind_heights_d.begin());
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_u_values.begin(), m_u_values.end(),
+        m_u_values_d.begin());
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_v_values.begin(), m_v_values.end(),
+        m_v_values_d.begin());
+    amrex::Gpu::copy(
+        amrex::Gpu::hostToDevice, m_w_values.begin(), m_w_values.end(),
+        m_w_values_d.begin());
+}
+
+WindSpongeForcing::~WindSpongeForcing() = default;
+
+void WindSpongeForcing::operator()(
+    const int lev,
+    const amrex::MFIter& mfi,
+    const amrex::Box& bx,
+    const FieldState fstate,
+    const amrex::Array4<amrex::Real>& src_term) const
+{
+    const auto& geom = m_mesh.Geom(lev);
+    const auto& dx = geom.CellSizeArray();
+    const auto& prob_lo = geom.ProbLoArray();
+    const auto& prob_hi = geom.ProbHiArray();
+    const auto& velocity =
+        m_velocity.state(field_impl::dof_state(fstate))(lev).const_array(mfi);
+    const amrex::Real sponge_start = m_meso_start;
+    const amrex::Real meso_timescale = m_time.delta_t();
+    const auto vsize = m_wind_heights_d.size();
+    const auto* wind_heights_d = m_wind_heights_d.data();
+    const auto* u_values_d = m_u_values_d.data();
+    const auto* v_values_d = m_v_values_d.data();
+    const auto* w_values_d = m_w_values_d.data();
+    const bool has_terrain = this->m_sim.repo().field_exists("terrain_height");
+    if (has_terrain) {
+        auto* const m_terrain_height =
+            &this->m_sim.repo().get_field("terrain_height");
+        const auto& terrain_height = (*m_terrain_height)(lev).const_array(mfi);
+        amrex::ParallelFor(
+            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real z = std::max(
+                    prob_lo[2] + (k + 0.5) * dx[2] - terrain_height(i, j, k),
+                    0.5 * dx[2]);
+                const amrex::Real zi = std::max(
+                    (z - sponge_start) / (prob_hi[2] - sponge_start), 0.0);
+                amrex::Real ref_windx = velocity(i, j, k, 0);
+                amrex::Real ref_windy = velocity(i, j, k, 1);
+                amrex::Real ref_windz = velocity(i, j, k, 2);
+                if (zi > 0) {
+                    ref_windx = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, u_values_d, z)
+                                    : velocity(i, j, k, 0);
+                    ref_windy = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, v_values_d, z)
+                                    : velocity(i, j, k, 1);
+                    ref_windz = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, w_values_d, z)
+                                    : velocity(i, j, k, 2);
+                }
+                src_term(i, j, k, 0) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 0) - ref_windx);
+                src_term(i, j, k, 1) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 1) - ref_windy);
+                src_term(i, j, k, 2) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 2) - ref_windz);
+            });
+    } else {
+        amrex::ParallelFor(
+            bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
+                const amrex::Real zi = std::max(
+                    (z - sponge_start) / (prob_hi[2] - sponge_start), 0.0);
+                amrex::Real ref_windx = velocity(i, j, k, 0);
+                amrex::Real ref_windy = velocity(i, j, k, 1);
+                amrex::Real ref_windz = velocity(i, j, k, 2);
+                if (zi > 0) {
+                    ref_windx = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, u_values_d, z)
+                                    : velocity(i, j, k, 0);
+                    ref_windy = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, v_values_d, z)
+                                    : velocity(i, j, k, 1);
+                    ref_windz = (vsize > 0)
+                                    ? interp::linear(
+                                          wind_heights_d,
+                                          wind_heights_d + vsize, w_values_d, z)
+                                    : velocity(i, j, k, 2);
+                }
+                src_term(i, j, k, 0) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 0) - ref_windx);
+                src_term(i, j, k, 1) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 1) - ref_windy);
+                src_term(i, j, k, 2) -=
+                    1.0 / meso_timescale * (velocity(i, j, k, 2) - ref_windz);
+            });
+    }
+}
+
+} // namespace amr_wind::pde::icns

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

@@ -30,7 +30,20 @@ private:
     amrex::Vector<amrex::Real> m_theta_values;
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_heights_d;
     amrex::Gpu::DeviceVector<amrex::Real> m_theta_values_d;
-    amrex::Real m_sponge_start{600};
+    amrex::Real m_meso_start{600};
+    amrex::Real m_meso_timescale{30};
+    bool m_horizontal_sponge{false};
+    amrex::Real m_sponge_strength{1.0};
+    amrex::Real m_sponge_density{1.0};
+    amrex::Real m_sponge_distance_west{-1000};
+    amrex::Real m_sponge_distance_east{1000};
+    amrex::Real m_sponge_distance_south{-1000};
+    amrex::Real m_sponge_distance_north{1000};
+    int m_sponge_west{0};
+    int m_sponge_east{1};
+    int m_sponge_south{0};
+    int m_sponge_north{1};
+    const CFDSim& m_sim;
 };
 
 } // namespace amr_wind::pde::temperature