IncfloDenFill/TracFill/TempFill: the outflow part of a direction_dependent face is never written (filcc skips BCType 7, functors are inflow-only) and the in/out test uses the constant input velocity, not the probtype profile

[WeiqunZhang]

Severity: medium · Category: correctness
Locations: src/prob/prob_bc.H:381-395, src/prob/prob_bc.H:396-426, src/prob/prob_bc.H:511-555, src/prob/prob_bc.H:639-682, src/prob/prob_bc.H:50-63 (IncfloVelFill profiles), src/boundary_conditions/boundary_conditions.cpp:108-113
Based on commit 46de3367 (line numbers refer to that tree).

The defect

init_bcs gives density, every tracer and temperature BCType::direction_dependent on a
direction_dependent face (boundary_conditions.cpp:111-113). AMReX's FilccCell has no case for that
type (AMReX_FilCC_C.H handle foextrap/hoextrap/hoextrapcc/reflect), and GpuBndryFuncFab simply calls
FilccCell followed by the user functor on every ghost cell, so the functor owns the whole fill.

IncfloVelFill does that correctly: inflow (norm_vel >= 0 on a low face) writes the prescribed values,
outflow copies the first interior cell (prob_bc.H:94-97 and the five twins). The three scalar functors
only have the inflow branch, e.g. prob_bc.H:381-387

