Reply To: Error When Using Checkpoint for Multi-Component
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Error When Using Checkpoint for Multi-Component › Reply To: Error When Using Checkpoint for Multi-Component
November 25, 2024 at 12:41 pm
#9559
yueyq
Participant
Hi Adrian,
Here is the code for checkpoint saving:
#include "olb2D.h"
#include "olb2D.hh"
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<FORCE,EXTERNAL_FORCE,STATISTIC>;
using COUPLING = ShanChenForcedPostProcessor<interaction::PsiEqualsRho>;
// Parameters for the simulation setup
const int nx = 800;
const int ny = 400;
const int maxIter = 40000;
// Stores geometry information in form of material numbers
void prepareGeometry( SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,"prepareGeometry" );
clout << "Prepare Geometry ..." << std::endl;
// Sets material number for fluid and boundary
superGeometry.rename( 0,1 );
Vector<T,2> origin1( -.5, -.5 );
Vector<T,2> origin2( -.5,ny/2. );
Vector<T,2> origin3( -.5, ny-1.5 );
Vector<T,2> extend1( nx+2, 1. );
Vector<T,2> extend2( nx+2., ny/2.+2. );
IndicatorCuboid2D<T> bottom( extend1, origin1 );
IndicatorCuboid2D<T> upper( extend2, origin2 );
IndicatorCuboid2D<T> top( extend1, origin3 );
superGeometry.rename( 1,2,upper );
superGeometry.rename( 1,3,bottom );
superGeometry.rename( 2,4,top );
// Removes all not needed boundary voxels outside the surface
//superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << "Prepare Geometry ... OK" << std::endl;
}
template <typename SuperLatticeCoupling>
void prepareLattice(SuperLattice<T,DESCRIPTOR>& sLatticeOne,
SuperLattice<T,DESCRIPTOR>& sLatticeTwo,
SuperLatticeCoupling& coupling,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout( std::cout,"prepareLattice" );
clout << "Prepare Lattice ..." << std::endl;
const T omega1 = 1.0;
const T omega2 = 1.0;
AnalyticalConst2D<T,T> rho0(0.0);
AnalyticalConst2D<T,T> rho1(1.0);
// The setup is: periodicity along horizontal direction, bounce-back on top
// and bottom. The upper half is initially filled with fluid 1 + random noise,
// and the lower half with fluid 2. Only fluid 1 experiences a forces,
// directed downwards.
using BulkDynamics = ForcedBGKdynamics<T,DESCRIPTOR,momenta::ExternalVelocityTuple>;
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 1);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 2);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 3);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 4);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 1);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 2);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 3);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 4);
setBounceBackBoundary(sLatticeOne, superGeometry, 3);
setBounceBackBoundary(sLatticeOne, superGeometry, 4);
setBounceBackBoundary(sLatticeTwo, superGeometry, 3);
setBounceBackBoundary(sLatticeTwo, superGeometry, 4);
sLatticeOne.setParameter<descriptors::OMEGA>(omega1);
sLatticeTwo.setParameter<descriptors::OMEGA>(omega2);
// A bounce-back node with fictitious density 0, which is experienced by the partner fluid
sLatticeOne.defineRho(superGeometry, 3, rho0);
sLatticeTwo.defineRho(superGeometry, 4, rho0);
// A bounce-back node with fictitious density 0, which is experienced by the partner fluid
sLatticeOne.defineRho(superGeometry, 4, rho1);
sLatticeTwo.defineRho(superGeometry, 3, rho1);
coupling.template setParameter<COUPLING::RHO0>({1, 1});
coupling.template setParameter<COUPLING::G>(T(3.));
coupling.template setParameter<COUPLING::OMEGA_A>(omega1);
coupling.template setParameter<COUPLING::OMEGA_B>(omega2);
sLatticeOne.addPostProcessor<stage::PreCoupling>(meta::id<RhoStatistics>());
sLatticeTwo.addPostProcessor<stage::PreCoupling>(meta::id<RhoStatistics>());
{
auto& communicator = sLatticeOne.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(1);
communicator.requestField<STATISTIC>();
communicator.exchangeRequests();
}
{
auto& communicator = sLatticeTwo.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(1);
communicator.requestField<STATISTIC>();
communicator.exchangeRequests();
}
clout << "Prepare Lattice ... OK" << std::endl;
}
void setBoundaryValues( SuperLattice<T, DESCRIPTOR>& sLatticeOne,
SuperLattice<T, DESCRIPTOR>& sLatticeTwo,
T force, int iT, SuperGeometry<T,2>& superGeometry )
{
if ( iT==0 ) {
AnalyticalConst2D<T,T> noise( 4.e-2 );
std::vector<T> v( 2,T() );
AnalyticalConst2D<T,T> zeroV( v );
AnalyticalConst2D<T,T> zero( 1.e-6 );
AnalyticalLinear2D<T,T> one( 0.,-force*invCs2<T,DESCRIPTOR>(),0.98+force*ny*invCs2<T,DESCRIPTOR>() );
AnalyticalConst2D<T,T> onePlus( 0.98+force*ny/2.*invCs2<T,DESCRIPTOR>() );
AnalyticalRandom2D<T,T> random;
AnalyticalIdentity2D<T,T> randomOne( random*noise+one );
AnalyticalIdentity2D<T,T> randomPlus( random*noise+onePlus );
std::vector<T> F( 2,T() );
F[1] = -force;
AnalyticalConst2D<T,T> f( F );
// for each material set the defineRhou and the Equilibrium
sLatticeOne.defineRhoU( superGeometry, 1, zero, zeroV );
sLatticeOne.iniEquilibrium( superGeometry, 1, zero, zeroV );
sLatticeOne.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 1, f );
sLatticeTwo.defineRhoU( superGeometry, 1, randomPlus, zeroV );
sLatticeTwo.iniEquilibrium( superGeometry, 1, randomPlus, zeroV );
sLatticeOne.defineRhoU( superGeometry, 2, randomOne, zeroV );
sLatticeOne.iniEquilibrium( superGeometry, 2, randomOne, zeroV );
sLatticeOne.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 2, f );
sLatticeTwo.defineRhoU( superGeometry, 2, zero, zeroV );
sLatticeTwo.iniEquilibrium( superGeometry, 2, zero, zeroV );
// Make the lattice ready for simulation
sLatticeOne.initialize();
sLatticeTwo.initialize();
}
}
void getResults( SuperLattice<T, DESCRIPTOR>& sLatticeTwo,
SuperLattice<T, DESCRIPTOR>& sLatticeOne, int iT,
SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer )
{
OstreamManager clout( std::cout,"getResults" );
SuperVTMwriter2D<T> vtmWriter( "rayleighTaylor2dsLatticeOne" );
const int vtkIter = 1000;
const int statIter = 1000;
if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLatticeOne, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLatticeOne );
SuperLatticeRank2D<T, DESCRIPTOR> rank( sLatticeOne );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}
// Get statistics
if ( iT%statIter==0 && iT > 0 ) {
// Timer console output
timer.update( iT );
timer.printStep();
clout << "averageRhoFluidOne=" << sLatticeOne.getStatistics().getAverageRho();
clout << "; averageRhoFluidTwo=" << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
}
// Writes the VTK files
if ( iT % vtkIter==0 ) {
sLatticeOne.setProcessingContext(ProcessingContext::Evaluation);
sLatticeTwo.setProcessingContext(ProcessingContext::Evaluation);
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity( sLatticeOne );
SuperLatticeDensity2D<T, DESCRIPTOR> density( sLatticeOne );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( density );
vtmWriter.write( iT );
BlockReduction2D2D<T> planeReduction( density, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
sLatticeOne.save( "rayleighTaylorOne" );
sLatticeTwo.save( "rayleighTaylorTwo" );
}
}
int main( int argc, char *argv[] )
{
// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( "./tmp/" );
OstreamManager clout( std::cout,"main" );
T force = 30./( T )ny/( T )ny;
// === 2nd Step: Prepare Geometry ===
// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cGeometry( 0, 0, 1, nx, ny, singleton::mpi().getSize() );
#else
CuboidGeometry2D<T> cGeometry( 0, 0, 1, nx, ny, 1 );
#endif
cGeometry.setPeriodicity( true, false );
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
SuperGeometry<T,2> superGeometry( cGeometry, loadBalancer, 2 );
prepareGeometry( superGeometry );
// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLatticeOne( superGeometry );
SuperLattice<T, DESCRIPTOR> sLatticeTwo( superGeometry );
SuperLatticeCoupling coupling(
ShanChenForcedPostProcessor<interaction::PsiEqualsRho>{},
names::A{}, sLatticeOne,
names::B{}, sLatticeTwo);
prepareLattice( sLatticeOne, sLatticeTwo, coupling, superGeometry );
// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << "starting simulation..." << std::endl;
util::Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();
//sLatticeOne.load( "rayleighTaylorOne" );
//sLatticeTwo.load( "rayleighTaylorTwo" );
for ( iT=0; iT<maxIter; ++iT ) {
// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( sLatticeOne, sLatticeTwo, force, iT, superGeometry );
// === 6th Step: Collide and Stream Execution ===
sLatticeOne.collideAndStream();
sLatticeTwo.collideAndStream();
sLatticeOne.executePostProcessors(stage::PreCoupling());
sLatticeTwo.executePostProcessors(stage::PreCoupling());
coupling.execute();
// === 7th Step: Computation and Output of the Results ===
getResults( sLatticeTwo, sLatticeOne, iT, superGeometry, timer );
}
timer.stop();
timer.printSummary();
}
Then, when it shows [Timer] step=3000; percent=7.5; passedTime=87.541; remTime=1079.67; MLUPs=10.9901 on the terminal, I run:
#include "olb2D.h"
#include "olb2D.hh"
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<FORCE,EXTERNAL_FORCE,STATISTIC>;
using COUPLING = ShanChenForcedPostProcessor<interaction::PsiEqualsRho>;
// Parameters for the simulation setup
const int nx = 800;
const int ny = 400;
const int maxIter = 40000;
// Stores geometry information in form of material numbers
void prepareGeometry( SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,"prepareGeometry" );
clout << "Prepare Geometry ..." << std::endl;
// Sets material number for fluid and boundary
superGeometry.rename( 0,1 );
Vector<T,2> origin1( -.5, -.5 );
Vector<T,2> origin2( -.5,ny/2. );
Vector<T,2> origin3( -.5, ny-1.5 );
Vector<T,2> extend1( nx+2, 1. );
Vector<T,2> extend2( nx+2., ny/2.+2. );
IndicatorCuboid2D<T> bottom( extend1, origin1 );
IndicatorCuboid2D<T> upper( extend2, origin2 );
IndicatorCuboid2D<T> top( extend1, origin3 );
superGeometry.rename( 1,2,upper );
superGeometry.rename( 1,3,bottom );
superGeometry.rename( 2,4,top );
// Removes all not needed boundary voxels outside the surface
//superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << "Prepare Geometry ... OK" << std::endl;
}
template <typename SuperLatticeCoupling>
void prepareLattice(SuperLattice<T,DESCRIPTOR>& sLatticeOne,
SuperLattice<T,DESCRIPTOR>& sLatticeTwo,
SuperLatticeCoupling& coupling,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout( std::cout,"prepareLattice" );
clout << "Prepare Lattice ..." << std::endl;
const T omega1 = 1.0;
const T omega2 = 1.0;
AnalyticalConst2D<T,T> rho0(0.0);
AnalyticalConst2D<T,T> rho1(1.0);
// The setup is: periodicity along horizontal direction, bounce-back on top
// and bottom. The upper half is initially filled with fluid 1 + random noise,
// and the lower half with fluid 2. Only fluid 1 experiences a forces,
// directed downwards.
using BulkDynamics = ForcedBGKdynamics<T,DESCRIPTOR,momenta::ExternalVelocityTuple>;
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 1);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 2);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 3);
sLatticeOne.defineDynamics<BulkDynamics>(superGeometry, 4);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 1);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 2);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 3);
sLatticeTwo.defineDynamics<BulkDynamics>(superGeometry, 4);
setBounceBackBoundary(sLatticeOne, superGeometry, 3);
setBounceBackBoundary(sLatticeOne, superGeometry, 4);
setBounceBackBoundary(sLatticeTwo, superGeometry, 3);
setBounceBackBoundary(sLatticeTwo, superGeometry, 4);
sLatticeOne.setParameter<descriptors::OMEGA>(omega1);
sLatticeTwo.setParameter<descriptors::OMEGA>(omega2);
// A bounce-back node with fictitious density 0, which is experienced by the partner fluid
sLatticeOne.defineRho(superGeometry, 3, rho0);
sLatticeTwo.defineRho(superGeometry, 4, rho0);
// A bounce-back node with fictitious density 0, which is experienced by the partner fluid
sLatticeOne.defineRho(superGeometry, 4, rho1);
sLatticeTwo.defineRho(superGeometry, 3, rho1);
coupling.template setParameter<COUPLING::RHO0>({1, 1});
coupling.template setParameter<COUPLING::G>(T(3.));
coupling.template setParameter<COUPLING::OMEGA_A>(omega1);
coupling.template setParameter<COUPLING::OMEGA_B>(omega2);
sLatticeOne.addPostProcessor<stage::PreCoupling>(meta::id<RhoStatistics>());
sLatticeTwo.addPostProcessor<stage::PreCoupling>(meta::id<RhoStatistics>());
{
auto& communicator = sLatticeOne.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(1);
communicator.requestField<STATISTIC>();
communicator.exchangeRequests();
}
{
auto& communicator = sLatticeTwo.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(1);
communicator.requestField<STATISTIC>();
communicator.exchangeRequests();
}
clout << "Prepare Lattice ... OK" << std::endl;
}
void setBoundaryValues( SuperLattice<T, DESCRIPTOR>& sLatticeOne,
SuperLattice<T, DESCRIPTOR>& sLatticeTwo,
T force, int iT, SuperGeometry<T,2>& superGeometry )
{
if ( iT==0 ) {
AnalyticalConst2D<T,T> noise( 4.e-2 );
std::vector<T> v( 2,T() );
AnalyticalConst2D<T,T> zeroV( v );
AnalyticalConst2D<T,T> zero( 1.e-6 );
AnalyticalLinear2D<T,T> one( 0.,-force*invCs2<T,DESCRIPTOR>(),0.98+force*ny*invCs2<T,DESCRIPTOR>() );
AnalyticalConst2D<T,T> onePlus( 0.98+force*ny/2.*invCs2<T,DESCRIPTOR>() );
AnalyticalRandom2D<T,T> random;
AnalyticalIdentity2D<T,T> randomOne( random*noise+one );
AnalyticalIdentity2D<T,T> randomPlus( random*noise+onePlus );
std::vector<T> F( 2,T() );
F[1] = -force;
AnalyticalConst2D<T,T> f( F );
// for each material set the defineRhou and the Equilibrium
sLatticeOne.defineRhoU( superGeometry, 1, zero, zeroV );
sLatticeOne.iniEquilibrium( superGeometry, 1, zero, zeroV );
sLatticeOne.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 1, f );
sLatticeTwo.defineRhoU( superGeometry, 1, randomPlus, zeroV );
sLatticeTwo.iniEquilibrium( superGeometry, 1, randomPlus, zeroV );
sLatticeOne.defineRhoU( superGeometry, 2, randomOne, zeroV );
sLatticeOne.iniEquilibrium( superGeometry, 2, randomOne, zeroV );
sLatticeOne.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 2, f );
sLatticeTwo.defineRhoU( superGeometry, 2, zero, zeroV );
sLatticeTwo.iniEquilibrium( superGeometry, 2, zero, zeroV );
// Make the lattice ready for simulation
sLatticeOne.initialize();
sLatticeTwo.initialize();
}
}
void getResults( SuperLattice<T, DESCRIPTOR>& sLatticeTwo,
SuperLattice<T, DESCRIPTOR>& sLatticeOne, int iT,
SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer )
{
OstreamManager clout( std::cout,"getResults" );
SuperVTMwriter2D<T> vtmWriter( "rayleighTaylor2dsLatticeOne" );
const int vtkIter = 1000;
const int statIter = 1000;
if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLatticeOne, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLatticeOne );
SuperLatticeRank2D<T, DESCRIPTOR> rank( sLatticeOne );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}
// Get statistics
if ( iT%statIter==0 && iT > 0 ) {
// Timer console output
timer.update( iT );
timer.printStep();
clout << "averageRhoFluidOne=" << sLatticeOne.getStatistics().getAverageRho();
clout << "; averageRhoFluidTwo=" << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
}
// Writes the VTK files
if ( iT % vtkIter==0 ) {
sLatticeOne.setProcessingContext(ProcessingContext::Evaluation);
sLatticeTwo.setProcessingContext(ProcessingContext::Evaluation);
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity( sLatticeOne );
SuperLatticeDensity2D<T, DESCRIPTOR> density( sLatticeOne );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( density );
vtmWriter.write( iT );
BlockReduction2D2D<T> planeReduction( density, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
sLatticeOne.save( "rayleighTaylorOne" );
sLatticeTwo.save( "rayleighTaylorTwo" );
}
}
int main( int argc, char *argv[] )
{
// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( "./tmp/" );
OstreamManager clout( std::cout,"main" );
T force = 30./( T )ny/( T )ny;
// === 2nd Step: Prepare Geometry ===
// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cGeometry( 0, 0, 1, nx, ny, singleton::mpi().getSize() );
#else
CuboidGeometry2D<T> cGeometry( 0, 0, 1, nx, ny, 1 );
#endif
cGeometry.setPeriodicity( true, false );
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
SuperGeometry<T,2> superGeometry( cGeometry, loadBalancer, 2 );
prepareGeometry( superGeometry );
// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLatticeOne( superGeometry );
SuperLattice<T, DESCRIPTOR> sLatticeTwo( superGeometry );
SuperLatticeCoupling coupling(
ShanChenForcedPostProcessor<interaction::PsiEqualsRho>{},
names::A{}, sLatticeOne,
names::B{}, sLatticeTwo);
prepareLattice( sLatticeOne, sLatticeTwo, coupling, superGeometry );
// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << "starting simulation..." << std::endl;
util::Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();
sLatticeOne.load( "rayleighTaylorOne" );
sLatticeTwo.load( "rayleighTaylorTwo" );
for ( iT=3000; iT<maxIter; ++iT ) {
// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( sLatticeOne, sLatticeTwo, force, iT, superGeometry );
// === 6th Step: Collide and Stream Execution ===
sLatticeOne.collideAndStream();
sLatticeTwo.collideAndStream();
sLatticeOne.executePostProcessors(stage::PreCoupling());
sLatticeTwo.executePostProcessors(stage::PreCoupling());
coupling.execute();
// === 7th Step: Computation and Output of the Results ===
getResults( sLatticeTwo, sLatticeOne, iT, superGeometry, timer );
}
timer.stop();
timer.printSummary();
}
Also, I wonder how to use the checkpoint for resolved-particle cases.
Thank you!
Best,
yueyq