Reply To: Convective heat transfer_Segmentation fault for OpenLB 1.6

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