diff --git a/Src/AmrCore/AMReX_Interp_3D_C.H b/Src/AmrCore/AMReX_Interp_3D_C.H
index b29ba4aca40..382962d11b2 100644
--- a/Src/AmrCore/AMReX_Interp_3D_C.H
+++ b/Src/AmrCore/AMReX_Interp_3D_C.H
@@ -294,13 +294,13 @@ facediv_int (int ci, int cj, int ck, int nf,
                                   + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
                                   +                dy3/(8*dz*xspys)*(w100+w112-w102-w110);
 
-    fine[1](fi+1, fj+1, fk  , nf) = Real(0.5)*(v001+v021)
+    fine[1](fi+1, fj+1, fk  , nf) = Real(0.5)*(v100+v120)
                                   + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u000+u210-u010-u200)
                                   +                dy3/(8*dx*yspzs)*(u001+u211-u011-u201)
                                   + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w100+w112-w102-w110)
                                   +                dy3/(8*dz*xspys)*(w000+w012-w002-w010);
 
-    fine[1](fi  , fj+1, fk+1, nf) = Real(0.5)*(v100+v120)
+    fine[1](fi  , fj+1, fk+1, nf) = Real(0.5)*(v001+v021)
                                   + dy*(2*dz*dz+dy*dy)/(8*dx*yspzs)*(u001+u211-u011-u201)
                                   +                dy3/(8*dx*yspzs)*(u000+u210-u010-u200)
                                   + dy*(2*dx*dx+dy*dy)/(8*dz*xspys)*(w000+w012-w002-w010)
@@ -319,17 +319,17 @@ facediv_int (int ci, int cj, int ck, int nf,
                                   +                dz3/(8*dy*zspxs)*(v100+v121-v101-v120);
 
     fine[2](fi  , fj+1, fk+1, nf) = Real(0.5)*(w010+w012)
-                                  + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
-                                  +                dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
-                                  + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
-                                  +                dz3/(8*dy*zspxs)*(v000+v021-v001-v020);
-
-    fine[2](fi+1, fj  , fk+1, nf) = Real(0.5)*(w100+w102)
                                   + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
                                   +                dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
                                   + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v000+v021-v001-v020)
                                   +                dz3/(8*dy*zspxs)*(v100+v121-v101-v120);
 
+    fine[2](fi+1, fj  , fk+1, nf) = Real(0.5)*(w100+w102)
+                                  + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u000+u201-u001-u200)
+                                  +                dz3/(8*dx*yspzs)*(u010+u211-u011-u210)
+                                  + dz*(2*dx*dx+dz*dz)/(8*dy*zspxs)*(v100+v121-v101-v120)
+                                  +                dz3/(8*dy*zspxs)*(v000+v021-v001-v020);
+
     fine[2](fi+1, fj+1, fk+1, nf) = Real(0.5)*(w110+w112)
                                   + dz*(2*dy*dy+dz*dz)/(8*dx*yspzs)*(u010+u211-u011-u210)
                                   +                dz3/(8*dx*yspzs)*(u000+u201-u001-u200)
diff --git a/Src/AmrCore/AMReX_Interpolater.H b/Src/AmrCore/AMReX_Interpolater.H
index c13fb283a35..6333112aae8 100644
--- a/Src/AmrCore/AMReX_Interpolater.H
+++ b/Src/AmrCore/AMReX_Interpolater.H
@@ -565,6 +565,9 @@ public:
 * of the divergence of the underlying crse cell. All fine cells overlying
 * a given coarse cell will have the same divergence, even when the coarse
 * grid divergence is spatially varying.
+* Based on Vanella et. al. (doi:10.1016/j.jcp.2010.05.003, section 3.2),
+* but solves the interior closure problem using least squares with an
+* initial guess equal to the average of fine face values across the cell.
 */
 class FaceDivFree
     :
