diff --git a/src/field/AddToStorage.h b/src/field/AddToStorage.h
index add0eb16d82095f1134245b36b8750759d45683a..4999b174639176547a878c3b15eadab44d688d2d 100644
--- a/src/field/AddToStorage.h
+++ b/src/field/AddToStorage.h
@@ -102,13 +102,15 @@ namespace internal {
 template< typename GhostLayerField_T, typename BlockStorage_T, class Enable = void >
 struct AddToStorage
 {
+   using Value_T = typename GhostLayerField_T::value_type;
    static BlockDataID add( const shared_ptr< BlockStorage_T > & blocks, const std::string & identifier,
                            const typename GhostLayerField_T::value_type & initValue, const Layout layout, const uint_t nrOfGhostLayers,
                            const bool /*alwaysInitialize*/, const std::function< void ( GhostLayerField_T * field, IBlock * const block ) > & initFunction,
                            const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
-                           const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = defaultSize )
+                           const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = defaultSize,
+                           const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr)
    {
-      auto dataHandling = walberla::make_shared< field::AlwaysInitializeBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize );
+      auto dataHandling = walberla::make_shared< field::AlwaysInitializeBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize, alloc );
       dataHandling->addInitializationFunction( initFunction );
       return blocks->addBlockData( dataHandling, identifier, requiredSelectors, incompatibleSelectors );
    }
@@ -119,20 +121,22 @@ struct AddToStorage< GhostLayerField_T, BlockStorage_T,
                      typename std::enable_if< ( std::is_integral< typename GhostLayerField_T::value_type >::value || std::is_floating_point< typename GhostLayerField_T::value_type >::value ) &&
 	                                            ! std::is_same< GhostLayerField_T, FlagField< typename GhostLayerField_T::value_type > >::value >::type >
 {
+   using Value_T = typename GhostLayerField_T::value_type;
    static BlockDataID add( const shared_ptr< BlockStorage_T > & blocks, const std::string & identifier,
                            const typename GhostLayerField_T::value_type & initValue, const Layout layout, const uint_t nrOfGhostLayers,
                            const bool alwaysInitialize, const std::function< void ( GhostLayerField_T * field, IBlock * const block ) > & initFunction,
                            const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
-                           const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = defaultSize )
+                           const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = defaultSize,
+                           const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr)
    {
       if( alwaysInitialize )
       {
-         auto dataHandling = walberla::make_shared< field::AlwaysInitializeBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize );
+         auto dataHandling = walberla::make_shared< field::AlwaysInitializeBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize, alloc );
          dataHandling->addInitializationFunction( initFunction );
          return blocks->addBlockData( dataHandling, identifier, requiredSelectors, incompatibleSelectors );
       }
 
-      auto dataHandling = walberla::make_shared< field::DefaultBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize );
+      auto dataHandling = walberla::make_shared< field::DefaultBlockDataHandling< GhostLayerField_T > >( blocks, nrOfGhostLayers, initValue, layout, calculateSize, alloc );
       dataHandling->addInitializationFunction( initFunction );
       return blocks->addBlockData( dataHandling, identifier, requiredSelectors, incompatibleSelectors );
    }
