diff --git a/Docs/sphinx_documentation/source/InputsAlgorithm.rst b/Docs/sphinx_documentation/source/InputsAlgorithm.rst
index e7ca962f..29ee6c54 100644
--- a/Docs/sphinx_documentation/source/InputsAlgorithm.rst
+++ b/Docs/sphinx_documentation/source/InputsAlgorithm.rst
@@ -26,4 +26,7 @@ The following inputs must be preceded by "incflo."
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
 |  mu_s                |  scalar diffusivity                                                   |  Real(s)    |  0.0         |
 +----------------------+-----------------------------------------------------------------------+-------------+--------------+
+|  use_tensor_solve    |  In velocity solve, use multicomponent :math:`\nabla \cdot \tau`      |  bool       |  true        |
+|                      |  otherwise use separate solves for each velocity component            |             |              |
++----------------------+-----------------------------------------------------------------------+-------------+--------------+
 
diff --git a/Docs/sphinx_documentation/source/InputsProblemDefinition.rst b/Docs/sphinx_documentation/source/InputsProblemDefinition.rst
index f2f44125..477c99cf 100644
--- a/Docs/sphinx_documentation/source/InputsProblemDefinition.rst
+++ b/Docs/sphinx_documentation/source/InputsProblemDefinition.rst
@@ -50,6 +50,7 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo",
 |                    | * 'po'  or 'pressure_outflow'                                             |             |           |
 |                    | * 'mi'  or 'mass_inflow'                                                  |             |           |
 |                    | * 'nsw' or 'no_slip_wall'                                                 |             |           |
+|                    | * 'mixed'                                                                 |             |           |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
 | pressure           | Sets boundary pressure for pressure inflows, outflows and mass inflows    |    Real     |  None     |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
@@ -60,5 +61,14 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo",
 | tracer             | Sets boundary tracer for mass inflows                                     |    Real     |  0.0      |
 +--------------------+---------------------------------------------------------------------------+-------------+-----------+
 
+   The 'mixed' boundary type allows for inflow and outflow on the same domain face.
+   The implementation requires that there is EB separating the inflow and outflow regions,
+   and only currently treats constant viscosity (i.e. requires the inputs file contains  ``incflo.use_tensor_solve = false``).
+   It does allow for non-constant scalar diffusivity. To create a new problem setup with mixed BCs, one must additionally
+   specify the Dirchlet and Neumann areas in ``prob/prob_bc.cpp`` In this file there are functions to set the needed BC info for the diffusion solver, the MAC projection, and with the same function, the nodal projection and advection. Comments in the code provide a detailed explaination for each instance. For addditional details on mixed BCs, also see AMReX-Hydro's documentation (:ref:`bcs`).
 
+..
+   To create a new setup:
+   initial conditions go in prob/prob_init_fluid.cpp
+   inflow boundary conditions are set in prob/prob_bc.H -- the BCType for mixed bc is set at foextrap here to allow the outflow to get appropriately filled, which would already have happened before reaching these functions.
 
diff --git a/src/boundary_conditions/boundary_conditions.cpp b/src/boundary_conditions/boundary_conditions.cpp
index b7c34ac2..c413fd49 100644
--- a/src/boundary_conditions/boundary_conditions.cpp
+++ b/src/boundary_conditions/boundary_conditions.cpp
@@ -91,6 +91,39 @@ void incflo::init_bcs ()
             // tracer has Dirichlet bcs
             pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
         }
+        else if (bc_type == "mixed" )
+        {
+            amrex::Print() << bcid << " set to mixed inflow outflow.\n";
+            m_has_mixedBC = true;
+#ifdef AMREX_USE_EB
+            // ReadParameters() already called
+            if (m_advection_type != "Godunov") { amrex::Abort("mixed BCs require Godunov"); }
+
+            ParmParse ipp("incflo");
+            std::string eb_geom = "null";
+            ipp.query("geometry", eb_geom);
+            eb_geom = amrex::toLower(eb_geom);
+            if (eb_geom == "null" || eb_geom == "all_regular")
+#endif
+            {
+                Abort("For now, mixed BCs must be separated by an EB");
+            }
+            Warning("Using BC type mixed requires that the Dirichlet and Neumann regions are separated by EB.");
+
+            m_bc_type[ori] = BC::mixed;
+
+            pp.get("pressure", m_bc_pressure[ori]);
+
+            std::vector<Real> v;
+            if (pp.queryarr("velocity", v, 0, AMREX_SPACEDIM)) {
+               for (int i=0; i<AMREX_SPACEDIM; i++){
+                   m_bc_velocity[ori][i] = v[i];
+               }
+            }
+
+            pp.query("density", m_bc_density[ori]);
+            pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
+        }
         else
         {
             m_bc_type[ori] = BC::undefined;
@@ -146,7 +179,9 @@ void incflo::init_bcs ()
             Orientation::Side side = ori.faceDir();
             auto const bct = m_bc_type[ori];
             if (bct == BC::pressure_inflow ||
-                bct == BC::pressure_outflow)
+                bct == BC::pressure_outflow ||
+                bct == BC::mixed ) // BC_mask will handle Dirichlet part.
+
             {
                 if (side == Orientation::low) {
                     AMREX_D_TERM(m_bcrec_velocity[0].setLo(dir, BCType::foextrap);,
@@ -220,7 +255,8 @@ void incflo::init_bcs ()
             auto const bct = m_bc_type[ori];
             if (bct == BC::pressure_inflow  ||
                 bct == BC::pressure_outflow ||
-                bct == BC::no_slip_wall)
+                bct == BC::no_slip_wall ||
+                bct == BC::mixed)
             {
                 if (side == Orientation::low) {
                     m_bcrec_density[0].setLo(dir, BCType::foextrap);
@@ -271,7 +307,8 @@ void incflo::init_bcs ()
             Orientation::Side side = ori.faceDir();
             auto const bct = m_bc_type[ori];
             if (bct == BC::pressure_inflow  ||
-                bct == BC::pressure_outflow)
+                bct == BC::pressure_outflow ||
+                bct == BC::mixed)
             {
                 if (side == Orientation::low) {
                     for (auto& b : m_bcrec_tracer) b.setLo(dir, BCType::foextrap);
diff --git a/src/boundary_conditions/incflo_set_bcs.cpp b/src/boundary_conditions/incflo_set_bcs.cpp
index a2ef9713..eec157a4 100644
--- a/src/boundary_conditions/incflo_set_bcs.cpp
+++ b/src/boundary_conditions/incflo_set_bcs.cpp
@@ -6,6 +6,177 @@
 
 using namespace amrex;
 
+// Make the iMultiFab needed by the Nodal Projection to support mixed BCs
+iMultiFab
+incflo::make_nodalBC_mask(int lev)
+{
+    // We use the solver's overset mask to implement mixed BCs. The mask is nodal, and
+    // so exists on the boundary face, where 0 indicates a Dirichlet boundary condition
+    // (i,e, no solve needed for that point because we know the solution there already).
+    // Because it's an overset mask and not a BC, we must fill the entire MF.
+    iMultiFab new_mask(amrex::convert(grids[lev], IntVect::TheNodeVector()), dmap[lev], 1, 0);
+    new_mask = 1;
+    int inflow = 1;
+    int outflow = 0;
+
+    Geometry const& gm = Geom(lev);
+    Box const& domain = gm.Domain();
+    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+        Orientation olo(dir,Orientation::low);
+        Orientation ohi(dir,Orientation::high);
+        if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) {
+            Box dlo = (m_bc_type[olo] == BC::mixed) ? surroundingNodes(bdryLo(domain,dir)) : Box().setType(IndexType::TheNodeType());
+            Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryHi(domain,dir)) : Box().setType(IndexType::TheNodeType());
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            for (MFIter mfi(new_mask); mfi.isValid(); ++mfi) {
+                Box blo = mfi.validbox() & dlo;
+                Box bhi = mfi.validbox() & dhi;
+                Array4<int> const& mask_arr = new_mask.array(mfi);
+                if (blo.ok()) {
+                    prob_set_BC_MF(olo, blo, mask_arr, lev, inflow, outflow, "projection");
+                }
+                if (bhi.ok()) {
+                    prob_set_BC_MF(ohi, bhi, mask_arr, lev, inflow, outflow, "projection");
+                }
+            }
+        }
+    }
+
+    return new_mask;
+}
+
+// Make a position dependent BC MultiFab for the advection routines to support
+// mixed BCs.
+std::unique_ptr<iMultiFab>
+incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
+                   std::string const& field)
+{
+    auto ncomp = static_cast<int>(bcs.size());
+    // The advection routines expect that the BC type is stored in the first ghost
+    // cell of a cell-centered MF.
+    std::unique_ptr<iMultiFab> BC_MF(new iMultiFab(grids[lev], dmap[lev], ncomp, 1));
+    // initialize to 0 == BCType::int_dir, not sure this is needed really
+    *BC_MF = 0;
+    // For advection, we use FOEXTRAP for outflow regions per incflo::init_bcs()
+    int inflow = BCType::ext_dir;
+    int outflow = BCType::foextrap;
+
+    Box const& domain = geom[lev].Domain();
+    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+        Orientation olo(dir,Orientation::low);
+        Orientation ohi(dir,Orientation::high);
+
+        Box dlo = adjCellLo(domain,dir);
+        Box dhi = adjCellHi(domain,dir);
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+        for (MFIter mfi(*BC_MF); mfi.isValid(); ++mfi) {
+            Box blo = mfi.growntilebox() & dlo;
+            Box bhi = mfi.growntilebox() & dhi;
+            Array4<int> const& bc_arr = BC_MF->array(mfi);
+            BCRec const* const bc_ptr = bcs.data();
+            if (m_bc_type[olo] == BC::mixed) {
+                prob_set_BC_MF(olo, blo, bc_arr, lev, inflow, outflow, field);
+            } else {
+                amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    for (int n = 0; n < ncomp; n++) {
+                        bc_arr(i,j,k,n) = bc_ptr[n].lo(dir);
+                    }
+                });
+            }
+            if (m_bc_type[ohi] == BC::mixed) {
+                prob_set_BC_MF(ohi, bhi, bc_arr, lev, inflow, outflow, field);
+            } else {
+                amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+                {
+                    for (int n = 0; n < ncomp; n++) {
+                        bc_arr(i,j,k,n) = bc_ptr[n].hi(dir);
+                    }
+                });
+            }
+        }
+    }
+
+    return BC_MF;
+}
+
+// Create the MutliFabs needed by the linear solver to define Robin BCs, which incflo
+// uses to support mixed BCs in the MAC Projection and diffusion.
+Vector<MultiFab>
+incflo::make_robinBC_MFs(int lev, MultiFab* state)
+{
+    // MLMG defines the Robin BC with
+    //     a u + b du/dn = f
+    // The values are stored in the first ghost cell, although the BC is considered on
+    // the face. Only the ghost cells of robin BC arrays are used in MLMG.
+    int nghost = 1;
+
+    Vector<MultiFab> robin(3);
+    for (int n = 0; n < 3; n++) {
+        robin[n].define(grids[lev], dmap[lev], 1, nghost);
+    }
+
+    MultiFab& robin_a = robin[0];
+    MultiFab& robin_b = robin[1];
+    MultiFab& robin_f = robin[2];
+
+    Box const& domain = Geom(lev).Domain();
+    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
+        Orientation olo(dir,Orientation::low);
+        Orientation ohi(dir,Orientation::high);
+
+        // Only need to fill Robin BC sides, MLMG will check for Robin BC first
+        if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) {
+            Box dlo = (m_bc_type[olo] == BC::mixed) ? adjCellLo(domain,dir) : Box();
+            Box dhi = (m_bc_type[ohi] == BC::mixed) ? adjCellHi(domain,dir) : Box();
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            // FIXME - do we want tiling here...
+            for (MFIter mfi(robin_a,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+                Box const& gbx = amrex::grow(mfi.validbox(),nghost);
+                Box blo = gbx & dlo;
+                Box bhi = gbx & dhi;
+                Array4<Real> const& a_arr     = robin_a.array(mfi);
+                Array4<Real> const& b_arr     = robin_b.array(mfi);
+                Array4<Real> const& f_arr     = robin_f.array(mfi);
+
+                if (state) {
+                    // state has filled boundry, so can use its ghost cells for the BC values
+                    Array4<Real const> const& bcv = state->const_array(mfi);
+                    // Diffusion, solving for velocity, density, tracer
+                    //
+                    // Robin BC:   a u + b du/dn = f  -- inflow,  Dirichlet a=1, b=0, f=bcv
+                    //                                -- outflow, Neumann   a=0, b=1, f=0
+                    if (blo.ok()) {
+                        prob_set_diffusion_robinBCs(olo, blo, a_arr, b_arr, f_arr, bcv, lev);
+                    }
+                    if (bhi.ok()) {
+                        prob_set_diffusion_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, bcv, lev);
+                    }
+                } else {
+                    // MAC, solving for phi (~pressure)
+                    //
+                    // Robin BC:   a u + b du/dn = f  -- inflow,  Neumann   a=0, b=1, f=0
+                    //                                -- outflow, Dirichlet a=1, b=0, f=0
+                    if (blo.ok()) {
+                        prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev);
+                    }
+                    if (bhi.ok()) {
+                        prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev);
+                    }
+                }
+            }
+        }
+    }
+
+    return robin;
+}
+
 #ifdef AMREX_USE_EB
 void
 incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)
diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp
index 9332c4a5..31680fde 100644
--- a/src/convection/incflo_compute_MAC_projected_velocities.cpp
+++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp
@@ -110,6 +110,15 @@ incflo::compute_MAC_projected_velocities (
             macproj->initProjector(lp_info, GetVecOfArrOfConstPtrs(inv_rho));
         }
         macproj->setDomainBC(get_mac_projection_bc(Orientation::low), get_mac_projection_bc(Orientation::high));
+
+        if ( m_has_mixedBC ) {
+            for (int lev = 0; lev <= finest_level; ++lev)
+            {
+                auto const robin = make_robinBC_MFs(lev);
+                macproj->setLevelBC(lev, nullptr,
+                                    &robin[0], &robin[1], &robin[2]);
+            }
+        }
     } else {
 #ifndef AMREX_USE_EB
         if (m_constant_density) {
@@ -117,7 +126,7 @@ incflo::compute_MAC_projected_velocities (
         } else
 #endif
         {
-            macproj->updateBeta(GetVecOfArrOfConstPtrs(inv_rho));
+            macproj->updateCoeffs(GetVecOfArrOfConstPtrs(inv_rho));
         }
     }
 
@@ -138,6 +147,19 @@ incflo::compute_MAC_projected_velocities (
             l_advection_type = "Godunov";
         }
 
+        std::unique_ptr<iMultiFab> BC_MF;
+        if (m_has_mixedBC) {
+            // Create a MF to hold the BCType info. Note that this is different than the
+            // bcs for the MAC projection because the MAC operates on phi, this is velocity.
+            //
+            // TODO? Could consider creating an incflo member variable to save the BC_MF
+            //     amrex::Vector<std::unique_ptr<amrex::iMultiFab> > m_BC_MF;
+            // Could stack components as vel, density, tracer, and then could use for
+            // scalars' advective step as well. But not sure it really matters one way or
+            // the other.
+            BC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
+        }
+
         // Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
         //    returned from this call are on face CENTROIDS
         HydroUtils::ExtrapVelToFaces(*vel[lev], *vel_forces[lev],
@@ -149,7 +171,7 @@ incflo::compute_MAC_projected_velocities (
                                       m_eb_flow.enabled ? get_velocity_eb()[lev] : nullptr,
 #endif
                                       m_godunov_ppm, m_godunov_use_forces_in_trans,
-                                      l_advection_type);
+                                      l_advection_type, PPM::default_limiter, BC_MF.get());
     }
 
     Vector<Array<MultiFab*,AMREX_SPACEDIM> > mac_vec(finest_level+1);
@@ -190,7 +212,6 @@ incflo::compute_MAC_projected_velocities (
     } else {
         macproj->project(m_mac_mg_rtol,m_mac_mg_atol);
     }
-
     // Note that the macproj->project call above ensures that the MAC velocities are averaged down --
     //      we don't need to do that again here
 }
diff --git a/src/convection/incflo_compute_advection_term.cpp b/src/convection/incflo_compute_advection_term.cpp
index 1e52b2f6..ef29eeed 100644
--- a/src/convection/incflo_compute_advection_term.cpp
+++ b/src/convection/incflo_compute_advection_term.cpp
@@ -271,6 +271,24 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
         // ************************************************************************
         // Compute advective fluxes
         // ************************************************************************
