Reply To: Convective heat transfer_Segmentation fault for OpenLB 1.6
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Convective heat transfer_Segmentation fault for OpenLB 1.6 › Reply To: Convective heat transfer_Segmentation fault for OpenLB 1.6
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();
}