@@ -152,7 +156,7 @@ BlockDataID addToStorage( const shared_ptr< BlockStorage_T > & blocks,
                           const std::function< void ( GhostLayerField_T * field, IBlock * const block ) > & initFunction =
                              std::function< void ( GhostLayerField_T * field, IBlock * const block ) >(),
                           const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(),
-                          const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                          const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet())
 {
    return internal::AddToStorage< GhostLayerField_T, BlockStorage_T >::add( blocks, identifier, initValue, layout, nrOfGhostLayers,
                                                                             alwaysInitialize, initFunction, requiredSelectors, incompatibleSelectors );
@@ -179,7 +183,7 @@ BlockDataID addToStorage( const shared_ptr< BlockStorage_T > & blocks,
 template< typename GhostLayerField_T, typename BlockStorage_T >
 BlockDataID addToStorage( const shared_ptr< BlockStorage_T > & blocks,
                           const std::string & identifier,
-                          const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize,
+                          const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize,
                           const typename GhostLayerField_T::value_type & initValue = typename GhostLayerField_T::value_type(),
                           const Layout layout = zyxf,
                           const uint_t nrOfGhostLayers = uint_t(1),
@@ -187,15 +191,15 @@ BlockDataID addToStorage( const shared_ptr< BlockStorage_T > & blocks,
                           const std::function< void ( GhostLayerField_T * field, IBlock * const block ) > & initFunction =
                           std::function< void ( GhostLayerField_T * field, IBlock * const block ) >(),
                           const Set<SUID> & requiredSelectors = Set<SUID>::emptySet(),
-                          const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                          const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet(),
+                          const shared_ptr< field::FieldAllocator<typename GhostLayerField_T::value_type> > alloc = nullptr)
 {
    return internal::AddToStorage< GhostLayerField_T, BlockStorage_T >::add( blocks, identifier, initValue, layout, nrOfGhostLayers,
                                                                             alwaysInitialize, initFunction, requiredSelectors,
-                                                                            incompatibleSelectors, calculateSize );
+                                                                            incompatibleSelectors, calculateSize, alloc );
 }
 
 
-
 template< typename GhostLayerField_T, typename BlockStorage_T >
 BlockDataID addToStorage( const shared_ptr< BlockStorage_T > & blocks,
                           const std::string & identifier,
diff --git a/src/field/blockforest/BlockDataHandling.h b/src/field/blockforest/BlockDataHandling.h
index 2939352343ca8b3b6d3ed2b98bfa111f9b9e75f5..a156e08b787a15b35797f6942720d93f625ae928 100644
--- a/src/field/blockforest/BlockDataHandling.h
+++ b/src/field/blockforest/BlockDataHandling.h
@@ -333,55 +333,65 @@ namespace internal
 
 template< typename GhostLayerField_T >
 inline GhostLayerField_T * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl,
-                                     const typename GhostLayerField_T::value_type & v, Layout l )
+                                     const typename GhostLayerField_T::value_type & v, Layout l,
+                                     const shared_ptr< field::FieldAllocator<typename GhostLayerField_T::value_type> > & alloc=nullptr)
 {
-   return new GhostLayerField_T(x,y,z,gl,v,l);
+   return new GhostLayerField_T(x,y,z,gl,v,l,alloc);
 }
 template<>
-inline FlagField<uint8_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint8_t &, Layout )
+inline FlagField<uint8_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint8_t &, Layout,
+                                      const shared_ptr< field::FieldAllocator<uint8_t> > & alloc)
 {
-   return new FlagField<uint8_t>(x,y,z,gl);
+   return new FlagField<uint8_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint16_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint16_t &, Layout )
+inline FlagField<uint16_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint16_t &, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint16_t> > & alloc)
 {
-   return new FlagField<uint16_t>(x,y,z,gl);
+   return new FlagField<uint16_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint32_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint32_t &, Layout )
+inline FlagField<uint32_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint32_t &, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint32_t> > & alloc)
 {
-   return new FlagField<uint32_t>(x,y,z,gl);
+   return new FlagField<uint32_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint64_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint64_t &, Layout )
+inline FlagField<uint64_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, const uint64_t &, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint64_t> > & alloc)
 {
-   return new FlagField<uint64_t>(x,y,z,gl);
+   return new FlagField<uint64_t>(x,y,z,gl,alloc);
 }
 
 template< typename GhostLayerField_T >
-inline GhostLayerField_T * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout l )
+inline GhostLayerField_T * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout l,
+                                     const shared_ptr< field::FieldAllocator<typename GhostLayerField_T::value_type> > & alloc=nullptr)
 {
-   return new GhostLayerField_T(x,y,z,gl,l);
+   return new GhostLayerField_T(x,y,z,gl,l, alloc);
 }
 template<>
