Reply To: Problem with defining linear temprature in RayleighBernard2D example

[navidkyo]

Thank you Mathias for your help. I used the following code for defining the linear function and it seems the problem is solved: (The rest f the code is the same as above)

//*****************************************************************************************
AnalyticalLinear2D<T,T> T_ini((T)0, (converter.getLatticeTemperature(Tcold), converter.getLatticeTemperature(Thot))/ly,converter.getLatticeTemperature(Thot) );
//*****************************************************************************************

However, as I checked the initial velocity at t=0, I got a linear symmetrical distribution for the initial velocity as well, and I wondered why this happened. Is it due to velocity shift by forcing term? or I did I miss a point here?

I attached my code:

//*****************************************************************************************

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

using namespace olb; //OK
using namespace olb::descriptors; //OK
using namespace olb::graphics; //OK

typedef double T; //OK

typedef D2Q9<FORCE> NSDESCRIPTOR; //OK
typedef D2Q5<VELOCITY> TDESCRIPTOR; //OK

// Parameters for the simulation setup

// const int SavingIter = 10.; ///use din line 316
const int N = 60; // resolution of the model //OK
const T L = 0.1/N; // latticeL //OK
const T lx = 6. ; // length of the channel //OK
const T ly = 4.+L; // height of the channel //OK
// const T Ra = 1e4; // Rayleigh number2000
// const T Pr = 0.71; // Prandtl number
const T Re = 76; //OK
const T maxPhysT = 500; // max. simulation time in s, SI unit //OK
const T epsilon = 1.e-5; // precision of the convergence (residuum) //NOT IN JUNWEI’S CODE*******************************

const T Thot = 5.; // temperature of the lower wall in Kelvin //OK
const T Tcold = 1.; // temperature of the fluid in Kelvin //OK
const T Tave = (Thot+Tcold)/2.0; // temperature of the fluid in Kelvin //OK
const T Tperturb = 1./5. * Tcold + 4./5. * Thot; // temperature of the perturbation //OK

const T centerCylinderX = 2.5; //OK
const T centerCylinderY = ly/2.0 +L/2.; //OK
const T radiusCylinder = 0.05; //OK

/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry<T,2>& superGeometry,
ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter)
{

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

superGeometry.rename(0,2); //OK
superGeometry.rename(2,1,{0,1}); //original: superGeometry.rename(2,1,0,1); //OK

std::vector<T> extend( 2, T(0) ); //OK
std::vector<T> origin( 2, T(0) ); //OK

origin[1] = -converter.getPhysLength(1); //OK
extend[0] = lx ; //OK
extend[1] = 2.*converter.getPhysLength(1); //OK
IndicatorCuboid2D<T> bottom(extend, origin); //OK

origin[1] = ly-converter.getPhysLength(1); //OK
IndicatorCuboid2D<T> top(extend, origin); //OK

superGeometry.rename(2,2,1,bottom);
superGeometry.rename(2,3,1,top);

// origin[0] = lx/2.;
// origin[1] = converter.getPhysLength(1);
// extend[0] = converter.getPhysLength(1);
// extend[1] = converter.getPhysLength(1);
// IndicatorCuboid2D<T> perturbation(extend, origin);

//************************ Inflow ******************** added by me
origin[0] = 0.0; //OK
origin[1] = 0.0; //OK
extend[0] = converter.getPhysLength(1.); //OK
extend[1] = ly; //OK
IndicatorCuboid2D<T> inflow(extend, origin); //OK
superGeometry.rename( 1,5,inflow ); //OK

//************************ outflow ******************** added by me
origin[0] = lx-converter.getPhysLength(1); //OK
extend[0] = 2.*converter.getPhysLength(1); //OK
IndicatorCuboid2D<T> outflow(extend, origin); //OK
superGeometry.rename( 1,6,1,outflow ); //????????????????

//************************ Cylinder ******************** added by me

Vector<T,2> center( centerCylinderX,centerCylinderY ); //OK
IndicatorCircle2D<T> circle( center, radiusCylinder ); //OK
superGeometry.rename( 1,4,circle ); //OK

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

superGeometry.print(); //OK

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

void prepareLattice( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,2>& superGeometry )
{

OstreamManager clout(std::cout,”prepareLattice”); //OK

T Tomega = converter.getLatticeThermalRelaxationFrequency(); //OK
T NSomega = converter.getLatticeRelaxationFrequency(); //ADDED BY ME //OK

/// define lattice Dynamics
clout << “defining dynamics” << std::endl; //OK

ADlattice.defineDynamics<NoDynamics>(superGeometry, 0); //OK
ADlattice.defineDynamics<AdvectionDiffusionBGKdynamics>(superGeometry.getMaterialIndicator({1, 2, 3,5,6})); //OK
ADlattice.defineDynamics<BounceBack>(superGeometry, 4); //OK

NSlattice.defineDynamics<NoDynamics>(superGeometry, 0); //OK
NSlattice.defineDynamics<ForcedBGKdynamics>(superGeometry, 1); //OK
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 2); //OK
// setSlipBoundary(NSlattice, superGeometry, 2); //OK
// NSlattice.defineDynamics<NoDynamics>(superGeometry, 3); //OK
// setSlipBoundary(NSlattice, superGeometry, 3); //OK
// NSlattice.defineDynamics<BounceBack>(superGeometry, 4); //OK

NSlattice.defineDynamics<BounceBack>(superGeometry, 2); //OK
NSlattice.defineDynamics<BounceBack>(superGeometry, 3); //OK

// Define Boundaries
setInterpolatedVelocityBoundary(NSlattice, NSomega, superGeometry, 5); //OK
setInterpolatedVelocityBoundary(NSlattice, NSomega, superGeometry, 6); //OK
// setLocalPressureBoundary<T,NSDESCRIPTOR>(NSlattice, NSomega, superGeometry, 6);

setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2); //OK
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 3); //OK

