Improve Stability Checker

apps/tutorials/lbm/01_BasicLBM.cpp

@@ -93,8 +93,9 @@ int main( int argc, char ** argv )
    timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID ) ), "LB stream & collide" );
 
    // LBM stability check
+   auto checkFunction = [](PdfField_T::value_type value) {return math::finite( value );};
    timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >( walberlaEnv.config(), blocks, pdfFieldId,
-                                                                                                             flagFieldId, fluidFlagUID ) ),
+                                                                                                             flagFieldId, fluidFlagUID, checkFunction ) ),
                                   "LBM stability check" );
 
    // log remaining time

apps/tutorials/lbm/01_BasicLBM.dox

@@ -261,10 +261,15 @@ The StabilityChecker instance can be controlled via the configuration file,
 for more information see \ref docStabilityChecker.
 Since field::makeStabilityChecker returns a shared pointer, we use makeSharedFunctor in order to wrap the shared pointer into an object
 that can be passed to the time loop.
+Note that NaNs are not defined if waLBerla is build using FASTMATH. For this case the field::StabilityChecker accepts
+a checkFunction as input argument. The checkFunction receives a value of the type the field::StabilityChecker is applied on
+and returns a bool. This function is applied on each value on each cell. If no checkFunction is provided a default is used
+which is exactly the one shown in the code below.
 
 \code
+auto checkFunction = [](PdfField_T::value_type value) {return math::finite( value );};
 timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< PdfField_T, FlagField_T >(
-                                     walberlaEnv.config(), blocks, pdfFieldId, flagFieldId, fluidFlagUID ) ),
+                               walberlaEnv.config(), blocks, pdfFieldId, flagFieldId, fluidFlagUID, checkFunction ) ),
                                "LBM stability check" );
 \endcode
 

src/field/EvaluationFilter.h

@@ -74,13 +74,13 @@ public:
    void operator()( const IBlock & block )
    {
       flagField_ = block.template getData< const FlagField_T >( flagFieldId_ );
-      WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField_ )
       evaluationMask_ = flagField_->getMask( cellsToEvaluate_ );
    }
 
    bool operator()( const cell_idx_t x, const cell_idx_t y, const cell_idx_t z ) const
    {
-      WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
+      WALBERLA_ASSERT_NOT_NULLPTR( flagField_ )
       return flagField_->isPartOfMaskSet( x, y, z, evaluationMask_ );
    }
 

src/waLBerlaDefinitions.in.h

@@ -16,6 +16,7 @@
 
 // Debugging options
 #cmakedefine WALBERLA_ENABLE_GUI
+#cmakedefine WALBERLA_BUILD_WITH_FASTMATH
 
 
 // External libraries

tests/field/CMakeLists.txt

@@ -28,6 +28,9 @@ waLBerla_execute_test( NAME FieldTiming  )
 waLBerla_compile_test( FILES FlagFieldTest.cpp)
 waLBerla_execute_test( NAME FlagFieldTest )
 
+waLBerla_compile_test( FILES StabilityCheckerTest.cpp DEPENDS blockforest field timeloop )
+waLBerla_execute_test( NAME StabilityCheckerTest )
+
 waLBerla_compile_test( FILES interpolators/InterpolationTest.cpp)
 waLBerla_execute_test( NAME InterpolationTest )
 

tests/field/StabilityCheckerTest.cpp

+//======================================================================================================================
+//
+//  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 StabilityCheckerTest.cpp
+//! \ingroup field
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/all.h"
+#include "core/all.h"
+#include "domain_decomposition/all.h"
+#include "field/all.h"
+#include "timeloop/all.h"
+
+
+
+namespace walberla {
+
+using Field_T = GhostLayerField<real_t, 1>;
+
+class TestSweep
+{
+ public:
+   TestSweep(BlockDataID fieldID) : fieldID_(fieldID) {}
+
+   void operator()(IBlock* const block)
+   {
+      Field_T* field = block->getData< Field_T >(fieldID_);
+
+      WALBERLA_FOR_ALL_CELLS(fieldIt, field, { *fieldIt += Field_T::value_type(1); }) // WALBERLA_FOR_ALL_CELLS
+   }
+
+ private:
+   BlockDataID fieldID_;
+};
+
+int main( int argc, char ** argv )
+{
+   debug::enterTestMode();
+   walberla::Environment walberlaEnv( argc, argv );
+
+   auto blocks = blockforest::createUniformBlockGrid( 1, 1, 1,
+                                                      4, 4, 4,
+                                                      1.0);
+
+   BlockDataID fieldID = field::addToStorage<Field_T>( blocks, "Field", Field_T::value_type(0));
+   SweepTimeloop timeloop(blocks->getBlockStorage(), uint_c(2));
+
+   timeloop.add() << Sweep(TestSweep(fieldID), "Test Sweep");
+
+   // LBM stability check
+   auto checkFunction = [](Field_T::value_type value) {return value < math::abs(Field_T::value_type(5));};
+   timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< Field_T >( blocks, fieldID, uint_c(1), checkFunction) ),"Stability check" );
+   timeloop.run();
+
+   return EXIT_SUCCESS;
+}
+}
+
+int main( int argc, char ** argv )
+{
+   walberla::main(argc, argv);
+}