diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
index b871ee5305d..d354c484ff0 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensorOp.cpp
@@ -420,6 +420,7 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
     const Box& domain = geom.Domain();
     const auto dlo = amrex::lbound(domain);
     const auto dhi = amrex::ubound(domain);
+    const GpuArray<int,AMREX_SPACEDIM>& is_periodic = geom.isPeriodicArray();
     const auto& bcondloc = *m_bcondloc[amrlev][mglev];
 
     Array<MultiFab,AMREX_SPACEDIM> const& etamf = m_b_coeffs[amrlev][mglev];
@@ -494,15 +495,15 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
 	      AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM
 	      ( xbx, txbx,
 		{
-		  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct);
+		  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct,is_periodic[0]);
 		}
 		, ybx, tybx,
 		{
-		  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct);
+		  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct,is_periodic[1]);
 		}
 		, zbx, tzbx,
 		{
-		  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct);
+		  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct,is_periodic[2]);
 		}
 	      );
 	  }
@@ -516,15 +517,15 @@ MLEBTensorOp::compCrossTerms(int amrlev, int mglev, MultiFab const& mf) const
 	    AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM
 	    ( xbx, txbx,
 	      {
-		mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv,dlo,dhi,bct);
+		mlebtensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,apx,flag,dxinv,dlo,dhi,bct,is_periodic[0]);
 	      }
 	      , ybx, tybx,
 	      {
-		mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv,dlo,dhi,bct);
+		mlebtensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,apy,flag,dxinv,dlo,dhi,bct,is_periodic[1]);
 	      }
 	      , zbx, tzbx,
 	      {
-		mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv,dlo,dhi,bct);
+		mlebtensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,apz,flag,dxinv,dlo,dhi,bct,is_periodic[2]);
 	      }
 	      );
 	  }
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
index 419330db8bb..cd8fc9aa355 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_2D_K.H
@@ -21,7 +21,8 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                                 Array4<EBCellFlag const> const& flag,
                                 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                                 const Dim3& dlo, const Dim3& dhi,
-                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                                const int& is_periodic) noexcept
 {
     const Real dyi = dxinv[1];
     const auto lo = amrex::lbound(box);
@@ -52,7 +53,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                 Real dvdy = (0.5*dyi) * ((vel(i  ,jhip,0,1)-vel(i  ,jhim,0,1))*whi
                                         +(vel(i-1,jlop,0,1)-vel(i-1,jlom,0,1))*wlo);
 
-                if (i==dhi.x+1){
+                if (i==dhi.x+1 and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudy = (dyi) * ((vel(i-1,jlop,0,0)-vel(i-1,jlom,0,0))*wlo);
                     }
@@ -66,7 +67,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                         dvdy = (dyi) * ((vel(i  ,jhip,0,1)-vel(i  ,jhim,0,1))*whi);
                     }
                 }
-                if (i==dlo.x){
+                if (i==dlo.x and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudy = (dyi) * ((vel(i  ,jhip,0,0)-vel(i  ,jhim,0,0))*whi);
                     }
@@ -101,7 +102,8 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                                 Array4<EBCellFlag const> const& flag,
                                 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                                 const Dim3& dlo, const Dim3& dhi,
-                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                                const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const auto lo = amrex::lbound(box);
@@ -131,7 +133,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                 Real dvdx = (0.5*dxi) * ((vel(ihip,j  ,0,1)-vel(ihim,j  ,0,1))*whi
                                         +(vel(ilop,j-1,0,1)-vel(ilom,j-1,0,1))*wlo);
 
-                if (j==dhi.y+1){
+                if (j==dhi.y+1 and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudx = (dxi) * ((vel(ilop,j-1,0,0)-vel(ilom,j-1,0,0))*wlo);
                     }
@@ -145,7 +147,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                         dvdx = (dxi) * ((vel(ihip,j  ,0,1)-vel(ihim,j  ,0,1))*whi);
                     }
                 }
