diff --git a/src/pde/sweeps/Multigrid.impl.h b/src/pde/sweeps/Multigrid.impl.h
index e4b5decf62bcd597fc7ef11eea707d1b3f85f42a..da7c0ab6615781b492adcb4e6a9ce55c804e0d6f 100644
--- a/src/pde/sweeps/Multigrid.impl.h
+++ b/src/pde/sweeps/Multigrid.impl.h
@@ -26,18 +26,7 @@
 
 #include "stencil/D3Q7.h"
 
-#ifdef _MSC_VER
-#  pragma warning(push)
-// disable warning for multi_array: "declaration of 'extents' hides global declaration"
-#  pragma warning( disable : 4459 )
-#endif //_MSC_VER
-
-#include <boost/multi_array.hpp>
-
-#ifdef _MSC_VER
-#  pragma warning(pop)
-#endif //_MSC_VER
-
+#include <array>
 
 namespace walberla {
 namespace pde {
@@ -195,32 +184,34 @@ void CoarsenStencilFieldsGCA< stencil::D3Q7 >::operator()( const std::vector<Blo
          StencilField_T * coarse = block->getData< StencilField_T >( stencilFieldId[lvl] );
 
 
-         typedef boost::multi_array<real_t, 3> Array3D;
-         Array3D p(boost::extents[7][7][7]);
-         Array3D r(boost::extents[2][2][2]);
+         typedef std::array<std::array<std::array<real_t,7>,7>,7> Array3D_7;
+         typedef std::array<std::array<std::array<real_t,2>,2>,2> Array3D_2;
+
+         Array3D_2 r;
+         Array3D_7 p;
 
          // Set to restriction weights from constant restrict operator
-          for(Array3D::index z=0; z<2; ++z) {
-            for(Array3D::index y=0; y<2; ++y) {
-               for(Array3D::index x=0; x<2; ++x) {
+          for(auto z = uint_t(0); z < uint_t(2); ++z) {
+            for(auto y = uint_t(0); y < uint_t(2); ++y) {
+               for(auto x = uint_t(0); x < uint_t(2); ++x) {
                   r[x][y][z] = real_t(1);
                }
             }
          }
 
          // Init to 0 //
-         for(Array3D::index k=0; k<7; ++k) {
-            for(Array3D::index j=0; j<7; ++j) {
-               for(Array3D::index i=0; i<7; ++i) {
+         for(auto k = uint_t(0); k < uint_t(7); ++k) {
+            for(auto j = uint_t(0); j < uint_t(7); ++j) {
+               for(auto i = uint_t(0); i < uint_t(7); ++i) {
                   p[i][j][k] = real_t(0);
                }
             }
          }
 
          // Set center to prolongation weights, including overrelaxation factor (latter therefore no longer needed in ProlongateAndCorrect)
-         for(Array3D::index k=0; k<2; ++k) {
-            for(Array3D::index j=0; j<2; ++j) {
-               for(Array3D::index i=0; i<2; ++i) {
+         for(auto k = uint_t(0); k < uint_t(2); ++k) {
+            for(auto j = uint_t(0); j < uint_t(2); ++j) {
+               for(auto i = uint_t(0); i < uint_t(2); ++i) {
                   p[i+2][j+2][k+2] = real_t(0.125)/overrelaxFact_;   // Factor 0.125 such that same prolongation operator for DCA and GCA
                }
             }
@@ -229,35 +220,35 @@ void CoarsenStencilFieldsGCA< stencil::D3Q7 >::operator()( const std::vector<Blo
 
          WALBERLA_FOR_ALL_CELLS_XYZ(coarse,
 
-            Array3D ap(boost::extents[7][7][7]);
+            Array3D_7 ap;
 
             const cell_idx_t fx = 2*x;
             const cell_idx_t fy = 2*y;
             const cell_idx_t fz = 2*z;
 
-            for(Array3D::index k=0; k<7; ++k)
-               for(Array3D::index j=0; j<7; ++j)
-                  for(Array3D::index i=0; i<7; ++i)
+            for(auto k = uint_t(0); k < uint_t(7); ++k)
+               for(auto j = uint_t(0); j < uint_t(7); ++j)
+                  for(auto i = uint_t(0); i < uint_t(7); ++i)
                      ap[i][j][k] = real_t(0);
 
 
             // Tested for spatially varying stencils! //
-            for(Array3D::index k=1; k<5; ++k)
-               for(Array3D::index j=1; j<5; ++j)
-                  for(Array3D::index i=1; i<5; ++i) {
+            for(auto k = uint_t(1); k < uint_t(5); ++k)
+               for(auto j = uint_t(1); j < uint_t(5); ++j)
+                  for(auto i = uint_t(1); i < uint_t(5); ++i){
                      ap[i][j][k] = real_t(0);
                      for(auto d = stencil::D3Q7::begin(); d != stencil::D3Q7::end(); ++d ){
-                        ap[i][j][k] += p[ i+d.cx() ] [ j+d.cy() ] [k+d.cz() ] * fine->get( fx+cell_idx_c(i%2), fy+cell_idx_c(j%2), fz+cell_idx_c(k%2), d.toIdx() ); // contains elements of one row of coarse grid matrix
+                        ap[i][j][k] += p[ uint_c(cell_idx_c(i)+d.cx()) ] [ uint_c(cell_idx_c(j)+d.cy()) ] [uint_c(cell_idx_c(k)+d.cz()) ] * fine->get( fx+cell_idx_c(i%2), fy+cell_idx_c(j%2), fz+cell_idx_c(k%2), d.toIdx() ); // contains elements of one row of coarse grid matrix
                      }
                   }
 
             // Checked, correct! //
             for(auto d = stencil::D3Q7::begin(); d != stencil::D3Q7::end(); ++d ){
                real_t sum = 0;
-               for(Array3D::index k=0; k<2; ++k)
-                  for(Array3D::index j=0; j<2; ++j)
-                     for(Array3D::index i=0; i<2; ++i) {
-                        sum += ap[i+2-2*d.cx()] [j+2-2*d.cy()] [k+2-2*d.cz()] * r[i][j][k];
+               for(auto k = uint_t(0); k < uint_t(2); ++k)
+                  for(auto j = uint_t(0); j < uint_t(2); ++j)
+                     for(auto i = uint_t(0); i < uint_t(2); ++i) {
+                        sum += ap[uint_c(cell_idx_c(i)+2-2*d.cx())] [uint_c(cell_idx_c(j)+2-2*d.cy())] [uint_c(cell_idx_c(k)+2-2*d.cz())] * r[i][j][k];
                         // either i+0 or i+4    // either j+0 or j+4    // either k+0 or k+4       // always 1 here
                      }
 
diff --git a/tests/geometry/VoxelFileTest.cpp b/tests/geometry/VoxelFileTest.cpp
index d8ab6ddb679cff372a72fea03d056b9d0e2bd523..8d3e8ace32aa3959f7597c8c85e91bba19a2b607 100644
--- a/tests/geometry/VoxelFileTest.cpp
+++ b/tests/geometry/VoxelFileTest.cpp
@@ -20,19 +20,6 @@
 //======================================================================================================================
 
 
-
-#ifdef _MSC_VER
-#  pragma warning(push)
-// disable warning boost/multi_array/base.hpp(475): warning C4189: 'bound_adjustment' : local variable is initialized but not referenced
-#  pragma warning( disable : 4189 )
-#endif //_MSC_VER
-
-#include <boost/multi_array/base.hpp>
-
-#ifdef _MSC_VER
-#  pragma warning(pop)
-#endif //_MSC_VER
-
 #include "geometry/structured/VoxelFileReader.h"
 #include "core/Abort.h"
 #include "core/DataTypes.h"
@@ -44,21 +31,8 @@
 
 #include <random>
 
-#ifdef _MSC_VER
-#  pragma warning(push)
-// disable warning for multi_array: "declaration of 'extents' hides global declaration"
-#  pragma warning( disable : 4459 )
-#endif //_MSC_VER
-
-#include <boost/multi_array.hpp>
-
-#ifdef _MSC_VER
-#  pragma warning(pop)
-#endif //_MSC_VER
-
 typedef std::mersenne_twister_engine< walberla::uint32_t, 32, 351, 175, 19, 0xccab8ee7, 11, 0xffffffff, 7, 0x31b6ab00, 15, 0xffe50000, 17, 0xa37d3c92 > mt11213b;
 
-
 /// randomize the memory underlying the vector up the maximal size (== capacity)
 template<typename T>
 void randomizeVector( std::vector<T> & v )
@@ -102,48 +76,41 @@ void randomizeVector( std::vector<char> & v )
 }
 
 template<typename T>
-void makeRandomMultiArray( boost::multi_array<T, 3> & ma)
+void makeRandomMultiArray( std::vector<T> & ma)
 {
    static_assert(sizeof(T) > sizeof(char), "cannot use char");
 
    mt11213b rng;
    std::uniform_int_distribution<T> dist( std::numeric_limits<T>::min(), std::numeric_limits<T>::max() );
 
-   for(T* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
+   for(auto it = ma.begin(); it != ma.end(); ++it)
       *it = dist(rng);
 }
 
 template<>
-void makeRandomMultiArray( boost::multi_array<unsigned char, 3> & ma)
-{
+void makeRandomMultiArray( std::vector<unsigned char> & ma) {
    mt11213b rng;
-   std::uniform_int_distribution<walberla::int16_t> dist( std::numeric_limits<unsigned char>::min(), std::numeric_limits<unsigned char>::max() );
+   std::uniform_int_distribution<walberla::int16_t> dist(std::numeric_limits<unsigned char>::min(), std::numeric_limits<unsigned char>::max());
 
-   for(unsigned char* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
-      *it = static_cast<unsigned char>( dist(rng) );
+   for (auto it = ma.begin(); it != ma.end(); ++it)
+      *it = static_cast<unsigned char>( dist(rng));
 }
 
 template<>
-void makeRandomMultiArray( boost::multi_array<char, 3> & ma)
+void makeRandomMultiArray( std::vector<char> & ma)
 {
    mt11213b rng;
    std::uniform_int_distribution<int16_t> dist( std::numeric_limits<char>::min(), std::numeric_limits<char>::max() );
-
-   for(char* it = ma.data(); it != ma.data() + ma.num_elements(); ++it)
+   for (auto it = ma.begin(); it != ma.end(); ++it)
       *it = static_cast<char>( dist(rng) );
 }
 
 template<typename T>
 void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t zSize);
 
-template<typename T>
-void makeRandomMultiArray( boost::multi_array<T, 3> & ma);
-
 void modifyHeader(std::string inputFilename, std::string outputFilename,
                   size_t xSize, size_t ySize, size_t zSize);
 
-template<typename T>
-void makeRandomMultiArray( boost::multi_array<T, 3> & ma);
 
 int main(int argc, char** argv)
 {
@@ -281,13 +248,12 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
 
    data.clear();
 
-   using bindex = boost::multi_array_types::index;
+   std::vector<T> reference(zSize * ySize * xSize);
 
-   boost::multi_array<T, 3> reference(boost::extents[walberla::numeric_cast<bindex>(zSize)][walberla::numeric_cast<bindex>(ySize)][walberla::numeric_cast<bindex>(xSize)]);
    makeRandomMultiArray(reference);
 
-   data.resize(reference.num_elements());
-   std::copy( reference.data(), reference.data() + reference.num_elements(), data.begin() );
+   data.resize(reference.size());
+   std::copy( reference.begin(), reference.end(), data.begin() );
 
    geometryFile.open(filename);
    WALBERLA_CHECK( geometryFile.isOpen() );
@@ -302,8 +268,8 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
    randomizeVector(data);
    geometryFile.read(aabb, data);
 
-   WALBERLA_CHECK_EQUAL(reference.num_elements(), data.size());
-   WALBERLA_CHECK( std::equal(reference.data(), reference.data() + reference.num_elements(), data.begin()) );
+   WALBERLA_CHECK_EQUAL(reference.size(), data.size());
+   WALBERLA_CHECK( std::equal(reference.begin(), reference.end(), data.begin()) );
 
    randomizeVector(data);
 
@@ -317,8 +283,8 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
 
    geometryFile.read(aabb, data);
 
-   WALBERLA_CHECK_EQUAL(reference.num_elements(), data.size());
-   WALBERLA_CHECK( std::equal(reference.data(), reference.data() + reference.num_elements(), data.begin()) );
+   WALBERLA_CHECK_EQUAL(reference.size(), data.size());
+   WALBERLA_CHECK( std::equal(reference.begin(), reference.end(), data.begin()) );
 
    std::vector<size_t> blockSizesX;
    blockSizesX.push_back(std::max(xSize / 2, size_t(1)));
@@ -359,17 +325,16 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
                      geometryFile.read(blockAABB, data);
                      WALBERLA_CHECK_EQUAL( data.size(), blockAABB.numCells() );
 
-                     typename boost::multi_array<T, 3>::template array_view<3>::type myview =
-                        reference[ boost::indices[boost::multi_array_types::index_range(blockAABB.zMin(), blockAABB.zMax() + 1)]
-                                                 [boost::multi_array_types::index_range(blockAABB.yMin(), blockAABB.yMax() + 1)]
-                                                 [boost::multi_array_types::index_range(blockAABB.xMin(), blockAABB.xMax() + 1)] ];
+                     auto zIndexOffset = walberla::numeric_cast<size_t>(blockAABB.zMin());
+                     auto yIndexOffset = walberla::numeric_cast<size_t>(blockAABB.yMin());
+                     auto xIndexOffset = walberla::numeric_cast<size_t>(blockAABB.xMin());
 
                      size_t vectorIdx = 0;
                      for(size_t z = 0; z < blockAABB.zSize(); ++z)
                         for(size_t y = 0; y < blockAABB.ySize(); ++y)
                            for(size_t x = 0; x < blockAABB.xSize(); ++x)
                            {
-                              WALBERLA_CHECK_EQUAL(data[vectorIdx], myview[walberla::numeric_cast<bindex>(z)][walberla::numeric_cast<bindex>(y)][walberla::numeric_cast<bindex>(x)]);
+                              WALBERLA_CHECK_EQUAL(data[vectorIdx], reference[(zIndexOffset+z)*ySize*xSize + (yIndexOffset+y)*xSize + (xIndexOffset+x)]);
                               ++vectorIdx;
                            }
                   }
@@ -408,10 +373,9 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
                      WALBERLA_CHECK_LESS(blockAABB.yMax(), ySize);
                      WALBERLA_CHECK_LESS(blockAABB.zMax(), zSize);
 
-                     typename boost::multi_array<T, 3>::template array_view<3>::type myview =
-                        reference[ boost::indices[boost::multi_array_types::index_range(blockAABB.zMin(), blockAABB.zMax() + 1)]
-                                                 [boost::multi_array_types::index_range(blockAABB.yMin(), blockAABB.yMax() + 1)]
-                                                 [boost::multi_array_types::index_range(blockAABB.xMin(), blockAABB.xMax() + 1)] ];
+                     auto zIndexOffset = walberla::numeric_cast<size_t>(blockAABB.zMin());
+                     auto yIndexOffset = walberla::numeric_cast<size_t>(blockAABB.yMin());
+                     auto xIndexOffset = walberla::numeric_cast<size_t>(blockAABB.xMin());
 
                      data.resize(blockAABB.numCells());
                      size_t vectorIdx = 0;
@@ -419,7 +383,7 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
                         for(size_t y = 0; y < blockAABB.ySize(); ++y)
                            for(size_t x = 0; x < blockAABB.xSize(); ++x)
                            {
-                              data[vectorIdx] = myview[walberla::numeric_cast<bindex>(z)][walberla::numeric_cast<bindex>(y)][walberla::numeric_cast<bindex>(x)];
+                              data[vectorIdx] = reference[(zIndexOffset+z)*ySize*xSize + (yIndexOffset+y)*xSize + (xIndexOffset+x)];
                               ++vectorIdx;
                            }
 
@@ -428,8 +392,8 @@ void runTests(const std::string & filename, size_t xSize, size_t ySize, size_t z
                }
             }
             geometryFile.read(aabb, data);
-            WALBERLA_CHECK_EQUAL(reference.num_elements(), data.size());
-            WALBERLA_CHECK( std::equal(reference.data(), reference.data() + reference.num_elements(), data.begin()) );
+            WALBERLA_CHECK_EQUAL(reference.size(), data.size());
+            WALBERLA_CHECK( std::equal(reference.begin(), reference.end(), data.begin()) );
          }
       }
    }