Damping refactor and enhancement

src/equation_systems/CMakeLists.txt

@@ -11,3 +11,4 @@ add_subdirectory(sdr)
 add_subdirectory(levelset)
 add_subdirectory(vof)
 add_subdirectory(passive_scalar)
+add_subdirectory(source_terms)

src/equation_systems/source_terms/CMakeLists.txt

@@ -0,0 +1,3 @@
+target_sources(${kynema_sgf_lib_name} PRIVATE
+  DampingLayerSource.cpp
+  )

src/equation_systems/source_terms/DampingLayerSource.H

@@ -0,0 +1,153 @@
+#ifndef DAMPINGLAYERSOURCE_H
+#define DAMPINGLAYERSOURCE_H
+
+#include <string>
+
+#include "src/CFDSim.H"
+#include "src/core/FieldRepo.H"
+#include "src/equation_systems/density/DensitySource.H"
+#include "src/equation_systems/sdr/SDRSource.H"
+#include "src/equation_systems/SourceTerm.H"
+#include "src/equation_systems/icns/MomentumSource.H"
+#include "src/equation_systems/temperature/TemperatureSource.H"
+#include "src/equation_systems/tke/TKESource.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_REAL.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf::pde {
+
+template <typename T>
+struct DampingLayerSourceTraits
+{
+    static constexpr bool is_specialized = false;
+    static constexpr const char* field_name = "";
+};
+
+template <>
+struct DampingLayerSourceTraits<MomentumSource>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "velocity";
+};
+
+template <>
+struct DampingLayerSourceTraits<TemperatureSource>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "temperature";
+};
+
+template <>
+struct DampingLayerSourceTraits<DensitySource>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "density";
+};
+
+template <>
+struct DampingLayerSourceTraits<TKESource>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "tke";
+};
+
+template <>
+struct DampingLayerSourceTraits<SDRSource>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "sdr";
+};
+
+template <>
+struct DampingLayerSourceTraits<SourceTerm>
+{
+    static constexpr bool is_specialized = true;
+    static constexpr const char* field_name = "";
+};
+
+template <typename T>
+class DampingLayerSource : public T::template Register<DampingLayerSource<T>>
+{
+public:
+    static_assert(
+        DampingLayerSourceTraits<T>::is_specialized,
+        "DampingLayerSourceTraits<T> specialization missing for this source "
+        "type. Add a specialization with field_name.");
+
+    static std::string identifier() { return "DampingLayerSource"; }
+
+    static constexpr const char* damped_field_name()
+    {
+        return DampingLayerSourceTraits<T>::field_name;
+    }
+
+    explicit DampingLayerSource(const CFDSim& sim)
+        : m_repo(sim.repo())
+        , m_damped_field_name(damped_field_name())
+        , m_damped_field(sim.repo().get_field(damped_field_name()))
+    {
+
+        // Insert gets and queries about target field profile, i.e., what to
+        // damp towards
+        amrex::ParmParse pp(identifier() + "." + m_damped_field_name);
+    }
+
+    ~DampingLayerSource() override = default;
+
+    void operator()(
+        int lev, FieldState fstate, amrex::MultiFab& src_term) const override
+    {
+
+        /*const auto& dx = geom.CellSizeArray();
+        const auto& prob_lo = geom.ProbLoArray();
+        const auto& prob_hi = geom.ProbHiArray();*/
+
+        Field* damping_layer_ptr{nullptr};
+
+        auto const& src_arrs = src_term.arrays();
+        auto const& damped_field_arrs =
+            m_damped_field.state(field_impl::dof_state(fstate))(lev)
+                .const_arrays();
+
+        for (int bc_idx = 0; bc_idx < 6; ++bc_idx) {
+            const auto damping_coefficient_field_name =
+                "damping_layer_" + m_damped_field_name + "_" +
+                m_bc_names[bc_idx];
+            if (m_repo.field_exists(damping_coefficient_field_name)) {
+                // Form damping layer field name
+                damping_layer_ptr =
+                    &m_repo.get_field(damping_coefficient_field_name);
+                auto& damping_mfab = (*damping_layer_ptr)(lev);
+                auto const& damping_arrs = damping_mfab.const_arrays();
+
+                amrex::ParallelFor(
+                    damping_mfab, amrex::IntVect(0), m_damped_field.num_comp(),
+                    [=] AMREX_GPU_DEVICE(
+                        int nbx, int i, int j, int k, int n) noexcept {
+                        // Modify src term based on current value of damped
+                        // field, damping coefficient field, target value
+
+                        // Placeholder for now!!!!
+                        const auto target_value = 0.0_rt;
+
+                        src_arrs[nbx](i, j, k, n) +=
+                            damping_arrs[nbx](i, j, k) *
+                            (target_value - damped_field_arrs[nbx](i, j, k, n));
+                    });
+            }
+        }
+    }
+
+private:
+    FieldRepo& m_repo;
+    std::string m_damped_field_name;
+    Field& m_damped_field;
+    const amrex::Vector<std::string> m_bc_names = {"xlo", "xhi", "ylo",
+                                                   "yhi", "zlo", "zhi"};
+};
+
+} // namespace kynema_sgf::pde
+
+#endif /* DAMPINGLAYERSOURCE_H */