Vector<T,2> center(centerCylinderX,centerCylinderY );
IndicatorCircle2D<T> circle ( center, radiusCylinder );
setBouzidiZeroVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 4, circle);

// Vector<T,2> center( centerCylinderX,centerCylinderY );
// IndicatorCircle2D<T> circle ( center, radiusCylinder );
// NSlattice.defineDynamics<NoDynamics<T,NSDESCRIPTOR>>(superGeometry, 7);
// setBouzidiZeroVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 7, circle);
// ADlattice.defineDynamics<BounceBack<T,TDESCRIPTOR>>(superGeometry, 7);

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

/// define initial conditions
AnalyticalConst2D<T,T> rho(1.); //OK
AnalyticalConst2D<T,T> u0(0.0, 0.0); //OK
AnalyticalConst2D<T,T> T_cold(converter.getLatticeTemperature(Tcold)); //OK
AnalyticalConst2D<T,T> T_hot(converter.getLatticeTemperature(Thot)); //OK
AnalyticalConst2D<T,T> T_perturb(converter.getLatticeTemperature(Tperturb)); //OK
// AnalyticalLinear2D<T,T> T_ini( 0, -0.25, 0.25); //OK&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

AnalyticalLinear2D<T,T> T_ini((T)0, (converter.getLatticeTemperature(Tcold)-converter.getLatticeTemperature(Thot))/4.,converter.getLatticeTemperature(Thot) );

AnalyticalConst2D<T,T> T_ave((converter.getLatticeTemperature(Thot)+converter.getLatticeTemperature(Tcold))/2); //OK

/// for each material set Rho, U and the Equilibrium

// NSlattice.defineU(superGeometry,5,u0); //OK
// NSlattice.defineU(superGeometry,6,u0); //OK

// ADlattice.defineRho(superGeometry, 1, T_ini);
// ADlattice.iniEquilibrium(superGeometry, 1, T_ini, u0);
// ADlattice.defineRho(superGeometry, 2, T_hot);
// ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
// ADlattice.defineRho(superGeometry, 3, T_cold);
// ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

// ADlattice.defineRho(superGeometry, 5, T_ini);
// ADlattice.iniEquilibrium(superGeometry, 5, T_ini, u0);
// ADlattice.defineRho(superGeometry, 6, T_ini);
// ADlattice.iniEquilibrium(superGeometry, 6, T_ini, u0);

ADlattice.defineRho(superGeometry, 1, T_ini);
ADlattice.iniEquilibrium(superGeometry, 1, T_ini, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u0);

ADlattice.defineRho(superGeometry, 5, T_ini);
ADlattice.iniEquilibrium(superGeometry, 5, T_ini, u0);
ADlattice.defineRho(superGeometry, 6, T_ini);
ADlattice.iniEquilibrium(superGeometry, 6, T_ini, u0);