else if ( (i < domain_box.smallEnd(0)) &&
          ( (bc.lo(0) == amrex::BCType::ext_dir) ||
            (bc.lo(0) == amrex::BCType::direction_dependent &&
             bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
{
    rho(i,j,k) = bcv_den[...x low...];
}

and nothing else for direction_dependent. Two consequences:

1. Where the prescribed normal velocity points out of the domain (xlo.velocity = -1 0 0, or the high
   side with a positive value), the density/tracer/temperature ghost cells outside the domain are never
   written by anyone: not by FillPatch (physical ghosts are not covered by any grid), not by filcc, not
   by the functor. They keep whatever the MultiFab held (uninitialised at start, stale later).
2. The direction test uses the raw input value bcv_vel[face][dir], but IncfloVelFill decides
   inflow/outflow with the probtype-modified profile (prob_bc.H:50-75: 42 -> time, 31/311, 43 ->
   6y(1-y)-1, 41, 32/322, 33/333). For probtype 43 the velocity fill makes the middle of x-lo inflow and
   the outer bands outflow, while the scalar fills see bcv_vel = 0 >= 0 and put the inflow Dirichlet
   value on the whole face, including the outflow bands.

Why it matters

The scalar ghost cells on a direction_dependent face are consumed by (a) the Godunov/PLM slope stencils
next to the boundary (hydro_godunov_plm.H:77 treats only dd && umac >= 0 as one-sided; otherwise the
ghost value enters the slope), (b) the implicit diffusion solve, where get_diffuse_scalar_bc maps
direction_dependent to LinOpBCType::Dirichlet, so the ghost value is the boundary value, and
(c) EB_interp_CellCentroid_to_FaceCentroid of eta with the tracer BCRec. Case 1 is an uninitialised read
that becomes a wrong Dirichlet value; case 2 gives the inflow density/tracer/temperature as Dirichlet
value on an outflow segment (wrong flux and wrong diffusion BC there).

How to reach it

Any direction_dependent face with incflo.constant_density = false, incflo.advect_tracer = 1 or
incflo.use_temperature = 1. Concretely: test_no_eb_2d/benchmark.inout (probtype 43) with
incflo.advect_tracer = 1 xlo.tracer = 1 xhi.tracer = 0 shows case 2; the same input with
xlo.velocity = -1 0 0 (uniform prescribed outflow, no probtype profile) shows case 1. DIM 2 or 3,
EB or not, CPU or GPU. (The inout benchmark currently also aborts at start-up for the reason given in
issue 010; the ghost fill defect is independent of that.)

Suggested fix

Factor the prescribed normal velocity (input value plus probtype profiles) into one AMREX_GPU_HOST_DEVICE
helper incflo_bc_normal_velocity(probtype, face, i, j, k, domain, time, bcv_vel) in prob_bc.H, use it in
all four functors, and give the scalar functors the same "outflow: copy the first interior cell" branch
that IncfloVelFill has. The Burggraf case (probtype 16) is deliberately not in the helper: it is a moving
wall, not an inflow, and is handled by issue 011. The diff below leaves the HIGH J block of IncfloVelFill
untouched so that it does not overlap with 011; the velocity LOW I / HIGH I / LOW J / LOW K blocks are
switched to the helper (pure refactor, same values).

Diff checked with git apply --check (audit/notes/T4-scratch/012.diff, 475 lines); it is reproduced in
full below.

--- a/src/prob/prob_bc.H
+++ b/src/prob/prob_bc.H
@@ -5,6 +5,95 @@
 #include <AMReX_Geometry.H>
 #include <AMReX_PhysBCFunct.H>
 
+// Normal velocity prescribed on a domain face, including the probtype-specific
+// inflow profiles.  Every fill functor uses this so that the direction_dependent
+// inflow/outflow decision is the same for velocity, density, tracer and
+// temperature.  (Probtype 16 is a moving wall, not an inflow, and is handled in
+// IncfloVelFill directly.)
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real incflo_bc_normal_velocity (int probtype, amrex::Orientation face,
+                                       int i, int j, int k,
+                                       amrex::Box const& domain_box, amrex::Real time,
+                                       amrex::GpuArray<amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>, AMREX_SPACEDIM*2> const& bcv_vel)
+{
+    amrex::ignore_unused(i,j,k,time);
+    const int dir = face.coordDir();
+    amrex::Real norm_vel = bcv_vel[face][dir];
+
+    if (dir == 0 && face.isLow())
+    {
+        if (42 == probtype)
+        {
+            norm_vel = time;
+        }
+        else if (31 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1));
+            norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y);
+        }
+        else if (43 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1);
+        }
+#if (AMREX_SPACEDIM == 3)
+        else if (311 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
+            norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z);
+        }
+        else if (41 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
+            norm_vel = amrex::Real(0.5) * z;
+        }
+#endif
+    }
+    else if (dir == 0 && face.isHigh())
+    {
+        if (42 == probtype)
+        {
+            norm_vel = time;
+        }
+        else if (43 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1);
+        }
+    }
+    else if (dir == 1 && face.isLow())
+    {
+#if (AMREX_SPACEDIM == 3)
+        if (32 == probtype)
+        {
+            amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2));
+            norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z);
+        }
+#endif
+        if (322 == probtype)
+        {
+            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+            norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x);
+        }
+    }
+#if (AMREX_SPACEDIM == 3)
+    else if (dir == 2 && face.isLow())
+    {
+        if (33 == probtype)
+        {
+            amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
+            norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x);
+        }
+        else if (333 == probtype)
+        {
+            amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
+            norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y);
+        }
+    }
+#endif
+    return norm_vel;
+}
+
 struct IncfloVelFill
 {
     int probtype;
@@ -44,35 +133,9 @@
             if (i < domain_box.smallEnd(0))
             {
                 int dir = 0;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (42 == probtype)
-                {
-                    norm_vel = time;
-                }
-                else if (31 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1)/domain_box.length(1));
-                    norm_vel = amrex::Real(6.) * y * (amrex::Real(1)-y);
-                }
-                else if (43 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel = amrex::Real(6) * y * (amrex::Real(1)-y) - amrex::Real(1);
-                }
-#if (AMREX_SPACEDIM == 3)
-                else if (311 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
-                    norm_vel = amrex::Real(6.) * z * (amrex::Real(1)-z);
-                }
-                else if (41 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1)/domain_box.length(2));
-                    norm_vel = amrex::Real(0.5) * z;
-                }
-#endif
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::x,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1101 == probtype)
@@ -103,18 +166,9 @@
             if (i > domain_box.bigEnd(0))
             {
                 int dir = 0;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (42 == probtype)
-                {
-                    norm_vel = time;
-                }
-                else if (43 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel = amrex::Real(6) * y * (amrex::Real(1.0)-y) - amrex::Real(1);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::x,amrex::Orientation::high),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1101 == probtype)
@@ -145,21 +199,9 @@
             if (j < domain_box.smallEnd(1))
             {
                 int dir = 1;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-#if (AMREX_SPACEDIM == 3)
-                if (32 == probtype)
-                {
-                    amrex::Real z = amrex::Real(k+0.5)*(amrex::Real(1.0)/domain_box.length(2));
-                    norm_vel *= amrex::Real(6.) * z * (amrex::Real(1.0)-z);
-                }
-#endif
-                if (322 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel *= amrex::Real(6.) * x * (amrex::Real(1.0)-x);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::y,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 if ( (bc.lo(dir) == amrex::BCType::ext_dir) ||
                      (bc.lo(dir) == amrex::BCType::direction_dependent && norm_vel >= 0.) )
@@ -225,19 +267,9 @@
             if (k < domain_box.smallEnd(2))
             {
                 int dir = 2;
-                amrex::Real norm_vel = bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][dir];
-
-                // This may modify the normal velocity for specific problems
-                if (33 == probtype)
-                {
-                    amrex::Real x = amrex::Real(i+0.5)*(amrex::Real(1.0)/domain_box.length(0));
-                    norm_vel *= amrex::Real(6.0) * x * (amrex::Real(1.0)-x);
-                }
-                else if (333 == probtype)
-                {
-                    amrex::Real y = amrex::Real(j+0.5)*(amrex::Real(1.0)/domain_box.length(1));
-                    norm_vel *= amrex::Real(6.0) * y * (amrex::Real(1.0)-y);
-                }
+                amrex::Real norm_vel = incflo_bc_normal_velocity(probtype,
+                    amrex::Orientation(amrex::Direction::z,amrex::Orientation::low),
+                    i, j, k, domain_box, time, bcv_vel);
 
                 // This is a special case -- all the logic is contained here
                 if (1100 == probtype)
@@ -314,7 +346,7 @@
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& rho,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -381,48 +413,74 @@
         else if ( (i < domain_box.smallEnd(0)) &&
                   ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                     (bc.lo(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
         }
         else if ( (i > domain_box.bigEnd(0)) &&
                   ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                     (bc.hi(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)];
         }
+        // outflow part of a direction_dependent face: extrapolate the interior cell
+        // (filcc does nothing for BCType::direction_dependent)
+        else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(domain_box.smallEnd(0),j,k);
+        }
+        else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(domain_box.bigEnd(0),j,k);
+        }
 
         if ( (j < domain_box.smallEnd(1)) &&
               ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                 (bc.lo(1) == amrex::BCType::direction_dependent &&
-                 bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                 incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)];
         }
         else if ( (j > domain_box.bigEnd(1)) &&
                   ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                     (bc.hi(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)];
         }
+        else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,domain_box.smallEnd(1),k);
+        }
+        else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,domain_box.bigEnd(1),k);
+        }
 #if (AMREX_SPACEDIM == 3)
          if ( (k < domain_box.smallEnd(2)) &&
               ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                 (bc.lo(2) == amrex::BCType::direction_dependent &&
-                 bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                 incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)];
         }
