[Fany]

Hello Adrian,

Many thanks for your help. The code works now.

Have a nice day!

[Fany]

Hello Adrian,

I am sorry for the wrong version of code. The full code (mpic++) is seen below:

#include “olb2D.h”
#include “olb2D.hh”

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace std;

using T = FLOATING_POINT_TYPE;

using NSDESCRIPTOR = D2Q9<FORCE>;
using BulkDynamics = BGKdynamics<T,NSDESCRIPTOR>;
using ForcedBulkDynamics = ForcedBGKdynamics<T,NSDESCRIPTOR>;

using TDESCRIPTOR = D2Q5<VELOCITY>;

//simulation parameters
const T lx1 = 0.003 ; // length of the channel
const T ly1 = 0.0008; // width of the channel
const T lx2 = 0.002 ; // length of the channel
const T ly2 = 0.0006 ; // width of the channel
const T theta1 = 0.0; // bifurcation angle at 1 level
const T theta2 = std::acos(ly1/2/ly2); // bifurcation angle at 2 level

const int N = 35; // resolution of the model
const int M = 10; // time discretization refinement
// const bool bouzidiOn = false; // choice of boundary condition
const T maxPhysT = 0.8; // max. simulation time in s, SI unit
const T Re = 30.; // Reynolds number
const T Pr = 5.4; // Prandtl number
const T epsilon = 1.e-5; // precision of the convergence (residuum)
const T L= 0.0006 ; //charPhysLength

const T Thot = 313.15; // temperature of the lower wall in Kelvin
const T Tcold = 293.15; // temperature of the fluid in Kelvin
const T physU = Re*8.05e-7/L; // physical velocity ;Re*8.05e-7/0.02246

// Stores data from stl file in geometry in form of material numbers
void prepareGeometry( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter, SuperGeometry<T,2>& superGeometry, auto & indicator1)
// , IndicatorF2D<T>& indicator2
{

OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename( 0,2, indicator1); //
superGeometry.rename( 2, 1,{1,1} );
superGeometry.clean();

// creating the inflow
std::vector<T> inlet(2,T());
inlet[0] = -lx1/2+0.5*converter.getPhysLength(1);
inlet[1] = 0.0;
IndicatorCuboid2D<T> inflow(converter.getPhysLength(1), ly1, inlet, theta1);

// creating the outflow0
std::vector<T> outlet0(2,T());
outlet0[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2)-converter.getPhysLength(1);
outlet0[1] = (lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2);
std::vector<T> extend1( 2, T(0) );
extend1[0] = converter.getPhysLength(1);
extend1[1] = ly2/std::cos(theta2);
IndicatorCuboid2D<T> outflow0(extend1, outlet0, theta1);
// auto outflow0=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly2, outlet0, theta2) ;

// creating the outflow1
std::vector<T> outlet1(2,T());
outlet1[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2)-converter.getPhysLength(1);
outlet1[1] = -(lx2)*std::sin(theta2)+ly2/std::cos(theta2)*std::pow(std::sin(theta2),2)-ly2/std::cos(theta2);
std::vector<T> extend2( 2, T(0) );
extend2[0] = converter.getPhysLength(1);
extend2[1] = ly2/std::cos(theta2);
IndicatorCuboid2D<T> outflow1(extend2, outlet1, theta1);
// auto outflow1=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly2, outlet1, theta2) ;

superGeometry.rename( 2,3,1, inflow );
superGeometry.rename( 2,4,1, outflow0 );
superGeometry.rename( 2,5,1, outflow1 );

// 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;

}

// Set up the geometry of the simulation
void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

// // material=0 –> do nothing
// ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

// material=1 –> bulk dynamics
//lattice.defineDynamics( superGeometry,1,&bulkDynamics );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
NSlattice.defineDynamics<BulkDynamics>(superGeometry, 1);

// material=2 –> bounceBack dynamics
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
setBounceBackBoundary(NSlattice, superGeometry, 2);
// setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 2);

// material=3 –> bulk dynamics + velocity (inflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 3 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 3);

// material=4,5 –> bulk dynamics + pressure (outflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 4 );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 5 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 4 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 5 );

// Boundary condition
setAdvectionDiffusionTemperatureBoundary(ADlattice, Tomega, superGeometry, 2);

setInterpolatedVelocityBoundary(NSlattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary(NSlattice, omega, superGeometry, 5);
setInterpolatedPressureBoundary(NSlattice, omega, superGeometry, 4);

// AnalyticalConst2D<T,T> rhoF( 1 );
// std::vector<T> velocity( 3,T() );
// AnalyticalConst2D<T,T> uF( velocity );
AnalyticalConst2D<T,T> rhoF(1.);
AnalyticalConst2D<T,T> uF( converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst2D<T,T> u0(0.0, 0.0);
AnalyticalConst2D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst2D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, u0);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, u0);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
NSlattice.setParameter<descriptors::OMEGA>(omega);

// Lattice initialize
NSlattice.initialize();
ADlattice.initialize();

clout << “Prepare Lattice … OK” << std::endl;
}