-                if (j==dlo.y){
+                if (j==dlo.y and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudx = (dxi) * ((vel(ihip,j  ,0,0)-vel(ihim,j  ,0,0))*whi);
                     }
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H
index 6afba949892..65b0f45e760 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLEBTensor_3D_K.H
@@ -21,7 +21,8 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                                 Array4<EBCellFlag const> const& flag,
                                 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                                 const Dim3& dlo, const Dim3& dhi,
-                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                                const int& is_periodic) noexcept
 {
     const Real dyi = dxinv[1];
     const Real dzi = dxinv[2];
@@ -55,7 +56,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                     Real dvdy = (0.5*dyi) * ((vel(i  ,jhip,k,1)-vel(i  ,jhim,k,1))*whi
                                             +(vel(i-1,jlop,k,1)-vel(i-1,jlom,k,1))*wlo);
 
-                    if (i==dhi.x+1){
+                    if (i==dhi.x+1 and !is_periodic){
                         if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                             dudy = (dyi) * ((vel(i-1,jlop,k,0)-vel(i-1,jlom,k,0))*wlo);
                         }
@@ -69,7 +70,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                             dvdy = (dyi) * ((vel(i  ,jhip,k,1)-vel(i  ,jhim,k,1))*whi);
                         }
                     }
-                    if (i==dlo.x){
+                    if (i==dlo.x and !is_periodic){
                         if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                             dudy = (dyi) * ((vel(i  ,jhip,k,0)-vel(i  ,jhim,k,0))*whi);
                         }
@@ -97,7 +98,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                     Real dwdz = (0.5*dzi) * ((vel(i  ,j,khip,2)-vel(i  ,j,khim,2))*whi
                                             +(vel(i-1,j,klop,2)-vel(i-1,j,klom,2))*wlo);
 
-                    if (i==dhi.x+1){
+                    if (i==dhi.x+1 and !is_periodic){
                         if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                             dudz = (dzi) * ((vel(i-1,j,klop,0)-vel(i-1,j,klom,0))*wlo);
                         }
@@ -111,7 +112,7 @@ void mlebtensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                             dwdz = (dzi) * ((vel(i  ,j,khip,2)-vel(i  ,j,khim,2))*whi);
                         }
                     }
-                    if (i==dlo.x){
+                    if (i==dlo.x and !is_periodic){
                         if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                             dudz = (dzi) * ((vel(i  ,j,khip,0)-vel(i  ,j,khim,0))*whi);
                         }
@@ -148,7 +149,8 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                                 Array4<EBCellFlag const> const& flag,
                                 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                                 const Dim3& dlo, const Dim3& dhi,
-                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                                const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const Real dzi = dxinv[2];
@@ -181,7 +183,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                     Real dvdx = (0.5*dxi) * ((vel(ihip,j  ,k,1)-vel(ihim,j  ,k,1))*whi
                                             +(vel(ilop,j-1,k,1)-vel(ilom,j-1,k,1))*wlo);
 
-                    if (j==dhi.y+1){ 
+                    if (j==dhi.y+1 and !is_periodic){ 
                         if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){
                             dudx = (dxi) * ((vel(ilop,j-1,k,0)-vel(ilom,j-1,k,0))*wlo);
                         }
@@ -195,7 +197,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                             dvdx = (dxi) * ((vel(ihip,j  ,k,1)-vel(ihim,j  ,k,1))*whi);
                         }
                     }
-                    if (j==dlo.y){
+                    if (j==dlo.y and !is_periodic){
                         if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){
                             dudx = (dxi) * ((vel(ihip,j  ,k,0)-vel(ihim,j  ,k,0))*whi);
                         }
@@ -223,7 +225,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                     Real dwdz = (0.5*dzi) * ((vel(i,j  ,khip,2)-vel(i,j  ,khim,2))*whi
                                             +(vel(i,j-1,klop,2)-vel(i,j-1,klom,2))*wlo);
 
-                    if (j==dhi.y+1){
+                    if (j==dhi.y+1 and !is_periodic){
                         if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){
                             dvdz = (dzi) * ((vel(i,j-1,klop,1)-vel(i,j-1,klom,1))*wlo);
                         }
@@ -237,7 +239,7 @@ void mlebtensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                             dwdz = (dzi) * ((vel(i,j  ,khip,2)-vel(i,j  ,khim,2))*whi);
                         }
                     }
