diff --git a/IO/Reader.cpp b/IO/Reader.cpp
index c9595a8b..dedb9b82 100644
--- a/IO/Reader.cpp
+++ b/IO/Reader.cpp
@@ -44,7 +44,6 @@ std::vector<std::string> IO::readTimesteps( const std::string& filename )
     FILE *fid= fopen(filename.c_str(),"rb");
     if ( fid==NULL )
         ERROR("Error opening file");
-    auto pos = std::min(filename.find_last_of(47),filename.find_last_of(90));
     std::vector<std::string> timesteps;
     char buf[1000];
     while (fgets(buf,sizeof(buf),fid) != NULL) {
@@ -160,7 +159,6 @@ std::shared_ptr<IO::Mesh> IO::getMesh( const std::string& path, const std::strin
 #ifdef USE_SILO
         const DatabaseEntry& database = meshDatabase.domains[domain];
         std::string filename = path + "/" + timestep + "/" + database.file;
-        int rank = std::stoi(database.file.substr(0,database.file.find(".silo")).c_str());
         auto fid = silo::open( filename, silo::READ );
         if ( meshDatabase.meshClass=="PointList" ) {
             Array<double> coords = silo::readPointMesh<double>( fid, database.name );
@@ -262,7 +260,6 @@ std::shared_ptr<IO::Variable> IO::getVariable( const std::string& path, const st
         const auto& database = meshDatabase.domains[domain];
         auto variableDatabase = meshDatabase.getVariableDatabase( variable );
         std::string filename = path + "/" + timestep + "/" + database.file;
-        int rank = std::stoi(database.file.substr(0,database.file.find(".silo")).c_str());
         auto fid = silo::open( filename, silo::READ );
         var.reset( new Variable( variableDatabase.dim, variableDatabase.type, variable ) );
         if ( meshDatabase.meshClass=="PointList" ) {
diff --git a/IO/silo.hpp b/IO/silo.hpp
index d1ecee0c..312f32d8 100644
--- a/IO/silo.hpp
+++ b/IO/silo.hpp
@@ -322,9 +322,7 @@ void readTriMesh( DBfile* fid, const std::string& meshname, Array<TYPE>& coords,
     }
     auto zones = mesh->zones;
     int N_zones = zones->nzones;
-    int ndim_zones = zones->ndims;
     ASSERT( zones->nshapes==1 );
-    int shape_type = zones->shapetype[0];
     int shapesize = zones->shapesize[0];
     tri.resize(N_zones,shapesize);
     for (int i=0; i<N_zones; i++) {
diff --git a/analysis/distance.cpp b/analysis/distance.cpp
index c5aa4d1e..f539aebf 100644
--- a/analysis/distance.cpp
+++ b/analysis/distance.cpp
@@ -1,60 +1,12 @@
 #include "analysis/distance.h"
 
 
-// Check if we need to recompute distance after updating ghose values
-static bool checkUpdate( const Array<Vec> &d, double dx, double dy, double dz )
-{
-    auto s = d.size();
-    bool test[3] = { false, false, false };
-    // Check x-direction
-    Vec v1, v2;
-    for (size_t k=1; k<s[2]-1; k++) {
-        for (size_t j=1; j<s[1]-1; j++) {
-            v1 = d(1,j,k);
-            v2 = d(0,j,k);
-            v2.x += dx;
-            test[0] = test[0] || v2 < v1;
-            v1 = d(s[0]-2,j,k);
-            v2 = d(s[0]-1,j,k);
-            v2.x -= dx;
-            test[0] = test[0] || v2 < v1;
-        }
-    }
-    // Check y-direction
-    for (size_t k=1; k<s[2]-1; k++) {
-        for (size_t i=1; i<s[0]-1; i++) {
-            v1 = d(i,1,k);
-            v2 = d(i,0,k);
-            v2.y += dy;
-            test[1] = test[1] || v2 < v1;
-            v1 = d(i,s[1]-2,k);
-            v2 = d(i,s[1]-1,k);
-            v2.y -= dy;
-            test[1] = test[1] || v2 < v1;
-        }
-    }
-    // Check z-direction
-    for (size_t j=1; j<s[1]-1; j++) {
-        for (size_t i=1; i<s[0]-1; i++) {
-            v1 = d(i,j,1);
-            v2 = d(i,j,0);
-            v2.z += dz;
-            test[2] = test[2] || v2 < v1;
-            v1 = d(i,j,s[2]-2);
-            v2 = d(i,j,s[2]-1);
-            v2.z -= dz;
-            test[2] = test[2] || v2 < v1;
-        }
-    }
-    return test[0] || test[1] || test[2];
-}
-
-
 /******************************************************************
 * A fast distance calculation                                     *
 ******************************************************************/
 template<class TYPE>
-void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, const std::array<bool,3>& periodic )
+void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
+    const std::array<bool,3>& periodic, const std::array<double,3>& dx )
 {
     ASSERT( Distance.size() == ID.size() );
     std::array<int,3> n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 };
@@ -64,15 +16,121 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, c
     for (size_t i=0; i<ID.length(); i++)
         id(i) = ID(i) == 0 ? -1:1;
     fillData.fill( id );
-    CalcVecDist( vecDist, id, Dm, periodic );
+    CalcVecDist( vecDist, id, Dm, periodic, dx );
     for (size_t i=0; i<Distance.length(); i++)
         Distance(i) = id(i)*vecDist(i).norm();
 }
-void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm, const std::array<bool,3>& periodic )
+
+
+/******************************************************************
+* Vector-based distance calculation                               *
+* Initialize cells adjacent to boundaries                         *
+******************************************************************/
+static void calcVecInitialize( Array<Vec> &d, const Array<int> &ID, double dx, double dy, double dz )
+{
+    d.fill( Vec( 1e50, 1e50, 1e50 ) );
+    const double dx0 = 0.5*dx;
+    const double dy0 = 0.5*dy;
+    const double dz0 = 0.5*dz;
+    //const double dxy0 = 0.25*sqrt( dx*dx + dy*dy );
+    //const double dxz0 = 0.25*sqrt( dx*dx + dz*dz );
+    //const double dyz0 = 0.25*sqrt( dy*dy + dz*dz );
+    //const double dxyz0 = sqrt( dx*dx + dy*dy + dz*dz );
+    int Nx = d.size(0);
+    int Ny = d.size(1);
+    int Nz = d.size(2);
+    for (int k=1; k<Nz-1; k++) {
+        for (int j=1; j<Ny-1; j++) {
+            for (int i=1; i<Nx-1; i++) {
+                int id = ID(i,j,k);
+                bool x[2] = { id != ID(i-1,j,k), id != ID(i+1,j,k) };
+                bool y[2] = { id != ID(i,j-1,k), id != ID(i,j+1,k) };
+                bool z[2] = { id != ID(i,j,k-1), id != ID(i,j,k+1) };
+                if ( x[0] )  d(i,j,k) = Vec( dx0, 0, 0 );
+                if ( x[1] )  d(i,j,k) = Vec( -dx0, 0, 0 );
+                if ( y[0] )  d(i,j,k) = Vec( 0, dy0, 0 );
+                if ( y[1] )  d(i,j,k) = Vec( 0, -dy0, 0 );
+                if ( z[0] )  d(i,j,k) = Vec( 0, 0, dz0 );
+                if ( z[1] )  d(i,j,k) = Vec( 0, 0, -dz0 );
+                /*if ( x[0] && y[0] )  d(i,j,k) = Vec( dxy0, dxy0, 0 );
+                if ( x[0] && y[1] )  d(i,j,k) = Vec( dxy0, -dxy0, 0 );
+                if ( x[1] && y[0] )  d(i,j,k) = Vec( -dxy0, dxy0, 0 );
+                if ( x[1] && y[1] )  d(i,j,k) = Vec( -dxy0, -dxy0, 0 );
+                if ( x[0] && z[0] )  d(i,j,k) = Vec( dxz0, 0, dxz0 ); 
+                if ( x[0] && z[1] )  d(i,j,k) = Vec( dxz0, 0, -dxz0 );
+                if ( x[1] && z[0] )  d(i,j,k) = Vec( -dxz0, 0, dxz0 );
+                if ( x[1] && z[1] )  d(i,j,k) = Vec( -dxz0, 0, -dxz0 );
+                if ( y[0] && z[0] )  d(i,j,k) = Vec( 0, dyz0, dyz0 );
+                if ( y[0] && z[1] )  d(i,j,k) = Vec( 0, dyz0, -dyz0 );
+                if ( y[1] && z[0] )  d(i,j,k) = Vec( 0, -dyz0, dyz0 );
+                if ( y[1] && z[1] )  d(i,j,k) = Vec( 0, -dyz0, -dyz0 );*/
+            }
+        }
+    }
+
+}
+
+
+/******************************************************************
+* Vector-based distance calculation                               *
+* Update interior cells                                           *
+******************************************************************/
+static double calcVecUpdateInterior( Array<Vec> &d, double dx, double dy, double dz )
+{
+    double err = 0;
+    int Nx = d.size(0);
+    int Ny = d.size(1);
+    int Nz = d.size(2);
+    // Propagate (+,+,+)
+    for (int k=1; k<Nz; k++) {
+        for (int j=1; j<Ny; j++) {
+            for (int i=1; i<Nx; i++) {
+                auto vx = d(i-1,j,k);
+                auto vy = d(i,j-1,k);
+                auto vz = d(i,j,k-1);
+                vx.x += dx;
+                vy.y += dy;
+                vz.z += dz;
+                auto v = std::min( std::min(vx,vy), vz );
+                double d1 = v.norm2();
+                double d2 = d(i,j,k).norm2();
+                if ( d1 < d2 ) {
+                    d(i,j,k) = v;
+                    err = std::max( err, sqrt(d2)-sqrt(d1) );
+                }
+            }
+        }
+    }
+    // Propagate (-,-,-)
+    for (int k=Nz-2; k>=0; k--) {
+        for (int j=Ny-2; j>=0; j--) {
+            for (int i=Nx-2; i>=0; i--) {
+                auto vx = d(i+1,j,k);
+                auto vy = d(i,j+1,k);
+                auto vz = d(i,j,k+1);
+                vx.x -= dx;
+                vy.y -= dy;
+                vz.z -= dz;
+                auto v = std::min( std::min(vx,vy), vz );
+                double d1 = v.norm2();
+                double d2 = d(i,j,k).norm2();
+                if ( d1 < d2 ) {
+                    d(i,j,k) = v;
+                    err = std::max( err, sqrt(d2)-sqrt(d1) );
+                }
+            }
+        }
+    }
+    return err;
+}
+
+
+/******************************************************************
+* Vector-based distance calculation                               *
+******************************************************************/
+void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
+    const std::array<bool,3>& periodic, const std::array<double,3>& dx )
 {
-    const double dx = 1.0;
-    const double dy = 1.0;
-    const double dz = 1.0;
     std::array<int,3> N = { Dm.Nx, Dm.Ny, Dm.Nz };
     std::array<int,3> n = { Dm.Nx-2, Dm.Ny-2, Dm.Nz-2 };
     // Create ID with ghosts
@@ -102,112 +160,30 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm, const
     fillDataID.fill( ID );
     // Create communicator for distance
     fillHalo<Vec> fillData( Dm.Comm, Dm.rank_info, n, {1,1,1}, 50, 1, {true,false,false}, periodic );
-    // Initialize the vector distance
-    d.fill( Vec( 1e50, 1e50, 1e50 ) );
-    const double dx0 = 0.5*dx;
-    const double dy0 = 0.5*dy;
-    const double dz0 = 0.5*dz;
-    //const double dxy0 = 0.25*sqrt( dx*dx + dy*dy );
-    //const double dxz0 = 0.25*sqrt( dx*dx + dz*dz );
-    //const double dyz0 = 0.25*sqrt( dy*dy + dz*dz );
-    //const double dxyz0 = sqrt( dx*dx + dy*dy + dz*dz );
-    for (int k=1; k<Dm.Nz-1; k++) {
-        for (int j=1; j<Dm.Ny-1; j++) {
-            for (int i=1; i<Dm.Nx-1; i++) {
-                int id = ID(i,j,k);
-                bool x[2] = { id != ID(i-1,j,k), id != ID(i+1,j,k) };
-                bool y[2] = { id != ID(i,j-1,k), id != ID(i,j+1,k) };
-                bool z[2] = { id != ID(i,j,k-1), id != ID(i,j,k+1) };
-                if ( x[0] )  d(i,j,k) = Vec( dx0, 0, 0 );
-                if ( x[1] )  d(i,j,k) = Vec( -dx0, 0, 0 );
-                if ( y[0] )  d(i,j,k) = Vec( 0, dy0, 0 );
-                if ( y[1] )  d(i,j,k) = Vec( 0, -dy0, 0 );
-                if ( z[0] )  d(i,j,k) = Vec( 0, 0, dz0 );
-                if ( z[1] )  d(i,j,k) = Vec( 0, 0, -dz0 );
-                /*if ( x[0] && y[0] )  d(i,j,k) = Vec( dxy0, dxy0, 0 );
-                if ( x[0] && y[1] )  d(i,j,k) = Vec( dxy0, -dxy0, 0 );
-                if ( x[1] && y[0] )  d(i,j,k) = Vec( -dxy0, dxy0, 0 );
-                if ( x[1] && y[1] )  d(i,j,k) = Vec( -dxy0, -dxy0, 0 );
-                if ( x[0] && z[0] )  d(i,j,k) = Vec( dxz0, 0, dxz0 ); 
-                if ( x[0] && z[1] )  d(i,j,k) = Vec( dxz0, 0, -dxz0 );
-                if ( x[1] && z[0] )  d(i,j,k) = Vec( -dxz0, 0, dxz0 );
-                if ( x[1] && z[1] )  d(i,j,k) = Vec( -dxz0, 0, -dxz0 );
-                if ( y[0] && z[0] )  d(i,j,k) = Vec( 0, dyz0, dyz0 );
-                if ( y[0] && z[1] )  d(i,j,k) = Vec( 0, dyz0, -dyz0 );
-                if ( y[1] && z[0] )  d(i,j,k) = Vec( 0, -dyz0, dyz0 );
-                if ( y[1] && z[1] )  d(i,j,k) = Vec( 0, -dyz0, -dyz0 );*/
-            }
-        }
+    // Calculate the local distances
+    calcVecInitialize( d, ID, dx[0], dx[1], dx[2] );
+    double err = 1e100;
+    double tol = 0.5 * std::min( std::min(dx[0],dx[1]), dx[2] );
+    for (int it=0; it<=50 && err>tol; it++) {
+        err = calcVecUpdateInterior( d, dx[0], dx[1], dx[2] );
     }
-    int N_it = Dm.nprocx() + Dm.nprocy() + Dm.nprocz() + 3;
+    // Calculate the global distances
+    int N_it = Dm.nprocx() + Dm.nprocy() + Dm.nprocz() + 100;
     for ( int it=0; it<N_it; it++ ) {
-        //if ( Dm.rank() == 0 )
-        //    printf("Runing iteration %i\n",it+1);
-        // Propagate +/- x-direction
-        for (int k=0; k<Dm.Nz; k++) {
-            for (int j=0; j<Dm.Ny; j++) {
-                for (int i=1; i<Dm.Nx; i++) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i-1,j,k);
-                    v2.x += dx;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-                for (int i=Dm.Nx-2; i>=0; i--) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i+1,j,k);
-                    v2.x -= dx;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-            }
-        }
-        // Propagate +/- y-direction
-        for (int k=0; k<Dm.Nz; k++) {
-            for (int i=0; i<Dm.Nx; i++) {
-                for (int j=1; j<Dm.Ny; j++) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i,j-1,k);
-                    v2.y += dy;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-                for (int j=Dm.Ny-2; j>=0; j--) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i,j+1,k);
-                    v2.y -= dy;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-            }
-        }
-        // Propagate +/- z-direction
-        for (int j=0; j<Dm.Ny; j++) {
-            for (int i=0; i<Dm.Nx; i++) {
-                for (int k=1; k<Dm.Nz; k++) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i,j,k-1);
-                    v2.z += dz;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-                for (int k=Dm.Nz-2; k>=0; k--) {
-                    auto v1 = d(i,j,k);
-                    auto v2 = d(i,j,k+1);
-                    v2.z -= dz;
-                    if ( v2 < v1 )
-                        d(i,j,k) = v2;
-                }
-            }
-        }
+        // Update ghosts
         fillData.fill( d );
