Hi All,
I switched to Palabos for porous media simulation. Using the example provided by Wim, I created my own code. I am trying to simulate pressure driven flow with no flow boundary conditions on the side walls and bounceback dynamics when particle nodes. It compiles well but I get nan as my average energy and density. Below is a copy of the code. Help will be appreciated.
#include “palabos3D.h”
#ifndef PLB_PRECOMPILED // Unless precompiled version is used,
#include “palabos3D.hh” // include full template code
#endif
#include
#include
#include
#include
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor
T pgradient = 1.e-6;
T tol = 1.e-6;
T val1 = 0;
T simhasconverged(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, IncomprFlowParam const& parameters)
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
T val2 = computeAverage( *computeVelocityNorm(lattice, Box3D(nx-1, nx-1, 0, ny-1, 0, nz-1)) );
T error = fabs( 1 - val1/val2);
val1 = val2;
return error;
}
T densitygradient(plint iX, IncomprFlowParam const& parameters){
T Lx = parameters.getNx()-1;
//T delp = deltaP/L * (T)iX * DESCRIPTOR::invCs2;
return 1. - pgradient * (T)iX * DESCRIPTOR::invCs2;
}
template
class Pressuregradient {
public:
Pressuregradient(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
T operator() (plint iX, plint iY, plint iZ) const {
return densitygradient(iX, parameters);
}
private:
IncomprFlowParam parameters;
};
template
class PressureandZeroVelocity {
public:
PressureandZeroVelocity(IncomprFlowParam parameters_)
: parameters(parameters_)
{ }
void operator() (plint iX, plint iY, plint iZ, T& rhox, Array<T,3>& vel) const {
rhox = densitygradient(iX, parameters);
vel[0] = T();
vel[1] = T();
vel[2] = T();
}
private:
IncomprFlowParam parameters;
};
void porousmediadomain( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>& boundaryCondition,
MultiScalarField3D& mediadomain)
{
const plint nx = lattice.getNx();
const plint ny = lattice.getNy();
const plint nz = lattice.getNz();
pcout << "Definition of inlet/outlet." << endl;
Box3D inlet(0,0, 1,ny-2, 1,nz-2);
// boundaryCondition->addPressureBoundary0N(inlet, lattice);
boundaryCondition.setPressureConditionOnBlockBoundaries (lattice, inlet );
setBoundaryDensity(lattice, inlet, 1.);
Box3D outlet(nx-1,nx-1, 1,ny-2, 1,nz-2);
// boundaryCondition->addPressureBoundary0P(outlet, lattice);
boundaryCondition.setPressureConditionOnBlockBoundaries (lattice, outlet );
setBoundaryDensity(lattice, outlet, 1. - ( pgradient * (T)(nx-1) * DESCRIPTOR<T>::invCs2 ) );
//set velocity boundaries on the other sides
//1.=>
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(0, nx-1, 0, ny-1, 0, 0) );
setBoundaryVelocity (
lattice, Box3D(0, nx-1, 0, ny-1, 0, 0),
Array<T,3>(0., 0., 0.) );
//2.=>
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(0, nx-1, 0, ny-1, nz-1, nz-1) );
setBoundaryVelocity (
lattice, Box3D(0, nx-1, 0, ny-1, nz-1, nz-1),
Array<T,3>(0., 0., 0.) );
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(1, nx-2, 0, 0, 1, nz-2) );
setBoundaryVelocity (
lattice, Box3D(1, nx-2, 0, 0, 1, nz-2),
Array<T,3>(0., 0., 0.) );
//4.=>
boundaryCondition.setVelocityConditionOnBlockBoundaries (
lattice, Box3D(1, nx-2, ny-1, ny-1, 1, nz-2) );
setBoundaryVelocity (
lattice, Box3D(1, nx-2, ny-1, ny-1, 1, nz-2),
Array<T,3>(0., 0., 0.) );
//set value of density on those boundaries
setBoundaryDensity (lattice, lattice.getBoundingBox(), Pressuregradient<T>(parameters) );
//set the density to be varying
//set the
pcout << "Definition of the geometry." << endl;
// Where "geometry" evaluates to 1, use bounce-back.
defineDynamics(lattice, mediadomain, new BounceBack<T,DESCRIPTOR>(), true);
pcout << "Initilization of rho and u." << endl;
initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureandZeroVelocity<T>(parameters) );
lattice.initialize();
}
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
IncomprFlowParam<T> parameters(
(T) 1e-2, // uMax
(T) 100., // Re
40, // N
1., // lx
1., // ly
1. // lz
);
const T logT = (T)1/(T)100;
const T imSave = (T)1/(T)40;
const T vtkSave = (T)1/40;
const T maxT = (T)10.1;
writeLogFile(parameters, "Porous media");
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();
const plint nz = parameters.getNz();
const T omega = parameters.getOmega();
pcout << "Creation of the lattice." << endl;
MultiBlockLattice3D<T, DESCRIPTOR> lattice (nx, ny, nz, new BGKdynamics<T, DESCRIPTOR>(omega));
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition
= createLocalBoundaryCondition3D<T,DESCRIPTOR>();
// Switch off periodicity.
lattice.periodicity().toggleAll(false);
pcout << "Reading the geometry file." << endl;
MultiScalarField3D<T> mediadomain(nx, ny, nz);
plb_ifstream fileinput("sim.txt");
if(!fileinput.is_open()) {
pcout << "Error: could not open geometry file " << endl;
return -1;
}
fileinput >> mediadomain;
porousmediadomain(lattice, parameters, *boundaryCondition, mediadomain);
// define convergence
pcout << "Simulation begins" << endl;
plint iT=0;
for(iT = 0; iT<parameters.nStep(maxT); ++iT){
//do {
if (iT%parameters.nStep(logT)==0) {
pcout << "step " << iT
<< "; t=" << iT*parameters.getDeltaT()
<< "; av energy="
<< getStoredAverageEnergy<T>(lattice)
<< "; av rho="
<< getStoredAverageDensity<T>(lattice) << endl << flush;
}
if (iT%parameters.nStep(imSave)==0 && iT>0) {
pcout << "Saving Gif ..." << endl << flush;
writeGifs(lattice, parameters, iT);
}
if (iT%parameters.nStep(vtkSave)==0 && iT>0) {
pcout << "Saving VTK file ..." << endl << flush;
writeVTK(lattice, parameters, iT);
}
lattice.collideAndStream();
if( simhasconverged(lattice, parameters) < tol)
break;
}
averagevelocity(lattice, parameters, iT);
delete boundaryCondition;
}