+        else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,j,domain_box.smallEnd(2));
+        }
 
         if ( (k > domain_box.bigEnd(2)) &&
                   ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                     (bc.hi(2) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             rho(i,j,k) = bcv_den[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)];
         }
+        else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+        {
+            rho(i,j,k) = rho(i,j,domain_box.bigEnd(2));
+        }
 #endif
     }
 };
@@ -442,7 +500,7 @@
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& tracer,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -511,47 +569,73 @@
             else if ( (i < domain_box.smallEnd(0)) &&
                       ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                         (bc.lo(0) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][n];
             }
             else if ( (i > domain_box.bigEnd(0)) &&
                       ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                         (bc.hi(0) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][n];
             }
+            // outflow part of a direction_dependent face: extrapolate the interior cell
+            // (filcc does nothing for BCType::direction_dependent)
+            else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(domain_box.smallEnd(0),j,k,n);
+            }
+            else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(domain_box.bigEnd(0),j,k,n);
+            }
 
             if ( (j < domain_box.smallEnd(1)) &&
                   ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                     (bc.lo(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][n];
             }
             else if ( (j > domain_box.bigEnd(1)) &&
                       ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                         (bc.hi(1) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][n];
             }
+            else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,domain_box.smallEnd(1),k,n);
+            }
+            else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,domain_box.bigEnd(1),k,n);
+            }
 #if (AMREX_SPACEDIM == 3)
             if ( (k < domain_box.smallEnd(2)) &&
                  ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                    (bc.lo(2) == amrex::BCType::direction_dependent &&
-                    bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                    incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][n];
             }
             else if ( (k > domain_box.bigEnd(2)) &&
                       ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                         (bc.hi(2) == amrex::BCType::direction_dependent &&
-                         bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                         incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
             {
                 tracer(i,j,k,n) = bcv_tra[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][n];
             }
+            else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,j,domain_box.smallEnd(2),n);
+            }
+            else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+            {
+                tracer(i,j,k,n) = tracer(i,j,domain_box.bigEnd(2),n);
+            }
 #endif
         }
     }
