Another way to fix the UDF MAC fills

amr-wind/core/FieldBCOps.H

@@ -11,6 +11,18 @@
 
 namespace amr_wind {
 
+namespace {
+template <int face_dir>
+struct FaceOp
+{
+    AMREX_GPU_DEVICE AMREX_FORCE_INLINE int operator()(const int idir)
+    {
+        return (idir == face_dir) ? 1 : 0;
+    }
+};
+
+} // namespace
+
 class Field;
 
 /** Interface to setup up boundary conditions on a field
@@ -42,6 +54,8 @@ public:
 struct FieldBCNoOp
 {
     using FunctorType = FieldBCNoOp;
+    template <int /* Index */>
+    using FunctorFType = FieldBCNoOp;
 
     AMREX_GPU_HOST
     constexpr FieldBCNoOp() = default;
@@ -50,6 +64,11 @@ struct FieldBCNoOp
     constexpr explicit FieldBCNoOp(const Field& /*unused*/) {}
 
     FunctorType operator()() const { return *this; }
+    template <typename T>
+    T opf() const
+    {
+        return *this;
+    }
 
     AMREX_GPU_DEVICE
     void operator()(
@@ -110,7 +129,7 @@ struct ConstDirichlet
 /** Sets constant values on specified boundaries
  *  \ingroup field_fillpatch
  */
-template <typename InflowOp, typename WallOp>
+template <typename InflowOp, typename WallOp, int face_dir = -1>
 struct DirichletOp
 {
     using InflowOpType = InflowOp;
@@ -168,10 +187,13 @@ struct DirichletOp
                     continue;
                 }
 
-                // Check if the point in question is on the boundary
+                // Check if the point in question is on the boundary, which
+                // depends on the location of the field variable
+                const int bigEnd_for_floc =
+                    domain_box.bigEnd(idir) + FaceOp<face_dir>()(idir);
                 const bool is_bndry =
                     (ori.isLow() ? (iv[idir] < domain_box.smallEnd(idir))
-                                 : (iv[idir] > domain_box.bigEnd(idir)));
+                                 : (iv[idir] > bigEnd_for_floc));
                 if (!is_bndry) {
                     continue;
                 }
@@ -204,10 +226,17 @@ struct DirichletOp
 struct FieldBCDirichlet
 {
     using FunctorType = DirichletOp<ConstDirichlet, ConstDirichlet>;
+    template <int Index>
+    using FunctorFType = DirichletOp<ConstDirichlet, ConstDirichlet, Index>;
 
     explicit FieldBCDirichlet(const Field& fld) : m_field(fld) {}
 
     inline FunctorType operator()() const { return FunctorType(m_field); }
+    template <typename T>
+    inline T opf() const
+    {
+        return T(m_field);
+    }
 
     const Field& m_field;
 };
@@ -218,6 +247,8 @@ struct BCOpCreator
     using InflowOpType = typename InflowOp::DeviceType;
     using WallOpType = typename WallOp::DeviceType;
     using FunctorType = DirichletOp<InflowOpType, WallOpType>;
+    template <int Index>
+    using FunctorFType = DirichletOp<InflowOpType, WallOpType, Index>;
 
     explicit BCOpCreator(const Field& fld)
         : m_field(fld), m_inflow_op(fld), m_wall_op(fld)
@@ -235,6 +266,14 @@ struct BCOpCreator
             m_wall_op.device_instance()};
     }
 
+    template <typename T>
+    inline T opf() const
+    {
+        return T{
+            m_field, m_inflow_op.device_instance(),
+            m_wall_op.device_instance()};
+    }
+
     const Field& m_field;
     InflowOp m_inflow_op;
     WallOp m_wall_op;

amr-wind/core/FieldFillPatchOps.H