-        bool test = checkUpdate( d, dx, dy, dz );
-        test = sumReduce( Dm.Comm, test );
-        if ( !test )
+        // Update distance
+        double err = calcVecUpdateInterior( d, dx[0], dx[1], dx[2] );
+        // Check if we are finished
+        err = maxReduce( Dm.Comm, err );
+        if ( err < tol )
             break;
     }
 }
-template void CalcDist<float>( Array<float>&, const Array<char>&, const Domain&, const std::array<bool,3>& );
-template void CalcDist<double>( Array<double>&, const Array<char>&, const Domain&, const std::array<bool,3>& );
+
+
+// Explicit instantiations
+template void CalcDist<float>( Array<float>&, const Array<char>&, const Domain&, const std::array<bool,3>&, const std::array<double,3>& );
+template void CalcDist<double>( Array<double>&, const Array<char>&, const Domain&, const std::array<bool,3>&, const std::array<double,3>& );
 
 
diff --git a/analysis/distance.h b/analysis/distance.h
index b1bb9af5..b84a93f8 100644
--- a/analysis/distance.h
+++ b/analysis/distance.h
@@ -11,6 +11,7 @@ struct Vec {
     inline Vec(): x(0), y(0), z(0) {}
     inline Vec( double x_, double y_, double z_ ): x(x_), y(y_), z(z_) {}
     inline double norm() const { return sqrt(x*x+y*y+z*z); }
+    inline double norm2() const { return x*x+y*y+z*z; }
 };
 inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.z < r.x*r.x+r.y*r.y+r.z*r.z; }
 