-inline FlagField<uint8_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout )
+inline FlagField<uint8_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout,
+                                      const shared_ptr< field::FieldAllocator<uint8_t> > & alloc)
 {
-   return new FlagField<uint8_t>(x,y,z,gl);
+   return new FlagField<uint8_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint16_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout )
+inline FlagField<uint16_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint16_t> > & alloc)
 {
-   return new FlagField<uint16_t>(x,y,z,gl);
+   return new FlagField<uint16_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint32_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout )
+inline FlagField<uint32_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint32_t> > & alloc)
 {
-   return new FlagField<uint32_t>(x,y,z,gl);
+   return new FlagField<uint32_t>(x,y,z,gl,alloc);
 }
 template<>
-inline FlagField<uint64_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout )
+inline FlagField<uint64_t> * allocate( const uint_t x, const uint_t y, const uint_t z, const uint_t gl, Layout,
+                                       const shared_ptr< field::FieldAllocator<uint64_t> > & alloc)
 {
-   return new FlagField<uint64_t>(x,y,z,gl);
+   return new FlagField<uint64_t>(x,y,z,gl,alloc);
 }
 
 inline Vector3< uint_t > defaultSize( const shared_ptr< StructuredBlockStorage > & blocks, IBlock * const block )
@@ -401,19 +411,22 @@ public:
    using Value_T = typename GhostLayerField_T::value_type;
 
    DefaultBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks,
-                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( uint_t(1) ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize )
+                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                             const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( uint_t(1) ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize ), alloc_(alloc)
    {}
 
    DefaultBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks, const uint_t nrOfGhostLayers,
-                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize )
+                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                             const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize ), alloc_(alloc)
    {}
 
    DefaultBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks, const uint_t nrOfGhostLayers,
                              const Value_T & initValue, const Layout layout = zyxf,
-                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_( initValue ), layout_( layout ), calculateSize_( calculateSize )
+                             const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                             const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_( initValue ), layout_( layout ), calculateSize_( calculateSize ), alloc_(alloc)
    {
       static_assert( !std::is_same< GhostLayerField_T, FlagField< Value_T > >::value,
                      "When using class FlagField, only constructors without the explicit specification of an initial value and the field layout are available!" );
@@ -427,16 +440,16 @@ protected:
       WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'DefaultBlockDataHandling' for a block storage object that doesn't exist anymore" );
       const Vector3< uint_t > size = calculateSize_( blocks, block );
       return internal::allocate< GhostLayerField_T >( size[0], size[1], size[2],
-                                                      nrOfGhostLayers_, initValue_, layout_ );
+                                                      nrOfGhostLayers_, initValue_, layout_, alloc_ );
    }
 
    GhostLayerField_T * reallocate( IBlock * const block ) override
    {
       auto blocks = blocks_.lock();
-      WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'DefaultBlockDataHandling' for a block storage object that doesn't exist anymore" );
+      WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'DefaultBlockDataHandling' for a block storage object that doesn't exist anymore" )
       const Vector3< uint_t > size = calculateSize_( blocks, block );
       return internal::allocate< GhostLayerField_T >( size[0], size[1], size[2],
-                                                      nrOfGhostLayers_, layout_ );
+                                                      nrOfGhostLayers_, layout_, alloc_ );
    }
 
 private:
@@ -447,6 +460,7 @@ private:
    Value_T initValue_;
    Layout  layout_;
    const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize_;
+   const shared_ptr< field::FieldAllocator<Value_T> > alloc_;
 
 }; // class DefaultBlockDataHandling
 
@@ -464,19 +478,22 @@ public:
    using InitializationFunction_T = std::function<void (GhostLayerField_T *, IBlock *const)>;
 
    AlwaysInitializeBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks,
-                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( uint_t(1) ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize )
+                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                                      const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( uint_t(1) ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize ), alloc_(alloc)
    {}
 
    AlwaysInitializeBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks, const uint_t nrOfGhostLayers,