void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,2>& superGeometry)
{
// / initial condition
OstreamManager clout( std::cout,”setBoundaryValues” );
// No of time steps for smooth start-up
int iTmaxStart = converter.getLatticeTime( maxPhysT*0.25 );
int iTupdate = 30;

if (iT%iTupdate == 0 && iT <= iTmaxStart) {

// T u_Re = converter.getLatticeVelocity( Re * converter.getPhysViscosity() / converter.getCharPhysLength() );
// Smooth start curve, sinus
SinusStartScale<T,int> nSinusStartScale(iTmaxStart, converter.getCharLatticeVelocity());

// Smooth start curve, polynomial
int iTvec[1] = { iT };
T frac[1] = {T()};
nSinusStartScale( frac,iTvec );

T maxVelocity = converter.getCharLatticeVelocity()*3./2.*frac[0];
T distance2Wall = converter.getConversionFactorLength()/2.;
Poiseuille2D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall );

NSlattice.defineU( superGeometry, 3, poiseuilleU );
NSlattice.setProcessingContext< Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
ProcessingContext::Simulation);
}

// if (iT==0){

// std::vector<T> velocity( 2,T() );
// velocity[0] = converter.getCharLatticeVelocity();
// AnalyticalConst2D<T,T> uF( velocity );
// NSlattice.defineU( superGeometry, 3, uF );
// NSlattice.setProcessingContext< Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
// }
}

// Output to console and files
void getResults( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,2>& superGeometry, Timer<T>& timer, bool converged)
{

OstreamManager clout( std::cout,”getResults” );

SuperVTMwriter2D<T> vtmWriter( “test2d” );
SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity( NSlattice, converter );
SuperLatticePhysPressure2D<T, NSDESCRIPTOR> pressure(NSlattice, converter);
SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
SuperLatticePhysHeatFlux2D<T,NSDESCRIPTOR,TDESCRIPTOR> HeatFlux1(ADlattice,converter);

vtmWriter.addFunctor( pressure );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( temperature );
vtmWriter.addFunctor( HeatFlux1 );

AnalyticalFfromSuperF2D<T> interpolation(velocity, true);

const int vtkIter = converter.getLatticeTime( 0.005 );
// const int statIter = converter.getLatticeTime( 10.0 );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank2D<T, NSDESCRIPTOR> rank(NSlattice);
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );

vtmWriter.createMasterFile();
}

// Writes the vtk files
if ( iT%vtkIter==0|| converged ) {

ADlattice.setProcessingContext(ProcessingContext::Evaluation);
NSlattice.setProcessingContext(ProcessingContext::Evaluation);

timer.update(iT);
timer.printStep();

/// NSLattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));
/// ADLattice statistics console output
ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

vtmWriter.write( iT );

//
{
SuperEuklidNorm2D<T, TDESCRIPTOR> normVel( temperature );
BlockReduction2D2D<T> planeReduction(normVel, 600, BlockDataSyncMode::ReduceOnly);

// // write output as ppm
// BlockGifWriter<T> gifWriter;
// gifWriter.write(planeReduction, Tcold*0.98, Thot*1.02, iT, “temperature”);

// write output as JPEG
heatmap::plotParam<T> jpeg_Tem;
jpeg_Tem.maxValue = Thot;
jpeg_Tem.minValue = Tcold;
heatmap::write(planeReduction, iT, jpeg_Tem);
}

{
SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel2( velocity );
BlockReduction2D2D<T> planeReduction2(normVel2, 600, BlockDataSyncMode::ReduceOnly);
// // write output as ppm
// BlockGifWriter<T> gifWriter2;
// gifWriter2.write( planeReduction2, iT, “velocity” );

// write output as JPEG
heatmap::plotParam<T> jpeg_Vel;
jpeg_Vel.maxValue = physU*1.1;
jpeg_Vel.minValue = 0.0;
heatmap::write(planeReduction2, iT, jpeg_Vel);
}

{
SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel3( pressure );
BlockReduction2D2D<T> planeReduction3( normVel3, 600, BlockDataSyncMode::ReduceOnly );
// // write output as ppm
// BlockGifWriter<T> gifWriter3;
// gifWriter3.write( planeReduction3, iT, “pressure” );

// write output as JPEG
heatmap::plotParam<T> jpeg_Pres;
jpeg_Pres.maxValue = physU*1.1;
jpeg_Pres.minValue = 0.0;
heatmap::write(planeReduction3, iT, jpeg_Pres);
}

}

// Writes output on the console

if ( NSlattice.getStatistics().getMaxU() > 0.8 ) {
clout << “PROBLEM uMax=” << NSlattice.getStatistics().getMaxU() << std::endl;
vtmWriter.write( iT );
std::exit( 0 );
}
}