-                    if (j==dlo.y){
+                    if (j==dlo.y and !is_periodic){
                         if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){
                             dvdz = (dzi) * ((vel(i,j  ,khip,1)-vel(i,j  ,khim,1))*whi);
                         }
@@ -274,7 +276,8 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                                 Array4<EBCellFlag const> const& flag,
                                 GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                                 const Dim3& dlo, const Dim3& dhi,
-                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                                Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                                const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const Real dyi = dxinv[1];
@@ -307,7 +310,7 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                     Real dwdx = (0.5*dxi) * ((vel(ihip,j,k  ,2)-vel(ihim,j,k  ,2))*whi
                                             +(vel(ilop,j,k-1,2)-vel(ilom,j,k-1,2))*wlo);
 
-                    if (k==dhi.z+1){
+                    if (k==dhi.z+1 and !is_periodic){
                         if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_NEUMANN){
                             dudx = (dxi) * ((vel(ilop,j,k-1,0)-vel(ilom,j,k-1,0))*wlo);
                         }
@@ -321,7 +324,7 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                             dwdx = (dxi) * ((vel(ihip,j,k  ,2)-vel(ihim,j,k  ,2))*whi);
                         }
                     }
-                    if (k==dlo.z){
+                    if (k==dlo.z and !is_periodic){
                         if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_NEUMANN){
                             dudx = (dxi) * ((vel(ihip,j,k  ,0)-vel(ihim,j,k  ,0))*whi);
                         }
@@ -349,7 +352,7 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                     Real dwdy = (0.5*dyi) * ((vel(i,jhip,k  ,2)-vel(i,jhim,k  ,2))*whi
                                             +(vel(i,jlop,k-1,2)-vel(i,jlom,k-1,2))*wlo);
 
-                    if (k==dhi.z+1){
+                    if (k==dhi.z+1 and !is_periodic){
                         if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_NEUMANN){
                             dvdy = (dyi) * ((vel(i,jlop,k-1,1)-vel(i,jlom,k-1,1))*wlo);
                         }
@@ -363,7 +366,7 @@ void mlebtensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                             dwdy = (dyi) * ((vel(i,jhip,k  ,2)-vel(i,jhim,k  ,2))*whi);
                         }
                     }
-                    if (k==dlo.z){
+                    if (k==dlo.z and !is_periodic){
                         if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_NEUMANN){
                             dvdy = (dyi) * ((vel(i,jhip,k  ,1)-vel(i,jhim,k  ,1))*whi);
                         }
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp
index 396069dae33..9e19196c2db 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp
+++ b/Src/LinearSolvers/MLMG/AMReX_MLTensorOp.cpp
@@ -215,6 +215,7 @@ MLTensorOp::apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc
     const Box& domain = geom.Domain();
     const auto dlo = amrex::lbound(domain);
     const auto dhi = amrex::ubound(domain);
+    const GpuArray<int,AMREX_SPACEDIM>& is_periodic = geom.isPeriodicArray();
     const auto& bcondloc = *m_bcondloc[amrlev][mglev];
 
     Array<MultiFab,AMREX_SPACEDIM> const& etamf = m_b_coeffs[amrlev][mglev];
@@ -262,15 +263,15 @@ MLTensorOp::apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc
             AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM
             ( xbx, txbx,
               {
-                  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct,is_periodic[0]);
               }
             , ybx, tybx,
               {
-                  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct,is_periodic[1]);
               }
             , zbx, tzbx,
               {
-                  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct,is_periodic[2]);
               }
             );
 