-                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize )
+                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                                      const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_(), layout_( zyxf ), calculateSize_( calculateSize ), alloc_(alloc)
    {}
 
    AlwaysInitializeBlockDataHandling( const weak_ptr< StructuredBlockStorage > & blocks, const uint_t nrOfGhostLayers,
                                       const Value_T & initValue, const Layout layout,
-                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize = internal::defaultSize ) :
-      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_( initValue ), layout_( layout ), calculateSize_( calculateSize )
+                                      const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) >& calculateSize = internal::defaultSize,
+                                      const shared_ptr< field::FieldAllocator<Value_T> > alloc = nullptr) :
+      blocks_( blocks ), nrOfGhostLayers_( nrOfGhostLayers ), initValue_( initValue ), layout_( layout ), calculateSize_( calculateSize ), alloc_(alloc)
    {
       static_assert( ! std::is_same< GhostLayerField_T, FlagField< Value_T > >::value,
                      "When using class FlagField, only constructors without the explicit specification of an initial value and the field layout are available!" );
@@ -487,10 +504,10 @@ public:
    GhostLayerField_T * initialize( IBlock * const block ) override
    {
       auto blocks = blocks_.lock();
-      WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'AlwaysInitializeBlockDataHandling' for a block storage object that doesn't exist anymore" );
+      WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'AlwaysInitializeBlockDataHandling' for a block storage object that doesn't exist anymore" )
       Vector3<uint_t> size = calculateSize_( blocks, block );
       GhostLayerField_T * field = internal::allocate< GhostLayerField_T >( size[0], size[1], size[2],
-                                                                           nrOfGhostLayers_, initValue_, layout_ );
+                                                                           nrOfGhostLayers_, initValue_, layout_, alloc_ );
       if( initFunction_ )
          initFunction_( field, block );
 
@@ -505,6 +522,7 @@ private:
    Value_T initValue_;
    Layout  layout_;
    const std::function< Vector3< uint_t > ( const shared_ptr< StructuredBlockStorage > &, IBlock * const ) > calculateSize_;
+   const shared_ptr< field::FieldAllocator<Value_T> > alloc_;
 
    InitializationFunction_T initFunction_;
 
diff --git a/src/lbm/field/AddToStorage.h b/src/lbm/field/AddToStorage.h
index cd3319a8b0d3305217523ed11ff54bbd0e82a9ee..f36da19b132fa3ca2c14f25243389b33cccf4881 100644
--- a/src/lbm/field/AddToStorage.h
+++ b/src/lbm/field/AddToStorage.h
@@ -45,10 +45,10 @@ public:
 
    PdfFieldHandling( const weak_ptr< StructuredBlockStorage > & blocks, const LatticeModel_T & latticeModel,
                      const bool _initialize, const Vector3<real_t> & initialVelocity, const real_t initialDensity,
-                     const uint_t nrOfGhostLayers, const field::Layout & layout ) :
+                     const uint_t nrOfGhostLayers, const field::Layout & layout, const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr ) :
       blocks_( blocks ), latticeModel_( latticeModel ),
       initialize_( _initialize ), initialVelocity_( initialVelocity ), initialDensity_( initialDensity ),
-      nrOfGhostLayers_( nrOfGhostLayers ), layout_( layout ) {}
+      nrOfGhostLayers_( nrOfGhostLayers ), layout_( layout ), alloc_( alloc ){}
 
    inline void serialize( IBlock * const block, const BlockDataID & id, mpi::SendBuffer & buffer ) override
    {
@@ -136,7 +136,7 @@ private:
       latticeModel_.configure( *block, *blocks );
 
       return new PdfField_T( blocks->getNumberOfXCells( *block ), blocks->getNumberOfYCells( *block ), blocks->getNumberOfZCells( *block ),
-                             latticeModel_, _initialize, initialVelocity_, initialDensity, nrOfGhostLayers_, layout_ );
+                             latticeModel_, _initialize, initialVelocity_, initialDensity, nrOfGhostLayers_, layout_, alloc_ );
    }
 
    weak_ptr< StructuredBlockStorage > blocks_;
@@ -147,6 +147,7 @@ private:
    real_t            initialDensity_;
    uint_t            nrOfGhostLayers_;
    field::Layout     layout_;