int main( int argc, char* argv[] )
{

// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );
// display messages from every single mpi process
//clout.setMultiOutput(true);

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> converter(
(T) L/N, // physDeltaX: spacing between two lattice cells in __m__
(T) L/(M), // physDeltaT: time step in __s__
(T) L, // charPhysLength: reference length of simulation geometry
(T) 0.045, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 8.05e-7, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 995.7, // physDensity: physical density in __kg / m^3__
(T) 0.618, // physThermalConductivity: physical density in __W/m/K__
(T) Pr * 0.618 / 8.05e-7 / 995.7 , // physSpecificHeatCapacity, J/kg/K
(T) 0.00425, // physThermalExpansionCoefficient K^-1
(T) Tcold, // charPhysLowTemperature
(T) Thot // charPhysHighTemperature
);
// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
// converter.write(“test2D”);

// === 2nd Step: Prepare Geometry ===

// // creating the inflow
// Vector<T,2> inlet(-lx1+converter.getPhysLength(1),0.0);
// // IndicatorCuboid2D<T> inflow(converter.getPhysLength(1), ly1, inlet, theta1);
// // std::shared_ptr<IndicatorF2D<T> > inflow= IndicatorCuboid2D<T> (converter.getPhysLength(1), ly1, inlet, theta1);
// auto inflow=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly1, inlet, theta1) ;

// creating channel1 at the first level
std::vector<T> origin(2,T());
std::shared_ptr<IndicatorF2D<T>> channel1=std::make_shared<IndicatorCuboid2D<T>> (lx1, ly1, origin, theta1) ;

// creating channel2 at the second level
std::vector<T> center2(2,T());
center2[0] = lx1/2+lx2*std::cos(theta2)/2-ly2*std::sin(theta2)/2 ;
center2[1] = lx2*std::sin(theta2)/2+ly2*std::cos(theta2)/2 ;
std::shared_ptr<IndicatorF2D<T>> channel2=std::make_shared<IndicatorCuboid2D<T>> (lx2, ly2, center2, theta2) ;

// // creating the outflow0
// Vector<T,2> outlet0(lx1/2+(lx2-converter.getPhysLength(1)/2)*cos(theta2)+ly2*sin(theta2)/2, (lx2-converter.getPhysLength(1)/2)*sin(theta2)+ly2*cos(theta2)/2);
// // IndicatorCuboid2D<T> outflow0(converter.getPhysLength(1), ly2, outlet0, theta2);
// auto outflow0=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly2, outlet0, theta2) ;

// creating channel3 at the second level
std::vector<T> center3(2,T()) ;
center3[0] = lx1/2+lx2*std::cos(theta2)/2-ly2*std::sin(theta2)/2 ;
center3[1] = -lx2*std::sin(theta2)/2-ly2*std::cos(theta2)/2 ;
std::shared_ptr<IndicatorF2D<T>> channel3=std::make_shared<IndicatorCuboid2D<T>> (lx2, ly2, center3, -theta2) ;

// creating the triangle1
std::vector<T> PointA1(2,T());
PointA1[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2);
PointA1[1] = (lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2) ;
std::vector<T> PointB1(2,T());
PointB1[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2);
PointB1[1] = (lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2) + ly2/std::cos(theta2);
std::vector<T> PointC1(2,T());
PointC1[0] = lx1/2 + lx2*std::cos(theta2);
PointC1[1] = lx2*std::sin(theta2);
std::shared_ptr<IndicatorF2D<T>> triangle1=std::make_shared<IndicatorTriangle2D<T>> (PointA1, PointB1, PointC1) ;

// creating the triangle1
std::vector<T> PointA2(2,T());
PointA2[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2);
PointA2[1] = -((lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2)) ;
std::vector<T> PointB2(2,T());
PointB2[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2);
PointB2[1] = -((lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2) + ly2/std::cos(theta2));
std::vector<T> PointC2(2,T());
PointC2[0] = lx1/2 + lx2*std::cos(theta2);
PointC2[1] = -lx2*std::sin(theta2);
std::shared_ptr<IndicatorF2D<T>> triangle2=std::make_shared<IndicatorTriangle2D<T>> (PointA2, PointB2, PointC2) ;

//Addition of the single geometry to overall geometry
auto bionics = channel1 + channel2 + channel3 – triangle1 – triangle2 ;
IndicatorLayer2D<T> extendedDomain( *bionics, converter.getConversionFactorLength() );

/// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = 2*singleton::mpi().getSize();
#else
const int noOfCuboids = 2;
#endif
CuboidGeometry2D<T> cuboidGeometry(extendedDomain, converter.getPhysDeltaX(), noOfCuboids);

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
SuperGeometry<T,2> superGeometry( cuboidGeometry, loadBalancer, 2);

prepareGeometry(converter, superGeometry, *bionics); //, extendedDomain

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

// ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(
// converter.getLatticeRelaxationFrequency(),
// instances::getBulkMomenta<T,NSDESCRIPTOR>());

// AdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> TbulkDynamics (
// converter.getLatticeThermalRelaxationFrequency(),
// instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>());

// T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
// converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

Timer<T> timer1( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer1.start();

prepareLattice(converter,NSlattice, ADlattice, superGeometry );

timer1.stop();
timer1.printSummary();

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer.start();

util::ValueTracer<T> converge(converter.getLatticeTime(20.),epsilon);
for ( std::size_t iT = 0; iT <= converter.getLatticeTime( maxPhysT ); iT++ ) {

if (converge.hasConverged()) {
clout << “Simulation converged.” << endl;

getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

clout << “Time ” << iT << “.” << std::endl;

break;
}

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

// === 6th Step: Collide and Stream Execution ===

NSlattice.collideAndStream();
ADlattice.collideAndStream();

NSlattice.executeCoupling();

// === 7th Step: Computation and Output of the Results ===
getResults( converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

converge.takeValue(ADlattice.getStatistics().getAverageEnergy(),true);

}

timer.stop();
timer.printSummary();
}

[Fany]

Hello Adrian,

Thanks for your information. I referred to the example “poiseuille2d” to declare the “descriptors::FORCE” at beginning. But the same errors occurred. Could you please have a look at the whole code as follows. Many thanks.

#include “olb2D.h”
#include “olb2D.hh”

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace std;

using T = FLOATING_POINT_TYPE;

using NSDESCRIPTOR = D2Q9<FORCE>;
using BulkDynamics = BGKdynamics<T,NSDESCRIPTOR>;
using ForcedBulkDynamics = ForcedBGKdynamics<T,NSDESCRIPTOR>;
using TDESCRIPTOR = D2Q5<VELOCITY>;

// Stores data from stl file in geometry in form of material numbers
void prepareGeometry( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter, SuperGeometry<T,2>& superGeometry, auto & indicator1)
// , IndicatorF2D<T>& indicator2
{

OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename( 0,2, indicator1); //
superGeometry.rename( 2, 1,{1,1} );
superGeometry.clean();

// creating the inflow
std::vector<T> inlet(2,T());
inlet[0] = -lx1/2+0.5*converter.getPhysLength(1);
inlet[1] = 0.0;
IndicatorCuboid2D<T> inflow(converter.getPhysLength(1), ly1, inlet, theta1);

// creating the outflow0
std::vector<T> outlet0(2,T());
outlet0[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2)-converter.getPhysLength(1);
outlet0[1] = (lx2)*std::sin(theta2)-ly2/std::cos(theta2)*std::pow(std::sin(theta2),2);
std::vector<T> extend1( 2, T(0) );
extend1[0] = converter.getPhysLength(1);
extend1[1] = ly2/std::cos(theta2);
IndicatorCuboid2D<T> outflow0(extend1, outlet0, theta1);
// auto outflow0=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly2, outlet0, theta2) ;

// creating the outflow1
std::vector<T> outlet1(2,T());
outlet1[0] = lx1/2+(lx2)*std::cos(theta2)-ly2*std::sin(theta2)-converter.getPhysLength(1);
outlet1[1] = -(lx2)*std::sin(theta2)+ly2/std::cos(theta2)*std::pow(std::sin(theta2),2)-ly2/std::cos(theta2);
std::vector<T> extend2( 2, T(0) );
extend2[0] = converter.getPhysLength(1);
extend2[1] = ly2/std::cos(theta2);
IndicatorCuboid2D<T> outflow1(extend2, outlet1, theta1);
// auto outflow1=std::make_shared<IndicatorCuboid2D<T>> (converter.getPhysLength(1), ly2, outlet1, theta2) ;

superGeometry.rename( 2,3,1, inflow );
superGeometry.rename( 2,4,1, outflow0 );
superGeometry.rename( 2,5,1, outflow1 );

// 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;

}

// Set up the geometry of the simulation
void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

// // material=0 –> do nothing
// ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

// material=1 –> bulk dynamics
//lattice.defineDynamics( superGeometry,1,&bulkDynamics );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
NSlattice.defineDynamics<BulkDynamics>(superGeometry, 1);

// material=2 –> bounceBack dynamics
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
setBounceBackBoundary(NSlattice, superGeometry, 2);
// setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 2);

// material=3 –> bulk dynamics + velocity (inflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 3 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 3);

// material=4,5 –> bulk dynamics + pressure (outflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 4 );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 5 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 4 );
NSlattice.defineDynamics<BulkDynamics>( superGeometry, 5 );

// Boundary condition
setAdvectionDiffusionTemperatureBoundary(ADlattice, Tomega, superGeometry, 2);

setInterpolatedVelocityBoundary(NSlattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary(NSlattice, omega, superGeometry, 5);
setInterpolatedPressureBoundary(NSlattice, omega, superGeometry, 4);

// AnalyticalConst2D<T,T> rhoF( 1 );
// std::vector<T> velocity( 3,T() );
// AnalyticalConst2D<T,T> uF( velocity );
AnalyticalConst2D<T,T> rhoF(1.);
AnalyticalConst2D<T,T> uF( converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst2D<T,T> u0(0.0, 0.0);
AnalyticalConst2D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst2D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, u0);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, u0);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
NSlattice.setParameter<descriptors::OMEGA>(omega);

// Lattice initialize
NSlattice.initialize();
ADlattice.initialize();

clout << “Prepare Lattice … OK” << std::endl;
}

void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,2>& superGeometry)
{
// / initial condition
OstreamManager clout( std::cout,”setBoundaryValues” );
// No of time steps for smooth start-up
int iTmaxStart = converter.getLatticeTime( maxPhysT*0.25 );
int iTupdate = 30;

if (iT%iTupdate == 0 && iT <= iTmaxStart) {

// T u_Re = converter.getLatticeVelocity( Re * converter.getPhysViscosity() / converter.getCharPhysLength() );
// Smooth start curve, sinus
SinusStartScale<T,int> nSinusStartScale(iTmaxStart, converter.getCharLatticeVelocity());

// Smooth start curve, polynomial
int iTvec[1] = { iT };
T frac[1] = {T()};
nSinusStartScale( frac,iTvec );

T maxVelocity = converter.getCharLatticeVelocity()*3./2.*frac[0];
T distance2Wall = converter.getConversionFactorLength()/2.;
Poiseuille2D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall );

NSlattice.defineU( superGeometry, 3, poiseuilleU );
NSlattice.setProcessingContext< Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
ProcessingContext::Simulation);
}