@@ -24,7 +25,8 @@ inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.
  * @param[in] periodic      Directions that are periodic
  */
 template<class TYPE>
-void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, const std::array<bool,3>& periodic = {true,true,true} );
+void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
+    const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );
 
 /*!
  * @brief  Calculate the distance using a simple method
@@ -34,7 +36,7 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm, c
  * @param[in] Dm            Domain information
  * @param[in] periodic      Directions that are periodic
  */
-void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm, const std::array<bool,3>& periodic );
-
+void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
+    const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );
 
 #endif
diff --git a/common/MPI_Helpers.h b/common/MPI_Helpers.h
index 5a02397c..0cae743b 100644
--- a/common/MPI_Helpers.h
+++ b/common/MPI_Helpers.h
@@ -170,6 +170,12 @@ void unpack( std::set<TYPE>& data, const char *buffer );
 
 
 // Helper functions
+inline double sumReduce( MPI_Comm comm, double x )
+{
+    double y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_DOUBLE,MPI_SUM,comm);
+    return y;
+}
 inline float sumReduce( MPI_Comm comm, float x )
 {
     float y = 0;
@@ -199,12 +205,24 @@ inline std::vector<int> sumReduce( MPI_Comm comm, const std::vector<int>& x )
 	MPI_Allreduce(x.data(),y.data(),x.size(),MPI_INT,MPI_SUM,comm);
     return y;
 }