diff --git a/Tests/DivFreePatch/main.cpp b/Tests/DivFreePatch/main.cpp
index 7358946c434..857461a1f64 100644
--- a/Tests/DivFreePatch/main.cpp
+++ b/Tests/DivFreePatch/main.cpp
@@ -94,11 +94,14 @@ void CoarsenToFine(MultiFab& div_refined_coarse,
 
 
 Real MFdiff(const MultiFab& lhs, const MultiFab& rhs,
-            int strt_comp, int num_comp, int nghost, const std::string name = "")
+            int strt_comp, int num_comp, int nghost, const std::string name = "",
+            bool relative = false)
 {
     MultiFab temp(lhs.boxArray(), lhs.DistributionMap(), lhs.nComp(), nghost);
     Copy(temp, lhs, strt_comp, strt_comp, num_comp, nghost);
     temp.minus(rhs, strt_comp, num_comp, nghost);
+    if (relative) {
+        temp.divide(rhs, strt_comp, num_comp, nghost); }
 
     if (name != "")
         { amrex::VisMF::Write(temp, std::string("pltfiles/" + name)); }
@@ -301,7 +304,7 @@ void main_main ()
 //  Setup initial value on the coarse faces.
     for (int i=0; i<AMREX_SPACEDIM; ++i)
     {
-        setupMF(c_mf_faces[i], 1);
+        setupMF(c_mf_faces[i], i);
     }
 
     amrex::UtilCreateDirectoryDestructive("pltfiles");
@@ -403,16 +406,21 @@ void main_main ()
                   amrex::VisMF::Write(f_mf_faces_wg[1], std::string("pltfiles/fwgy"));,
                   amrex::VisMF::Write(f_mf_faces_wg[2], std::string("pltfiles/fwgz"));  );
 
-    amrex::Print() << " Max InterpFromCoarse divergence error: "
-                   << MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diff") << '\n';
+    amrex::Print() << " Max InterpFromCoarse divergence error: absolute         relative\n "
+                   << "                                       "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diff")
+                   << "  "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
+                   << '\n';
 
 // ***************************************************************
 
     // Change coarse values, save current fine values for checking
     //   the final result.
+    amrex::Print()<<" Cyclically permute the velocity components. Should not change div\n";
     for (int i=0; i<AMREX_SPACEDIM; ++i)
     {
-        setupMF(c_mf_faces[i], 2, BoxArray(amrex::coarsen(f_geom.Domain(), ratio).convert(c_mf_faces[i].ixType())));
+        setupMF(c_mf_faces[i], (i+1)%3, BoxArray(amrex::coarsen(f_geom.Domain(), ratio).convert(c_mf_faces[i].ixType())));
         Copy(f_mf_copy[i], f_mf_faces[i], 0, 0, 1, 0);
     }
 
@@ -424,8 +432,13 @@ void main_main ()
     CoarsenToFine(div_refined_coarse, div_coarse, c_geom, f_geom_all, ratio);
     amrex::VisMF::Write(div_coarse, std::string("pltfiles/coarsetofineB"));
 
-    amrex::Print() << " Checking new adjustment hasn't changed solution on fine region: "
-                   << MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP") << '\n';
+    amrex::Print() << " Change to divergence on fine region:  absolute         relative\n "
+                   << "                                       "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f)
+                   << "  "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
+                   << '\n';
+
 
 
 // ***************************************************************
@@ -485,20 +498,23 @@ void main_main ()
 // ================================================
 
     // Checking fine valid cells are identical.
-    Real max_diff = 0;
     for (int i=0; i<AMREX_SPACEDIM; ++i)
     {
-        Real max_i = std::abs( MFdiff(f_mf_copy[i], f_mf_faces[i], 0, 1, 0) );
-        max_diff = (max_diff > max_i) ? max_diff : max_i;
+        amrex::Print() << " Fine valid region maximum change, comp "<<i<<": "
+                       <<MFdiff(f_mf_copy[i], f_mf_faces[i], 0, 1, 0)
+                       << '\n';
     }
-    amrex::Print() << " Fine values maximum change: " << max_diff << '\n';
 
     // Check fine divergence = coarse divergence in ghost cells.
     calcDiv(f_mf_faces, div_fine, f_geom.CellSizeArray());
     amrex::VisMF::Write(div_fine, std::string("pltfiles/fineFP"));
 
-    amrex::Print() << " Max FillPatchTwoLevels divergence error: "
-                   << MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP") << '\n';
+    amrex::Print() << " Max FillPatchTwoLevels divergence error:  absolute         relative\n "
+                   << "                                       "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "diffFP")
+                   << "  "
+                   <<MFdiff(div_fine, div_refined_coarse, 0, 1, nghost_f, "", true)
+                   << '\n';
 
     for (int i=0; i<AMREX_SPACEDIM; ++i)
     {