+        //
+        // Create BC MF
+        // FIXME -- would it be worth it to save this from vel extrap???
+        std::unique_ptr<iMultiFab> velBC_MF;
+        if (m_has_mixedBC) {
+            velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
+        }
+        std::unique_ptr<iMultiFab> densBC_MF;
+        if (m_has_mixedBC) {
+            densBC_MF = make_BC_MF(lev, m_bcrec_density_d, "density");
+        }
+        std::unique_ptr<iMultiFab> tracBC_MF;
+        if (m_advect_tracer  && (m_ntrac>0)) {
+            if (m_has_mixedBC) {
+                tracBC_MF = make_BC_MF(lev, m_bcrec_tracer_d, "tracer");
+            }
+        }
+
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
@@ -322,6 +340,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
             int face_comp = 0;
             int ncomp = AMREX_SPACEDIM;
             bool is_velocity = true;
+            Array4<int const> const& velbc_arr = velBC_MF ? (*velBC_MF).const_array(mfi)
+                                                          : Array4<int const>{};
             HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
                                                      (m_advect_momentum) ? rhovel[lev].array(mfi) : vel[lev]->const_array(mfi),
                                                      AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
@@ -350,8 +370,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 #endif
                                                      m_godunov_ppm, m_godunov_use_forces_in_trans,
                                                      is_velocity, fluxes_are_area_weighted,
-                                                     m_advection_type);
-
+                                                     m_advection_type, PPM::default_limiter, velbc_arr);
 
             // ************************************************************************
             // Density
@@ -361,6 +380,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                 face_comp = AMREX_SPACEDIM;
                 ncomp = 1;
                 is_velocity = false;
+                Array4<int const> const& densbc_arr = densBC_MF ? (*densBC_MF).const_array(mfi)
+                                                                : Array4<int const>{};
                 HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
                                                          density[lev]->const_array(mfi),
                                                          AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
@@ -384,7 +405,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 #endif
                                                          m_godunov_ppm, m_godunov_use_forces_in_trans,
                                                          is_velocity, fluxes_are_area_weighted,
-                                                         m_advection_type);
+                                                         m_advection_type, PPM::default_limiter, densbc_arr);
             }
 
             // ************************************************************************
@@ -414,7 +435,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
                    face_comp = AMREX_SPACEDIM+1;
                 ncomp = m_ntrac;
                 is_velocity = false;
-
+                Array4<int const> const& tracbc_arr = tracBC_MF ? (*tracBC_MF).const_array(mfi)
+                                                                : Array4<int const>{};
                 HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, ro_trac,
                                           AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
                                                        flux_y[lev].array(mfi,face_comp),
@@ -438,11 +460,12 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 #endif
                                           m_godunov_ppm, m_godunov_use_forces_in_trans,
                                           is_velocity, fluxes_are_area_weighted,
-                                          m_advection_type);
+                                          m_advection_type, PPM::default_limiter, tracbc_arr);
             }
         } // mfi
     } // lev
 
+
     // In order to enforce conservation across coarse-fine boundaries we must be sure to average down the fluxes
     //    before we use them.  Note we also need to average down the face states if we are going to do
     //    convective differencing
@@ -646,40 +669,39 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
         drdt_tmp.FillBoundary(geom[lev].periodicity());
         dtdt_tmp.FillBoundary(geom[lev].periodicity());
 
-//Was this OMP intentionally left off?
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
         for (MFIter mfi(*density[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
         {
-        Box const& bx = mfi.tilebox();
-
-        // velocity
-        auto const& bc_vel = get_velocity_bcrec_device_ptr();
-        redistribute_term(mfi, *conv_u[lev], dvdt_tmp,
-                  (m_advect_momentum) ? rhovel[lev] : *vel[lev],
-                  bc_vel, lev);
-
-        // density
-        if (!m_constant_density) {
-        auto const& bc_den = get_density_bcrec_device_ptr();
-        redistribute_term(mfi, *conv_r[lev], drdt_tmp,
-                  *density[lev], bc_den, lev);
-        } else {
-        auto const& drdt = conv_r[lev]->array(mfi);
-        amrex::ParallelFor(bx,
+            Box const& bx = mfi.tilebox();
+
+            // velocity
+            auto const& bc_vel = get_velocity_bcrec_device_ptr();
+            redistribute_term(mfi, *conv_u[lev], dvdt_tmp,
+                              (m_advect_momentum) ? rhovel[lev] : *vel[lev],
+                              bc_vel, lev);
+
+            // density
+            if (!m_constant_density) {
+                auto const& bc_den = get_density_bcrec_device_ptr();
+                redistribute_term(mfi, *conv_r[lev], drdt_tmp,
+                                  *density[lev], bc_den, lev);
+            } else {
+                auto const& drdt = conv_r[lev]->array(mfi);
+                amrex::ParallelFor(bx,
                 [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-        {
-            drdt(i,j,k) = 0.;
-        });
-        }
+                {
+                    drdt(i,j,k) = 0.;
+                });
+            }
 
-        if (m_advect_tracer) {
-        auto const& bc_tra = get_tracer_bcrec_device_ptr();
-        redistribute_term(mfi, *conv_t[lev], dtdt_tmp,
-                  rhotrac[lev], bc_tra, lev);
-        }
-    } // mfi
+            if (m_advect_tracer) {
+                auto const& bc_tra = get_tracer_bcrec_device_ptr();
+                redistribute_term(mfi, *conv_t[lev], dtdt_tmp,
+                                  rhotrac[lev], bc_tra, lev);
+            }
+        } // mfi
 #endif
     } // lev
 }
diff --git a/src/diffusion/DiffusionScalarOp.cpp b/src/diffusion/DiffusionScalarOp.cpp
index d8250ad8..f0517f72 100644
--- a/src/diffusion/DiffusionScalarOp.cpp
+++ b/src/diffusion/DiffusionScalarOp.cpp
@@ -53,6 +53,7 @@ DiffusionScalarOp::DiffusionScalarOp (incflo* a_incflo)
                                                      m_incflo->DistributionMap(0,finest_level),
                                                      info_apply, ebfact);
             m_eb_scal_apply_op->setMaxOrder(m_mg_maxorder);
+
             m_eb_scal_apply_op->setDomainBC(m_incflo->get_diffuse_scalar_bc(Orientation::low),
                                             m_incflo->get_diffuse_scalar_bc(Orientation::high));
         }
@@ -185,6 +186,11 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& tracer,
 #ifdef AMREX_USE_EB
         if (m_eb_scal_solve_op)
         {
+            if ( m_incflo->m_has_mixedBC && comp>0 ) {
+                // Must reset scalars (and Acoef, done below) to reuse solver with Robin BC
+                m_eb_scal_solve_op->setScalars(1.0, dt);
+            }
+
             for (int lev = 0; lev <= finest_level; ++lev) {
                 Array<MultiFab,AMREX_SPACEDIM> b = m_incflo->average_scalar_eta_to_faces(lev, comp, *eta[lev]);
                 m_eb_scal_solve_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
@@ -219,7 +225,20 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& tracer,
 
 #ifdef AMREX_USE_EB
             if (m_eb_scal_solve_op) {
-                m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]);
+                if ( m_incflo->m_has_mixedBC ) {
+                    if ( comp>0 ) {
+                        // Must reset Acoef to reuse solver with Robin BC
+                        m_eb_scal_solve_op->setACoeffs(lev, *density[lev]);
+                    }
+
+                    auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]);
+
+                    m_eb_scal_solve_op->setLevelBC(lev, &phi[lev],
+                                                   &robin[0], &robin[1], &robin[2]);
+                }
+                else {
+                    m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]);
+                }
 
                 // For when we use the stencil for centroid values
                 // m_eb_scal_solve_op->setPhiOnCentroid();
@@ -366,8 +385,15 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
                 // For when we use the stencil for centroid values
                 // m_eb_vel_solve_op->setPhiOnCentroid();
 
-                m_eb_vel_solve_op->setLevelBC(lev, &phi[lev]);
+                if ( m_incflo->m_has_mixedBC ) {
+                    auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]);
 
+                    m_eb_vel_solve_op->setLevelBC(lev, &phi[lev],
+                                                  &robin[0], &robin[1], &robin[2]);
+                }
+                else {
+                    m_eb_vel_solve_op->setLevelBC(lev, &phi[lev]);
+                }
             } else
 #endif
             {
@@ -446,15 +472,15 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
         // For when we use the stencil for centroid values
         // m_eb_scal_apply_op->setPhiOnCentroid();
 
-        // This should have no effect since the first scalar is 0
-        for (int lev = 0; lev <= finest_level; ++lev) {
-            m_eb_scal_apply_op->setACoeffs(lev, *a_density[lev]);
-        }
-
         // FIXME? Can we do the solve together now?
         for (int comp = 0; comp < m_incflo->m_ntrac; ++comp) {
             int eta_comp = comp;
 
+            if ( m_incflo->m_has_mixedBC && comp>0 ){
+                // Must reset scalars to use solver with Robin BC
+                m_eb_scal_apply_op->setScalars(0.0, -1.0);
+            }
+
             Vector<MultiFab> laps_comp;
             Vector<MultiFab> scalar_comp;
             for (int lev = 0; lev <= finest_level; ++lev) {
@@ -465,7 +491,17 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
                     b = m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]);
 
                 m_eb_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);
