OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Convective heat transfer_Segmentation fault for OpenLB 1.6

Viewing 7 posts - 1 through 7 (of 7 total)
• Author
Posts
• #7396
Fany
Participant

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,
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 );
NSlattice.defineDynamics<ForcedBGKdynamics>(superGeometry, 1);

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

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

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

// Boundary condition
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 );

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

// Lattice initialize
NSlattice.initialize();

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

#7397
Keymaster

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.

#7401
Fany
Participant

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,
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
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

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

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

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

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

// Boundary condition

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

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

// Lattice initialize
NSlattice.initialize();

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

void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
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,
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);

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 ) {

NSlattice.setProcessingContext(ProcessingContext::Evaluation);

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

/// NSLattice statistics console output
NSlattice.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();
// converter.write(“test2D”);

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

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 superGeometry

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

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

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

// converter.getLatticeThermalRelaxationFrequency(),

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

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

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 ===

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

NSlattice.collideAndStream();

NSlattice.executeCoupling();

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

}

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

#7411
Keymaster

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

#7412
Fany
Participant

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,
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
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

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

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

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

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

// Boundary condition

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

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

// Lattice initialize
NSlattice.initialize();

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

void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
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,
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);

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 ) {

NSlattice.setProcessingContext(ProcessingContext::Evaluation);

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

/// NSLattice statistics console output
NSlattice.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();
// 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 superGeometry

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

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

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

// converter.getLatticeThermalRelaxationFrequency(),

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

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

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 ===

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

NSlattice.collideAndStream();

NSlattice.executeCoupling();

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

}

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

#7413
Keymaster

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.

#7414
Fany
Participant