+   shared_ptr< field::FieldAllocator<real_t> > alloc_;
 
 }; // class PdfFieldHandling
 
@@ -159,25 +160,26 @@ BlockDataID addPdfFieldToStorage( const shared_ptr< BlockStorage_T > & blocks, c
                                   const LatticeModel_T & latticeModel,
                                   const field::Layout & layout = field::zyxf,
                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
-                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet(),
+                                  const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr)
 {
    return blocks->addBlockData( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                   blocks, latticeModel, true, Vector3<real_t>(0), real_t(1), uint_t(1), layout ),
+                                   blocks, latticeModel, true, Vector3<real_t>(0), real_t(1), uint_t(1), layout, alloc ),
                                 identifier, requiredSelectors, incompatibleSelectors );
 }
 
 
-
 template< typename LatticeModel_T, typename BlockStorage_T >
 BlockDataID addPdfFieldToStorage( const shared_ptr< BlockStorage_T > & blocks, const std::string & identifier,
                                   const LatticeModel_T & latticeModel,
                                   const uint_t ghostLayers,
                                   const field::Layout & layout = field::zyxf,
                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
-                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet(),
+                                  const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr)
 {
    return blocks->addBlockData( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                   blocks, latticeModel, true, Vector3<real_t>(0), real_t(1), ghostLayers, layout ),
+                                   blocks, latticeModel, true, Vector3<real_t>(0), real_t(1), ghostLayers, layout, alloc ),
                                 identifier, requiredSelectors, incompatibleSelectors );
 }
 
@@ -189,10 +191,11 @@ BlockDataID addPdfFieldToStorage( const shared_ptr< BlockStorage_T > & blocks, c
                                   const Vector3< real_t > & initialVelocity, const real_t initialDensity,
                                   const field::Layout & layout = field::zyxf,
                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
-                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet(),
+                                  const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr)
 {
    return blocks->addBlockData( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                   blocks, latticeModel, true, initialVelocity, initialDensity, uint_t(1), layout ),
+                                   blocks, latticeModel, true, initialVelocity, initialDensity, uint_t(1), layout, alloc ),
                                 identifier, requiredSelectors, incompatibleSelectors );
 }
 
@@ -205,10 +208,11 @@ BlockDataID addPdfFieldToStorage( const shared_ptr< BlockStorage_T > & blocks, c
                                   const uint_t ghostLayers,
                                   const field::Layout & layout = field::zyxf,
                                   const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
-                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+                                  const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet(),
+                                  const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr)
 {
    return blocks->addBlockData( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                   blocks, latticeModel, true, initialVelocity, initialDensity, ghostLayers, layout ),
+                                   blocks, latticeModel, true, initialVelocity, initialDensity, ghostLayers, layout, alloc ),
                                 identifier, requiredSelectors, incompatibleSelectors );
 }
 
@@ -223,36 +227,40 @@ struct PdfFieldCreator : public domain_decomposition::BlockDataCreator< lbm::Pdf
    PdfFieldCreator( const shared_ptr< StructuredBlockStorage > & blocks,
                     const std::string & identifier, const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
                     const LatticeModel_T & latticeModel,
-                    const field::Layout & layout = field::zyxf ) :
+                    const field::Layout & layout = field::zyxf,
+                    const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr) :
       domain_decomposition::BlockDataCreator< lbm::PdfField< LatticeModel_T > >( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                                                                    blocks, latticeModel, false, Vector3<real_t>(0), real_t(1), uint_t(1), layout ),
+                                                                                    blocks, latticeModel, false, Vector3<real_t>(0), real_t(1), uint_t(1), layout, alloc ),
                                                                                  identifier, requiredSelectors, incompatibleSelectors )
    {}
 
    PdfFieldCreator( const shared_ptr< StructuredBlockStorage > & blocks,
                     const std::string & identifier, const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
                     const LatticeModel_T & latticeModel, const uint_t ghostLayers,
-                    const field::Layout & layout = field::zyxf ) :
+                    const field::Layout & layout = field::zyxf,
+                    const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr) :
       domain_decomposition::BlockDataCreator< lbm::PdfField< LatticeModel_T > >( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                                                                    blocks, latticeModel, false, Vector3<real_t>(0), real_t(1), ghostLayers, layout ),