-                m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]);
+
+                if ( m_incflo->m_has_mixedBC ) {
+
+                    auto const robin = m_incflo->make_robinBC_MFs(lev, &scalar_comp[lev]);
+
+                    m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev],
+                                                   &robin[0], &robin[1], &robin[2]);
+                }
+                else {
+                    m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]);
+                }
             }
 
             MLMG mlmg(*m_eb_scal_apply_op);
@@ -474,12 +510,12 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
 
         for(int lev = 0; lev <= finest_level; lev++)
         {
-        amrex::single_level_redistribute(laps_tmp[lev],
-                         *a_laps[lev], 0, m_incflo->m_ntrac,
+            amrex::single_level_redistribute(laps_tmp[lev],
+                                             *a_laps[lev], 0, m_incflo->m_ntrac,
                                              m_incflo->Geom(lev));
-        // auto const& bc = m_incflo->get_tracer_bcrec_device_ptr();
+            // auto const& bc = m_incflo->get_tracer_bcrec_device_ptr();
             // m_incflo->redistribute_term(*a_laps[lev], laps_tmp[lev], *a_scalar[lev],
-        //                 bc, lev);
+            //                 bc, lev);
         }
     }
     else
@@ -558,11 +594,6 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
         // For when we use the stencil for centroid values
         // m_eb_vel_apply_op->setPhiOnCentroid();
 
-        // This should have no effect since the first scalar is 0
-        for (int lev = 0; lev <= finest_level; ++lev) {
-            m_eb_vel_apply_op->setACoeffs(lev, *a_density[lev]);
-        }
-
         int eta_comp = 0;
 
         for (int comp = 0; comp < a_divtau[0]->nComp(); ++comp)
@@ -570,6 +601,11 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
             Vector<MultiFab> divtau_single;
             Vector<MultiFab>    vel_single;
 
+            if ( m_incflo->m_has_mixedBC && comp>0 ){
+                // Must reset scalars to use solver with Robin BC
+                m_eb_vel_apply_op->setScalars(0.0, -1.0);
+            }
+
             // Because the different components may have different boundary conditions, we need to
             // reset these for each solve
             m_eb_vel_apply_op->setDomainBC(m_incflo->get_diffuse_velocity_bc(Orientation::low ,comp),
@@ -577,8 +613,17 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
 
             for (int lev = 0; lev <= finest_level; ++lev) {
                 divtau_single.emplace_back(divtau_tmp[lev],amrex::make_alias,comp,1);
-                   vel_single.emplace_back(       vel[lev],amrex::make_alias,comp,1);
-                m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]);
+                vel_single.emplace_back(       vel[lev],amrex::make_alias,comp,1);
+
+                if ( m_incflo->m_has_mixedBC ) {
+                    auto const robin = m_incflo->make_robinBC_MFs(lev, &vel_single[lev]);
+
+                    m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev],
+                                                  &robin[0], &robin[1], &robin[2]);
+                }
+                else {
+                    m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]);
+                }
 
                 Array<MultiFab,AMREX_SPACEDIM> b =
                     m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]);
@@ -592,12 +637,12 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
 
         for(int lev = 0; lev <= finest_level; lev++)
         {
-        amrex::single_level_redistribute(divtau_tmp[lev],
+            amrex::single_level_redistribute(divtau_tmp[lev],
                                              *a_divtau[lev], 0, a_divtau[lev]->nComp(),
                                              m_incflo->Geom(lev));
             // auto const& bc = m_incflo->get_velocity_bcrec_device_ptr();
             // m_incflo->redistribute_term(*a_divtau[lev], divtau_tmp[lev], *a_vel[lev],
-        //                 bc, lev);
+            //                 bc, lev);
         }
     }
     else
diff --git a/src/diffusion/incflo_diffusion.cpp b/src/diffusion/incflo_diffusion.cpp
index 08d37cf4..d171d545 100644
--- a/src/diffusion/incflo_diffusion.cpp
+++ b/src/diffusion/incflo_diffusion.cpp
@@ -193,6 +193,13 @@ incflo::get_diffuse_velocity_bc (Orientation::Side side, int comp) const noexcep
                 r[dir][dir] = LinOpBCType::Dirichlet;
                 break;
             }
+            case BC::mixed:
+            {
+                AMREX_D_TERM(r[0][dir] = LinOpBCType::Robin;,
+                             r[1][dir] = LinOpBCType::Robin;,
+                             r[2][dir] = LinOpBCType::Robin;);
+                break;
+            }
             default:
                 amrex::Abort("get_diffuse_tensor_bc: undefined BC type");
             };
@@ -237,6 +244,11 @@ incflo::get_diffuse_scalar_bc (Orientation::Side side) const noexcept
                 r[dir] = LinOpBCType::Dirichlet;
                 break;
             }
+            case BC::mixed:
+            {
+                r[dir] = LinOpBCType::Robin;
+                break;
+            }
             default:
                 amrex::Abort("get_diffuse_scalar_bc: undefined BC type");
             };
diff --git a/src/embedded_boundaries/eb_twocylinders.cpp b/src/embedded_boundaries/eb_twocylinders.cpp
index ff51a345..18d30a2d 100644
--- a/src/embedded_boundaries/eb_twocylinders.cpp
+++ b/src/embedded_boundaries/eb_twocylinders.cpp
@@ -22,10 +22,12 @@ void incflo::make_eb_twocylinders()
     Real radius2 = 0.5;
     Vector<Real> centervec1(3);
     Vector<Real> centervec2(3);
+    bool inside = false;
 
     // Get information from inputs file.                               *
     ParmParse pp("twocylinders");
 
+    pp.query("internal_flow", inside);
     pp.query("direction1", direction1);
     pp.query("direction2", direction2);
     pp.query("radius1", radius1);
@@ -59,7 +61,8 @@ void incflo::make_eb_twocylinders()
     // Build the implicit function as a union of two cylinders
     EB2::CylinderIF cyl1(radius1, direction1, center1, false);
     EB2::CylinderIF cyl2(radius2, direction2, center2, false);
-    auto twocylinders = EB2::makeUnion(cyl1, cyl2);
+    auto twocylinders = inside ? EB2::makeComplement(EB2::makeUnion(cyl1, cyl2))
+                               : EB2::makeUnion(cyl1, cyl2);
 
     // Generate GeometryShop
     auto gshop = EB2::makeShop(twocylinders);
diff --git a/src/embedded_boundaries/embedded_boundaries.cpp b/src/embedded_boundaries/embedded_boundaries.cpp
index 6c7d6bde..342e2634 100644
--- a/src/embedded_boundaries/embedded_boundaries.cpp
+++ b/src/embedded_boundaries/embedded_boundaries.cpp
@@ -89,7 +89,7 @@ void incflo::MakeEBGeometry()
                    << " Will build all regular geometry." << std::endl;
         make_eb_regular();
     }
-    amrex::Print() << "Done making the geometry ebfactory.\n" << std::endl;
+    amrex::Print() << "Done making the EB geometry index space.\n" << std::endl;
 
     if (m_write_geom_chk) {
        const auto& is = amrex::EB2::IndexSpace::top();
diff --git a/src/incflo.H b/src/incflo.H
index 93264290..ddd9934e 100644
--- a/src/incflo.H
+++ b/src/incflo.H
@@ -110,6 +110,16 @@ public:
     //
     ///////////////////////////////////////////////////////////////////////////
 
+    std::unique_ptr<amrex::iMultiFab> make_BC_MF (int lev,
+                                                  amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
+                                                  std::string const& field);
+    amrex::iMultiFab make_nodalBC_mask (int lev);
+    amrex::Vector<amrex::MultiFab> make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr);
+    // void make_ccBC_mask (int lev, const amrex::BoxArray& ba,
+    //                      const amrex::DistributionMapping& dm);
+    // void make_nodalBC_mask (int lev, const amrex::BoxArray& ba,
+    //                         const amrex::DistributionMapping& dm);
+
 #ifdef AMREX_USE_EB
     void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost);
     void set_eb_density (int lev, amrex::Real time, amrex::MultiFab& eb_density, int nghost);
@@ -210,6 +220,20 @@ public:
     ///////////////////////////////////////////////////////////////////////////
 
     void prob_init_fluid (int lev);