// if (iT==0){

// std::vector<T> velocity( 2,T() );
// velocity[0] = converter.getCharLatticeVelocity();
// AnalyticalConst2D<T,T> uF( velocity );
// NSlattice.defineU( superGeometry, 3, uF );
// NSlattice.setProcessingContext< Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
// }
}

// Output to console and files
void getResults( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,2>& superGeometry, Timer<T>& timer, bool converged)
{

OstreamManager clout( std::cout,”getResults” );

SuperVTMwriter2D<T> vtmWriter( “test2d” );
SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity( NSlattice, converter );
SuperLatticePhysPressure2D<T, NSDESCRIPTOR> pressure(NSlattice, converter);
SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
SuperLatticePhysHeatFlux2D<T,NSDESCRIPTOR,TDESCRIPTOR> HeatFlux1(ADlattice,converter);

vtmWriter.addFunctor( pressure );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( temperature );
vtmWriter.addFunctor( HeatFlux1 );

AnalyticalFfromSuperF2D<T> interpolation(velocity, true);

const int vtkIter = converter.getLatticeTime( 0.005 );
// const int statIter = converter.getLatticeTime( 10.0 );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank2D<T, NSDESCRIPTOR> rank(NSlattice);
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );

vtmWriter.createMasterFile();
}

// Writes the vtk files
if ( iT%vtkIter==0|| converged ) {

ADlattice.setProcessingContext(ProcessingContext::Evaluation);
NSlattice.setProcessingContext(ProcessingContext::Evaluation);

timer.update(iT);
timer.printStep();

/// NSLattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));
/// ADLattice statistics console output
ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

vtmWriter.write( iT );

//
{
SuperEuklidNorm2D<T, TDESCRIPTOR> normVel( temperature );
BlockReduction2D2D<T> planeReduction(normVel, 600, BlockDataSyncMode::ReduceOnly);

// // write output as ppm
// BlockGifWriter<T> gifWriter;
// gifWriter.write(planeReduction, Tcold*0.98, Thot*1.02, iT, “temperature”);

// write output as JPEG
heatmap::plotParam<T> jpeg_Tem;
jpeg_Tem.maxValue = Thot;
jpeg_Tem.minValue = Tcold;
heatmap::write(planeReduction, iT, jpeg_Tem);
}

{
SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel2( velocity );
BlockReduction2D2D<T> planeReduction2(normVel2, 600, BlockDataSyncMode::ReduceOnly);
// // write output as ppm
// BlockGifWriter<T> gifWriter2;
// gifWriter2.write( planeReduction2, iT, “velocity” );

// write output as JPEG
heatmap::plotParam<T> jpeg_Vel;
jpeg_Vel.maxValue = physU*1.1;
jpeg_Vel.minValue = 0.0;
heatmap::write(planeReduction2, iT, jpeg_Vel);
}

{
SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel3( pressure );
BlockReduction2D2D<T> planeReduction3( normVel3, 600, BlockDataSyncMode::ReduceOnly );
// // write output as ppm
// BlockGifWriter<T> gifWriter3;
// gifWriter3.write( planeReduction3, iT, “pressure” );

// write output as JPEG
heatmap::plotParam<T> jpeg_Pres;
jpeg_Pres.maxValue = physU*1.1;
jpeg_Pres.minValue = 0.0;
heatmap::write(planeReduction3, iT, jpeg_Pres);
}

}

// Writes output on the console

if ( NSlattice.getStatistics().getMaxU() > 0.8 ) {
clout << “PROBLEM uMax=” << NSlattice.getStatistics().getMaxU() << std::endl;
vtmWriter.write( iT );
std::exit( 0 );
}
}

int main( int argc, char* argv[] )
{

// === 1st Step: Initialization ===
olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );
// display messages from every single mpi process
//clout.setMultiOutput(true);

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> converter(
(T) L/N, // physDeltaX: spacing between two lattice cells in __m__
(T) L/(M), // physDeltaT: time step in __s__
(T) L, // charPhysLength: reference length of simulation geometry
(T) 0.045, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 8.05e-7, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 995.7, // physDensity: physical density in __kg / m^3__
(T) 0.618, // physThermalConductivity: physical density in __W/m/K__
(T) Pr * 0.618 / 8.05e-7 / 995.7 , // physSpecificHeatCapacity, J/kg/K
(T) 0.00425, // physThermalExpansionCoefficient K^-1
(T) Tcold, // charPhysLowTemperature
(T) Thot // charPhysHighTemperature
);
// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
// converter.write(“test2D”);

// === 2nd Step: Prepare Geometry ===

STLreader<T> stlReader( “test3d.stl”, converter.getConversionFactorLength() );
IndicatorLayer2D<T> extendedDomain( *bionics, converter.getConversionFactorLength() );

/// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = 2*singleton::mpi().getSize();
#else
const int noOfCuboids = 2;
#endif
CuboidGeometry2D<T> cuboidGeometry(*bionics, converter.getPhysDeltaX(), noOfCuboids);

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
SuperGeometry<T,2> superGeometry( cuboidGeometry, loadBalancer, 2);

prepareGeometry(converter, superGeometry, *bionics); //, extendedDomain

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

// ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(
// converter.getLatticeRelaxationFrequency(),
// instances::getBulkMomenta<T,NSDESCRIPTOR>());

// AdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> TbulkDynamics (
// converter.getLatticeThermalRelaxationFrequency(),
// instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>());

// T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
// converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

Timer<T> timer1( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer1.start();

prepareLattice(converter,NSlattice, ADlattice, superGeometry );

timer1.stop();
timer1.printSummary();

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer.start();

util::ValueTracer<T> converge(converter.getLatticeTime(20.),epsilon);
for ( std::size_t iT = 0; iT <= converter.getLatticeTime( maxPhysT ); iT++ ) {

if (converge.hasConverged()) {
clout << “Simulation converged.” << endl;

getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

clout << “Time ” << iT << “.” << std::endl;

break;
}

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

// === 6th Step: Collide and Stream Execution ===

NSlattice.collideAndStream();
ADlattice.collideAndStream();

NSlattice.executeCoupling();

// === 7th Step: Computation and Output of the Results ===
getResults( converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

converge.takeValue(ADlattice.getStatistics().getAverageEnergy(),true);

}

timer.stop();
timer.printSummary();
}

[Fany]

Dear wf_guoyl,

Your code seemed to simulate single component multiphase flow. I am simulating the flow boiling based on Julius Weinmiller’s work. But it always had a problem of “density is NAN”(“uMax=0.523235; avEnergy=-nan; avRho=nan”). I found uMax sharply increased from 0.007 to 0.53 after the endTimeInit. Have you encounterred this problem?
I would like to ask how you set the MaxVelocity. should it be same with the inlet velocity? I appreciate very much if you can share you code to me.

Best regards,
Fany

[Fany]

Hello Julius,

I would like to ask if you make your code of flow boiling in a channel compatible in OpenLB 1.4. I am doing this topic, but I have no idea of how to implement it. Could you please release your new compatible version if you made it. Many thanks.

Best regards
Fany

[Fany]

Hi Fedor,

Many thanks for your help. I added
“setAdvectionDiffusionConvectionBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry.getMaterialIndicator({3, 4, 5}))”.

but there was still the segmentation fault:
”’
[prepareGeometry] Prepare Geometry … OK
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
[prepareLattice] Prepare Lattice …
————————————————————————–
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
————————————————————————–
————————————————————————–
mpiexec noticed that process rank 0 with PID 55086 on node chem001 exited on signal 11 (Segmentation fault).
”’

Do I need to set other BCs except the following BCs in my case?
”’
// Boundary condition
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 2);

setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({4, 5}));

setAdvectionDiffusionConvectionBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry.getMaterialIndicator({3, 4, 5}));

AnalyticalConst3D<T,T> rhoF(1.);
AnalyticalConst3D<T,T> uF( 0.0, converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst3D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, uF);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, uF);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

NSlattice.defineU( superGeometry, 3, uF );

ADlattice.setParameter<descriptors::OMEGA>(converter.getLatticeThermalRelaxationFrequency());
NSlattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());

// Lattice initialize
NSlattice.initialize();
ADlattice.initialize();

”’

BR，
Fan

[Fany]

Hi Fedor，

I am sorry that I did not see the periodic baoundaries in the rayleighBernard3d.

In my case, I need to simulate the temperature and velocity distribution when the water flows the aorta structure (heating wall). The material 1 is water flowing from inlet to outlet as cooling source. The material 3 is inlet, and 4, 5 are the outlets.

So I have just set the BC of wall (heating source) at constant temperature by the following functor:
”’
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2)
”’
I did not set the temperature BC for the inlet and outlet.

and set the BC of inlet velocity and outlet pressure :
”’
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({4, 5}));
”’

and then set the initial conditions:
”’
AnalyticalConst3D<T,T> rhoF(1.);
AnalyticalConst3D<T,T> uF( 0.0, converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst3D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, uF);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, uF);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

NSlattice.defineU( superGeometry, 3, uF );
”’

All above BC settings can make my case work well in the OpenLB V1.4 but fail in V1.5.

I am not sure if you mean I need to set the Temperature BC at one inlet and two outlets.

BR,
Fany

[Fany]

Hi Fedor,

Many thanks for your reply. I referred to the examples “aorta3d” and “rayleighBenard3d” . It seems that there was just “setAdvectionDiffusionTemperatureBoundary” for heating source needing to be consider. Other materials are just needing to set the initial conditions based on the setup. What BC functions are needed for the Temperature lattice in my case (forced convective heat transfer)?

Best regards,
Fany

[Fany]

Hi Fedor,

Thanks for your reply. There were BC of inlet velocity and outlet pressure, as the above mentioned code:
”’
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({4, 5}));
NSlattice.defineU( superGeometry, 3, uF ); # initial conditions
”’
By the way, the difference of my code between V1.4 (BulkDynamics) and v1.5 (ForcedBGKdynamics) was just different Dynamics. I am not know if there are any changes in the new version. The new guidance do not present it.

Best regards,
Fany

[Fany]

Hi Adrian,

The examples “aorta3d” and “rayleighBenard3d” worked in our cluster.
I was wanting to added the heat transfer into the case “aorta3d” (forced convective heat transfer). The new code can run in the OpenLB 1.4 but not 1.5.
The Segmentation fault always occurred in the step “preparelattice” as following:
”’
[prepareLattice] Prepare Lattice …
————————————————————————–
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
————————————————————————–
————————————————————————–
mpiexec noticed that process rank 1 with PID 43519 on node fluid028 exited on signal 11 (Segmentation fault).
————————————————————————–
2 total processes killed (some possibly by mpiexec during cleanup)
”’

Here is my preparelattice code:
void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
STLreader<T>& stlReader,
SuperGeometry<T,3>& superGeometry )
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

// material=0 –> do nothing
ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

// material=1 –> bulk dynamics
//lattice.defineDynamics( superGeometry,1,&bulkDynamics );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
NSlattice.defineDynamics<ForcedBGKdynamics>(superGeometry, 1);

// material=2 –> bounceBack dynamics
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
NSlattice.defineDynamics<BounceBack>(superGeometry, 2);
// setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 2);

// material=3 –> bulk dynamics + velocity (inflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 3);
NSlattice.defineDynamics<ForcedBGKdynamics>( superGeometry, 3);
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);

// material=4,5 –> bulk dynamics + pressure (outflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry.getMaterialIndicator({4, 5}));
NSlattice.defineDynamics<ForcedBGKdynamics>( superGeometry.getMaterialIndicator({4, 5}));

setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({4, 5}));

// AnalyticalConst3D<T,T> rhoF( 1 );
// std::vector<T> velocity( 3,T() );
// AnalyticalConst3D<T,T> uF( velocity );

AnalyticalConst3D<T,T> rhoF(1.);
AnalyticalConst3D<T,T> uF( 0.0, converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst3D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, uF);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, uF);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

NSlattice.defineU( superGeometry, 3, uF );

ADlattice.setParameter<descriptors::OMEGA>(converter.getLatticeThermalRelaxationFrequency());
NSlattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());

// Lattice initialize
NSlattice.initialize();
ADlattice.initialize();

clout << “Prepare Lattice … OK” << std::endl;
}

I would like to know if the new version transfer the “BulkDynamics” in NSDESCRIPTOR to DESCRIPTOR. In last version, I can use this function from NSDESCRIPTOR. but in version 1.5, there is just ForcedBGKdynamics in NSDESCRIPTOR.

Best regards,
Fany

[Fany]

Hi Adrian,
Thanks for your reply. The example rayleighBenard3d in OpenLB v1.5 can run on this cluster. I am checking the problem.

[Fany]

Hi Adrian,

I ran the code example aorta3d in the OpenLB v1.5 by mpi code. The error occurred as following.
The sequential mode was fine. I would ask what problem it could be. Many thanks.

”’
[main] starting simulation…
mlx5: csk090.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a008778 123c0bd3
mlx5: csk090.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a008775 0c44a0d2
mlx5: csk090.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a00bbdc 0d1d45d2
mlx5: csk090.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a008776 109fcfd2
”’
Here is my setup on cluster:
#!/bin/bash
#SBATCH –time=2:00:00
#SBATCH –nodes=1
#SBATCH –ntasks=8
module load gcc/9.2.0
module load paraview
module load openmpi
module load make
module load gnuplot
make clean
make cleanbuild
make
mpiexec ./aorta3d # mpiexec

[Fany]

Here is the whole error log of MPI mode

[prepareGeometry] Prepare Geometry … OK
mlx5: csk062.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a00db26 096ddcd3
mlx5: csk062.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a00db24 098809d3
[csk062:07276] *** An error occurred in MPI_Wait
[csk062:07276] *** reported by process [815464449,5]
[csk062:07276] *** on communicator MPI COMMUNICATOR 7 DUP FROM 0
[csk062:07276] *** MPI_ERR_INTERN: internal error
[csk062:07276] *** MPI_ERRORS_ARE_FATAL (processes in this communicator will now abort,
[csk062:07276] *** and potentially your MPI job)
[csk062.cluster:07267] 5 more processes have sent help message help-mpi-errors.txt / mpi_errors_are_fatal
[csk062.cluster:07267] Set MCA parameter “orte_base_help_aggregate” to 0 to see all help / error messages

[Fany]

Hi Adrian,
Many thanks for your recommendation. I changed it, but there was a problem (MPI mode) as following:
[prepareGeometry] Prepare Geometry …
mlx5: csk060.cluster: got completion with error:
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 1f006802 0a009a4b 002ffcd2

What’s more, when I ran it by serial mode (MPI mode off), there was still a problem:
[prepareLattice] Prepare Lattice …
————————————————————————–
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
————————————————————————–
————————————————————————–
mpiexec noticed that process rank 0 with PID 26350 on node csk060 exited on signal 11 (Segmentation fault).
————————————————————————–

Here is my code of prepareGeometry, prepareLattice and setboundary:

void prepareGeometry( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter, IndicatorF3D<T>& indicator,
STLreader<T>& stlReader, SuperGeometry<T,3>& superGeometry )
{

OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename( 0,2,indicator );

superGeometry.rename( 2,1,stlReader );
superGeometry.clean();

// Set material number for inflow
IndicatorCircle3D<T> inflow( 0.218125,0.249987,0.0234818, 0., 1.,0., 0.0112342 );
IndicatorCylinder3D<T> layerInflow( inflow, 2.*converter.getConversionFactorLength() );
superGeometry.rename( 2,3,1,layerInflow );

// Set material number for outflow0
//IndicatorCircle3D<T> outflow0(0.2053696,0.0900099,0.0346537, 2.5522,5.0294,-1.5237, 0.0054686 );
IndicatorCircle3D<T> outflow0( 0.2053696,0.0900099,0.0346537, 0.,-1.,0., 0.0054686 );
IndicatorCylinder3D<T> layerOutflow0( outflow0, 2.*converter.getConversionFactorLength() );
superGeometry.rename( 2,4,1,layerOutflow0 );

// Set material number for outflow1
//IndicatorCircle3D<T> outflow1(0.2388403,0.0900099,0.0343228, -1.5129,5.1039,-2.8431, 0.0058006 );
IndicatorCircle3D<T> outflow1( 0.2388403,0.0900099,0.0343228, 0.,-1.,0., 0.0058006 );
IndicatorCylinder3D<T> layerOutflow1( outflow1, 2.*converter.getConversionFactorLength() );
superGeometry.rename( 2,5,1,layerOutflow1 );

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean( 3 );

superGeometry.checkForErrors();

superGeometry.print();
clout << “Prepare Geometry … OK” << std::endl;
}

// Set up the geometry of the simulation
void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,

STLreader<T>& stlReader,
SuperGeometry<T,3>& superGeometry )
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

// material=0 –> do nothing
ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

// material=1 –> bulk dynamics
//lattice.defineDynamics( superGeometry,1,&bulkDynamics );
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 1);
NSlattice.defineDynamics<ForcedBGKdynamics>(superGeometry, 1);