+                                                                                    blocks, latticeModel, false, Vector3<real_t>(0), real_t(1), ghostLayers, layout, alloc ),
                                                                                  identifier, requiredSelectors, incompatibleSelectors )
    {}
 
    PdfFieldCreator( const shared_ptr< StructuredBlockStorage > & blocks,
                     const std::string & identifier, const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
                     const LatticeModel_T & latticeModel, const Vector3< real_t > & initialVelocity, const real_t initialDensity,
-                    const field::Layout & layout = field::zyxf ) :
+                    const field::Layout & layout = field::zyxf,
+                    const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr) :
       domain_decomposition::BlockDataCreator< lbm::PdfField< LatticeModel_T > >( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                                                                    blocks, latticeModel, true, initialVelocity, initialDensity, uint_t(1), layout ),
+                                                                                    blocks, latticeModel, true, initialVelocity, initialDensity, uint_t(1), layout, alloc ),
                                                                                  identifier, requiredSelectors, incompatibleSelectors )
    {}
 
    PdfFieldCreator( const shared_ptr< StructuredBlockStorage > & blocks,
                     const std::string & identifier, const Set<SUID> & requiredSelectors, const Set<SUID> & incompatibleSelectors,
                     const LatticeModel_T & latticeModel, const Vector3< real_t > & initialVelocity, const real_t initialDensity, const uint_t ghostLayers,
-                    const field::Layout & layout = field::zyxf ) :
+                    const field::Layout & layout = field::zyxf,
+                    const shared_ptr< field::FieldAllocator<real_t> > alloc = nullptr) :
       domain_decomposition::BlockDataCreator< lbm::PdfField< LatticeModel_T > >( make_shared< internal::PdfFieldHandling< LatticeModel_T > >(
-                                                                                    blocks, latticeModel, true, initialVelocity, initialDensity, ghostLayers, layout ),
+                                                                                    blocks, latticeModel, true, initialVelocity, initialDensity, ghostLayers, layout, alloc ),
                                                                                  identifier, requiredSelectors, incompatibleSelectors )
    {}
 };
