Reply To: Shan-Chen
Dear Luiz,
In our current investigation of multiphase flow systems, a critical issue persists regarding solid material numbers. The following code integrates the example of cylinder2d and microFluidcs2d. While the cylindrical region has been specifically assigned material number to zero through phase tagging to represent solid boundaries, subsequent visualization through ParaView reveals liquid density(rho) inside cylindrical region.The results is wrong obviously.So, What fundamental adjustments would ensure the correct running of this code.
The codes are as follows:
#include “olb2D.h”
#include “olb2D.hh”
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using T = FLOATING_POINT_TYPE;
typedef D2Q9<CHEM_POTENTIAL,FORCE> DESCRIPTOR;
const int N = 100;
const T nx = 800e-6;
const T ny = 100e-6;
const T dx = ny/N;
const T alpha = 1.;
const T kappa1 = 0.005;
const T kappa2 = 0.005;
const T gama = 10.;
const T h1 = 0.0001448;
const T h2 = -0.0001448;
const T inletVelocity = 0.01;
const T outletDensity = 1.;
const int maxIter = 100000;
const int vtkIter = 200;
const int statIter = 200;
void prepareGeometry( SuperGeometry<T,2>& superGeometry,
std::shared_ptr<IndicatorF2D<T>> section1
)
{
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;
superGeometry.rename(0,2,section1);
superGeometry.rename(2,1,{1,1});
IndicatorCuboid2D<T> inlet(dx,ny,{0.0,ny/T(2)});
IndicatorCuboid2D<T> outlet(dx,ny,{nx,ny/T(2)});
superGeometry.rename( 2,3,1,inlet );
superGeometry.rename( 2,4,1,outlet );
IndicatorCircle2D<T> circle({100.0*dx,ny/T(2)}, 12.5*dx);
superGeometry.rename( 1,0,circle );
superGeometry.clean();
//superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << “Prepare Geometry … OK” << std::endl;
}
void prepareLattice( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2,
UnitConverter<T, DESCRIPTOR>& converter,
SuperGeometry<T,2>& superGeometry)
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;
T omega = converter.getLatticeRelaxationFrequency();
clout << “Prepare Lattice: Define lattice dynamics …” << std::endl;
sLattice1.defineDynamics<ForcedBGKdynamics>(superGeometry, 1);
sLattice2.defineDynamics<FreeEnergyBGKdynamics>(superGeometry, 1);
clout << “Prepare Lattice: Add wall boundary …” << std::endl;
setFreeEnergyWallBoundary(sLattice1, superGeometry, 2, alpha, kappa1, kappa2, h1, h2, 1);
setFreeEnergyWallBoundary(sLattice2, superGeometry, 2, alpha, kappa1, kappa2, h1, h2, 2);
clout << “Prepare Lattice: Add inlet boundaries …” << std::endl;
auto inletIndicator = superGeometry.getMaterialIndicator(3);
setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice1, omega, inletIndicator, “velocity”, 1);
setFreeEnergyInletBoundary<T,DESCRIPTOR>(sLattice2, omega, inletIndicator, “velocity”, 2);
clout << “Prepare Lattice: Add outlet boundary …” << std::endl;
auto outletIndicator = superGeometry.getMaterialIndicator(4);
setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice1, omega, outletIndicator, “density”,1);
setFreeEnergyOutletBoundary<T,DESCRIPTOR>(sLattice2, omega, outletIndicator, “density”,2);
clout << “Prepare Lattice: Bulk initial conditions …” << std::endl;
std::vector<T> v( 2,T() );
AnalyticalConst2D<T,T> zeroVelocity( v );
AnalyticalConst2D<T,T> one( 1.0 );
AnalyticalConst2D<T,T> zero( 0.0 );
IndicatorCuboid2D<T> ind1(nx-dx,ny,{dx+(nx-dx)/T(2),ny/T(2)});
SmoothIndicatorCuboid2D<T,T> section1( ind1, 0 );
AnalyticalIdentity2D<T,T> c1( section1 );
AnalyticalIdentity2D<T,T> rho( one ); // c2 + c1
AnalyticalIdentity2D<T,T> phi( one – c1 -c1 );// c2 – c1
auto allIndicator = superGeometry.getMaterialIndicator({1, 2, 3, 4});
sLattice1.iniEquilibrium( allIndicator, rho, zeroVelocity );
sLattice2.iniEquilibrium( allIndicator, phi, zeroVelocity );
clout << “Prepare Lattice: Inlet boundary conditions …” << std::endl;
AnalyticalConst2D<T,T> inletU(inletVelocity);
sLattice1.defineU( inletIndicator, inletU );
sLattice2.defineRho( inletIndicator, phi );
clout << “Prepare Lattice: Outlet initial / Boundary conditions …” << std::endl;
AnalyticalConst2D<T,T> rhoOutlet( outletDensity );
AnalyticalIdentity2D<T,T> phiOutlet( zero );
sLattice1.defineRho( outletIndicator, rhoOutlet );
sLattice2.defineRho( outletIndicator, phiOutlet );
sLattice1.setParameter<descriptors::OMEGA>( omega );
sLattice2.setParameter<descriptors::OMEGA>( omega );
sLattice2.setParameter<collision::FreeEnergy::GAMMA>( gama );
clout << “Prepare Lattice: Initialise lattices …” << std::endl;sLattice1.initialize();
sLattice1.initialize();
sLattice2.initialize();
clout << “Prepare Lattice: Communicate …” << std::endl;
sLattice1.communicate();
sLattice2.communicate();
clout << “Prepare Lattice … OK” << std::endl;
}
void prepareCoupling( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2,
SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,”prepareCoupling” );
clout << “Add lattice coupling” << std::endl;
// 主体耦合
FreeEnergyChemicalPotentialGenerator2D<T,DESCRIPTOR> coupling2( alpha, kappa1, kappa2 ); //化学势耦合→coupling2
FreeEnergyForceGenerator2D<T,DESCRIPTOR> coupling3; //外力耦合→coupling3
// 进出口耦合
FreeEnergyDensityOutletGenerator2D<T,DESCRIPTOR> coupling1( outletDensity ); //密度出口耦合→coupling1
FreeEnergyInletOutletGenerator2D<T,DESCRIPTOR> coupling4; //进出口耦合→coupling4
sLattice1.addLatticeCoupling( superGeometry, 4, coupling1, sLattice2 );//密度出口耦合
sLattice1.addLatticeCoupling( superGeometry, 1, coupling2, sLattice2 );//化学势耦合
sLattice2.addLatticeCoupling( superGeometry, 1, coupling3, sLattice1 );//外力耦合
sLattice2.addLatticeCoupling( superGeometry, 3, coupling4, sLattice1 );//进出口耦合
sLattice2.addLatticeCoupling( superGeometry, 4, coupling4, sLattice1 );//进出口耦合
{
auto& communicator = sLattice1.getCommunicator(stage::PostCoupling());
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestOverlap(sLattice1.getOverlap());
communicator.exchangeRequests();
}
{
auto& communicator = sLattice2.getCommunicator(stage::PreCoupling());
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestOverlap(sLattice2.getOverlap());
communicator.exchangeRequests();
}
clout << “Add lattice coupling … OK!” << std::endl;
}
void getResults( SuperLattice<T, DESCRIPTOR>& sLattice1,
SuperLattice<T, DESCRIPTOR>& sLattice2, int iT,
SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer,
UnitConverter<T, DESCRIPTOR> converter )
{
OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter2D<T> vtmWriter( “FreeEnergy” );
if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice1, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice1 );
SuperLatticeRank2D<T, DESCRIPTOR> rank( sLattice1 );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}
// Get statistics
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();
sLattice1.getStatistics().print( iT, converter.getPhysTime(iT) );
sLattice2.getStatistics().print( iT, converter.getPhysTime(iT) );
clout << “averageRhoFluidOne=” << sLattice1.getStatistics().getAverageRho();
clout << “; averageRhoFluidTwo=” << sLattice2.getStatistics().getAverageRho() << std::endl;
}
// Writes the VTK files
if ( iT%vtkIter==0 ) {
sLattice1.setProcessingContext(ProcessingContext::Evaluation);
sLattice2.setProcessingContext(ProcessingContext::Evaluation);
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity( sLattice1 );
SuperLatticeDensity2D<T, DESCRIPTOR> rho( sLattice1 );
rho.getName() = “rho”;
SuperLatticeDensity2D<T, DESCRIPTOR> phi( sLattice2 );
phi.getName() = “phi”;
AnalyticalConst2D<T,T> half_( 0.5 );
SuperLatticeFfromAnalyticalF2D<T, DESCRIPTOR> half(half_, sLattice1);
SuperIdentity2D<T,T> c1 (half*(rho+phi));
c1.getName() = “density-fluid-1”;
SuperIdentity2D<T,T> c2 (half*(rho-phi));
c2.getName() = “density-fluid-2”;
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( rho );
vtmWriter.addFunctor( phi );
vtmWriter.addFunctor( c1 );
vtmWriter.addFunctor( c2 );
vtmWriter.write( iT );
// BlockReduction2D2D interpolates the data of a SuperF2D functor in a given resolution.
// This is primarily used for exporting GIF images via BlockGifWriter.
// BlockReduction2D2D在给定的分辨率内插一个SuperF2D函子的数据。这主要用于通过BlockGifWriter导出GIF图像。
// olb::BlockReduction2D2D<T>::BlockReduction2D2D(FunctorPtr<SuperF2D<T>>&&f,int resolution = 600,
// BlockDataSyncMode mode=BlockDataSyncMode::ReduceAndBcast)
// BlockReduction2D2D所需参数：FunctorPtr<SuperF2D<T>>&&f, resolution, ReduceAndBcast/ReduceOnly/None
BlockReduction2D2D<T> planeReduction( c1, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
}
}
int main( int argc, char *argv[] )
{
// === 第一步：初始化 ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );
UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> converter(
(T) N, // resolution
(T) 1, // lattice relaxation time (tau)
(T) ny, // charPhysLength: reference length of simulation geometry
(T) 0.1, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 1e-5, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1000. // physDensity: physical density in __kg / m^3__
);
// Prints the converter log as console output
converter.print();
converter.write(“FreeEnergy”);
// === 第二步：准备几何模型 ===
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 1;
#endif
// setup section1
Vector<T,2> extendChannel( nx, ny );
Vector<T,2> originChannel( 0, 0 );
std::shared_ptr<IndicatorF2D<T>> section1 = std::make_shared<IndicatorCuboid2D<T>>( extendChannel, originChannel );
CuboidGeometry2D<T> cGeometry( *(section1), converter.getConversionFactorLength(), noOfCuboids );
// 设置模型周期性边界
//cGeometry.setPeriodicity( false, true );
// Instantiation of loadbalancer
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
loadBalancer.print();
// Instantiation of superGeometry
SuperGeometry<T,2> superGeometry( cGeometry,loadBalancer );
prepareGeometry( superGeometry, section1 );
// === 第三步：准备格子 ===
SuperLattice<T, DESCRIPTOR> sLattice1( superGeometry );
SuperLattice<T, DESCRIPTOR> sLattice2( superGeometry );
// 准备格子和边界条件
prepareLattice( sLattice1, sLattice2, converter, superGeometry);
prepareCoupling( sLattice1, sLattice2, superGeometry);
// === 第四步：主循环 ===
int iT = 0;
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();
for ( iT=0; iT<=maxIter; ++iT ) {
// === 第五步：定义初始化条件及边界条件 ===
//setBoundaryValues( sLatticeOne, sLatticeTwo, force, iT, superGeometry );
// === 第六步：执行碰撞和迁移步骤 ===
sLattice1.collideAndStream();
sLattice2.collideAndStream();
// 执行格子间耦合
sLattice1.executeCoupling();
sLattice2.executeCoupling();
// === 第七步：计算并输出结果 ===
getResults( sLattice1, sLattice2, iT, superGeometry, timer, converter );
}
timer.stop();
timer.printSummary();
}
Thank you for your sincere help.
Best regards,
Bobbie