+    void prob_set_BC_MF (amrex::Orientation const& ori, amrex::Box const& bx,
+                         amrex::Array4<int> const& mask, int lev,
+                         int inflow_val, int outflow_val, std::string const& field);
+    void prob_set_MAC_robinBCs (amrex::Orientation const& ori, amrex::Box const& bx,
+                                amrex::Array4<amrex::Real> const& robin_a,
+                                amrex::Array4<amrex::Real> const& robin_b,
+                                amrex::Array4<amrex::Real> const& robin_f,
+                                int lev);
+    void prob_set_diffusion_robinBCs (amrex::Orientation const& ori, amrex::Box const& bx,
+                                      amrex::Array4<amrex::Real> const& robin_a,
+                                      amrex::Array4<amrex::Real> const& robin_b,
+                                      amrex::Array4<amrex::Real> const& robin_f,
+                                      amrex::Array4<amrex::Real const> const& bcval,
+                                      int lev);
 
 #include "incflo_prob_I.H"
 #include "incflo_prob_usr_I.H"
@@ -563,9 +587,11 @@ private:
 
     enum struct BC {
         pressure_inflow, pressure_outflow, mass_inflow, no_slip_wall, slip_wall,
-        periodic, undefined
+        periodic, mixed, undefined
     };
 
+    bool m_has_mixedBC = false;
+
     amrex::GpuArray<BC                         , AMREX_SPACEDIM*2> m_bc_type;
     amrex::GpuArray<amrex::Real                , AMREX_SPACEDIM*2> m_bc_pressure;
     amrex::GpuArray<amrex::Real                , AMREX_SPACEDIM*2> m_bc_density;
diff --git a/src/incflo.cpp b/src/incflo.cpp
index bb691d1c..9d7e72e6 100644
--- a/src/incflo.cpp
+++ b/src/incflo.cpp
@@ -59,7 +59,7 @@ void incflo::InitData ()
         InitialRedistribution();
 #endif
 
-        if (m_do_initial_proj) {
+    if (m_do_initial_proj) {
             InitialProjection();
         }
 
@@ -211,11 +211,11 @@ void incflo::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_gr
 #endif
 
     m_leveldata[lev] = std::make_unique<LevelData>(grids[lev], dmap[lev], *m_factory[lev],
-                                         m_ntrac, nghost_state(),
-                                         m_advection_type,
-                                         m_diff_type==DiffusionType::Implicit,
-                                           use_tensor_correction,
-                                         m_advect_tracer);
+                                                   m_ntrac, nghost_state(),
+                                                   m_advection_type,
+                                                   m_diff_type==DiffusionType::Implicit,
+                                                   use_tensor_correction,
+                                                   m_advect_tracer);
 
     m_t_new[lev] = time;
     m_t_old[lev] = time - Real(1.e200);
@@ -223,6 +223,7 @@ void incflo::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_gr
     if (m_restart_file.empty()) {
         prob_init_fluid(lev);
     }
+    //make_mixedBC_mask(lev, grids[lev], dmap[lev]);
 
 #ifdef AMREX_USE_EB
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,finest_level),
diff --git a/src/incflo_apply_predictor.cpp b/src/incflo_apply_predictor.cpp
index ee57f33a..0583e9a5 100644
--- a/src/incflo_apply_predictor.cpp
+++ b/src/incflo_apply_predictor.cpp
@@ -137,7 +137,7 @@ void incflo::ApplyPredictor (bool incremental_projection)
     // *************************************************************************************
     // Compute explicit viscous term
     // *************************************************************************************
-    if (need_divtau() || use_tensor_correction)
+    if (need_divtau() || use_tensor_correction )
     {
         compute_divtau(get_divtau_old(),get_velocity_old_const(),
                        get_density_old_const(),GetVecOfConstPtrs(vel_eta));
diff --git a/src/incflo_redistribute.cpp b/src/incflo_redistribute.cpp
index e243614c..37f7f827 100644
--- a/src/incflo_redistribute.cpp
+++ b/src/incflo_redistribute.cpp
@@ -26,7 +26,7 @@ incflo::redistribute_term ( MultiFab& result,
 #endif
     for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
-    redistribute_term(mfi, result, result_tmp, state, bc, lev);
+        redistribute_term(mfi, result, result_tmp, state, bc, lev);
     }
 }
 
@@ -55,51 +55,51 @@ incflo::redistribute_term ( MFIter const& mfi,
 
     if (!regular && !covered)
     {
-    auto const& vfrac = ebfact.getVolFrac().const_array(mfi);
-    auto const& ccc   = ebfact.getCentroid().const_array(mfi);
-    AMREX_D_TERM(auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);,
-             auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);,
-             auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi););
-    AMREX_D_TERM(auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);,
-             auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);,
-             auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi););
-
-    Box gbx = bx;
-
-    if (m_redistribution_type == "StateRedist") {
-        gbx.grow(3);
-    } else if (m_redistribution_type == "FluxRedist") {
-        gbx.grow(2);
-    }
-
-    FArrayBox scratch_fab(gbx,ncomp);
-    Array4<Real> scratch = scratch_fab.array();
-    Elixir eli_scratch = scratch_fab.elixir();
-
-    // This is scratch space if calling StateRedistribute
-    //  but is used as the weights (here set to 1) if calling
-    //  FluxRedistribute
-    amrex::ParallelFor(Box(scratch),
+        auto const& vfrac = ebfact.getVolFrac().const_array(mfi);
+        auto const& ccc   = ebfact.getCentroid().const_array(mfi);
+        AMREX_D_TERM(auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);,
+                     auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);,
+                     auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi););
+        AMREX_D_TERM(auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);,
+                     auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);,
+                     auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi););
+
+        Box gbx = bx;
+
+        if (m_redistribution_type == "StateRedist") {
+            gbx.grow(3);
+        } else if (m_redistribution_type == "FluxRedist") {
+            gbx.grow(2);
+        }
+
+        FArrayBox scratch_fab(gbx,ncomp);
+        Array4<Real> scratch = scratch_fab.array();
+        Elixir eli_scratch = scratch_fab.elixir();
+
+        // This is scratch space if calling StateRedistribute
+        //  but is used as the weights (here set to 1) if calling
+        //  FluxRedistribute
+        amrex::ParallelFor(Box(scratch),
         [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
-    {
-        scratch(i,j,k) = 1.;
-    });
+        {
+            scratch(i,j,k) = 1.;
+        });
 
-    Array4<Real const> state_arr = state.const_array(mfi);
-    // State redist acts on a state. What would that be for the diffusive term??
+        Array4<Real const> state_arr = state.const_array(mfi);
+        // State redist acts on a state. What would that be for the diffusive term??
         ApplyRedistribution( bx, ncomp, out, in, state_arr,
-                          scratch, flag,
-                            AMREX_D_DECL(apx, apy, apz), vfrac,
-                      AMREX_D_DECL(fcx, fcy, fcz), ccc,
-                  bc, geom[lev], m_dt, m_redistribution_type);
+                             scratch, flag,
+                             AMREX_D_DECL(apx, apy, apz), vfrac,
+                             AMREX_D_DECL(fcx, fcy, fcz), ccc,
+                             bc, geom[lev], m_dt, m_redistribution_type);
     }
     else
     {
-    amrex::ParallelFor(bx, ncomp,
+        amrex::ParallelFor(bx, ncomp,
         [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
-    {
-        out(i,j,k,n) = in(i,j,k,n);
-    });
+        {
+            out(i,j,k,n) = in(i,j,k,n);
+        });
     }
 }
 #endif
diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp
index 867c7afb..802de0dd 100644
--- a/src/incflo_regrid.cpp
+++ b/src/incflo_regrid.cpp
@@ -43,6 +43,8 @@ void incflo::MakeNewLevelFromCoarse (int lev,
     m_leveldata[lev] = std::move(new_leveldata);
     m_factory[lev] = std::move(new_fact);
 
+    //make_mixedBC_mask(lev, ba, dm);
+
     m_diffusion_tensor_op.reset();
     m_diffusion_scalar_op.reset();
 
@@ -60,7 +62,7 @@ void incflo::MakeNewLevelFromCoarse (int lev,
 // fill with existing fine and coarse data.
 // overrides the pure virtual function in AmrCore
 void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
-             const DistributionMapping& dm)
+                          const DistributionMapping& dm)
 {
     BL_PROFILE("incflo::RemakeLevel()");
 
@@ -94,6 +96,8 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
     m_leveldata[lev] = std::move(new_leveldata);
     m_factory[lev] = std::move(new_fact);
 
+    //make_mixedBC_mask(lev, ba, dm);
+
     m_diffusion_tensor_op.reset();
     m_diffusion_scalar_op.reset();
 
diff --git a/src/prob/CMakeLists.txt b/src/prob/CMakeLists.txt
index 0b8791aa..0eae391a 100644
--- a/src/prob/CMakeLists.txt
+++ b/src/prob/CMakeLists.txt
@@ -5,6 +5,7 @@ target_sources(incflo
    incflo_prob_I.H
    incflo_prob_usr_I.H
    prob_bc.H
+   prob_bc.cpp
    prob_init_fluid.cpp
    prob_init_fluid_usr.cpp
    )
diff --git a/src/prob/Make.package b/src/prob/Make.package
index a922cb47..2885ea73 100644
--- a/src/prob/Make.package
+++ b/src/prob/Make.package
@@ -1,5 +1,6 @@
 CEXE_sources += prob_init_fluid.cpp
 CEXE_sources += prob_init_fluid_usr.cpp
+CEXE_sources += prob_bc.cpp
 
 CEXE_headers += prob_bc.H
 CEXE_headers += incflo_prob_I.H
diff --git a/src/prob/incflo_prob_I.H b/src/prob/incflo_prob_I.H
index e9784122..d72e2b7f 100644
--- a/src/prob/incflo_prob_I.H
+++ b/src/prob/incflo_prob_I.H
@@ -146,4 +146,13 @@
                                amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
                                amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
                                amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi);
+
+    void init_jump (amrex::Box const& vbx, amrex::Box const& gbx,
+                    amrex::Array4<amrex::Real> const& vel,
+                    amrex::Array4<amrex::Real> const& density,
+                    amrex::Array4<amrex::Real> const& tracer,
+                    amrex::Box const& domain,
+                    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
+                    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
+                    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi) const;
 #endif
diff --git a/src/prob/prob_bc.H b/src/prob/prob_bc.H
index 7049f2a5..8d54d20f 100644
--- a/src/prob/prob_bc.H
+++ b/src/prob/prob_bc.H
@@ -24,6 +24,8 @@ struct IncfloVelFill
                      const int orig_comp) const
     {
         // do something for external Dirichlet (amrex::BCType::ext_dir)
+        // and for the Dirichlet part of mixed BCs - which gets FOEXTRAP so the
+        // Neumann part is already filled before arriving here during FillPatch.
         const int i = iv[0];
         const int j = iv[1];
 #if (AMREX_SPACEDIM == 3)
@@ -38,7 +40,36 @@ struct IncfloVelFill
         {
             const amrex::BCRec& bc = bcr[bcomp+nc];
 
-            if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
+            if (1101 == probtype && i < domain_box.smallEnd(0))
+            {
+                int direction = 1;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (j > half_num_cells) {
+                    // Here is the dirichlet portion
+                    vel(i,j,k,dcomp+nc) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][orig_comp+nc];
+                }
+            }
+            else if (1101 == probtype && i > domain_box.bigEnd(0))
+            {
+                int direction = 1;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (j <= half_num_cells) {
+                    // Here we take minus the inflow BC
+                    vel(i,j,k,dcomp+nc) = -bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][orig_comp+nc];
+                }
+            }
+#if (AMREX_SPACEDIM == 3)
+            else if (1102 == probtype && j > domain_box.bigEnd(1))
+            {
+                int direction = 2;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (k <= half_num_cells) {
+                    // Here we take minus the inflow BC specified in inputs file
+                    vel(i,j,k,dcomp+nc) = -bcv[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][orig_comp+nc];
+                }
+            }
+#endif
+            else if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
             {
                 vel(i,j,k,dcomp+nc) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][orig_comp+nc];
 
@@ -160,7 +191,35 @@ struct IncfloDenFill
         const Box& domain_box = geom.Domain();
         const BCRec& bc = bcr[bcomp];
 
-        if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
+        if (1101 == probtype && i < domain_box.smallEnd(0))
+        {
+            // this probtype uses mixed BC on x faces
+            int direction = 1;
+            int half_num_cells = domain_box.length(direction) / 2;
+            if (j > half_num_cells) {
+                // Here is the dirichlet portion; we take the inflow BC specified in inputs file
+                rho(i,j,k) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
+            }
+        }
+        else if (1101 == probtype && i > domain_box.bigEnd(0))
+        {
+            int direction = 1;
+            int half_num_cells = domain_box.length(direction) / 2;
+            if (j <= half_num_cells) {
+                rho(i,j,k) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)];
+            }
+        }
+#if (AMREX_SPACEDIM == 3)
+        else if (1102 == probtype && j > domain_box.bigEnd(1))
+        {
+            int direction = 2;
+            int half_num_cells = domain_box.length(direction) / 2;
+            if (k <= half_num_cells) {
+                rho(i,j,k) = bcv[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)];
+            }
+        }
+#endif
+        else if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
         {
             rho(i,j,k) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)];
         }
@@ -224,7 +283,34 @@ struct IncfloTracFill
         {
             const BCRec& bc = bcr[bcomp+n];
 
-            if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
+            if (1101 == probtype && i < domain_box.smallEnd(0))
+            {
+                int direction = 1;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (j > half_num_cells) {
+                    // Here we take the inflow BC specified in inputs file
+                    tracer(i,j,k,n) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][n];
+                }
+            }
+            else if (1101 == probtype && i > domain_box.bigEnd(0))
+            {
+                int direction = 1;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (j <= half_num_cells) {
+                    tracer(i,j,k,n) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)][n];
+                }
+            }
+#if (AMREX_SPACEDIM == 3)
+            else if (1102 == probtype && j > domain_box.bigEnd(1))
+            {
+                int direction = 2;
+                int half_num_cells = domain_box.length(direction) / 2;
+                if (k <= half_num_cells) {
+                    tracer(i,j,k,n) = bcv[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)][n];
+                }
+            }
+#endif
+            else if (bc.lo(0) == amrex::BCType::ext_dir && i < domain_box.smallEnd(0))
             {
                 tracer(i,j,k,n) = bcv[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)][n];
             }
