Skip to content

Poiseuille Profile Test – Diverge

OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Poiseuille Profile Test – Diverge

Viewing 2 posts - 1 through 2 (of 2 total)
  • Author
    Posts
  • July 4, 2016 at 3:37 pm #1843
    Alejandro_Claro
    Member

    Hi everyone,rnrnI am trying to implement the Poiseuille2D example with the next changes:rn

      rn1. The dimensions have been set to 0.2 x 1.0 (height x long).rn2. The fluid been simulated is the water at 20°C -> molecular viscosity: 1×10⁻⁶ m²/s; density: 998.2 kg/m³.rn3. The grid spacing is 0.001 m and the velocity in lattice units is 0.001.rn4. I have set a pressure level of 101325 Pa in order to indicate atmospheric condition.rn5. The Poiseuille profile is implemented in order to obtain the maximum velocity close to the top wall.rn

    rnrnI simulate the fluid flow at Re = 500 with the ForcedBGKdynamics class. The inlet condition use the Skodors algorithm with the regularized scheme of Latt. The outlet is defined as a outflow condition using the convection boundary method. The top wall is defined as a freeSlip conditon by means of a specular reflection. The bottom wall is simulated with the Skodors algorithm by setting a velocity of zero. rnrnNormally, this testing problem should not diverge but I at some point of the simulation my case “”explose””. Here is the link for seen the images that I get from the simulation: rn

    rnrnI copy my code here in order to present you all the option that I have applied to this case:rn

    Code:
    rn/* bgkPoiseuille2d.cpp:rn * This example examines a 2D Poseuille flow with a velocityrn * or pressure boundary at the inlet/outlet.rn * Computation of error norms via functors is shown as well.rn */rnrnrn#include “”olb2D.h””rn#ifndef OLB_PRECOMPILED // Unless precompiled version is used,rn#include “”olb2D.hh”” // include full template codern#endifrn#include <vector>rn#include <cmath>rn#include <iostream>rn#include <iomanip>rn#include <fstream>rnrnusing namespace olb;rnusing namespace olb::descriptors;rnusing namespace olb::graphics;rnusing namespace std;rnrn////////////////////////////////////////////////////////////////////////rnstruct ForcedShearSmagorinsky2dDescriptor {rn static const int numScalars = 3;rn static const int numSpecies = 2;rn static const int avShearIsAt = 0;rn static const int sizeOfAvShear = 1;rn static const int forceBeginsAt = 1;rn static const int sizeOfForce = 2;rn };rn rn struct ForcedShearSmagorinsky2dDescriptorBase {rn typedef ForcedShearSmagorinsky2dDescriptor ExternalField;rn };rn rn template <typename T> struct ForcedShearSmagorinskyD2Q9Descriptorrn : public D2Q9DescriptorBase<T>, public ForcedShearSmagorinsky2dDescriptorBase {rn };rn////////////////////////////////////////////////////////////////////////rnrntypedef enum {velocity, pressure} BoundaryType;rnrntypedef double T;rn//#define DESCRIPTOR D2Q9Descriptorrn//#define DESCRIPTOR ForcedD2Q9Descriptorrn#define DESCRIPTOR ForcedShearSmagorinskyD2Q9Descriptorrnrn////// turbulence model ////////////////////////////////////////////////rn#define Laminairern//#define ShearSmagorinskyWallrn////////////////////////////////////////////////////////////////////////rnrn// Parameters for the simulation setuprn//const T lx = 2.; // length of the channelrnconst T lx = 1.; // length of the channelrn//const T ly = 1.; // height of the channelrnconst T ly = 0.2; // height of the channelrnconst int N = 1000; // resolution of the modelrnconst int M = N; // time discretization refinementrn//const T Re = 10.; // Reynolds numberrnconst T Re = 8400.; // Reynolds numberrnconst T Re_initial = 500.; // Reynolds number initialrnconst T maxPhysT = 50.; // max. simulation time in s, SI unitrnrnrn/// Stores geometry information in form of material numbersrnvoid prepareGeometry(LBconverter<T> const& converter,rn SuperGeometry2D<T>& superGeometry) {rnrn OstreamManager clout(std::cout,””prepareGeometry””);rn clout << “”Prepare Geometry …”” << std::endl;rnrn superGeometry.rename(0,2);rnrn superGeometry.rename(2,1,1,1);rnrn std::vector<T> extend(2,T());rn std::vector<T> origin(2,T());rnrn /// Set material number for inflowrn extend[1] = ly;rn IndicatorCuboid2D<T> inflow(extend, origin);rn superGeometry.rename(2,3,1,inflow);rnrn /// Set material number for outflowrn origin[0] = lx;rn IndicatorCuboid2D<T> outflow(extend, origin);rn superGeometry.rename(2,4,1,outflow);rn rn /// Set material number for top wallrn origin[0] = (T)0.;rn origin[1] = (T)ly-0.1*converter.getLatticeL();rn extend[0] = lx;rn extend[1] = (T)0.2*converter.getLatticeL();rn IndicatorCuboid2D<T> TopWall(extend, origin);rn superGeometry.rename(2,5,TopWall);rn rn /// Removes all not needed boundary voxels outside the surfacern superGeometry.clean();rn /// Removes all not needed boundary voxels inside the surfacern superGeometry.innerClean();rn superGeometry.checkForErrors();rnrn superGeometry.print();rnrn clout << “”Prepare Geometry … OK”” << std::endl;rn}rnrn/// Set up the geometry of the simulationrnvoid prepareLattice(LBconverter<T> const& converter,rn SuperLattice2D<T, DESCRIPTOR>& sLattice,rn Dynamics<T, DESCRIPTOR>& bulkDynamics,rn sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,rn BoundaryType inflowBoundary, BoundaryType outflowBoundary,rn SuperGeometry2D<T>& superGeometry) {rnrn OstreamManager clout(std::cout,””prepareLattice””);rn clout << “”Prepare Lattice …”” << endl;rnrn T omega = converter.getOmega();rnrn /// Material=0 –>do nothingrn sLattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());rnrn /// Material=1 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 1, &bulkDynamics);rnrn /// Material=2 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 2, &bulkDynamics);rnrn /// Material=3 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 3, &bulkDynamics);rnrn /// Material=4 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 4, &bulkDynamics);rn rn /// Material=4 –>bulk dynamicsrn sLattice.defineDynamics(superGeometry, 5, &bulkDynamics);rnrn /// Setting of the boundary conditionsrn sBoundaryCondition.addVelocityBoundary(superGeometry, 2, omega);rn //sBoundaryCondition.addVelocityBoundary(superGeometry, 5, omega);rn sBoundaryCondition.addSlipBoundary(superGeometry, 5);rnrn if (inflowBoundary==velocity) {rn sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);rn }rn else {rn //sBoundaryCondition.addPressureBoundary(superGeometry, 3, omega);rn sBoundaryCondition.addConvectionBoundary(superGeometry, 3, omega); rn }rnrn if (outflowBoundary==velocity) {rn sBoundaryCondition.addVelocityBoundary(superGeometry, 4, omega);rn }rn else {rn //sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);rn sBoundaryCondition.addConvectionBoundary(superGeometry, 4, omega);rn }rnrn /// Initial conditionsrn /*rn T Lx = converter.numCells(lx);rn T Ly = converter.numCells(ly);rn T p0 =8.*converter.getLatticeNu()*converter.getLatticeU()*Lx/(Ly*Ly);rn AnalyticalLinear2D<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2 , 0 , p0*DESCRIPTOR<T>::invCs2+1 );rn */rn AnalyticalConst2D<T,T> rho(1.);rn AnalyticalConst2D<T,T> u_zero(0.);rn rn T maxVelocity = converter.getLatticeU()*(Re_initial/Re);rn T distance2Wall = converter.getLatticeL();rn //////////////////////////////////////////////////////////////////////rn std::vector<T> axisPoint = superGeometry.getStatistics().getMaxPhysR(3);rn axisPoint[1] += T(distance2Wall); rn std::vector<int> discreteNormal = superGeometry.getStatistics().computeDiscreteNormal(3);rn std::vector<T> axisDirection;rn axisDirection.push_back((T)(discreteNormal[0]));rn axisDirection.push_back((T)(discreteNormal[1]));rn T radius = T(distance2Wall);rn for (int iD = 0; iD < 2; iD++) {rn radius += ((superGeometry.getStatistics().getPhysRadius(3)[iD])*(T)2.);rn }rn Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);rn //////////////////////////////////////////////////////////////////////rn //Poiseuille2D<T> u(superGeometry, 3, maxVelocity, distance2Wall);rnrn // Initialize all values of distribution functions to their local equilibriumrn sLattice.defineRhoU(superGeometry, 1, rho, u);rn sLattice.iniEquilibrium(superGeometry, 1, rho, u);rn //sLattice.defineRhoU(superGeometry, 1, rho, u_zero);rn //sLattice.iniEquilibrium(superGeometry, 1, rho, u_zero);rn rn sLattice.defineRhoU(superGeometry, 2, rho, u);rn sLattice.iniEquilibrium(superGeometry, 2, rho, u);rn //sLattice.defineRhoU(superGeometry, 2, rho, u_zero);rn //sLattice.iniEquilibrium(superGeometry, 2, rho, u_zero);rn rn sLattice.defineRhoU(superGeometry, 3, rho, u);rn sLattice.iniEquilibrium(superGeometry, 3, rho, u);rn rn sLattice.defineRhoU(superGeometry, 4, rho, u);rn sLattice.iniEquilibrium(superGeometry, 4, rho, u);rn //sLattice.defineRhoU(superGeometry, 4, rho, u_zero);rn //sLattice.iniEquilibrium(superGeometry, 4, rho, u_zero);rn rn //////////////////////////////////////////////////////////////////////rn // Initialize avShearrn std::vector<T> avShearVec(1,T());rn AnalyticalConst2D<T,T> avShear(avShearVec);rn sLattice.defineExternalField(superGeometry, 1,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 2,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 3,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 4,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn sLattice.defineExternalField(superGeometry, 5,rn DESCRIPTOR<T>::ExternalField::avShearIsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfAvShear, avShear );rn //////////////////////////////////////////////////////////////////////rn rn /// Initial conditions force termrn std::vector<T> extend(2,T());rn std::vector<T> origin(2,T());rn extend[0] = lx;rn extend[1] = ly;rn IndicatorCuboid2D<T> geo(extend, origin);rnrn //AnalyticalConst2D<T,T> force_x(-9.81/converter.physMasslessForce()*sin(atan(-1./1700.)));rn //AnalyticalConst2D<T,T> force_y(-9.81/converter.physMasslessForce()*cos(atan(-1./1700.)));rn //AnalyticalConst2D<T,T> force_x(-9.81/converter.physMasslessForce()*sin(atan(0)));rn //AnalyticalConst2D<T,T> force_y(-9.81/converter.physMasslessForce()*cos(atan(0)));rn AnalyticalConst2D<T,T> force_x(0.);rn AnalyticalConst2D<T,T> force_y(0.);rn AnalyticalComposed2D<T,T> force(force_x,force_y);rnrn // Initialize forcern /*rn sLattice.defineExternalField(superGeometry, 1,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.defineExternalField(superGeometry, 2,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.defineExternalField(superGeometry, 3,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn sLattice.defineExternalField(superGeometry, 4,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn */rn sLattice.defineExternalField(superGeometry, geo,rn DESCRIPTOR<T>::ExternalField::forceBeginsAt,rn DESCRIPTOR<T>::ExternalField::sizeOfForce, force );rn rn //////////////////////////////////////////////////////////////////////rn rn /// Make the lattice ready for simulationrn sLattice.initialize();rnrn clout << “”Prepare Lattice … OK”” << std::endl;rn}rnrn/// Compute error normsrnvoid error(SuperGeometry2D<T>& superGeometry,rn SuperLattice2D<T, DESCRIPTOR>& sLattice,rn LBconverter<T> const& converter,rn Dynamics<T, DESCRIPTOR>& bulkDynamics) {rnrn OstreamManager clout(std::cout,””error””);rnrn int tmp[]={int()};rn T result[2]={T(),T()}, result_tmp[2]={T(),T()};rn T result1;rnrn const T maxVelocity = converter.physVelocity(converter.getLatticeU());rn T distance2Wall = converter.getLatticeL();rn Poiseuille2D<T> uSol(superGeometry, 3, maxVelocity, distance2Wall);rn SuperLatticePhysVelocity2D<T,DESCRIPTOR> u(sLattice,converter);rn SuperLatticeFfromAnalyticalF2D<T,DESCRIPTOR> uSolLattice(uSol,sLattice,superGeometry);rnrn T Lx = converter.numCells(lx);rn T Ly = converter.numCells(ly);rn T p0 =8.*converter.getLatticeNu()*converter.getLatticeU()*Lx/(Ly*Ly);rn AnalyticalLinear2D<T,T> pressureSol(-converter.physPressureFromRho(p0*DESCRIPTOR<T>::invCs2+1)/lx , 0 , converter.physPressureFromRho(p0*DESCRIPTOR<T>::invCs2+1) );rn SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(sLattice,converter);rn SuperLatticeFfromAnalyticalF2D<T,DESCRIPTOR> pressureSolLattice(pressureSol,sLattice,superGeometry);rnrn PoiseuilleStrainRate2D<T,T> sSol(converter, ly);rn SuperLatticePhysStrainRate2D<T,DESCRIPTOR> s(sLattice,converter);rn SuperLatticeFfromAnalyticalF2D<T,DESCRIPTOR> sSolLattice(sSol,sLattice,superGeometry);rnrnrn // velocity errorrn SuperL1Norm2D<T,DESCRIPTOR> uL1Norm(uSolLattice-u,superGeometry,1);rn SuperL1Norm2D<T,DESCRIPTOR> uSolL1Norm(uSolLattice,superGeometry,1);rn uL1Norm(result,tmp);rn uSolL1Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”velocity-L1-error(abs)=”” << result[0] << “”; velocity-L1-error(rel)=”” << result1 << std::endl;rnrn SuperL2Norm2D<T,DESCRIPTOR> uL2Norm(uSolLattice-u,superGeometry,1);rn SuperL2Norm2D<T,DESCRIPTOR> uSolL2Norm(uSolLattice,superGeometry,1);rn uL2Norm(result,tmp);rn uSolL2Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”velocity-L2-error(abs)=”” << result[0] << “”; velocity-L2-error(rel)=”” << result1 << std::endl;rnrn SuperLinfNorm2D<T,DESCRIPTOR> uLinfNorm(uSolLattice-u,superGeometry,1);rn uLinfNorm(&result1,tmp);rn clout << “”velocity-Linf-error(abs)=”” << result1 << std::endl;rnrn // density errorrn SuperL1Norm2D<T,DESCRIPTOR> pressureL1Norm(pressureSolLattice-pressure,superGeometry,1);rn SuperL1Norm2D<T,DESCRIPTOR> pressureSolL1Norm(pressureSolLattice,superGeometry,1);rn pressureL1Norm(result,tmp);rn pressureSolL1Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”pressure-L1-error(abs)=”” << result[0] << “”; pressure-L1-error(rel)=”” << result1 << std::endl;rnrn SuperL2Norm2D<T,DESCRIPTOR> pressureL2Norm(pressureSolLattice-pressure,superGeometry,1);rn SuperL2Norm2D<T,DESCRIPTOR> pressureSolL2Norm(pressureSolLattice,superGeometry,1);rn pressureL2Norm(result,tmp);rn pressureSolL2Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”pressure-L2-error(abs)=”” << result[0] << “”; pressure-L2-error(rel)=”” << result1 << std::endl;rnrn SuperLinfNorm2D<T,DESCRIPTOR> pressureLinfNorm(pressureSolLattice-pressure,superGeometry,1);rn pressureLinfNorm(&result1,tmp);rn clout << “”pressure-Linf-error(abs)=”” << result1 << std::endl;rnrn // strain rate errorrn SuperL1Norm2D<T,DESCRIPTOR> sL1Norm(sSolLattice-s,superGeometry,1);rn SuperL1Norm2D<T,DESCRIPTOR> sSolL1Norm(sSolLattice,superGeometry,1);rn sL1Norm(result,tmp);rn sSolL1Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”strainRate-L1-error(abs)=”” << result[0] << “”; strainRate-L1-error(rel)=”” << result1 << std::endl;rnrn SuperL2Norm2D<T,DESCRIPTOR> sL2Norm(sSolLattice-s,superGeometry,1);rn SuperL2Norm2D<T,DESCRIPTOR> sSolL2Norm(sSolLattice,superGeometry,1);rn sL2Norm(result,tmp);rn sSolL2Norm(result_tmp,tmp);rn result1=result[0]/result_tmp[0];rn clout << “”strainRate-L2-error(abs)=”” << result[0] << “”; strainRate-L2-error(rel)=”” << result1 << std::endl;rnrn SuperLinfNorm2D<T,DESCRIPTOR> sLinfNorm(sSolLattice-s,superGeometry,1);rn sLinfNorm(&result1,tmp);rn clout << “”strainRate-Linf-error(abs)=”” << result1 << std::endl;rn}rnrn/// Output to console and filesrnvoid getResults(SuperLattice2D<T,DESCRIPTOR>& sLattice, Dynamics<T, DESCRIPTOR>& bulkDynamics,rn LBconverter<T>& converter, int iT,rn SuperGeometry2D<T>& superGeometry, Timer<T>& timer) {rnrn OstreamManager clout(std::cout,””getResults””);rnrn SuperVTKwriter2D<T> vtkWriter(“”Channel2d””);rn SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);rn SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);rn vtkWriter.addFunctor( velocity );rn vtkWriter.addFunctor( pressure );rnrn const int vtkIter = converter.numTimeSteps(maxPhysT/100.);rn const int statIter = converter.numTimeSteps(maxPhysT/100.);rn const int saveIter = converter.numTimeSteps(maxPhysT/100.);rnrn if (iT==0) {rn /// Writes the converter log filern writeLogFile(converter, “”Channel2d””);rnrn /// Writes the geometry, cuboid no. and rank no. as vti file for visualizationrn SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, superGeometry);rn SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLattice);rn SuperLatticeRank2D<T, DESCRIPTOR> rank(sLattice);rn superGeometry.rename(0,2);rn vtkWriter.write(geometry);rn vtkWriter.write(cuboid);rn vtkWriter.write(rank);rnrn vtkWriter.createMasterFile();rn }rnrn /// Writes the vtk files and profile text filern if (iT%vtkIter==0) {rn vtkWriter.write(iT);rnrn SuperEuklidNorm2D<T, DESCRIPTOR> normVel(velocity);rn BlockLatticeReduction2D<T, DESCRIPTOR> planeReduction(normVel);rn BlockGifWriter<T> gifWriter;rn gifWriter.write(planeReduction, iT, “”vel””); // scaledrnrn ofstream *ofile = 0;rn if (singleton::mpi().isMainProcessor()) {rn ofile = new ofstream((singleton::directories().getLogOutDir()+””centerVel.dat””).c_str());rn }rn T Ly = converter.numCells(ly);rn for (int iY=0; iY<=Ly; ++iY) {rn /*rn T dx = converter.getDeltaX();rn const T maxVelocity = converter.physVelocity(converter.getLatticeU());rn T point[2]={T(),T()};rn point[0] = lx/2.;rn point[1] = (T)iY/Ly;rn T distance2Wall = converter.getLatticeL();rn Poiseuille2D<T> uSol(superGeometry, 3, maxVelocity, distance2Wall);rn T analytical[2] = {T(),T()};rn uSol(analytical,point);rn SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);rn AnalyticalFfromSuperLatticeF2D<T, DESCRIPTOR> intpolateVelocity(velocity, true);rn T numerical[2] = {T(),T()};rn intpolateVelocity(numerical,point);rn if (singleton::mpi().isMainProcessor()) {rn *ofile << iY*dx << “” “” << analytical[0]rn << “” “” << numerical[0] << “”n””;rn }rn */rn T dx = converter.getDeltaX();rn T point[2]={T(),T()};rn point[0] = lx/2.;rn point[1] = (T)iY/Ly;rn SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);rn AnalyticalFfromSuperLatticeF2D<T, DESCRIPTOR> intpolateVelocity(velocity, true);rn T numerical[2] = {T(),T()};rn intpolateVelocity(numerical,point);rn if (singleton::mpi().isMainProcessor()) {rn *ofile << iY*dx << “” “” << numerical[0] << “”n””;rn }rn }rn delete ofile;rn }rnrn /// Writes output on the consolern if (iT%statIter==0) {rn /// Timer console outputrn timer.update(iT);rn timer.printStep();rnrn /// Lattice statistics console outputrn sLattice.getStatistics().print(iT,converter.physTime(iT));rn rn /*rn /// Error normsrn error(superGeometry, sLattice, converter, bulkDynamics);rn */ rn }rn rn /// Saves lattice datarn if (iT%saveIter==0 && iT>0) {rn clout << “”Checkpointing the system at t=”” << iT << endl;rn //sLattice.save(“”Channel2d.checkpoint””);rn std::string file_name = createFileName(“”Channel2d””,iT);rn file_name = file_name + “”.checkpoint””;rn sLattice.save(file_name);rn }rn rn /// Load lattice datarn //if (iT==0) {rn // // The data can be reloaded usingrn // sLattice.load(“”channel2dtest.checkpoint””);rn //}rn}rnrnint main(int argc, char* argv[]) {rnrn /// === 1st Step: Initialization ===rn olbInit(&argc, &argv);rn singleton::directories().setOutputDir(“”./tmp/””);rn OstreamManager clout(std::cout,””main””);rn // display messages from every single mpi processrn //clout.setMultiOutput(true);rnrn const BoundaryType inflowBoundary = velocity;rn const BoundaryType outflowBoundary = pressure;rn rn /*rn LBconverter<T> converter(rn (int) 2, // dimrn 1./N, // latticeL_rn 1./M, // latticeU_rn (T) 1./Re, // charNu_rn (T) 1. // charL_ = 1,rn );rn */rn rn LBconverter<T> converter(rn (int) 2, // dimrn 1./N, // latticeLrn 1./M, // latticeUrn (T) 0.000001, // charNurn (T) 0.4, // charLrn (T) Re*0.000001/0.4, // charUrn (T) 998.2, // charRhorn (T) 101325.0 // pressureLevelrn );rn rn converter.print();rnrn /// === 2nd Step: Prepare Geometry ===rn std::vector<T> extend(2,T());rn extend[0] = lx;rn extend[1] = ly;rn std::vector<T> origin(2,T());rn IndicatorCuboid2D<T> cuboid(extend, origin);rnrn /// Instantiation of a cuboidGeometry with weightsrn#ifdef PARALLEL_MODE_MPIrn const int noOfCuboids = singleton::mpi().getSize();rn#elsern const int noOfCuboids = 7;rn#endifrn CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getLatticeL(), noOfCuboids);rnrn /// Instantiation of a loadBalancerrn HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);rnrn /// Instantiation of a superGeometryrn SuperGeometry2D<T> superGeometry(cuboidGeometry, loadBalancer, 2);rnrn prepareGeometry(converter, superGeometry);rnrn /// === 3rd Step: Prepare Lattice ===rn SuperLattice2D<T, DESCRIPTOR> sLattice(superGeometry);rn rn Dynamics<T, DESCRIPTOR>* bulkDynamics;rnrn#if defined(Laminaire)rn bulkDynamics = new ForcedBGKdynamics<T, DESCRIPTOR> (converter.getOmega(), rn instances::getBulkMomenta<T,DESCRIPTOR>() );rn rn#else //ShearSmagorinskyWallrn bulkDynamics = new ShearSmagorinskyForcedBGKdynamics<T, DESCRIPTOR> (converter.getOmega(), rn instances::getBulkMomenta<T,DESCRIPTOR>(), 0.1, converter.getLatticeL(), rn converter.physTime() );rn#endifrnrn // choose between local and non-local boundary conditionrn sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryCondition(sLattice);rn createInterpBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);rn // createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryCondition);rnrn prepareLattice(converter, sLattice, *bulkDynamics, sBoundaryCondition, inflowBoundary, outflowBoundary, superGeometry);rnrn /// === 4th Step: Main Loop with Timer ===rn clout << “”starting simulation…”” << endl;rn Timer<T> timer(converter.numTimeSteps(maxPhysT), superGeometry.getStatistics().getNvoxel() );rn timer.start();rnrn for (int iT = 0; iT < converter.numTimeSteps(maxPhysT); ++iT) {rnrn /// === 5th Step: Definition of Initial and Boundary Conditions ===rn // in this application no boundary conditions have to be adjustedrnrn /// === 6th Step: Collide and Stream Execution ===rn sLattice.collideAndStream();rnrn /// === 7th Step: Computation and Output of the Results ===rn getResults(sLattice, *bulkDynamics, converter, iT, superGeometry, timer);rn }rnrn timer.stop();rn timer.printSummary();rn}rn

    rnrnDoes anyone know the reasons which explain this behavior? rnrnBest regards,rnrnAlejandrorn

    July 20, 2016 at 10:15 am #2405
    mathias
    Keymaster

    Dear Alejandro,rnrn There are several issues to be discussed for that topic. Further, we will have a spring school in March 2017 for LBM and OpenLB, where we can help to get started. Please contact me to proceed.rnrnBestrnMathias

  • Author
    Posts
Viewing 2 posts - 1 through 2 (of 2 total)
  • You must be logged in to reply to this topic.