if ( bouzidiOn ) {
// // material=2 –> no dynamics + bouzidi zero velocity
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);

NSlattice.defineDynamics<NoDynamics>(superGeometry, 2);
setBouzidiZeroVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 2, stlReader);

// // material=3 –> no dynamics + bouzidi velocity (inflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 3);
NSlattice.defineDynamics<NoDynamics>( superGeometry, 3);
setBouzidiVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 3, stlReader);

}
else {
// material=2 –> bounceBack dynamics
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry, 2);
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
NSlattice.defineDynamics<BounceBack>(superGeometry, 2);
// setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 2);

// material=3 –> bulk dynamics + velocity (inflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry, 3);
NSlattice.defineDynamics<ForcedBGKdynamics>( superGeometry, 3);
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);

}

// material=4,5 –> bulk dynamics + pressure (outflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry.getMaterialIndicator({4, 5}));
NSlattice.defineDynamics<ForcedBGKdynamics>( superGeometry.getMaterialIndicator({4, 5}));

setInterpolatedPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({4, 5}));

// AnalyticalConst3D<T,T> rhoF( 1 );
// std::vector<T> velocity( 3,T() );
// AnalyticalConst3D<T,T> uF( velocity );

ADlattice.setParameter<descriptors::OMEGA>(converter.getLatticeThermalRelaxationFrequency());
NSlattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());

clout << “Prepare Lattice … OK” << std::endl;
}

void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,3>& superGeometry)
{

if (iT == 0) {

// T u_Re = converter.getLatticeVelocity( Re * converter.getPhysViscosity() / converter.getCharPhysLength() );

AnalyticalConst3D<T,T> rhoF(1.);
AnalyticalConst3D<T,T> uF( 0.0, converter.getLatticeVelocity(Re * converter.getPhysViscosity() / converter.getCharPhysLength()), 0.0);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst3D<T,T> T_hot(converter.getLatticeTemperature(Thot));

// Initialize all values of distribution functions to their local equilibrium
NSlattice.defineRhoU(superGeometry, 3, rhoF, uF);
NSlattice.iniEquilibrium(superGeometry, 3, rhoF, uF);
NSlattice.defineRhoU( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );
NSlattice.iniEquilibrium( superGeometry.getMaterialIndicator({1, 4, 5}),rhoF, u0 );

ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({ 4, 5}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({ 4, 5}), T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

NSlattice.defineU( superGeometry, 3, uF );

// Lattice initialize
NSlattice.initialize();
ADlattice.initialize();

}

}

[Fany]

Here is parts of my code:

#include “olb3D.h”
#include “olb3D.hh”

// #ifndef OLB_PRECOMPILED // Unless precompiled version is used,
// #include “olb3D.hh” // include full template code;
// #endif
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace std;

typedef double T;

typedef D3Q19<FORCE> NSDESCRIPTOR;
typedef D3Q7<VELOCITY> TDESCRIPTOR;

//simulation parameters
const int N = 35; // resolution of the model
const int M = 100; // time discretization refinement
const bool bouzidiOn = false; // choice of boundary condition
const T maxPhysT = 50.; // max. simulation time in s, SI unit

// === 2nd Step: Prepare Geometry ===

// Instantiation of the STLreader class
// file name, voxel size in meter, stl unit in meter, outer voxel no., inner voxel no.
STLreader<T> stlReader( “test3d.stl”, converter.getConversionFactorLength(), 0.001, 0, true );
IndicatorLayer3D<T> extendedDomain( stlReader, converter.getConversionFactorLength() );

// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids =std::min( 16*N,2*singleton::mpi().getSize() ) ; //2*singleton::mpi().getSize()
#else
const int noOfCuboids = 2;
#endif
CuboidGeometry3D<T> cuboidGeometry( extendedDomain, converter.getConversionFactorLength(), noOfCuboids );

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
SuperGeometry<T,3> superGeometry( cuboidGeometry, loadBalancer, 2);

prepareGeometry( converter, extendedDomain, stlReader, superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

// ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(
// converter.getLatticeRelaxationFrequency(),
// instances::getBulkMomenta<T,NSDESCRIPTOR>());

// AdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> TbulkDynamics (
// converter.getLatticeThermalRelaxationFrequency(),
// instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>());

// std::vector<T> dir{0.0, 1.0, 0.0};

// T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
// converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

///T boussinesqForcePrefactor=0;

//NavierStokesAdvectionDiffusionCouplingGenerator3D<T,NSDESCRIPTOR> coupling(0, converter.getLatticeLength(0.2388403+0.0058006 ), 0, converter.getLatticeLength(0.249987+0.0112342),
//0, converter.getLatticeLength(0.02246+0.0054686), boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir);

//NSlattice.addLatticeCoupling(coupling, ADlattice);

Timer<T> timer1( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer1.start();

prepareLattice(converter,
NSlattice, ADlattice,

stlReader,
superGeometry );

timer1.stop();
timer1.printSummary();

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer.start();

for ( std::size_t iT = 0; iT <= converter.getLatticeTime( maxPhysT ); iT++ ) {

// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

// === 6th Step: Collide and Stream Execution ===

NSlattice.collideAndStream();
ADlattice.collideAndStream();

//NSlattice.executeCoupling();

// === 7th Step: Computation and Output of the Results ===
getResults( converter, NSlattice, ADlattice, iT, superGeometry, timer, stlReader);

}

timer.stop();
timer.printSummary();
}

[Author]

Posts