Re: Pressure Boundary Condition
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Pressure Boundary Condition › Re: Pressure Boundary Condition
Okay, so here is one of the combination I’ve been trying recently with which I simply get a lot of “nan” results. The bold lines are the ones related to the pressure boundary condition I am trying to include and without them the simulation works.
#include “olb2D.hh”
#include <cstdlib>
#include <iostream>
using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace std;
typedef double T;
#define DESCRIPTOR ShanChenDynOmegaForcedD2Q9Descriptor
// Parameters for the simulation setup
const int nx = 100;
const int ny = 250;
const T centerDropletX = (nx-1)/2;
const T centerDropletY = 20.;
const T radiusDroplet = 10.;
const T wallThickness =4.;
const int maxIter = 6900;
/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry2D<T>& superGeometry) {
OstreamManager clout(std::cout,”prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;
// Sets material number for fluid and boundary
superGeometry.rename(0,1);
std::vector<T> origin1(2,T());
origin1[0] = 0;
origin1[1] = 0;
std::vector<T> origin2(2,T());
origin2[0] = 0;
origin2[1] = ny-1;
std::vector<T> origin3(2,T());
origin3[0] = 0;
origin3[1] = wallThickness;
std::vector<T> origin4(2,T());
origin4[0] = nx-wallThickness-1;
origin4[1] = wallThickness;
std::vector<T> extend1(2,T());
extend1[0] = nx;
extend1[1] = wallThickness;
std::vector<T> extend2(2,T());
extend2[0] = wallThickness;
extend2[1] = ny-wallThickness;
std::vector<T>extend3(2,T());
extend3[0] = nx;
extend3[1] = 1;
Vector<T,2> centerDroplet(centerDropletX,centerDropletY);
IndicatorCuboid2D<T> bottom(extend1, origin1);
IndicatorCuboid2D<T> top(extend3, origin2);
IndicatorCuboid2D<T> leftwall(extend2, origin3);
IndicatorCuboid2D<T> rightwall(extend2, origin4);
IndicatorCircle2D<T> droplet(centerDroplet, radiusDroplet);
superGeometry.rename(1,2,droplet);
superGeometry.rename(1,3,bottom);
superGeometry.rename(1,3,leftwall);
superGeometry.rename(1,3,rightwall);
superGeometry.rename(1,4,top);
/// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << “Prepare Geometry … OK” << std::endl;
}
void prepareLattice(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
Dynamics<T, DESCRIPTOR>& bulkDynamics2,
Dynamics<T, DESCRIPTOR>& bounceBackRho0,
Dynamics<T, DESCRIPTOR>& bounceBackRho1,
[b] sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeOne,
sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeOne,
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeTwo,
sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeTwo,[/b]
SuperGeometry2D<T>& superGeometry)
{
OstreamManager clout(std::cout,”prepareLattice”);
clout << “Prepare Lattice …” << std::endl;
[b] const T omega1 = 1.0;
const T omega2 = 1.0;[/b]
/// define lattice Dynamics
sLatticeOne.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 3, &bounceBackRho1);
sLatticeOne.defineDynamics(superGeometry, 4, &bulkDynamics1);
[b] sBoundaryConditionLatticeOne.addPressureBoundary(superGeometry, 4, omega1);[/b]
sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 3, &bounceBackRho1);
sLatticeTwo.defineDynamics(superGeometry, 4, &bulkDynamics2);
[b]
sBoundaryConditionLatticeTwo.addPressureBoundary(superGeometry, 4, omega2);[/b]
clout << “Prepare Lattice … OK” << std::endl;
}
void setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
T force, int iT, SuperGeometry2D<T>& superGeometry) {
if(iT==0) {
AnalyticalConst2D<T,T> noise(4.e-2);
std::vector<T> v(2,T());
AnalyticalConst2D<T,T> zeroV(v);
AnalyticalConst2D<T,T> zero(0.);
AnalyticalLinear2D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.98+force*ny*DESCRIPTOR<T>::invCs2);
AnalyticalConst2D<T,T> onePlus(0.98+force*ny/2.*DESCRIPTOR<T>::invCs2);
AnalyticalRandom2D<T,T> random;
AnalyticalIdentity2D<T,T> randomOne(random*noise+one);
AnalyticalIdentity2D<T,T> randomPlus(random*noise+onePlus);
std::vector<T> F(2,T());
F[1] = -force;
AnalyticalConst2D<T,T> f(F);
/// for each material set the defineRhou and the Equilibrium
sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 1,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );
sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
sLatticeOne.defineExternalField(superGeometry, 2,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );
sLatticeOne.defineRhoU(superGeometry, 4, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 4, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 4,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );
sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);
sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);
sLatticeTwo.defineRhoU(superGeometry, 4, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 4, randomPlus, zeroV);
/// Make the lattice ready for simulation
sLatticeOne.initialize();
sLatticeTwo.initialize();
}
}
void getResults(SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
SuperLattice2D<T, DESCRIPTOR>& sLatticeOne, int iT,
SuperGeometry2D<T>& superGeometry, Timer<T>& timer) {
OstreamManager clout(std::cout,”getResults”);
SuperVTKwriter2D<T> vtkWriter(“rayleighTaylor2dsLatticeOne”);
const int vtkIter = 100;
const int statIter = 10;
if (iT==0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLatticeOne, superGeometry);
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLatticeOne);
SuperLatticeRank2D<T, DESCRIPTOR> rank(sLatticeOne);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}
// Get statistics
if (iT%statIter==0 && iT > 0) {
// Timer console output
timer.update( iT );
timer.printStep();
clout << “averageRhoFluidOne=” << sLatticeOne.getStatistics().getAverageRho();
clout << “; averageRhoFluidTwo=” << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
}
/// Writes the VTK files
if (iT%vtkIter==0) {
clout << “Writing VTK …” << std::endl;
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity(sLatticeOne);
SuperLatticeDensity2D<T, DESCRIPTOR> density(sLatticeOne);
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( density );
vtkWriter.write(iT);
BlockLatticeReduction2D<T, DESCRIPTOR> planeReduction(density);
BlockGifWriter<T> gifWriter;
gifWriter.write(planeReduction, iT, “density”);
clout << “Writing VTK … OK” << std::endl;
}
}
int main(int argc, char *argv[]) {
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout,”main”);
const T omega1 = 1.0;
const T omega2 = 1.0;
const T G = 3.;
T force = -30./(T)ny/(T)ny;
/// === 2nd Step: Prepare Geometry ===
/// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, singleton::mpi().getSize());
#else
CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, 1);
#endif
cGeometry.setPeriodicity(false, false);
HeuristicLoadBalancer<T> loadBalancer(cGeometry);
SuperGeometry2D<T> superGeometry(cGeometry, loadBalancer, 2);
prepareGeometry(superGeometry);
/// === 3rd Step: Prepare Lattice ===
SuperLattice2D<T, DESCRIPTOR> sLatticeOne(superGeometry);
SuperLattice2D<T, DESCRIPTOR> sLatticeTwo(superGeometry);
ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
omega1, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );
ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics2 (
omega2, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );
// A bounce-back node with fictitious density 1,
// which is experienced by the partner fluid
BounceBack<T, DESCRIPTOR> bounceBackRho1( 1. );
// A bounce-back node with fictitious density 0,
// which is experienced by the partner fluid
BounceBack<T, DESCRIPTOR> bounceBackRho0( 0. );
std::vector<T> rho0;
rho0.push_back(1);
rho0.push_back(1);
PsiEqualsRho<T,T> interactionPotential;
ShanChenForcedGenerator2D<T,DESCRIPTOR> coupling(G,rho0,interactionPotential);
[b] // Pressure boundary
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeOne(sLatticeOne);
createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeOne);
sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeOne(sLatticeOne);
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeTwo(sLatticeTwo);
createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeTwo);
sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeTwo(sLatticeTwo);[/b]
sLatticeOne.addLatticeCoupling(superGeometry, 1, coupling, sLatticeTwo);
sLatticeOne.addLatticeCoupling(superGeometry, 2, coupling, sLatticeTwo);
prepareLattice(sLatticeOne, sLatticeTwo, bulkDynamics1, bulkDynamics2,
bounceBackRho0, bounceBackRho1,[b] sBoundaryConditionLatticeOne, sOffBoundaryConditionLatticeOne, sBoundaryConditionLatticeTwo, sOffBoundaryConditionLatticeTwo,[/b] superGeometry);
/// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << “starting simulation…” << endl;
Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();
for (iT=0; iT<maxIter; ++iT) {
/// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(sLatticeOne, sLatticeTwo, force, iT, superGeometry);
/// === 6th Step: Collide and Stream Execution ===
sLatticeOne.collideAndStream();
sLatticeTwo.collideAndStream();
sLatticeOne.communicate();
sLatticeTwo.communicate();
sLatticeOne.executeCoupling();
/// === 7th Step: Computation and Output of the Results ===
getResults(sLatticeTwo, sLatticeOne, iT, superGeometry, timer);
}
timer.stop();
timer.printSummary();
}