diff --git a/src/prob/prob_bc.cpp b/src/prob/prob_bc.cpp
new file mode 100644
index 00000000..cc3fe345
--- /dev/null
+++ b/src/prob/prob_bc.cpp
@@ -0,0 +1,298 @@
+#include <incflo.H>
+
+using namespace amrex;
+
+// Specify mixed BCs for both the nodal projection and advection. The necessary value
+// to indicate an in/outflow BC differs between advection and the NodalProj, so we pass it.
+// Note that the advection BCs are on the velocity and scalars, whereas for the Nodal
+// Projection we're creating an overset mask that operates on the solve for phi (~pressure)
+void incflo::prob_set_BC_MF (Orientation const& ori, Box const& bx,
+                             Array4<int> const& mask, int lev,
+                             int inflow_val, int outflow_val,
+                             std::string const& field)
+{
+    if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype)
+    {
+        int direction = 0;
+        if (1101 == m_probtype) {
+            direction = 1;
+        }
+        else if (1102 == m_probtype) {
+            direction = 2;
+        }
+        Box const& domain = geom[lev].Domain();
+        int half_num_cells  = domain.length(direction) / 2;
+
+        // for this problem, bcs are same for all fields, only ncomp varies
+        int ncomp =  field == "velocity" ? AMREX_SPACEDIM : 1;
+
+        Orientation::Side side = ori.faceDir();
+        if (side == Orientation::low) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                for ( int n = 0; n < ncomp; n++ ){
+                    if (direction == 0) {
+                        if (i <= half_num_cells) {
+                            mask(i,j,k,n) = outflow_val; // outflow on bottom
+                        } else {
+                            mask(i,j,k,n) = inflow_val; // inflow on top
+                        }
+                    }
+                    else if (direction == 1) {
+                        if (j <= half_num_cells) {
+                            mask(i,j,k,n) = outflow_val;
+                        } else {
+                            mask(i,j,k,n) = inflow_val;
+                        }
+                    }
+                    else if (direction == 2) {
+                        if (k <= half_num_cells) {
+                            mask(i,j,k,n) = outflow_val;
+                        } else {
+                            mask(i,j,k,n) = inflow_val;
+                        }
+                    }
+                }
+            });
+        } else {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                for ( int n = 0; n < ncomp; n++ ){
+                    if (direction == 0) {
+                        if (i > half_num_cells) {
+                            mask(i,j,k,n) = outflow_val; // outflow on top
+                        } else {
+                            mask(i,j,k,n) = inflow_val;  // inflow on bottom
+                        }
+                    }
+                    else if (direction == 1) {
+                        if (j > half_num_cells) {
+                            mask(i,j,k,n) = outflow_val;
+                        } else {
+                            mask(i,j,k,n) = inflow_val;
+                        }
+                    }
+                    else if (direction == 2) {
+                        if (k > half_num_cells) {
+                            mask(i,j,k,n) = outflow_val;
+                        } else {
+                            mask(i,j,k,n) = inflow_val;
+                        }
+                    }
+                }
+            });
+        }
+    }
+    else
+    {
+        Abort("incflo::prob_set_BC_MF: No masking function for probtype "
+              +std::to_string(m_probtype));
+    }
+}
+
+void incflo::prob_set_MAC_robinBCs (Orientation const& ori, Box const& bx,
+                                    Array4<Real> const& robin_a,
+                                    Array4<Real> const& robin_b,
+                                    Array4<Real> const& robin_f,
+                                    int lev)
+{
+    // For the MAC Projection, the BC is on phi (~pressure), so always homogeneous
+    // Robin BC:   a u + b du/dn = f  -- inflow,  Neumann   a=0, b=1, f=0
+    //                                -- outflow, Dirichlet a=1, b=0, f=0
+
+    if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype)
+    {
+        int direction = 0;
+        if (1101 == m_probtype) {
+            direction = 1;
+        }
+        else if (1102 == m_probtype) {
+            direction = 2;
+        }
+        Box const& domain = geom[lev].Domain();
+        int half_num_cells  = domain.length(direction) / 2;
+
+        Orientation::Side side = ori.faceDir();
+        if (side == Orientation::low) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                robin_f(i,j,k) = 0.;
+
+                if (direction == 0) {
+                    if (i <= half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0; // outflow on bottom
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.; // inflow on top
+                    }
+                }
+                else if (direction == 1) {
+                    if (j <= half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0;
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.;
+                    }
+                }
+                else if (direction == 2) {
+                    if (k <= half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0;
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.;
+                    }
+                }
+            });
+        } else {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                robin_f(i,j,k) = 0.;
+
+                if (direction == 0) {
+                    if (i > half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0; // outflow on top
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.;  // inflow on bottom
+                    }
+                }
+                else if (direction == 1) {
+                    if (j > half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0;
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.;
+                    }
+                }
+                else if (direction == 2) {
+                    if (k > half_num_cells) {
+                        robin_a(i,j,k) = 1.;
+                        robin_b(i,j,k) = 0;
+                    } else {
+                        robin_a(i,j,k) = 0.;
+                        robin_b(i,j,k) = 1.;
+                    }
+                }
+            });
+        }
+    }
+    else
+    {
+        Abort("incflo::prob_set_MAC_robinBCs: No masking function for probtype "
+              +std::to_string(m_probtype));
+    }
+}
+
+void incflo::prob_set_diffusion_robinBCs (Orientation const& ori, Box const& bx,
+                                          Array4<Real> const& robin_a,
+                                          Array4<Real> const& robin_b,
+                                          Array4<Real> const& robin_f,
+                                          Array4<Real const> const& bcval,
+                                          int lev)
+{
+    // For diffusion, we also pass in the dirichlet bc (bcval)
+    // Robin BC:   a u + b du/dn = f  -- inflow,  Dirichlet a=1, b=0, f=bcval
+    //                                -- outflow, Neumann   a=0, b=1, f=0
+
+    if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype)
+    {
+        int direction = 0;
+        if (1101 == m_probtype) {
+            direction = 1;
+        }
+        else if (1102 == m_probtype) {
+            direction = 2;
+        }
+        Box const& domain = geom[lev].Domain();
+        int half_num_cells  = domain.length(direction) / 2;
+
+        Orientation::Side side = ori.faceDir();
+        if (side == Orientation::low) {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (direction == 0) {
+                    if (i <= half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.; // outflow on bottom
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.; // inflow on top
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+                else if (direction == 1) {
+                    if (j <= half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.;
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.;
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+                else if (direction == 2) {
+                    if (k <= half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.;
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.;
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+            });
+        } else {
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                robin_f(i,j,k,0) = 0.;
+
+                if (direction == 0) {
+                    if (i > half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.; // outflow on top
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.;  // inflow on bottom
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+                else if (direction == 1) {
+                    if (j > half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.;
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.;
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+                else if (direction == 2) {
+                    if (k > half_num_cells) {
+                        robin_a(i,j,k,0) = 0.;
+                        robin_b(i,j,k,0) = 1.;
+                        robin_f(i,j,k,0) = 0.;
+                    } else {
+                        robin_a(i,j,k,0) = 1.;
+                        robin_b(i,j,k,0) = 0.;
+                        robin_f(i,j,k,0) = bcval(i,j,k,0);
+                    }
+                }
+            });
+        }
+    }
+    else
+    {
+        Abort("incflo::prob_set_diffusion_robinBCs: No masking function for probtype "
+              +std::to_string(m_probtype));
+    }
+}
diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp
index 3ce284d4..b0bfc323 100644
--- a/src/prob/prob_init_fluid.cpp
+++ b/src/prob/prob_init_fluid.cpp
@@ -84,6 +84,14 @@ void incflo::prob_init_fluid (int lev)
                         ld.tracer.array(mfi),
                         domain, dx, problo, probhi);
         }
+        else if (1100 == m_probtype || 1101 == m_probtype || 1102 == m_probtype)
+        {
+            init_jump(vbx, gbx,
+                      ld.velocity.array(mfi),
+                      ld.density.array(mfi),
+                      ld.tracer.array(mfi),
+                      domain, dx, problo, probhi);
+        }
         else if (111 == m_probtype || 112 == m_probtype || 113 == m_probtype)
         {
             init_boussinesq_bubble(vbx, gbx,
@@ -607,6 +615,7 @@ void incflo::init_tuscan (Box const& vbx, Box const& /*gbx*/,
                           GpuArray<Real, AMREX_SPACEDIM> const& /*dx*/,
                           GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
                           GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/)
+
 {
     int half_num_cells = domain.length(AMREX_SPACEDIM-1) / 2;
     amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
@@ -619,7 +628,45 @@ void incflo::init_tuscan (Box const& vbx, Box const& /*gbx*/,
             tracer(i,j,k) = Real(0.0);
         } else {
             tracer(i,j,k) = Real(0.01);
-         }
+        }
+    });
+}
+
+void incflo::init_jump (Box const& vbx, Box const& /*gbx*/,
+                        Array4<Real> const& vel,
+                        Array4<Real> const& /*density*/,
+                        Array4<Real> const& /*tracer*/,
+                        Box const& domain,
+                        GpuArray<Real, AMREX_SPACEDIM> const& /*dx*/,
+                        GpuArray<Real, AMREX_SPACEDIM> const& /*problo*/,
+                        GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/) const
+{
+    int direction = 0;
+    if (1101 == m_probtype) {
+        direction = 1;
+    }
+    else if (1102 == m_probtype) {
+        direction = 2;
+    }
+
+    int half_num_cells = domain.length(direction) / 2;
+    amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+        if (direction == 0) {
+            if (i <= half_num_cells) {
+                vel(i,j,k,2) = -vel(i,j,k,2);
+            }
+        }
+        else if (direction == 1) {
+            if (j <= half_num_cells) {
+                vel(i,j,k,0) = -vel(i,j,k,0);
+            }
+        }
+        else if (direction == 2) {
+            if (k <= half_num_cells) {
+                vel(i,j,k,1) = -vel(i,j,k,1);
+            }
+        }
     });
 }
 
diff --git a/src/projection/incflo_apply_nodal_projection.cpp b/src/projection/incflo_apply_nodal_projection.cpp
index 71044a18..cb061be5 100644
--- a/src/projection/incflo_apply_nodal_projection.cpp
+++ b/src/projection/incflo_apply_nodal_projection.cpp
@@ -125,8 +125,6 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
                 (geom[lev], inflow_bcr, IncfloVelFill{m_probtype, m_bc_velocity});
             physbc(*vel[lev], 0, AMREX_SPACEDIM, nghost, time, 0);
 
-            // Note that we wouldn't need this if we trusted users to supply inflow data
-            // that obeys periodicity
             // We make sure to only fill "nghost" ghost cells so we don't accidentally
             // over-write good ghost cell values with unfilled ghost cell values
             vel[lev]->EnforcePeriodicity(0, AMREX_SPACEDIM, nghost, geom[lev].periodicity());
@@ -154,6 +152,17 @@ void incflo::ApplyNodalProjection (Vector<MultiFab const*> density,
           nodal_projector->getLinOp().setEBInflowVelocity(lev, *get_velocity_eb()[lev]);
        }
     }
