diff --git a/tests/lbm/Su.prm b/tests/lbm/Su.prm
index 470b5dea2ace6b7e7bd5c5ec15f76110dcf3b839..6c215a2cb4f92646f3e2030abe09d398246fbea4 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 cd63004be55b896fe696ab182503c0a2bb557783..cecac2e28b68c12455de2e398a8ed8d93b0c85f9 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;
}