@@ -421,6 +422,7 @@ MLTensorOp::compFlux (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& fluxes,
     const Box& domain = geom.Domain();
     const auto dlo = amrex::lbound(domain);
     const auto dhi = amrex::ubound(domain);
+    const GpuArray<int,AMREX_SPACEDIM>& is_periodic = geom.isPeriodicArray();
     const auto& bcondloc = *m_bcondloc[amrlev][mglev];
 
     Array<MultiFab,AMREX_SPACEDIM> const& etamf = m_b_coeffs[amrlev][mglev];
@@ -467,15 +469,15 @@ MLTensorOp::compFlux (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& fluxes,
             AMREX_LAUNCH_HOST_DEVICE_LAMBDA_DIM
             ( xbx, txbx,
               {
-                  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fx(txbx,fxfab,vfab,etaxfab,kapxfab,dxinv,dlo,dhi,bct,is_periodic[0]);
               }
             , ybx, tybx,
               {
-                  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fy(tybx,fyfab,vfab,etayfab,kapyfab,dxinv,dlo,dhi,bct,is_periodic[1]);
               }
             , zbx, tzbx,
               {
-                  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct);
+                  mltensor_cross_terms_fz(tzbx,fzfab,vfab,etazfab,kapzfab,dxinv,dlo,dhi,bct,is_periodic[2]);
               }
             );
 
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H
index 716af64a34e..6bb4b7215aa 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLTensor_2D_K.H
@@ -131,7 +131,8 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                               Array4<Real const> const& kapx,
                               GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                               const Dim3& dlo, const Dim3& dhi,
-                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                              const int& is_periodic) noexcept
 {
     const Real dyi = dxinv[1];
     const auto lo = amrex::lbound(box);
@@ -143,7 +144,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
         for (int i = lo.x; i <= hi.x; ++i) {
             Real dudy = (vel(i,j+1,0,0)+vel(i-1,j+1,0,0)-vel(i,j-1,0,0)-vel(i-1,j-1,0,0))*(0.25*dyi);
             Real dvdy = (vel(i,j+1,0,1)+vel(i-1,j+1,0,1)-vel(i,j-1,0,1)-vel(i-1,j-1,0,1))*(0.25*dyi);
-            if (i==dhi.x+1){
+            if (i==dhi.x+1 and !is_periodic){
                 if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                     dudy = (0.5*dyi) * (vel(i-1,j+1,0,0)-vel(i-1,j-1,0,0));
                 }
@@ -157,7 +158,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                     dvdy = (0.5*dyi) * (vel(i  ,j+1,0,1)-vel(i  ,j-1,0,1));
                 }
             }
-            if (i==dlo.x){
+            if (i==dlo.x and !is_periodic){
                 if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                     dudy = (0.5*dyi) * (vel(i  ,j+1,0,0)-vel(i  ,j-1,0,0));
                 }
@@ -188,7 +189,8 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                               Array4<Real const> const& kapy,
                               GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                               const Dim3& dlo, const Dim3& dhi,
-                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                              const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const auto lo = amrex::lbound(box);
@@ -200,7 +202,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
         for (int i = lo.x; i <= hi.x; ++i) {
             Real dudx = (vel(i+1,j,0,0)+vel(i+1,j-1,0,0)-vel(i-1,j,0,0)-vel(i-1,j-1,0,0))*(0.25*dxi);
             Real dvdx = (vel(i+1,j,0,1)+vel(i+1,j-1,0,1)-vel(i-1,j,0,1)-vel(i-1,j-1,0,1))*(0.25*dxi);
-            if (j==dhi.y+1){
+            if (j==dhi.y+1 and !is_periodic){
                 if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){
                     dudx = (0.5*dxi) * (vel(i+1,j-1,0,0)-vel(i-1,j-1,0,0));
                 }
@@ -214,7 +216,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                     dvdx = (0.5*dxi) * (vel(i+1,j  ,0,1)-vel(i-1,j  ,0,1));
                 }
             }
-            if (j==dlo.y){
+            if (j==dlo.y and !is_periodic){
                 if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){
                     dudx = (0.5*dxi) * (vel(i+1,j  ,0,0)-vel(i-1,j  ,0,0));
                 }
diff --git a/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H b/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H
index c10c1c3e7da..b2d23038d12 100644
--- a/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H
+++ b/Src/LinearSolvers/MLMG/AMReX_MLTensor_3D_K.H
@@ -1058,7 +1058,8 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                               Array4<Real const> const& kapx,
                               GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                               const Dim3& dlo, const Dim3& dhi,
-                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                              const int& is_periodic) noexcept
 {
     const Real dyi = dxinv[1];
     const Real dzi = dxinv[2];
@@ -1072,7 +1073,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
             for (int i = lo.x; i <= hi.x; ++i) {
                 Real dudy = (vel(i,j+1,k,0)+vel(i-1,j+1,k,0)-vel(i,j-1,k,0)-vel(i-1,j-1,k,0))*(0.25*dyi);
                 Real dvdy = (vel(i,j+1,k,1)+vel(i-1,j+1,k,1)-vel(i,j-1,k,1)-vel(i-1,j-1,k,1))*(0.25*dyi);
-                if (i==dhi.x+1){
+                if (i==dhi.x+1 and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudy = (0.5*dyi) * (vel(i-1,j+1,k,0)-vel(i-1,j-1,k,0));
                     }
@@ -1086,7 +1087,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                         dvdy = (0.5*dyi) * (vel(i  ,j+1,k,1)-vel(i  ,j-1,k,1));
                     }
                 }
-                if (i==dlo.x){
+                if (i==dlo.x and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudy = (0.5*dyi) * (vel(i  ,j+1,k,0)-vel(i  ,j-1,k,0));
                     }
@@ -1102,7 +1103,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                 }
                 Real dudz = (vel(i,j,k+1,0)+vel(i-1,j,k+1,0)-vel(i,j,k-1,0)-vel(i-1,j,k-1,0))*(0.25*dzi);
                 Real dwdz = (vel(i,j,k+1,2)+vel(i-1,j,k+1,2)-vel(i,j,k-1,2)-vel(i-1,j,k-1,2))*(0.25*dzi);
-                if (i==dhi.x+1){
+                if (i==dhi.x+1 and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudz = (0.5*dzi) * (vel(i-1,j,k+1,0)-vel(i-1,j,k-1,0));
                     }
@@ -1116,7 +1117,7 @@ void mltensor_cross_terms_fx (Box const& box, Array4<Real> const& fx,
                         dwdz = (0.5*dzi) * (vel(i  ,j,k+1,2)-vel(i  ,j,k-1,2));
                     }
                 }
-                if (i==dlo.x){
+                if (i==dlo.x and !is_periodic){
                     if (bct(Orientation(Direction::x,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudz = (0.5*dzi) * (vel(i  ,j,k+1,0)-vel(i  ,j,k-1,0));
                     }
@@ -1150,7 +1151,8 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                               Array4<Real const> const& kapy,
                               GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                               const Dim3& dlo, const Dim3& dhi,
-                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                              const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const Real dzi = dxinv[2];
@@ -1164,7 +1166,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
             for (int i = lo.x; i <= hi.x; ++i) {
                 Real dudx = (vel(i+1,j,k,0)+vel(i+1,j-1,k,0)-vel(i-1,j,k,0)-vel(i-1,j-1,k,0))*(0.25*dxi);
                 Real dvdx = (vel(i+1,j,k,1)+vel(i+1,j-1,k,1)-vel(i-1,j,k,1)-vel(i-1,j-1,k,1))*(0.25*dxi);
-                if (j==dhi.y+1){
+                if (j==dhi.y+1 and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudx = (0.5*dxi) * (vel(i+1,j-1,k,0)-vel(i-1,j-1,k,0));
                     }
@@ -1178,7 +1180,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                         dvdx = (0.5*dxi) * (vel(i+1,j  ,k,1)-vel(i-1,j  ,k,1));
                     }
                 }
-                if (j==dlo.y){
+                if (j==dlo.y and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudx = (0.5*dxi) * (vel(i+1,j  ,k,0)-vel(i-1,j  ,k,0));
                     }
@@ -1194,7 +1196,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                 }
                 Real dvdz = (vel(i,j,k+1,1)+vel(i,j-1,k+1,1)-vel(i,j,k-1,1)-vel(i,j-1,k-1,1))*(0.25*dzi);
                 Real dwdz = (vel(i,j,k+1,2)+vel(i,j-1,k+1,2)-vel(i,j,k-1,2)-vel(i,j-1,k-1,2))*(0.25*dzi);
-                if (j==dhi.y+1){
+                if (j==dhi.y+1 and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::high),1)==AMREX_LO_NEUMANN){
                         dvdz = (0.5*dzi) * (vel(i,j-1,k+1,1)-vel(i,j-1,k-1,1));
                     }
@@ -1208,7 +1210,7 @@ void mltensor_cross_terms_fy (Box const& box, Array4<Real> const& fy,
                         dwdz = (0.5*dzi) * (vel(i,j  ,k+1,2)-vel(i,j  ,k-1,2));
                     }
                 }
-                if (j==dlo.y){
+                if (j==dlo.y and !is_periodic){
                     if (bct(Orientation(Direction::y,Orientation::low),1)==AMREX_LO_NEUMANN){
                         dvdz = (0.5*dzi) * (vel(i,j  ,k+1,1)-vel(i,j  ,k-1,1));
                     }
@@ -1241,7 +1243,8 @@ void mltensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                               Array4<Real const> const& kapz,
                               GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
                               const Dim3& dlo, const Dim3& dhi,
-                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct) noexcept
+                              Array2D<BoundCond,0,2*AMREX_SPACEDIM-1,0,AMREX_SPACEDIM-1> const& bct,
+                              const int& is_periodic) noexcept
 {
     const Real dxi = dxinv[0];
     const Real dyi = dxinv[1];
@@ -1255,7 +1258,7 @@ void mltensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
             for (int i = lo.x; i <= hi.x; ++i) {
                 Real dudx = (vel(i+1,j,k,0)+vel(i+1,j,k-1,0)-vel(i-1,j,k,0)-vel(i-1,j,k-1,0))*(0.25*dxi);
                 Real dwdx = (vel(i+1,j,k,2)+vel(i+1,j,k-1,2)-vel(i-1,j,k,2)-vel(i-1,j,k-1,2))*(0.25*dxi);
-                if (k==dhi.z+1){
+                if (k==dhi.z+1 and !is_periodic){
                     if (bct(Orientation(Direction::z,Orientation::high),0)==AMREX_LO_NEUMANN){
                         dudx = (0.5*dxi) * (vel(i+1,j,k-1,0)-vel(i-1,j,k-1,0));
                     }
@@ -1269,7 +1272,7 @@ void mltensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                         dwdx = (0.5*dxi) * (vel(i+1,j,k  ,2)-vel(i-1,j,k  ,2));
                     }
                 }
-                if (k==dlo.z){
+                if (k==dlo.z and !is_periodic){
                     if (bct(Orientation(Direction::z,Orientation::low),0)==AMREX_LO_NEUMANN){
                         dudx = (0.5*dxi) * (vel(i+1,j,k  ,0)-vel(i-1,j,k  ,0));
                     }
@@ -1285,7 +1288,7 @@ void mltensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                 }
                 Real dvdy = (vel(i,j+1,k,1)+vel(i,j+1,k-1,1)-vel(i,j-1,k,1)-vel(i,j-1,k-1,1))*(0.25*dyi);
                 Real dwdy = (vel(i,j+1,k,2)+vel(i,j+1,k-1,2)-vel(i,j-1,k,2)-vel(i,j-1,k-1,2))*(0.25*dyi);
-                if (k==dhi.z+1){
+                if (k==dhi.z+1 and !is_periodic){
                     if (bct(Orientation(Direction::z,Orientation::high),1)==AMREX_LO_NEUMANN){
                         dvdy = (0.5*dyi) * (vel(i,j+1,k-1,1)-vel(i,j-1,k-1,1));
                     }
@@ -1299,7 +1302,7 @@ void mltensor_cross_terms_fz (Box const& box, Array4<Real> const& fz,
                         dwdy = (0.5*dyi) * (vel(i,j+1,k  ,2)-vel(i,j-1,k  ,2));
                     }
                 }
-                if (k==dlo.z){
+                if (k==dlo.z and !is_periodic){
                     if (bct(Orientation(Direction::z,Orientation::low),1)==AMREX_LO_NEUMANN){
                         dvdy = (0.5*dyi) * (vel(i,j+1,k  ,1)-vel(i,j-1,k  ,1));
                     }