[Bobbie]

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

[Bobbie]

The following complete simulation code combines the examples of cylinder2d and microFluidcs2d. However, despite explicitly setting the cylinder’s material identifier to 0 to designate it as a solid domain, the post-processing results demonstrate unexpected fluid-like behavior within the cylindrical region. Could you please advise on necessary code modifications to properly enforce solid phase characteristics?

#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();

}

[Bobbie]

I have verified that the material number is correctly assigned. Despite attempts to set the material number of the removed geometric model to 0 or treat it as a solid region, the results still indicate that it exhibits fluid-phase characteristics.

[Bobbie]

Yes, I’m truly grateful for your invaluable assistance, as it has provided me with valuable insights on how to enhance the efficiency of my calculations.

[Bobbie]

Yeah, Mathias. I found that information in Section 10.9 (Lesson 9: Run your Programs on a Parallel Machine) of the User Guide. After adhering to the instructions to modify the config.mk file, it seems that parallel computation across multiple CPU cores has been successfully implemented without encountering any errors or warnings. Thank you once again for your assistance.

[Bobbie]

Hello,Adrian. Thank you for your answer. Given that my current computer is equipped with multiple CPU cores and threads, I’d like to leverage CPU parallel computation to accelerate the simulation process. May I inquire if this approach aligns with the MPI (Message Passing Interface) you mentioned, as I’m not very familiar with these concepts?

[Bobbie]

Although the code is running, the calculation process is very long. Is this related to the change of value of noofCuboids?

[Bobbie]

Hello，Adrian. This code is used to simulate the multi-phase flow inside a porous media, so I combined the codes of porous media (like resolvedRock3d) with the example of multiphase flow (FreeEnergy multi-phase model). And then the above issue occurred at running time.
But when I changed the noofCuboids to 1, the code seemed to work, but the terminal still displayed the above issue(no discreteNormal is found).

[Bobbie]

Thank you for your attention. I have reviewed the geometrical model in Paraview, and the details are clearly visible. However, I am unsure which step is causing the above issue.

[Bobbie]

Dear Jan,
Thank you for sharing your method with great respect! I’ll give it a try and I hope it will succeed.
Best wishes,
Bobbie

[Bobbie]

Dear Nipin(Nipinl)
You are quite right and the case you mentioned seems incomplete. I have also appeared problems when I tried to make and run it. I’m also trying to improve it, but I haven’t made much progress yet. Perhaps the answer Jan gave is feasible, but I haven’t had a chance to try it. You can try it. And do you have any new findings that I hope we can continue to communicate.
Best regards,
Bobbie

[Bobbie]

Thank you for your help. Best wishes!

[Bobbie]

Thank you for your help again!

[Bobbie]

Thank you for your answers, Adrian. The geometric model was generated by SW and exported to stl format. The size of the model is 0.15 mm * 0.15 mm * 1 mm. So I followed the example and Users Guide to set the stl unit to 0.001 when importing the stl file,such as follows:

STLreader<T> stlReader( “multicylinder3d.stl”, converter.getConversionFactorLength(), 0.001 );
IndicatorLayer3D<T> extendedDomain( stlReader, converter.getConversionFactorLength() );

Should I adjust the value of 0.001 in above codes or make changes in other locations to prevent the occurrence of the aforementioned issues?

Thank you for your help again!
Best wishes!

• This reply was modified 10 months, 2 weeks ago by Bobbie.

[Bobbie]

OK, thank you for your help. Best wishes!

[Author]

Posts