First implimentation of stability checker for list

apps/benchmarks/ComplexGeometryList/listStabillityChecker.h

+//
+// Created by po60nani on 6/4/18.
+//
+
+#ifndef WALBERLA_LISTSTABILLITYCHECKER_H
+#define WALBERLA_LISTSTABILLITYCHECKER_H
+//======================================================================================================================
+//
+//  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 StabilityCheckerer.h
+//! \ingroup field
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/Set.h"
+#include "core/cell/CellSet.h"
+#include "core/debug/CheckFunctions.h"
+#include "core/logging/Logging.h"
+#include "core/math/Vector3.h"
+#include "core/mpi/Reduce.h"
+#include "core/uid/SUID.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/EvaluationFilter.h"
+#include "field/iterators/IteratorMacros.h"
+
+#include "vtk/BlockCellDataWriter.h"
+#include "vtk/DumpBlockStructureLevel.h"
+#include "vtk/DumpBlockStructureProcess.h"
+#include "vtk/VTKOutput.h"
+
+
+
+namespace walberla {
+    namespace field {
+
+        namespace internal {
+
+            const std::string stabilityCheckerVTKBase("vtk_out");
+            const std::string stabilityCheckerVTKFolder("output");
+            const std::string stabilityCheckerVTKIdentifier("error_field");
+
+            const bool stabilityCheckerVTKBinary(true);
+            const bool stabilityCheckerVTKLittleEndian(true);
+            const bool stabilityCheckerVTKMPIIO(true);
+            const bool stabilityCheckerVTKForcePVTU(false);
+
+            const std::string stabilityCheckerConfigBlock("StabilityChecker");
+
+            template<typename T>
+            inline bool stabilityCheckerIsFinite(const T &value) { return math::finite(value); }
+
+            template<>
+            inline bool stabilityCheckerIsFinite(const Vector3<real_t> &value) {
+                return math::finite(value[0]) && math::finite(value[1]) && math::finite(value[2]);
+            }
+
+        }
+
+        template<typename List_T, typename LatticeModel_T >
+        class StabilityChecker {
+
+        private:
+
+            typedef std::map< const IBlock *, std::map< Cell, std::set< cell_idx_t > > > BlockCellsMap;
+
+            ////////////////////////////////
+            //  VTK Output Helper Classes //
+            ////////////////////////////////
+
+            /// This cell filter selects only those cells in which at least one non-finite value (= infinite or NaN) was detected
+            class VTKCellFilter
+            {
+            public:
+
+                VTKCellFilter( BlockCellsMap & map ) : map_( map ) {}
+
+                void operator()( CellSet & filteredCells, const IBlock & block, const StructuredBlockStorage & storage, const uint_t ghostLayers ) const
+                {
+                    if( map_.find( &block ) != map_.end() )
+                    {
+                        const auto & cellMap = map_[ &block ];
+                        if( !cellMap.empty() )
+                        {
+                            const cell_idx_t gl    = cell_idx_c( ghostLayers );
+                            const cell_idx_t begin = cell_idx_c( -1 ) * gl;
+
+                            for( cell_idx_t z = begin; z < cell_idx_c( storage.getNumberOfZCells(block) ) + gl; ++z )
+                                for( cell_idx_t y = begin; y < cell_idx_c( storage.getNumberOfYCells(block) ) + gl; ++y )
+                                    for( cell_idx_t x = begin; x < cell_idx_c( storage.getNumberOfXCells(block) ) + gl; ++x )
+                                        if( cellMap.find( Cell(x,y,z) ) != cellMap.end() ) filteredCells.insert( x, y, z );
+                        }
+                    }
+                }
+
+            private:
+
+                BlockCellsMap & map_;
+            };
+
+            /// For each value of a cell, either '0' or '1' is stored in the VTK file.
+            /// If the value is non-finite (= infinite or NaN), '1' is written to file - otherwise '0'.
+        class FValueVTKWriter : public vtk::BlockCellDataWriter< uint8_t, List_T::Stencil::Q >
+            {
+            public:
+
+                FValueVTKWriter( BlockCellsMap & map, const std::string & id ) : vtk::BlockCellDataWriter< uint8_t, List_T::Stencil::Q >( id ), map_( map ) {}
+
+            protected:
+
+                void configure() {}
+
+                uint8_t evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t f )
+                {
+                    WALBERLA_ASSERT( map_.find( this->block_ ) != map_.end() );
+                    WALBERLA_ASSERT( map_[ this->block_ ].find( Cell(x,y,z) ) != map_[ this->block_ ].end() );
+
+                    return ( map_[ this->block_ ][ Cell(x,y,z) ].find( f ) != map_[ this->block_ ][ Cell(x,y,z) ].end() ) ? uint8_t(1) : uint8_t(0);
+                }
+
+            private:
+
+                BlockCellsMap & map_;
+            };
+
+            /// For each cell, the corresponding block local cell coordinates are stored in the VTK file.
+            class LocalCoordVTKWriter : public vtk::BlockCellDataWriter< cell_idx_t, 3 >
+            {
+            public:
+                LocalCoordVTKWriter( const std::string & id ) : vtk::BlockCellDataWriter< cell_idx_t, 3 >( id ) {}
+            protected:
+                void configure() {}
+                cell_idx_t evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t f )
+                {
+                    return ( f == cell_idx_t(0) ) ? x : ( ( f == cell_idx_t(1) ) ? y : z );
+                }
+            };
+
+            /// For each cell, the corresponding global cell coordinates are stored in the VTK file.
+            class GlobalCoordVTKWriter : public vtk::BlockCellDataWriter< cell_idx_t, 3 >
+            {
+            public:
+                GlobalCoordVTKWriter( const std::string & id ) : vtk::BlockCellDataWriter< cell_idx_t, 3 >( id ) {}
+            protected:
+                void configure() {}
+                cell_idx_t evaluate( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z, const cell_idx_t f )
+                {
+                    Cell cell(x,y,z);
+                    this->blockStorage_->transformBlockLocalToGlobalCell( cell, *(this->block_) );
+                    return ( f == cell_idx_t(0) ) ? cell.x(): ( ( f == cell_idx_t(1) ) ? cell.y() : cell.z() );
+                }
+            };
+
+
+
+        public:
+
+            StabilityChecker( const weak_ptr< StructuredBlockStorage > & blocks, const ConstBlockDataID & pdfListId,
+                              const uint_t checkFrequency,
+                              const bool outputToStream = true, const bool outputVTK = true,
+                              const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
+                              const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() ) :
+                    blocks_( blocks ), executionCounter_( uint_c(0) ), checkFrequency_( checkFrequency ),
+                    pdfListId_( pdfListId ), outputToStream_( outputToStream ), outputVTK_( outputVTK ),
+                    vtkBaseFolder_( internal::stabilityCheckerVTKBase ),
+                    vtkExecutionFolder_( internal::stabilityCheckerVTKFolder ),
+                    vtkIdentifier_( internal::stabilityCheckerVTKIdentifier ),
+                    vtkBinary_( internal::stabilityCheckerVTKBinary ),
+                    vtkLittleEndian_( internal::stabilityCheckerVTKLittleEndian ),
+                    vtkMPIIO_( internal::stabilityCheckerVTKMPIIO ),
+                    vtkForcePVTU_( internal::stabilityCheckerVTKForcePVTU ),
+                    requiredSelectors_(requiredSelectors), incompatibleSelectors_( incompatibleSelectors ) {}
+
+
+            void setVTKBaseFolder     ( const std::string & vtkBaseFolder      ) { vtkBaseFolder_      = vtkBaseFolder; }
+            void setVTKExecutionFolder( const std::string & vtkExecutionFolder ) { vtkExecutionFolder_ = vtkExecutionFolder; }
+            void setVTKIdentifier     ( const std::string & vtkIdentifier      ) { vtkIdentifier_      = vtkIdentifier; }
+
+            void setVTKBinary      ( const bool vtkBinary       ) { vtkBinary_       = vtkBinary; }
+            void setVTKLittleEndian( const bool vtkLittleEndian ) { vtkLittleEndian_ = vtkLittleEndian; }
+            void setVTKMPIIO       ( const bool vtkMPIIO        ) { vtkMPIIO_        = vtkMPIIO; }
+            void setVTKForcePVTU   ( const bool vtkForcePVTU    ) { vtkForcePVTU_    = vtkForcePVTU; }
+
+            void operator()();
+
+        private:
+
+            void checkBlock( const IBlock * const block );
+
+
+
+            weak_ptr< StructuredBlockStorage > blocks_;
+
+
+            uint_t executionCounter_;
+            uint_t checkFrequency_;
+
+            ConstBlockDataID pdfListId_;
+
+            BlockCellsMap failedCells_;
+
+            bool outputToStream_;
+            bool outputVTK_;
+
+            std::string vtkBaseFolder_;
+            std::string vtkExecutionFolder_;
+            std::string vtkIdentifier_;
+
+            bool vtkBinary_;
+            bool vtkLittleEndian_;
+            bool vtkMPIIO_;
+            bool vtkForcePVTU_;
+
+            Set<SUID> requiredSelectors_;
+            Set<SUID> incompatibleSelectors_;
+
+
+        }; // StabilityChecker
+
+
+        template< typename List_T, typename LatticeModel_T >
+        void StabilityChecker< List_T, LatticeModel_T >::operator()()
+        {
+            auto blocks = blocks_.lock();
+            WALBERLA_CHECK_NOT_NULLPTR( blocks, "Trying to access 'StabilityChecker' for a block storage object that doesn't exist anymore" );
+
+            for( auto block = blocks->begin( requiredSelectors_, incompatibleSelectors_ ); block != blocks->end(); ++block )
+                checkBlock( block.get() );
+
+            bool abort = !failedCells_.empty();
+            mpi::allReduceInplace( abort, mpi::LOGICAL_OR );
+
+            if( abort )
+            {
+                if( outputVTK_ )
+                {
+                    auto vtkWriter = vtk::createVTKOutput_BlockData( blocks, vtkIdentifier_, uint_t(1), uint_t(0), vtkForcePVTU_,
+                                                                     vtkBaseFolder_, vtkExecutionFolder_, false, vtkBinary_, vtkLittleEndian_, vtkMPIIO_ );
+
+                    vtkWriter->addCellInclusionFilter( VTKCellFilter( failedCells_ ) );
+
+                    vtkWriter->addCellDataWriter( walberla::make_shared< vtk::DumpBlockStructureProcess >( "process" ) );
+                    vtkWriter->addCellDataWriter( walberla::make_shared< vtk::DumpBlockStructureLevel >( "level" ) );
+                    vtkWriter->addCellDataWriter( walberla::make_shared< FValueVTKWriter >( boost::ref( failedCells_ ), "F" ) );
+                    vtkWriter->addCellDataWriter( walberla::make_shared< LocalCoordVTKWriter >( "blockLocalCoordinate" ) );
+                    vtkWriter->addCellDataWriter( walberla::make_shared< GlobalCoordVTKWriter >( "globalCoordinate" ) );
+
+                    vtkWriter->write();
+                }
+
+                WALBERLA_LOG_WARNING_ON_ROOT( "Field stability check failed - aborting program ..." );
+                WALBERLA_MPI_WORLD_BARRIER();
+
+                WALBERLA_ABORT_NO_DEBUG_INFO("");
+            }
+        }
+
+
+        template< typename List_T, typename LatticeModel_T >
+        void StabilityChecker< List_T, LatticeModel_T >::checkBlock( const IBlock * const block ) {
+            const List_T *lbmList = block->getData<List_T>(pdfListId_);
+            WALBERLA_CHECK_NOT_NULLPTR(lbmList);
+
+
+            for (uint_t iter = 0; iter != lbmList->numFluidCells(); ++iter) {
+                for (auto dirIt = LatticeModel_T::Stencil::begin(); dirIt != LatticeModel_T::Stencil::end(); ++dirIt) {
+
+                    auto cellPdf = lbmList->get(iter, *dirIt);
+
+                    if (math::isinf(cellPdf) || math::isnan(cellPdf)) {
+                        failedCells_[block][ lbmList->getCell(iter) ].insert(*dirIt);
+
+
+                        if( outputToStream_ ) {
+                            WALBERLA_LOG_WARNING ("#block = " << block << "\n #cell = " << iter << "\n cell position = " << lbmList->getCell(iter) << "\n #direction = " << *dirIt);
+                        }
+                    }
+                }
+            }
+        }
+
+
+
+
+///////////////////////////////////////////////////////////////
+// makeListStabilityChecker functions without configuration file //
+///////////////////////////////////////////////////////////////
+        template< typename List_T, typename LatticeModel_T >
+        shared_ptr< StabilityChecker< List_T, LatticeModel_T > > makeListStabilityChecker( const weak_ptr< StructuredBlockStorage > & blocks, const ConstBlockDataID & pdfListId,
+                                                                        const uint_t checkFrequency,
+                                                                        const bool outputToStream = true, const bool outputVTK = true,
+                                                                        const Set<SUID> & requiredSelectors     = Set<SUID>::emptySet(),
+                                                                        const Set<SUID> & incompatibleSelectors = Set<SUID>::emptySet() )
+        {
+            typedef StabilityChecker< List_T, LatticeModel_T > SC_T;
+            return shared_ptr< SC_T >( new SC_T( blocks, pdfListId, checkFrequency, outputToStream, outputVTK, requiredSelectors, incompatibleSelectors ) );
+        }
+
+
+
+    } // namespace field
+} // namespace walberla
+
+
+
+
+#endif //WALBERLA_LISTSTABILLITYCHECKER_H