diff --git a/tests/lbm/CMakeLists.txt b/tests/lbm/CMakeLists.txt
index a2ddd66905fda5cbc57601163746d1d96ab58421..87d27692a08414ffb67ab69f9c0d1f9ddea37d00 100644
--- a/tests/lbm/CMakeLists.txt
+++ b/tests/lbm/CMakeLists.txt
@@ -74,6 +74,10 @@ waLBerla_execute_test( NAME  SuViscoelasticityShortTest COMMAND $<TARGET_FILE:Su
 waLBerla_compile_test( FILES field/QCriterionTest.cpp DEPENDS field blockforest )
 waLBerla_execute_test( NAME QCriterionTest )
 
+add_subdirectory(field)
+waLBerla_compile_test( FILES field/AlignedPDFfieldTest.cpp DEPENDS field blockforest lbm )
+waLBerla_execute_test( NAME AlignedPDFfieldTest COMMAND $<TARGET_FILE:AlignedPDFfieldTest> ${CMAKE_CURRENT_SOURCE_DIR}/field/AlignedPDFfieldTest.prm )
+
 # Code Generation
 if( WALBERLA_BUILD_WITH_CODEGEN )
 add_subdirectory(codegen)
diff --git a/tests/lbm/field/AlignedPDFfieldTest.cpp b/tests/lbm/field/AlignedPDFfieldTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7c02cefa9afb1fbec5482ce27e6adc68185e054
--- /dev/null
+++ b/tests/lbm/field/AlignedPDFfieldTest.cpp
@@ -0,0 +1,80 @@
+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file AlignedPDFfieldTest.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+
+#include "core/all.h"
+
+#include "field/allocation/FieldAllocator.h"
+
+#include "lbm/all.h"
+
+namespace walberla
+{
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::SRT, false >;
+typedef lbm::PdfField< LatticeModel_T > PdfField_T;
+
+int main(int argc, char** argv)
+{
+   debug::enterTestMode();
+   walberla::Environment walberlaEnv(argc, argv);
+
+   auto blocks = blockforest::createUniformBlockGridFromConfig(walberlaEnv.config());
+
+   auto parameters    = walberlaEnv.config()->getOneBlock("Parameters");
+   const real_t omega = parameters.getParameter< real_t >("omega", real_c(1.4));
+
+   LatticeModel_T latticeModel = LatticeModel_T(omega);
+
+   shared_ptr< field::FieldAllocator< real_t > > alloc_32 = make_shared< field::AllocateAligned< real_t, 32 > >();
+   BlockDataID pdfFieldId_32 = lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, field::fzyx, Set<SUID>::emptySet(), Set<SUID>::emptySet(), alloc_32);
+
+   shared_ptr< field::FieldAllocator< real_t > > alloc_64 = make_shared< field::AllocateAligned< real_t, 64 > >();
+   BlockDataID pdfFieldId_64 = lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, field::fzyx, Set<SUID>::emptySet(), Set<SUID>::emptySet(), alloc_64);
+
+   shared_ptr< field::FieldAllocator< real_t > > alloc_128 = make_shared< field::AllocateAligned< real_t, 128 > >();
+   BlockDataID pdfFieldId_128 = lbm::addPdfFieldToStorage(blocks, "pdf field", latticeModel, field::fzyx, Set<SUID>::emptySet(), Set<SUID>::emptySet(), alloc_128);
+
+   for (auto& block : *blocks)
+   {
+      auto pdfField_32 = block.getData< PdfField_T >(pdfFieldId_32);
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_32->dataAt(0, 0, 0, 0) % 32, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_32->dataAt(0, 0, 0, 1) % 32, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_32->dataAt(0, 0, 1, 0) % 32, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_32->dataAt(0, 1, 0, 0) % 32, 0)
+
+      auto pdfField_64 = block.getData< PdfField_T >(pdfFieldId_64);
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_64->dataAt(0, 0, 0, 0) % 64, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_64->dataAt(0, 0, 0, 1) % 64, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_64->dataAt(0, 0, 1, 0) % 64, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_64->dataAt(0, 1, 0, 0) % 64, 0)
+
+      auto pdfField_128 = block.getData< PdfField_T >(pdfFieldId_128);
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_128->dataAt(0, 0, 0, 0) % 128, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_128->dataAt(0, 0, 0, 1) % 128, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_128->dataAt(0, 0, 1, 0) % 128, 0)
+      WALBERLA_CHECK_EQUAL((size_t) pdfField_128->dataAt(0, 1, 0, 0) % 128, 0)
+   }
+
+   return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::main(argc, argv); }
\ No newline at end of file
diff --git a/tests/lbm/field/AlignedPDFfieldTest.prm b/tests/lbm/field/AlignedPDFfieldTest.prm
new file mode 100644
index 0000000000000000000000000000000000000000..16e1bafce4c675bbcc8580c13f089ba46918f0e5
--- /dev/null
+++ b/tests/lbm/field/AlignedPDFfieldTest.prm
@@ -0,0 +1,12 @@
+
+Parameters
+{
+	omega           1.8;
+}
+
+DomainSetup
+{
+   blocks        <  1,   1,  1 >;
+   cellsPerBlock <  10, 10, 10 >;
+   periodic      <  0,   0,  0 >;
+}
diff --git a/tests/lbm/field/CMakeLists.txt b/tests/lbm/field/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7170fde6d819ad39d446a0ac82cd91e1b069b773
--- /dev/null
+++ b/tests/lbm/field/CMakeLists.txt
@@ -0,0 +1,7 @@
+#############################################################################################################################
+#
+# Tests for lbm module
+#
+#############################################################################################################################
+
+waLBerla_link_files_to_builddir( "*.prm" )
\ No newline at end of file