NSlattice.defineRhoU(superGeometry, 1, rho, u0); //OK
NSlattice.iniEquilibrium(superGeometry, 1, rho, u0); //OK
NSlattice.defineRhoU(superGeometry, 2, rho, u0); //OK
NSlattice.iniEquilibrium(superGeometry, 2, rho, u0); //OK
NSlattice.defineRhoU(superGeometry, 3, rho,u0); //OK
NSlattice.iniEquilibrium(superGeometry, 3, rho, u0); //OK
NSlattice.defineRhoU(superGeometry, 5, rho,u0); //OK
NSlattice.iniEquilibrium(superGeometry, 5, rho, u0); //OK
NSlattice.defineRhoU(superGeometry, 6, rho,u0); //OK
NSlattice.iniEquilibrium(superGeometry, 6, rho, u0); //OK

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

/// Make the lattice ready for simulation
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)
{
{

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

// No of time steps for smooth start-up
int iTmaxStart = converter.getLatticeTime( maxPhysT*0.4 ); //original: int iTmaxStart = converter.getLatticeTime( maxPhysT*0.4 );
int iTupdate = 20; ///frequency of printing max vel on the console

if ( iT%iTupdate==0 && iT<= iTmaxStart ) {
// Smooth start curve, sinus
// SinusStartScale<T,int> StartScale(iTmaxStart, T(1));

// Smooth start curve, polynomial
PolynomialStartScale<T,T> StartScale( iTmaxStart, T( 1 ) );

// Creates and sets the Poiseuille inflow profile using functors
T iTvec[1] = {T( iT )};
T frac[1] = {};
StartScale( frac,iTvec );
// std::vector<T> maxVelocity( 2,0 );
T maxVelocity = converter.getCharLatticeVelocity()*frac[0];
// T maxVelocity2 = 0.0*frac[0];
T distance2Wall = converter.getPhysLength(1)/2.;
Poiseuille2D<T> poiseuilleU( superGeometry, 5, maxVelocity, distance2Wall );
// AnalyticalConst2D<T,T> poiseuilleU(maxVelocity, 0.0);

NSlattice.defineU( superGeometry, 5, poiseuilleU );
NSlattice.defineU( superGeometry, 6, poiseuilleU );

NSlattice.setProcessingContext<Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
ProcessingContext::Simulation);
clout << “step=” << iT << “; maxVel=” << maxVelocity << “; (%)maxVel=” << maxVelocity/converter.getCharLatticeVelocity() << std::endl;

}
}

// OstreamManager clout( std::cout,”setBoundaryValues” );

// // No of time steps for smooth start-up
// int iTmaxStart = converter.getLatticeTime( maxPhysT*0.4 );
// int iTupdate = 5;

// if ( iT%iTupdate==0 && iT<= iTmaxStart ) {
// // Smooth start curve, sinus
// // SinusStartScale<T,int> StartScale(iTmaxStart, T(1));

// // Smooth start curve, polynomial
// PolynomialStartScale<T,T> StartScale( iTmaxStart, T( 1 ) );

// // Creates and sets the Poiseuille inflow profile using functors
// T iTvec[1] = {T( iT )};
// T frac[1] = {};
// StartScale( frac,iTvec );
// T maxVelocity = converter.getCharLatticeVelocity()/200.*frac[0]; //*************CHECK: SETS MAX INLET VELOCITY ORIGINAL: 3./20.
// //T maxVelocity = 0.005 ; IT DOES NOT WORK LIKE THIS!!!!! maxVelocity is an ARRAY
// T distance2Wall = converter.getPhysLength(1)/2.;
// //Poiseuille2D<T> poiseuilleU( superGeometry, 5, maxVelocity, distance2Wall );
// AnalyticalConst2D<T,T> poiseuilleU( maxVelocity );
// //AnalyticalConst2D<T,T> constantU(maxVelocity, 0.0);
// NSlattice.defineU( superGeometry, 5, poiseuilleU );
// clout << “step=” << iT << “; maxVel=” << maxVelocity << std::endl;

// }
}

