Issues with palabos

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.) );

//3.=>

	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;

}

Hi,

did you try to reduce a bit the reynolds number (or increase the size of your domain)?

Well, for porous media simulation, the Reynolds number does not influence the simulation. It is only used to fix the omega. It has no relationship to the pressure/density difference. I worked like this in OpenLB and it was good. But don’t know why palabos is going this way.

Yes I agree. But here from your parameters omega=1.94.

Which is a quite high value (close to 2). Maybe reducing it could help.