@@ -167,6 +167,8 @@
 {
 public:
     using Functor = typename BCOpCreator::FunctorType;
+    template <int Index>
+    using FunctorF = typename BCOpCreator::template FunctorFType<Index>;
     /**
      *  @param field Field whose patches are filled by this instance
      *  @param mesh  The mesh instance to determine amrex::Geometry at a level
@@ -270,13 +272,34 @@
     {
 
         if (lev == 0) {
-            amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> physbc(
-                m_mesh.Geom(lev), bcrec, bc_functor());
+            amrex::PhysBCFunct<amrex::GpuBndryFuncFab<FunctorF<0>>> physbcx(
+                m_mesh.Geom(lev), bcrec, bc_functor_f<FunctorF<0>>());
+            amrex::PhysBCFunct<amrex::GpuBndryFuncFab<FunctorF<1>>> physbcy(
+                m_mesh.Geom(lev), bcrec, bc_functor_f<FunctorF<1>>());
+            amrex::PhysBCFunct<amrex::GpuBndryFuncFab<FunctorF<2>>> physbcz(
+                m_mesh.Geom(lev), bcrec, bc_functor_f<FunctorF<2>>());
+            int i = 0;
+            std::cout << "sibling fields FPSL dir i = 0 (umac)\n";
+            amrex::FillPatchSingleLevel(
+                *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1,
+                m_mesh.Geom(lev), physbcx, i);
+            i = 1;
+            std::cout << "sibling fields FPSL dir i = 1 (vmac)\n";
+            amrex::FillPatchSingleLevel(
+                *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1,
+                m_mesh.Geom(lev), physbcy, i);
+            i = 2;
+            std::cout << "sibling fields FPSL dir i = 2 (wmac)\n";
+            amrex::FillPatchSingleLevel(
+                *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1,
+                m_mesh.Geom(lev), physbcz, i);
+            /*
             for (int i = 0; i < static_cast<int>(mfabs.size()); i++) {
                 amrex::FillPatchSingleLevel(
                     *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1,
                     m_mesh.Geom(lev), physbc, i);
             }
+            */
         } else {
             amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> cphysbc(
                 m_mesh.Geom(lev - 1), bcrec, bc_functor());
@@ -386,6 +409,11 @@
 
 protected:
     Functor bc_functor() { return m_op(); }
+    template <typename T>
+    T bc_functor_f()
+    {
+        return m_op.template opf<T>();
+    }
 
     void check_face_mapper()
     {

amr-wind/equation_systems/icns/icns_advection.H

@@ -301,7 +301,7 @@ struct AdvectionOp<ICNS, fvm::Godunov>
         // MAC projection
         m_macproj_op(fstate, dt);
 
-        // Fill mac velocities using velocity BCs
+        std::cout << "// Fill mac velocities using velocity BCs\n";
         if (fvm::Godunov::nghost_state > 0) {
             amrex::Array<Field*, AMREX_SPACEDIM> mac_vel = {
                 AMREX_D_DECL(&u_mac, &v_mac, &w_mac)};

amr-wind/physics/udfs/Rankine.H

@@ -45,6 +45,27 @@
                     vel_ref[2] + 0.0)};
 
             field(iv[0], iv[1], iv[2], dcomp + comp) = vel[orig_comp + comp];
+
+            amrex::IntVect ivhh(AMREX_D_DECL(41, 12, 1));
+            amrex::IntVect ivhi(AMREX_D_DECL(40, 12, 1));
+            amrex::IntVect ivlo(AMREX_D_DECL(0, 12, 1));
+            amrex::IntVect ivlm(AMREX_D_DECL(-1, 12, 1));
+            if (iv == ivlo)
+                amrex::Print()
+                    << "Rankine.H op: filling " << vel[orig_comp + comp]
+                    << " at i=0" << std::endl;
+            if (iv == ivlm)
+                amrex::Print()
+                    << "Rankine.H op: filling " << vel[orig_comp + comp]
+                    << " at i=-1" << std::endl;
+            if (iv == ivhi)
+                amrex::Print()
+                    << "Rankine.H op: filling " << vel[orig_comp + comp]
+                    << " at i=40" << std::endl;
+            if (iv == ivhh)
+                amrex::Print()
+                    << "Rankine.H op: filling " << vel[orig_comp + comp]
+                    << " at i=41" << std::endl;
         }
     };
     using DeviceType = DeviceOp;

test/test_files/rankine/rankine.inp

@@ -66,7 +66,11 @@ xlo.type = "mass_inflow"
 xlo.density = 1.0
 xlo.velocity.inflow_type = Rankine
 xlo.temperature = 300.0
-xhi.type = "pressure_outflow"
+xhi.type = "mass_inflow"
+xhi.density = 1.0
+xhi.velocity.inflow_type = Rankine
+xhi.temperature = 300.0
+#xhi.type = "pressure_outflow"
 ylo.type = "slip_wall"
 yhi.type =   "slip_wall"
 zlo.type =   "slip_wall"