Re: Pressure Boundary Condition

[ivan]

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.

Code:

#include “olb2D.h”
#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();
}