+inline double maxReduce( MPI_Comm comm, double x )
+{
+    double y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_DOUBLE,MPI_MAX,comm);
+    return y;
+}
 inline float maxReduce( MPI_Comm comm, float x )
 {
     float y = 0;
 	MPI_Allreduce(&x,&y,1,MPI_FLOAT,MPI_MAX,comm);
     return y;
 }
+inline int maxReduce( MPI_Comm comm, int x )
+{
+    int y = 0;
+	MPI_Allreduce(&x,&y,1,MPI_INT,MPI_MAX,comm);
+    return y;
+}
 
 
 #endif
diff --git a/tests/TestSegDist.cpp b/tests/TestSegDist.cpp
index f2b8b631..5866101d 100644
--- a/tests/TestSegDist.cpp
+++ b/tests/TestSegDist.cpp
@@ -17,11 +17,18 @@
 
 std::shared_ptr<Database> loadInputs( int nprocs )
 {
-    INSIST(nprocs==8, "TestSegDist: Number of MPI processes must be equal to 8");
+    std::vector<int> nproc;
+    if ( nprocs == 1 ) {
+        nproc = { 1, 1, 1 };
+    } else if ( nprocs == 8 ) {
+        nproc = { 2, 2, 2 };
+    } else {
+        ERROR("TestSegDist: Unsupported number of processors");
+    }
     auto db = std::make_shared<Database>( );
     db->putScalar<int>( "BC", 0 );
-    db->putVector<int>( "nproc", { 2, 2, 2 } );
-    db->putVector<int>( "n", { 100, 100, 100 } );
+    db->putVector<int>( "nproc", nproc );
+    db->putVector<int>( "n", { 200, 200, 200 } );
     db->putScalar<int>( "nspheres", 0 );
     db->putVector<double>( "L", { 1, 1, 1 } );
     return db;
@@ -63,35 +70,30 @@ int main(int argc, char **argv)
     int ny = Ny+2;
     int nz = Nz+2;
 
-    double BubbleRadius = 0.3*nx;
-
     // Initialize the bubble
-    double Cx = 1.0*nx;
-    double Cy = 1.0*ny;
-    double Cz = 1.0*nz;
+    double BubbleRadius = 0.15*Nx*Dm.nprocx();
+    double Cx = 0.40*Nx*Dm.nprocx();
+    double Cy = 0.45*Nx*Dm.nprocy();
+    double Cz = 0.50*Nx*Dm.nprocy();
 
     DoubleArray TrueDist(nx,ny,nz);
     Array<char> id(nx,ny,nz);
     id.fill(0);
 
-    for (int k=1;k<nz-1;k++){
-        for (int j=1;j<ny-1;j++){
-            for (int i=1;i<nx-1;i++){
-
+    for (int k=1; k<nz-1; k++) {
+        double z = k - 0.5 + Dm.kproc()*Nz;
+        for (int j=1; j<ny-1; j++) {
+            double y = j - 0.5 + Dm.jproc()*Ny;
+            for (int i=1; i<nx-1; i++) {
+                double x = i - 0.5 + Dm.iproc()*Nx;
                 // True signed distance
-                double x = (nx-2)*Dm.iproc()+i-1;
-                double y = (ny-2)*Dm.jproc()+j-1;
-                double z = (nz-2)*Dm.kproc()+k-1;
                 TrueDist(i,j,k) = sqrt((x-Cx)*(x-Cx)+(y-Cy)*(y-Cy)+(z-Cz)*(z-Cz)) - BubbleRadius;
-
                 // Initialize phase positions
                 if (TrueDist(i,j,k) < 0.0){
                     id(i,j,k) = 0;
                 } else{
                     id(i,j,k)=1;
                 }
-
-
             }
         }
     }
@@ -106,24 +108,18 @@ int main(int argc, char **argv)
     if (rank==0)
         printf("Total time: %f seconds \n",t2-t1);
 
-    double localError=0.0;
-    int localCount = 0;
-    for (int k=0;k<nz;k++){
-        for (int j=0;j<ny;j++){
-            for (int i=0;i<nx;i++){
-                if (fabs(TrueDist(i,j,k)) < 3.0){
-                    localError += (Distance(i,j,k)-TrueDist(i,j,k))*(Distance(i,j,k)-TrueDist(i,j,k));
-                    localCount++;
-                }
+    double err = 0.0;
+    for (int i=1; i<Nx-1; i++) {
+        for (int j=1; j<Ny-1; j++) {
+            for (int k=1; k<Nz-1; k++) {
+                err += (Distance(i,j,k)-TrueDist(i,j,k))*(Distance(i,j,k)-TrueDist(i,j,k));
             }
         }
     }
-    double globalError;
-    int globalCount;
-    MPI_Allreduce(&localError,&globalError,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    MPI_Allreduce(&localCount,&globalCount,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
-    double err2 = sqrt(globalError)/(double (globalCount));
-    if (rank==0) printf("Mean error %f \n", err2);
+    err = sumReduce( Dm.Comm, err );
+    err = sqrt( err / (nx*ny*nz*nprocs) );
+    if (rank==0)
+        printf("Mean error %0.4f \n", err);
 
     // Write the results to visit
     Array<int> ID0(id.size());