Flather BC

src/boundary_conditions/field_boundary_fill/BoundaryPlane.cpp

@@ -499,8 +499,7 @@ void BoundaryPlane::initialize_data()
         if (m_repo.field_exists(fname)) {
             auto& fld = m_repo.get_field(fname);
             if (m_io_mode == io_mode::input) {
-                fld.register_fill_patch_op<PlaneFillInflow>(
-                    m_mesh, m_time, *this);
+                fld.add_fill_patch_op<PlaneFillInflow>(m_mesh, m_time, *this);
             }
             m_fields.emplace_back(&fld);
         } else {

src/boundary_conditions/field_boundary_fill/CMakeLists.txt

@@ -4,6 +4,8 @@ target_sources(${kynema_sgf_lib_name}
   OceanWavesBoundary.cpp
   PlaneFillInflow.cpp
   BoundaryPlane.cpp
+  FillFlather.cpp
+  Flather.cpp
   FillMPL.cpp
   ModulatedPowerLaw.cpp
   )

src/boundary_conditions/field_boundary_fill/FillFlather.H

@@ -0,0 +1,60 @@
+#ifndef FILLFLATHER_H
+#define FILLFLATHER_H
+
+#include "src/core/FieldFillPatchOps.H"
+#include "src/core/FieldBCOps.H"
+#include "src/boundary_conditions/field_boundary_fill/Flather.H"
+
+namespace kynema_sgf {
+
+/** Fill patch operator that applies Flather boundary updates. */
+class FillFlather : public FieldFillPatchOps<FieldBCDirichlet>
+{
+public:
+    FillFlather(
+        Field& field,
+        const amrex::AmrCore& mesh,
+        const SimTime& time,
+        Flather& flather);
+
+    ~FillFlather() override;
+
+    void fillpatch(
+        int lev,
+        amrex::Real time,
+        amrex::MultiFab& mfab,
+        const amrex::IntVect& nghost,
+        FieldState fstate = FieldState::New) override;
+
+    void fillpatch_sibling_fields(
+        int lev,
+        amrex::Real time,
+        amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& mfabs,
+        amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& ffabs,
+        amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& cfabs,
+        const amrex::IntVect& nghost,
+        const amrex::Vector<amrex::BCRec>& bcrec,
+        const amrex::Vector<amrex::BCRec>& /* unused */,
+        FieldState fstate = FieldState::New) override;
+
+    void fillpatch_from_coarse(
+        int lev,
+        amrex::Real time,
+        amrex::MultiFab& mfab,
+        const amrex::IntVect& nghost,
+        FieldState fstate = FieldState::New) override;
+
+    void fillphysbc(
+        int lev,
+        amrex::Real time,
+        amrex::MultiFab& mfab,
+        const amrex::IntVect& nghost,
+        FieldState fstate = FieldState::New) override;
+
+protected:
+    Flather& m_flather;
+};
+
+} // namespace kynema_sgf
+
+#endif /* FILLFLATHER_H */

src/boundary_conditions/field_boundary_fill/FillFlather.cpp

@@ -0,0 +1,75 @@
+#include <utility>
+
+#include "src/boundary_conditions/field_boundary_fill/FillFlather.H"
+
+namespace kynema_sgf {
+
+FillFlather::FillFlather(
+    Field& field,
+    const amrex::AmrCore& mesh,
+    const SimTime& time,
+    Flather& flather)
+    : FieldFillPatchOps<FieldBCDirichlet>(
+          field, mesh, time, FieldInterpolator::CellConsLinear)
+    , m_flather(flather)
+{}
+
+FillFlather::~FillFlather() = default;
+
+void FillFlather::fillpatch(
+    const int lev,
+    const amrex::Real time,
+    amrex::MultiFab& mfab,
+    const amrex::IntVect& /* nghost */,
+    const FieldState fstate)
+{
+    if (m_field.base_name() == "velocity") {
+        m_flather.update_flather_variables(lev, fstate);
+        m_flather.set_velocity(lev, time, m_field, mfab);
+    }
+}
+
+void FillFlather::fillpatch_from_coarse(
+    const int lev,
+    const amrex::Real time,
+    amrex::MultiFab& mfab,
+    const amrex::IntVect& /* nghost */,
+    const FieldState fstate)
+{
+    if (m_field.base_name() == "velocity") {
+        m_flather.update_flather_variables(lev, fstate);
+        m_flather.set_velocity(lev, time, m_field, mfab);
+    }
+}
+
+void FillFlather::fillphysbc(
+    const int lev,
+    const amrex::Real time,
+    amrex::MultiFab& mfab,
+    const amrex::IntVect& /* nghost */,
+    const FieldState fstate)
+{
+    if (m_field.base_name() == "velocity") {
+        m_flather.update_flather_variables(lev, fstate);
+        m_flather.set_velocity(lev, time, m_field, mfab);
+    }
+}
+
+void FillFlather::fillpatch_sibling_fields(
+    const int lev,
+    const amrex::Real time,
+    amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& mfabs,
+    amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& /* ffabs */,
+    amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& /* cfabs */,
+    const amrex::IntVect& /* nghost */,
+    const amrex::Vector<amrex::BCRec>& /* bcrec */,
+    const amrex::Vector<amrex::BCRec>& /* unused */,
+    const FieldState fstate)
+{
+    m_flather.update_flather_variables(lev, fstate, true);
+    for (int i = 0; std::cmp_less(i, mfabs.size()); ++i) {
+        m_flather.set_velocity(lev, time, m_field, *mfabs[i], 0, i);
+    }
+}
+
+} // namespace kynema_sgf

src/boundary_conditions/field_boundary_fill/FillMPL.cpp

@@ -20,12 +20,9 @@ void FillMPL::fillpatch(
     const int lev,
     const amrex::Real time,
     amrex::MultiFab& mfab,
-    const amrex::IntVect& nghost,
-    const FieldState fstate)
+    const amrex::IntVect& /* nghost */,
+    const FieldState /* fstate */)
 {
-    FieldFillPatchOps<FieldBCDirichlet>::fillpatch(
-        lev, time, mfab, nghost, fstate);
-
     if (m_field.base_name() == "velocity") {
         m_abl_mpl.set_velocity(lev, time, m_field, mfab);
     } else if (m_field.base_name() == "temperature") {
@@ -37,12 +34,9 @@ void FillMPL::fillpatch_from_coarse(
     const int lev,
     const amrex::Real time,
     amrex::MultiFab& mfab,
-    const amrex::IntVect& nghost,
-    const FieldState fstate)
+    const amrex::IntVect& /* nghost */,
+    const FieldState /* fstate */)
 {
-    FieldFillPatchOps<FieldBCDirichlet>::fillpatch_from_coarse(
-        lev, time, mfab, nghost, fstate);
-
     if (m_field.base_name() == "velocity") {
         m_abl_mpl.set_velocity(lev, time, m_field, mfab);
     } else if (m_field.base_name() == "temperature") {
@@ -54,12 +48,9 @@ void FillMPL::fillphysbc(
     const int lev,
     const amrex::Real time,
     amrex::MultiFab& mfab,
-    const amrex::IntVect& nghost,
-    const FieldState fstate)
+    const amrex::IntVect& /* nghost */,
+    const FieldState /* fstate */)
 {
-    FieldFillPatchOps<FieldBCDirichlet>::fillphysbc(
-        lev, time, mfab, nghost, fstate);
-
     if (m_field.base_name() == "velocity") {
         m_abl_mpl.set_velocity(lev, time, m_field, mfab);
     } else if (m_field.base_name() == "temperature") {
@@ -78,39 +69,36 @@ void FillMPL::fillpatch_sibling_fields(
     const amrex::Vector<amrex::BCRec>& /* unused */,
     const FieldState fstate)
 {
-    if (m_field.base_name() == "velocity") {
-        // For an ABL MPL, we first just foextrap the mac velocities
-        amrex::Vector<amrex::BCRec> lbcrec(m_field.num_comp());
-        const auto& ibctype = m_field.bc_type();
-        for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
-            auto ori = oit();
-            const auto side = ori.faceDir();
-            const auto bct = ibctype[ori];
-            const int dir = ori.coordDir();
-            for (int i = 0; i < m_field.num_comp(); ++i) {
-                if ((bct == BC::mass_inflow) ||
-                    (bct == BC::mass_inflow_outflow)) {
-                    if (side == amrex::Orientation::low) {
-                        lbcrec[i].setLo(dir, amrex::BCType::foextrap);
-                    } else {
-                        lbcrec[i].setHi(dir, amrex::BCType::foextrap);
-                    }
+    // For an ABL MPL, we first just foextrap the mac velocities
+    amrex::Vector<amrex::BCRec> lbcrec(m_field.num_comp());
+    const auto& ibctype = m_field.bc_type();
+    for (amrex::OrientationIter oit; oit != nullptr; ++oit) {
+        auto ori = oit();
+        const auto side = ori.faceDir();
+        const auto bct = ibctype[ori];
+        const int dir = ori.coordDir();
+        for (int i = 0; i < m_field.num_comp(); ++i) {
+            if ((bct == BC::mass_inflow) || (bct == BC::mass_inflow_outflow)) {
+                if (side == amrex::Orientation::low) {
+                    lbcrec[i].setLo(dir, amrex::BCType::foextrap);
                 } else {
-                    if (side == amrex::Orientation::low) {
-                        lbcrec[i].setLo(dir, bcrec[i].lo(dir));
-                    } else {
-                        lbcrec[i].setHi(dir, bcrec[i].hi(dir));
-                    }
+                    lbcrec[i].setHi(dir, amrex::BCType::foextrap);
+                }
+            } else {
+                if (side == amrex::Orientation::low) {
+                    lbcrec[i].setLo(dir, bcrec[i].lo(dir));
+                } else {
+                    lbcrec[i].setHi(dir, bcrec[i].hi(dir));
                 }
             }
         }
+    }
 
-        FieldFillPatchOps<FieldBCDirichlet>::fillpatch_sibling_fields(
-            lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, lbcrec, fstate);
+    FieldFillPatchOps<FieldBCDirichlet>::fillpatch_sibling_fields(
+        lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, lbcrec, fstate);
 
-        for (int i = 0; std::cmp_less(i, mfabs.size()); i++) {
-            m_abl_mpl.set_velocity(lev, time, m_field, *mfabs[i], 0, i);
-        }
+    for (int i = 0; std::cmp_less(i, mfabs.size()); i++) {
+        m_abl_mpl.set_velocity(lev, time, m_field, *mfabs[i], 0, i);
     }
 }
 

src/boundary_conditions/field_boundary_fill/Flather.H

@@ -0,0 +1,107 @@
+#ifndef FLATHER_H
+#define FLATHER_H
+
+#include "src/core/Field.H"
+#include "src/core/IntField.H"
+#include "src/CFDSim.H"
+#include "src/boundary_conditions/field_boundary_fill/FieldBoundary.H"
+#include "src/utilities/MultiLevelVector.H"
+#include "src/utilities/trig_ops.H"
+#include "AMReX_REAL.H"
+#include <limits>
+
+using namespace amrex::literals;
+
+namespace kynema_sgf {
+
+/** Flather-style open boundary for velocity
+ *
+ *  This boundary reads current boundary ghost-cell values and adjacent interior
+ *  values, computes a blended target, and writes updated boundary values.
+ */
+class Flather : public FieldBoundary::Register<Flather>
+{
+public:
+    static std::string identifier() { return "Flather"; }
+
+    explicit Flather(CFDSim& sim);
+
+    void post_init_actions() override;
+
+    void pre_advance_work() override;
+
+    void post_advance_work() override {}
+
+    void set_velocity(
+        int lev,
+        amrex::Real time,
+        const Field& fld,
+        amrex::MultiFab& mfab,
+        int dcomp = 0,
+        int orig_comp = 0) const override;
+
+    void update_flather_variables(const int lev, const FieldState fstate)
+    {
+        compute_boundary_z_averages(lev, fstate, false);
+    }
+
+    void update_flather_variables(
+        const int lev, const FieldState fstate, const bool use_mac_fields)
+    {
+        compute_boundary_z_averages(lev, fstate, use_mac_fields);
+    }
+
+    void accumulate_boundary(
+        int current_level,
+        int idir,
+        int phase_switch,
+        bool is_low,
+        MultiLevelVector& out_uvec,
+        MultiLevelVector& out_hvec,
+        bool sample_boundary,
+        FieldState fstate,
+        bool use_mac_fields = false) const;
+
+private:
+    void compute_internal_z_averages();
+    void compute_boundary_z_averages(
+        int lev, FieldState fstate, bool use_mac_fields = false);
+
+    const CFDSim& m_sim;
+    const kynema_sgf::SimTime& m_time;
+    const FieldRepo& m_repo;
+    const amrex::AmrCore& m_mesh;
+    Field& m_velocity;
+    const Field& m_u_mac;
+    const Field& m_v_mac;
+    const Field& m_vof;
+    const IntField* m_terrain_blank{nullptr};
+    amrex::Vector<amrex::Real> m_gravity{0.0_rt, 0.0_rt, -9.81_rt};
+    amrex::Real m_rho1{1000.0_rt};
+    amrex::Real m_rho2{1.0_rt};
+
+    MultiLevelVector m_xlo_uliq{FieldLoc::CELL};
+    MultiLevelVector m_xhi_uliq{FieldLoc::CELL};
+    MultiLevelVector m_ylo_uliq{FieldLoc::CELL};
+    MultiLevelVector m_yhi_uliq{FieldLoc::CELL};
+    MultiLevelVector m_xlo_umix{FieldLoc::CELL};
+    MultiLevelVector m_xhi_umix{FieldLoc::CELL};
+    MultiLevelVector m_ylo_umix{FieldLoc::CELL};
+    MultiLevelVector m_yhi_umix{FieldLoc::CELL};
+    MultiLevelVector m_xlo_bnd_uvof{FieldLoc::CELL};
+    MultiLevelVector m_xhi_bnd_uvof{FieldLoc::CELL};
+    MultiLevelVector m_ylo_bnd_uvof{FieldLoc::CELL};
+    MultiLevelVector m_yhi_bnd_uvof{FieldLoc::CELL};
+    MultiLevelVector m_xlo_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_xhi_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_ylo_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_yhi_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_xlo_bnd_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_xhi_bnd_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_ylo_bnd_h_avg{FieldLoc::CELL};
+    MultiLevelVector m_yhi_bnd_h_avg{FieldLoc::CELL};
+};
+
+} // namespace kynema_sgf
+
+#endif /* FLATHER_H */