@@ -572,7 +656,7 @@
     AMREX_GPU_DEVICE
     void operator() (const amrex::IntVect& iv, amrex::Array4<amrex::Real> const& temperature,
                      const int /*dcomp*/, const int /*numcomp*/,
-                     amrex::GeometryData const& geom, const amrex::Real /*time*/,
+                     amrex::GeometryData const& geom, const amrex::Real time,
                      const amrex::BCRec* bcr, const int bcomp,
                      const int /*orig_comp*/) const
     {
@@ -639,47 +723,73 @@
         else if ( (i < domain_box.smallEnd(0)) &&
                   ( (bc.lo(0) == amrex::BCType::ext_dir) ||
                     (bc.lo(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][0] >= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
         }
         else if ( (i > domain_box.bigEnd(0)) &&
                   ( (bc.hi(0) == amrex::BCType::ext_dir) ||
                     (bc.hi(0) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][0] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::x,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)];
         }
+        // outflow part of a direction_dependent face: extrapolate the interior cell
+        // (filcc does nothing for BCType::direction_dependent)
+        else if ( (i < domain_box.smallEnd(0)) && (bc.lo(0) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(domain_box.smallEnd(0),j,k);
+        }
+        else if ( (i > domain_box.bigEnd(0)) && (bc.hi(0) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(domain_box.bigEnd(0),j,k);
+        }
 
         if ( (j < domain_box.smallEnd(1)) &&
              ( (bc.lo(1) == amrex::BCType::ext_dir) ||
                (bc.lo(1) == amrex::BCType::direction_dependent &&
-                bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)][1] >= 0.) ) )
+                incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)];
         }
         else if ( (j > domain_box.bigEnd(1)) &&
                   ( (bc.hi(1) == amrex::BCType::ext_dir) ||
                     (bc.hi(1) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][1] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::y,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)];
         }
+        else if ( (j < domain_box.smallEnd(1)) && (bc.lo(1) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,domain_box.smallEnd(1),k);
+        }
+        else if ( (j > domain_box.bigEnd(1)) && (bc.hi(1) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,domain_box.bigEnd(1),k);
+        }
 #if (AMREX_SPACEDIM == 3)
         if ( (k < domain_box.smallEnd(2)) &&
              ( (bc.lo(2) == amrex::BCType::ext_dir) ||
                (bc.lo(2) == amrex::BCType::direction_dependent &&
-                bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)][2] >= 0.) ) )
+                incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::low), i, j, k, domain_box, time, bcv_vel) >= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::low)];
         }
         else if ( (k > domain_box.bigEnd(2)) &&
                   ( (bc.hi(2) == amrex::BCType::ext_dir) ||
                     (bc.hi(2) == amrex::BCType::direction_dependent &&
-                     bcv_vel[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)][2] <= 0.) ) )
+                     incflo_bc_normal_velocity(probtype, amrex::Orientation(amrex::Direction::z,amrex::Orientation::high), i, j, k, domain_box, time, bcv_vel) <= 0.) ) )
         {
             temperature(i,j,k) = bcv_tem[amrex::Orientation(amrex::Direction::z,amrex::Orientation::high)];
         }
+        else if ( (k < domain_box.smallEnd(2)) && (bc.lo(2) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,j,domain_box.smallEnd(2));
+        }
+        else if ( (k > domain_box.bigEnd(2)) && (bc.hi(2) == amrex::BCType::direction_dependent) )
+        {
+            temperature(i,j,k) = temperature(i,j,domain_box.bigEnd(2));
+        }
 #endif
     }
 };

[asalmgren]

closed by #265