src/equation_systems/source_terms/DampingLayerSource.cpp

@@ -0,0 +1,12 @@
+#include "src/equation_systems/source_terms/DampingLayerSource.H"
+
+namespace kynema_sgf::pde {
+
+template class DampingLayerSource<MomentumSource>;
+template class DampingLayerSource<TemperatureSource>;
+template class DampingLayerSource<DensitySource>;
+template class DampingLayerSource<TKESource>;
+template class DampingLayerSource<SDRSource>;
+template class DampingLayerSource<SourceTerm>;
+
+} // namespace kynema_sgf::pde

src/physics/CMakeLists.txt

@@ -20,6 +20,7 @@ target_sources(${kynema_sgf_lib_name}
   ActuatorSourceTagging.cpp
   Intermittency.cpp
   ForestDrag.cpp
+  DampingLayer.cpp
   )
 
 add_subdirectory(multiphase)

src/physics/DampingLayer.H

@@ -0,0 +1,120 @@
+#ifndef DAMPINGLAYER_H
+#define DAMPINGLAYER_H
+
+#include <numbers>
+#include "src/core/Physics.H"
+#include "src/core/Field.H"
+#include "src/CFDSim.H"
+#include "AMReX_Gpu.H"
+
+namespace kynema_sgf::damping_layer {
+
+enum class BlendingFunctionType : std::uint8_t {
+    Linear = 0,
+    Quadratic = 1,
+    Exponential = 2,
+    Cosine = 3
+};
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BlendingFunctionType
+string_to_blending_function_type(const std::string& op_str)
+{
+    std::string op = amrex::toLower(op_str);
+    if (op == "linear") {
+        return BlendingFunctionType::Linear;
+    }
+    if (op == "quadratic") {
+        return BlendingFunctionType::Quadratic;
+    }
+    if (op == "exponential") {
+        return BlendingFunctionType::Exponential;
+    }
+    if (op == "cosine") {
+        return BlendingFunctionType::Cosine;
+    }
+    amrex::Abort(
+        "Invalid blending function type: " + op_str +
+        ". Must be one of: linear, quadratic, exponential, cosine.");
+    // Unreachable, but avoids compiler warning
+    return BlendingFunctionType::Cosine;
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real blending_function(
+    amrex::Real position_nondim, BlendingFunctionType function_type)
+{
+    switch (function_type) {
+    case BlendingFunctionType::Linear:
+        return 1.0_rt - position_nondim;
+    case BlendingFunctionType::Quadratic:
+        return 1.0_rt - position_nondim * position_nondim;
+    case BlendingFunctionType::Exponential: // exponential difference and ratio
+        return 1.0_rt - (std::expm1(std::pow(position_nondim, 3.5_rt)) /
+                         std::expm1(1.0_rt));
+    case BlendingFunctionType::Cosine:
+    default:
+        return 0.5_rt +
+               0.5_rt *
+                   std::cos(std::numbers::pi_v<amrex::Real> * position_nondim);
+    }
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real damping_calc(
+    amrex::Real distance_from_boundary,
+    amrex::Real layer_thickness,
+    amrex::Real blending_fraction,
+    BlendingFunctionType function_type)
+{
+    const amrex::Real length_full_damp =
+        layer_thickness * (1.0_rt - blending_fraction);
+    const amrex::Real blend_position_nondim =
+        (distance_from_boundary - length_full_damp) /
+        (layer_thickness * blending_fraction);
+    if (blend_position_nondim <= 0.0_rt) {
+        return 1.0_rt;
+    }
+    if (blend_position_nondim <= 1.0_rt) {
+        return blending_function(blend_position_nondim, function_type);
+    }
+    return 0.0_rt;
+}
+
+class DampingLayer : public Physics::Register<DampingLayer>
+{
+public:
+    static std::string identifier() { return "DampingLayer"; }
+
+    explicit DampingLayer(CFDSim& sim);
+
+    ~DampingLayer() override = default;
+
+    void initialize_fields(int level, const amrex::Geometry& geom) override;
+
+    void pre_init_actions() override {}
+
+    void post_init_actions() override {}
+
+    void post_regrid_actions() override;
+
+    void pre_advance_work() override {}
+
+    void post_advance_work() override {}
+
+private:
+    FieldRepo& m_repo;
+    const amrex::AmrCore& m_mesh;
+
+    const amrex::Vector<std::string> m_bc_names = {"xlo", "xhi", "ylo",
+                                                   "yhi", "zlo", "zhi"};
+
+    //! Quantities for constructing damping layers
+    amrex::Vector<std::string> m_field_names;
+    amrex::Gpu::DeviceVector<amrex::Array<amrex::Real, 6>> m_layers_thickness;
+    amrex::Gpu::DeviceVector<amrex::Array<amrex::Real, 6>>
+        m_layers_blending_fraction;
+    amrex::Gpu::DeviceVector<amrex::Array<amrex::Real, 4>> m_layers_min_height;
+    amrex::Gpu::DeviceVector<amrex::Array<BlendingFunctionType, 6>>
+        m_layers_blending_function_type;
+};
+} // namespace kynema_sgf::damping_layer
+
+#endif