Commit 993bce66 authored by Christoph Rettinger's avatar Christoph Rettinger

Removed boost multi_array where possible

parent 55db38ed
......@@ -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
}
......
......@@ -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()) );
}
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment