From 0def78855ddd0ec5055179b5034e097333963f28 Mon Sep 17 00:00:00 2001
From: Cameron Stewart <cstewart@icp.uni-stuttgart.de>
Date: Wed, 28 Nov 2018 10:55:41 +0100
Subject: [PATCH] test first draft

---
 tests/lbm/Su.prm                    |  28 +--
 tests/lbm/SuViscoelasticityTest.cpp | 342 ++++++++++++++--------------
 2 files changed, 185 insertions(+), 185 deletions(-)

diff --git a/tests/lbm/Su.prm b/tests/lbm/Su.prm
index 470b5dea2..6c215a2cb 100644
--- a/tests/lbm/Su.prm
+++ b/tests/lbm/Su.prm
@@ -1,35 +1,25 @@
 Parameters 
 {
 	eta_s           0.5;
-	timesteps       1000;
+	timesteps       3000;
 	force           0.00001;
-	lambda_p        30.0;
+	lambda_p        300.0;
 	eta_p           0.5;
     period          1;
-    baseFolder      vtk_out;
-    writePeriod     1;
-
-    checkpointFrequency 1000;
-	remainingTimeLoggerFrequency 6.0; // in seconds
+    L               11;
+    H               33;
+    uMax            1.28737;
+    t0              370;
+    t1              326;
 }
 
 DomainSetup
 {
-   blocks        <  3,3,1 >;
+   blocks        <  1,3,1 >;
    cartesianSetup false;
-   cellsPerBlock <  10,10, 1 >;
+   cellsPerBlock <  11,11, 1 >;
    periodic      <  1, 0, 0 >;
 }
 
-StabilityChecker
-{
-   checkFrequency 1;
-   streamOutput   false;
-   vtkOutput      true;
-}
 
-Boundaries 
-{
-	    Border { direction N,S;  walldistance -1;  NoSlip    {} }
-}
 
diff --git a/tests/lbm/SuViscoelasticityTest.cpp b/tests/lbm/SuViscoelasticityTest.cpp
index cd63004be..cecac2e28 100644
--- a/tests/lbm/SuViscoelasticityTest.cpp
+++ b/tests/lbm/SuViscoelasticityTest.cpp
@@ -66,138 +66,92 @@ using namespace walberla;
 // TYPEDEFS //
 //////////////
 
+
 typedef GhostLayerField< Vector3<real_t>, 1> ForceField_T;
-typedef lbm::D2Q9< lbm::collision_model::TRT, false, lbm::force_model::GuoField< ForceField_T > >  LatticeModel_T;
-typedef LatticeModel_T::CommunicationStencil    CommunicationStencil_T;
+
+typedef lbm::collision_model::TRT CollisionModel_T;
+typedef lbm::force_model::GuoField< ForceField_T > ForceModel_T;
+typedef lbm::D2Q9< CollisionModel_T, false, ForceModel_T >  LatticeModel_T;
+typedef LatticeModel_T::Stencil    Stencil_T;
+
 typedef GhostLayerField<Matrix3<real_t>, 1> StressField_T;
 typedef GhostLayerField< Vector3<real_t>, 1> VelocityField_T;
 typedef lbm::PdfField< LatticeModel_T >  PdfField_T;
 typedef walberla::uint8_t flag_t;
 typedef FlagField< flag_t > FlagField_T;
-typedef lbm::ExtendedBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
+const uint_t FieldGhostLayers = 2;
 
-class WriteCheckpoint {
-public:
-   WriteCheckpoint(Timeloop &timeloop, BlockStorage & blockStorage, uint_t checkpointFrequency, uint_t timesteps,
-                   BlockDataID pdfFieldId, BlockDataID stressId, std::string dir, std::string paramFile) : timeloop_(timeloop), blockStorage_(blockStorage),
-                                                                     checkpointFrequency_(checkpointFrequency), timesteps_(timesteps),
-                                                                     pdfFieldId_(pdfFieldId), stressId_(stressId), dir_(dir + "_chk"), paramFile_(paramFile) {}
-   void operator()(){
-      if(checkpointFrequency_ > 0 && timeloop_.getCurrentTimeStep() % checkpointFrequency_== 0 && timeloop_.getCurrentTimeStep() > 0){
-         Save();
-      }
-      if(timeloop_.getCurrentTimeStep() == timesteps_-1 && filesystem::is_directory(dir_)){
-         filesystem::remove_all(dir_);
-      }
-   }
+typedef lbm::NoSlip< LatticeModel_T, flag_t> NoSlip_T;
+typedef boost::tuples::tuple< NoSlip_T > BoundaryConditions_T;
+typedef BoundaryHandling< FlagField_T, Stencil_T, BoundaryConditions_T> BoundaryHandling_T;
 
-   bool canLoad()
-   {
-      if(!filesystem::is_directory(dir_)){
-         std::cout << "No Checkpoint Found" << std::endl;
-         return false;
-      }
+///////////
+// FLAGS //
+///////////
 
-      std::ifstream f1;
-      std::ifstream f2;
-      f1.open(paramFile_, std::ifstream::binary|std::ifstream::ate);
-      f2.open(dir_ + "/" + paramFile_, std::ifstream::binary|std::ifstream::ate);
+const FlagUID Fluid_Flag( "fluid" );
+const FlagUID NoSlip_Flag( "no slip" );
 
-      if (f1.fail() || f2.fail()) {
-         std::cout << "Failed to check parameter file" << std::endl;
-         return false; //file problem
-      }
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
 
-      if (f1.tellg() != f2.tellg()) {
-         std::cout << "Ignoring checkpoint for different size config file" << std::endl;
-         return false; //size mismatch
-      }
+class MyBoundaryHandling
+{
+public:
 
-      f1.seekg(0, std::ifstream::beg);
-      f2.seekg(0, std::ifstream::beg);
-      if(!std::equal(std::istreambuf_iterator<char>(f1.rdbuf()), std::istreambuf_iterator<char>(), std::istreambuf_iterator<char>(f2.rdbuf()))){
-         std::cout << "Ignoring checkpoint for different config file" << std::endl;
-         return false;
-      }
-      return true;
-   }
+   MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID ) : flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ) {}
 
-   void Load()
-   {
-      std::cout << "Loading Checkpoint" << std::endl;
-      field::readFromFile<PdfField_T>(dir_ + "/pdfField.chk", blockStorage_, pdfFieldId_);
-      field::readFromFile<StressField_T>(dir_ + "/stressField.chk", blockStorage_, stressId_);
-      std::ifstream f;
-      f.open(dir_ + "/chkpoint");
-      uint_t currentTime;
-      f >> currentTime;
-      timeloop_.setCurrentTimeStep(currentTime + 1);
-   }
+   BoundaryHandling_T * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const;
 
-   void Save(){
+private:
 
-      auto tmpdir = dir_ + "~";
-      if(filesystem::is_directory(tmpdir)){
-         filesystem::remove_all(tmpdir);
-      }
+   const BlockDataID flagFieldID_;
+   const BlockDataID pdfFieldID_;
 
-      filesystem::create_directory(tmpdir);
+}; // class MyBoundaryHandling
 
-      filesystem::copy(paramFile_, tmpdir + "/" + paramFile_);
-      field::writeToFile<PdfField_T>(tmpdir + "/pdfField.chk", blockStorage_, pdfFieldId_);
-      field::writeToFile<StressField_T>(tmpdir + "/stressField.chk", blockStorage_, stressId_);
-      std::ofstream f;
-      f.open(tmpdir + "/chkpoint");
-      f << timeloop_.getCurrentTimeStep();
-      f.close();
 
-      if(filesystem::is_directory(dir_)) {
-         filesystem::remove_all(dir_);
-      }
-      filesystem::rename(tmpdir, dir_);
-   }
+BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const StructuredBlockStorage * const storage ) const
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( block );
+   WALBERLA_ASSERT_NOT_NULLPTR( storage );
 
-private:
-   Timeloop & timeloop_;
-   BlockStorage & blockStorage_;
-   uint_t checkpointFrequency_, timesteps_;
-   BlockDataID pdfFieldId_, stressId_;
-   std::string dir_, paramFile_;
-};
+   FlagField_T * flagField       = block->getData< FlagField_T >( flagFieldID_ );
+   PdfField_T *  pdfField        = block->getData< PdfField_T > ( pdfFieldID_ );
 
-class WriteVelocity
-{
-public:
-   WriteVelocity(Timeloop & timeloop, uint_t writeFrequency, shared_ptr< StructuredBlockForest > blocks, BlockDataID pdfFieldId) : timeloop_(timeloop), writeFrequency_(writeFrequency), blocks_(blocks),
-   pdfFieldId_(pdfFieldId){}
+   const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
 
-   void operator()(){
-      if(writeFrequency_ > 0 && timeloop_.getCurrentTimeStep() % writeFrequency_==0 ){
+   BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
+                                                           boost::tuples::make_tuple(    NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ) ) );
 
-         std::ofstream f;
-         f.open("vel.dat", std::ios_base::app);
-         for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
-         {
-            PdfField_T * pdf = blockIt.get()->getData< PdfField_T >(pdfFieldId_);
+   const auto noSlip = flagField->getFlag(NoSlip_Flag);
 
-            if(blockIt.get()->getAABB().contains(15,15,0)){
-               std::cout<< pdf->getVelocity(5,5,0).length() << std::endl;
-               f << std::setprecision(10) << pdf->getVelocity(5,5,0).length() << ",";
-            }
+   CellInterval domainBB = storage->getDomainCellBB();
+   domainBB.xMin() -= cell_idx_c( 1 );
+   domainBB.xMax() += cell_idx_c( 1 );
 
-         }
-         f.close();
-      }
-   }
+   domainBB.yMin() -= cell_idx_c( 1 );
+   domainBB.yMax() += cell_idx_c( 1 );
 
-private:
-   Timeloop & timeloop_;
-   uint_t writeFrequency_;
-   shared_ptr< StructuredBlockForest > blocks_;
-   BlockDataID pdfFieldId_;
-};
+   // SOUTH
+   CellInterval south( domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMin(), domainBB.zMax() );
+   storage->transformGlobalToBlockLocalCellInterval( south, *block );
+   handling->forceBoundary( noSlip, south );
+
+   // NORTH
+   CellInterval north( domainBB.xMin(), domainBB.yMax(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax() );
+   storage->transformGlobalToBlockLocalCellInterval( north, *block );
+   handling->forceBoundary( noSlip, north );
 
+   handling->fillWithDomain( FieldGhostLayers );
+
+   return handling;
+}
 
+//////////////////////
+// Poiseuille Force //
+//////////////////////
 
 template < typename BoundaryHandling_T >
 class ConstantForce
@@ -228,96 +182,152 @@ private:
    real_t force_;
 };
 
+////////////////////
+// Take Test Data //
+////////////////////
+
+class TestData
+{
+public:
+   TestData(Timeloop & timeloop, shared_ptr< StructuredBlockForest > blocks, BlockDataID pdfFieldId, BlockDataID stressFieldId, uint_t timesteps, uint_t blockSize, real_t L, real_t H, real_t uExpected)
+            : timeloop_(timeloop), blocks_(blocks), pdfFieldId_(pdfFieldId), stressFieldId_(stressFieldId), timesteps_(timesteps), blockSize_(blockSize), t0_(0), t1_(0), L_(L), H_(H), uMax_(0.0), uPrev_(0.0),
+            uExpected_(uExpected), uSteady_(0.0) {}
+
+   void operator()() {
+      for (auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt) {
+         PdfField_T *pdf = blockIt.get()->getData<PdfField_T>(pdfFieldId_);
+
+         if (blockIt.get()->getAABB().contains(L_/2.0, H_/2.0, 0)) {
+            uCurr_ = pdf->getVelocity(int32_c(blockSize_/2), int32_c(blockSize_/2), 0)[0]/uExpected_;
+            tCurr_ = timeloop_.getCurrentTimeStep();
+            if (tCurr_ == timesteps_ - 1){
+               uSteady_ = uCurr_;
+            }
+            if (maxFlag_ == 0) {
+               if (uCurr_ >= uPrev_) {
+                  uMax_ = uCurr_;
+               } else {
+                  t0_ = tCurr_;
+                  maxFlag_ = 1;
+               }
+               uPrev_ = uCurr_;
+            }
+            else if (maxFlag_ == 1) {
+               if ((uCurr_ - 1.0) <= (uMax_ - 1.0) / std::exp(1)) {
+                  t1_ = tCurr_ - t0_;
+                  maxFlag_ = 2;
+               }
+            }
+         }
+      }
+   }
+
+   real_t getUSteady(){
+      return uSteady_;
+   }
+
+   real_t getUMax(){
+      return uMax_;
+   }
+
+   real_t getT0(){
+      return real_c(t0_);
+   }
+
+   real_t getT1(){
+      return real_c(t1_);
+   }
+
+private:
+   Timeloop & timeloop_;
+   shared_ptr< StructuredBlockForest > blocks_;
+   BlockDataID pdfFieldId_, stressFieldId_;
+   uint_t timesteps_, blockSize_, t0_, t1_, tCurr_;
+   uint_t maxFlag_ = 0;
+   real_t L_, H_, uMax_, uPrev_, uCurr_, uExpected_, uSteady_;
+};
+
+
+//////////
+// Main //
+//////////
+
 int main(int argc, char ** argv ){
 
    walberla::Environment env( argc, argv );
 
-   // create block storage from the config file
+   // read parameter
    shared_ptr<StructuredBlockForest> blocks = blockforest::createUniformBlockGridFromConfig( env.config() );
-
-   // read parameters
    auto parameters = env.config()->getOneBlock( "Parameters" );
 
    // extract some constants from the parameters
-   const std::string paramFile(argv[1]);
    const real_t          eta_s           = parameters.getParameter< real_t >         ( "eta_s",           real_c( 1.4 ) );
-   const real_t          force           = parameters.getParameter< real_t >         ( "force",           real_t( 1 )   );
+   const real_t          force           = parameters.getParameter< real_t >         ( "force",           real_c( 1 )   );
    const real_t          eta_p           = parameters.getParameter< real_t >         ( "eta_p",           real_c( 1.0 ) );
+   const real_t          lambda_p        = parameters.getParameter< real_t >         ( "lambda_p",        real_c( 10  ) );
    const uint_t          period          = parameters.getParameter< uint_t >         ( "period",          uint_c( 1   ) );
-   const uint_t          writePeriod     = parameters.getParameter< uint_t >         ( "writePeriod",     uint_c( 10  ) );
+   const real_t          L               = parameters.getParameter< real_t >         ( "L",               real_c( 10  ) );
+   const real_t          H               = parameters.getParameter< real_t >         ( "H",               real_c( 30  ) );
+   const uint_t          blockSize       = parameters.getParameter< uint_t >         ( "blockSize",       uint_c( 11  ) );
+   const uint_t          timesteps       = parameters.getParameter< uint_t >         ( "timesteps",       uint_c( 10  ) );
 
-   const uint_t checkpointFrequency = parameters.getParameter< uint_t > ("checkpointFrequency", uint_c(5000));
-   const real_t remainingTimeLoggerFrequency = parameters.getParameter< real_t >( "remainingTimeLoggerFrequency", real_c(3.0) ); // in seconds
-
-   // take parameters from command line args preferably
-   char* end;
-   const real_t lambda_p = ((argc > 2)?        atof(argv[2])              : parameters.getParameter<real_t>( "lambda_p", real_c(10)));
-   const std::string baseFolder = ((argc > 2)? std::string(argv[2])       : parameters.getParameter<std::string>( "baseFolder", std::string("vtk_out")));
-   const uint_t timesteps = ((argc > 3)?       strtoul(argv[3],&end,10)   : parameters.getParameter<uint_t>( "timesteps", uint_c(10) ));
+   // reference data
+   const real_t          uExpected = force*H*H/(8.0*(eta_s + eta_p));
+   const real_t          uMax      = parameters.getParameter< real_t > ("uMax");
+   const real_t          t0        = parameters.getParameter< real_t > ("t0");
+   const real_t          t1        = parameters.getParameter< real_t > ("t1");
 
    // create fields
-   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field", uint_c(2));
-   BlockDataID forceFieldId = field::addToStorage<ForceField_T>( blocks, "Force Field", Vector3<real_t>(0.0), field::zyxf, uint_c(2));
+   BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field", FieldGhostLayers);
+   BlockDataID forceFieldId = field::addToStorage<ForceField_T>( blocks, "Force Field", Vector3<real_t>(0.0), field::zyxf, FieldGhostLayers);
    LatticeModel_T latticeModel = LatticeModel_T(lbm::collision_model::TRT::constructWithMagicNumber( walberla::lbm::collision_model::omegaFromViscosity(eta_s)), lbm::force_model::GuoField<ForceField_T>( forceFieldId ) );
-   BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, Vector3<real_t>(), real_c(1), uint_c(2) );
-   BlockDataID stressId = walberla::field::addToStorage<StressField_T>( blocks, "Stress Field", Matrix3<real_t>(0.0), field::zyxf, uint_c(2));
-   BlockDataID stressOldId = walberla::field::addToStorage<StressField_T>( blocks, "Old Stress Field", Matrix3<real_t>(0.0), field::zyxf, uint_c(2));
-   BlockDataID velocityId = walberla::field::addToStorage<VelocityField_T> (blocks, "Velocity Field", Vector3<real_t>(0.0), field::zyxf, uint_c(1));
-
-   // create and initialize boundary handling
-   const FlagUID fluidFlagUID( "Fluid" );
-
-   auto boundariesConfig = env.config()->getOneBlock( "Boundaries" );
-
-   BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage( blocks, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID );
-   geometry::initBoundaryHandling<BHFactory::BoundaryHandling>( *blocks, boundaryHandlingId, boundariesConfig );
-   geometry::setNonBoundaryCellsToDomain<BHFactory::BoundaryHandling> ( *blocks, boundaryHandlingId );
+   BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( blocks, "pdf field", latticeModel, Vector3<real_t>(), real_c(1), FieldGhostLayers );
+   BlockDataID stressId = walberla::field::addToStorage<StressField_T>( blocks, "Stress Field", Matrix3<real_t>(0.0), field::zyxf, FieldGhostLayers);
+   BlockDataID stressOldId = walberla::field::addToStorage<StressField_T>( blocks, "Old Stress Field", Matrix3<real_t>(0.0), field::zyxf, FieldGhostLayers);
+   BlockDataID velocityId = walberla::field::addToStorage<VelocityField_T> (blocks, "Velocity Field", Vector3<real_t>(0.0), field::zyxf, FieldGhostLayers);
 
+   // add boundary handling
+   BlockDataID boundaryHandlingId = blocks->addStructuredBlockData< BoundaryHandling_T >( MyBoundaryHandling( flagFieldId, pdfFieldId ), "boundary handling" );
 
    // create time loop
    SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
 
    // create communication for PdfField
-   blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > communication( blocks );
+   blockforest::communication::UniformBufferedScheme< Stencil_T > communication( blocks );
    communication.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldId ) );
+   auto testData = make_shared< TestData >(TestData(timeloop, blocks, pdfFieldId, stressId, timesteps, blockSize, L, H, uExpected));
 
    timeloop.add() << BeforeFunction( communication, "communication" )
-                  << Sweep( BHFactory::BoundaryHandling::getBlockSweep( boundaryHandlingId ), "boundary handling" );
-   timeloop.add() << BeforeFunction( lbm::viscoelastic::Su<LatticeModel_T, BHFactory::BoundaryHandling>(blocks, forceFieldId, pdfFieldId, boundaryHandlingId, stressId, stressOldId, velocityId,
+                  << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingId ), "boundary handling" );
+   timeloop.add() << BeforeFunction( lbm::viscoelastic::Su<LatticeModel_T, BoundaryHandling_T>(blocks, forceFieldId, pdfFieldId, boundaryHandlingId, stressId, stressOldId, velocityId,
                                                                                                lambda_p, eta_p, period, true), "viscoelasticity")
-                  << Sweep( ConstantForce<BHFactory::BoundaryHandling>(forceFieldId, boundaryHandlingId, force),"Poiseuille Force");
-   timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, fluidFlagUID ) ),
-                            "LB stream & collide" );
-
-   // LBM stability check
-   timeloop.addFuncAfterTimeStep( makeSharedFunctor( field::makeStabilityChecker< ForceField_T, FlagField_T >( env.config(), blocks, forceFieldId,
-                                                                                                             flagFieldId, fluidFlagUID ) ),
-                                  "LBM stability check" );
-
-   // log remaining time
-   timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger( timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency ), "remaining time logger" );
-
-   // Add checkpointing
-   //auto checkpointing = WriteCheckpoint(timeloop, blocks->getBlockStorage(), checkpointFrequency, timesteps, pdfFieldId, stressId, baseFolder, paramFile);
-   //timeloop.addFuncAfterTimeStep(checkpointing);
-   //if(checkpointing.canLoad()){checkpointing.Load();}
+                  << Sweep( ConstantForce<BoundaryHandling_T>(forceFieldId, boundaryHandlingId, force),"Poiseuille Force");
+   timeloop.add() << Sweep( makeSharedSweep( lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldId, flagFieldId, Fluid_Flag ) ),
+                            "LB stream & collide" )
+                  << AfterFunction(makeSharedFunctor(testData), "test data");
 
-   // write midvelocity
-   //auto writeVelocity = WriteVelocity(timeloop, writePeriod, blocks, pdfFieldId);
-   //timeloop.addFuncAfterTimeStep(writeVelocity);
+   timeloop.run();
 
-   // add VTK outputs
-   timeloop.addFuncAfterTimeStep(field::createVTKOutput<StressField_T>(stressId, *blocks, "stress", writePeriod, uint_c(0), false, baseFolder, std::string("step"),
-                                                                       false, true, true, 0, Set<SUID>::emptySet(), Set<SUID>::emptySet(), true, timeloop.getCurrentTimeStep()));
-   timeloop.addFuncAfterTimeStep(field::createVTKOutput<VelocityField_T>(velocityId, *blocks, "velocity", writePeriod, uint_c(0), false, baseFolder, std::string("step"),
-                                                                       false, true, true, 0, Set<SUID>::emptySet(), Set<SUID>::emptySet(), true, timeloop.getCurrentTimeStep()));
 
-   auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", writePeriod, 0, false, baseFolder );
-   flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldId, "FlagField" ) );
-   timeloop.addFuncAfterTimeStep(vtk::writeFiles( flagFieldVTK ), "Flag field VTK");
+   // compare to reference data
+   real_t errSteady = abs(testData->getUSteady() - 1.0)/1.0;
+   real_t errMax = abs(testData->getUMax() - uMax)/uMax;
+   real_t errt0 = abs(testData->getT0() - t0)/t0;
+   real_t errt1 = abs(testData->getT1() - t1)/t1;
 
-   timeloop.run();
+   std::cout << "Relative Errors" << std::endl << std::endl;
+   std::cout << "Steady State Velocity: " << errSteady << std::endl;
+   std::cout << "Maximum Velocity: " << errMax << std::endl;
+   std::cout << "Time of Maximum: " << errt0 << std::endl;
+   std::cout << "Decay Time: " << errt1 << std::endl << std::endl;
 
+   if (errSteady < 0.01 && errMax < 0.01 && errt0 < 0.01 && errt1 < 0.01){
+      std::cout << "Success" << std::endl;
+      return EXIT_SUCCESS;
+   }
+   else {
+      std::cout << "Failure" << std::endl;
+      return EXIT_FAILURE;
+   }
 
-   return EXIT_SUCCESS;
 }
-- 
GitLab