void getResults(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,2>& superGeometry,
util::Timer<T>& timer,
bool converged)
{

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

SuperVTMwriter2D<T> vtkWriter(“rayleighBenard2d”);
SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity(NSlattice, converter);
SuperLatticePhysPressure2D<T, NSDESCRIPTOR> presure(NSlattice, converter);
SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
vtkWriter.addFunctor( presure );
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( temperature );

SuperLatticeVelocity2D<T, NSDESCRIPTOR> velocityLatticeUnit(NSlattice);

const int saveIter = converter.getLatticeTime(1.); //fffffffffffffffffffffffffffffff

SuperLatticeField2D<T,NSDESCRIPTOR,olb::descriptors::FORCE> bodyForce (NSlattice); //ADDED BY ME
vtkWriter.addFunctor( bodyForce ); //ADDED BY ME

if (iT == 0) {
/// Writes the converter log file
// writeLogFile(converter,”rayleighBenard2d”);

/// 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);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);

vtkWriter.createMasterFile();
}

/// Writes the VTK files and prints statistics
if (iT%saveIter == 0 || converged) {
/// Timer console output
timer.update(iT);
timer.printStep();

/// Lattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));

vtkWriter.write(iT);

BlockReduction2D2D<T> planeReduction(temperature, 600, BlockDataSyncMode::ReduceOnly);
BlockGifWriter<T> gifWriter;
gifWriter.write(planeReduction, Tcold-0.1, Thot+0.1, iT, “temperature”);

// &***************************** Drag, lift, pressure drop &**********************************************
static Gnuplot<T> gplot( “drag” );
const int vtkIter = converter.getLatticeTime( .3 );

// Writes the vtk files
if ( iT%vtkIter == 0 && iT > 0 ) {
vtkWriter.write( iT );
}

// write pdf at last time step
if ( iT == converter.getLatticeTime( maxPhysT )-1 ) {
// writes pdf
gplot.writePDF();
}

//******************************
AnalyticalFfromSuperF2D<T> intpolatePressure( presure, true );
SuperLatticePhysDrag2D<T,NSDESCRIPTOR> drag( NSlattice, superGeometry, 4, converter );

T point1[2] = {};
T point2[2] = {};

point1[0] = centerCylinderX – radiusCylinder;
point1[1] = centerCylinderY;

point2[0] = centerCylinderX + radiusCylinder;
point2[1] = centerCylinderY;

T p1, p2;
intpolatePressure( &p1,point1 );
intpolatePressure( &p2,point2 );

clout << “pressure1=” << p1;
clout << “; pressure2=” << p2;

T pressureDrop = p1-p2;
clout << “; pressureDrop=” << pressureDrop;

int input[3] = {};
T _drag[drag.getTargetDim()];
drag( _drag,input );
clout << “; drag=” << _drag[0] << “; lift=” << _drag[1] << std::endl;

// set data for gnuplot: input={xValue, yValue(s), names (optional), position of key (optional)}
gplot.setData( converter.getPhysTime( iT ), {_drag[0], _drag[1]}, {“drag(openLB)”, “lift(OpenLB)”}, “bottom right”, {‘l’,’l’} );
// writes a png in one file for every timestep, if the file is open it can be used as a “liveplot”
gplot.writePNG();

// every (iT%vtkIter) write an png of the plot
if ( iT%( vtkIter ) == 0 ) {
// writes pngs: input={name of the files (optional), x range for the plot (optional)}
gplot.writePNG( iT, maxPhysT );
}

}
}

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

/// === 1st Step: Initialization ===
OstreamManager clout(std::cout,”main”);
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> converter(
(T) 0.1/N, // physDeltaX
(T) 0.00043860, // physDeltaT = charLatticeVelocity / charPhysVelocity * physDeltaX
(T) 0.1, // charPhysLength
(T) Re*1.e-4/0.1, // charPhysVelocity
(T) 1.e-4, // physViscosity
(T) 1., // physDensity
(T) 0.03, // physThermalConductivity
(T) 210000, // physSpecificHeatCapacity
(T) 0.000143573, // physThermalExpansionCoefficient
(T) Tcold, // charPhysLowTemperature
(T) Thot // charPhysHighTemperature
);
converter.print();
converter.write(“olbUnitConverter”);

/// === 2nd Step: Prepare Geometry ===
std::vector<T> extend(2,T());
extend[0] = lx;
extend[1] = ly;
std::vector<T> origin(2,T());
IndicatorCuboid2D<T> cuboid(extend, origin);

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

cuboidGeometry.setPeriodicity(false, false);

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

SuperGeometry<T,2> superGeometry(cuboidGeometry, loadBalancer);

prepareGeometry(superGeometry, converter);

/// === 3rd Step: Prepare Lattice ===

SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
// This coupling must be necessarily be put on the Navier-Stokes lattice!!
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//

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

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

NavierStokesAdvectionDiffusionCouplingGenerator2D<T,NSDESCRIPTOR> coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), boussinesqForcePrefactor, (converter.getLatticeTemperature(Tcold)+converter.getLatticeTemperature(Thot))/2., 1., dir);

NSlattice.addLatticeCoupling(coupling, ADlattice);

prepareLattice(converter, NSlattice, ADlattice, superGeometry);

/// === 4th Step: Main Loop with Timer ===
util::Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );
timer.start();

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

if (converge.hasConverged()) {
clout << “Simulation converged.” << std::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 ===
ADlattice.collideAndStream();
NSlattice.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();
}

//*****************************************************************************************

Thank you very much