+
+    if(m_has_mixedBC)
+    {
+        // We use an overset mask to effectively apply the mixed BC
+        for(int lev = 0; lev <= finest_level; ++lev)
+        {
+            const auto mask = make_nodalBC_mask(lev);
+            nodal_projector->getLinOp().setOversetMask(lev, mask);
+        }
+    }
+
 #endif
 
     nodal_projector->project(m_nodal_mg_rtol, m_nodal_mg_atol);
diff --git a/src/projection/incflo_projection_bc.cpp b/src/projection/incflo_projection_bc.cpp
index de4bd2de..0bd24ccd 100644
--- a/src/projection/incflo_projection_bc.cpp
+++ b/src/projection/incflo_projection_bc.cpp
@@ -20,6 +20,7 @@ incflo::get_nodal_projection_bc (Orientation::Side side) const noexcept
                 break;
             }
             case BC::mass_inflow:
+            case BC::mixed:
             {
                 r[dir] = LinOpBCType::inflow;
                 break;
@@ -62,6 +63,11 @@ incflo::get_mac_projection_bc (Orientation::Side side) const noexcept
                 r[dir] = LinOpBCType::Neumann;
                 break;
             }
+            case BC::mixed:
+            {
+                r[dir] = LinOpBCType::Robin;
+                break;
+            }
             default:
                 amrex::Abort("get_mac_projection_bc: undefined BC type");
             };
diff --git a/src/setup/init.cpp b/src/setup/init.cpp
index 7385f1a5..562ebaca 100644
--- a/src/setup/init.cpp
+++ b/src/setup/init.cpp
@@ -375,21 +375,18 @@ incflo::InitialRedistribution ()
         // We also need any physical boundary conditions imposed if we are
         //    calling state redistribution (because that calls the slope routine)
 
-        EB_set_covered(ld.velocity, 0, AMREX_SPACEDIM, ld.velocity.nGrow(), 0.0);
         ld.velocity.FillBoundary(geom[lev].periodicity());
         MultiFab::Copy(ld.velocity_o, ld.velocity, 0, 0, AMREX_SPACEDIM, ld.velocity.nGrow());
         fillpatch_velocity(lev, m_t_new[lev], ld.velocity_o, 3);
 
         if (!m_constant_density)
         {
-            EB_set_covered(ld.density, 0, 1, ld.density.nGrow(), 0.0);
             ld.density.FillBoundary(geom[lev].periodicity());
             MultiFab::Copy(ld.density_o, ld.density, 0, 0, 1, ld.density.nGrow());
             fillpatch_density(lev, m_t_new[lev], ld.density_o, 3);
         }
         if (m_advect_tracer)
         {
-            EB_set_covered(ld.tracer, 0, m_ntrac, ld.tracer.nGrow(), 0.0);
             ld.tracer.FillBoundary(geom[lev].periodicity());
             MultiFab::Copy(ld.tracer_o, ld.tracer, 0, 0, 1, ld.tracer.nGrow());
             fillpatch_tracer(lev, m_t_new[lev], ld.tracer_o, 3);
diff --git a/test_2d/benchmark.split b/test_2d/benchmark.split
new file mode 100644
index 00000000..b3517c26
--- /dev/null
+++ b/test_2d/benchmark.split
@@ -0,0 +1,93 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+max_step                =  50           # Max number of time steps
+steady_state            =  0            # Steady-state solver? 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+incflo.cfl = 0.45
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+amr.plot_int            =   100         # Steps between plot files
+amr.plt_regtest         =   1
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity          =   0.  0.  0.  # Gravitational force (3D)
+incflo.ro_0             =   1.          # Reference density 
+
+incflo.fluid_model      =   "newtonian" # Fluid model (rheology)
+incflo.mu               =   0.01        # Dynamic viscosity coefficient
+incflo.use_tensor_solve =   false
+incflo.advect_momentum  =   true
+incflo.constant_density =   false
+
+incflo.advect_tracer    =   true
+incflo.ntrac            =   1           # Number of tracers
+incflo.mu_s             =   0.001       # Scalar diffusion coefficient
+
+incflo.advection_type        =  "Godunov" 
+incflo.diffusion_type        =  1       # Crank-Nicolson
+incflo.redistribution_type   =  "StateRedist" 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              =   32 32 8     # Grid cells at coarsest AMRlevel
+amr.max_level           =   1           # Max AMR level in hierarchy 
+amr.max_grid_size_x     =   1024
+amr.max_grid_size_y     =   1024
+amr.max_grid_size_z     =   1024
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =  -1. -1. -.25 # Lo corner coordinates
+geometry.prob_hi        =   1.  1.  .25 # Hi corner coordinates
+geometry.is_periodic    =   0   0   1   # Periodicity x y z (0/1)
+
+# Boundary conditions
+xlo.type                =   "mixed"
+xlo.velocity            =   1.0  0.0  0.0
+xlo.pressure            =   0.0
+xlo.density             =   2.0
+xlo.tracer              =   2.0
+
+xhi.type                =   "mixed"
+xhi.velocity            =   1.0  0.0  0.0
+xhi.pressure            =   0.0
+xhi.density             =   2.0
+xhi.tracer              =   2.0
+
+ylo.type                =   "nsw"
+yhi.type                =   "nsw"
+
+# Add box
+incflo.geometry = "box"
+box.internal_flow = false
+
+box.Lo = -0.4  -0.2  -1.1
+box.Hi =  1.1   0.2   1.1
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#           INITIAL CONDITIONS          #
+#.......................................#
+incflo.probtype = 1101
+
+incflo.ic_u             =   0.0         #
+incflo.ic_v             =   0.0         #
+incflo.ic_w             =   0.0         #
+incflo.ic_p             =   0.0         #
+incflo.ic_t             =   0.0         #
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              VERBOSITY                #
+#.......................................#
+incflo.verbose            =   1         # incflo itself
+mac_proj.verbose          =   1         # MAC Projector
+nodal_proj.verbose        =   1         # Nodal Projector
diff --git a/test_3d/benchmark.upipe b/test_3d/benchmark.upipe
new file mode 100644
index 00000000..a87b33d1
--- /dev/null
+++ b/test_3d/benchmark.upipe
@@ -0,0 +1,89 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+max_step                =  50           # Max number of time steps
+steady_state            =  0            # Steady-state solver? 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+incflo.cfl = 0.45
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+amr.plot_int            =   100         # Steps between plot files
+amr.plt_regtest         =   1
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity          =   0.  0.  1.  # Gravitational force (3D)
+incflo.ro_0             =   1.          # Reference density 
+
+incflo.fluid_model      =   "newtonian" # Fluid model (rheology)
+incflo.mu               =   0.001       # Dynamic viscosity coefficient
+incflo.use_tensor_solve =   false
+incflo.advect_momentum  =   false
+incflo.constant_density =   false       #
+
+incflo.advect_tracer    =   true
+incflo.ntrac            =   1           # Number of tracers
+incflo.mu_s             =   0.001       # Scalar diffusion coefficient
+
+incflo.advection_type        =  "Godunov" 
+incflo.diffusion_type        =  1       # Crank-Nicolson
+incflo.redistribution_type   =  "StateRedist" 
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              =   8 32 32     # Grid cells at coarsest AMRlevel
+amr.max_level           =   1           # Max AMR level in hierarchy 
+amr.max_grid_size_x     =   1024
+amr.max_grid_size_y     =   1024
+amr.max_grid_size_z     =   1024
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        = -.25  -1. -1. # Lo corner coordinates
+geometry.prob_hi        =  .25   1.  1. # Hi corner coordinates
+geometry.is_periodic    =    1   0   0  # Periodicity x y z (0/1)
+
+# Boundary conditions
+ylo.type                =   "nsw"
+
+yhi.type                =   "mixed"
+yhi.velocity            =   0.0  1.0  0.0
+yhi.pressure            =   0.0
+yhi.density             =   2.0
+yhi.tracer              =   2.0
+
+zlo.type                =   "nsw"
+zhi.type                =   "nsw"
+
+# Add box
+incflo.geometry = "box"
+box.internal_flow = false
+
+box.Lo = -1.1  -0.4  -0.2
+box.Hi =  1.1   1.1   0.2
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#           INITIAL CONDITIONS          #
+#.......................................#
+incflo.probtype = 1102
+
+incflo.ic_u             =   0.0         #
+incflo.ic_v             =   0.0         #
+incflo.ic_w             =   0.0         #
+incflo.ic_p             =   0.0         #
+incflo.ic_t             =   0.0         #
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              VERBOSITY                #
+#.......................................#
+incflo.verbose            =   1         # incflo itself
+mac_proj.verbose          =   1         # MAC Projector
+nodal_proj.verbose        =   1         # Nodal Projector