Convective heat transfer_Segmentation fault for OpenLB 1.6

[Fany]

Dear all,

I am working on the forced convective heat transfer in a 2D rectangular fractal channel. But I was always faced with the problem of Segmentation fault (11) in the Preparing Lattice. Is that caused by wrong set of boundary conditions? I would like to ask how to set the boundary conditions. The code is as follows:

”’
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=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);
setBounceBackBoundary(NSlattice, 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);

// material=4,5 –> bulk dynamics + pressure (outflow)
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>( superGeometry.getMaterialIndicator({4, 5}));
NSlattice.defineDynamics<ForcedBGKdynamics>( superGeometry.getMaterialIndicator({4, 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.);
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;
}
”’

[Adrian]

Your code excerpt looks fine. The only potential issue I can currently suspect is that you did not declare the descriptors::FORCE field required by the ForcedBGKdynamics dynamics in the NSDESCRIPTOR. This is not necessary and would work if you actually explicitly initialized the field in your code (which would transparently allocate it). However, if it is not declared and not initialized this could potentially cause an access error depending on the platform (which would of course still be a problem to be fixed with a better error message).

If this is not the cause I’ll need more details / a full code example.

[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();
}

[Adrian]

Thanks, but I don’t think that this is the full code. It fails to compile due to various missing variables such as lx1.

[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();
}

[Adrian]

The issue is that the geometry contains boundary features for which setAdvectionDiffusionTemperatureBoundary is not well defined, causing a faulty parameter definition bug. You can mitigate this by either removing line 96 of src/boundary/setAdvectionDiffusionTemperatureBoundary2D.hh or replacing it by

      if (dynamics) {
        dynamics->getParameters(block).template set<descriptors::OMEGA>(omega);
      }

although this separate parameter definition is a legacy artifact and unnecessary due to the following global parameter setter.

[Fany]

Hello Adrian,

Many thanks for your help. The code works now.

Have a nice day!

[Author]

Posts