# Problem coupling Multi-component Lattices

Hi there,
I recently started using Palabos to calculate the permeability of a fluid through porous medium. To do this I simply took the Permeability tutorial from the documentation and modified it slightly. I ended up getting reasonable results for single component/phase flow.
Then I modified the code to take 2 lattices (1 for each fluid) and attempted to couple them such that I could simulate multi-component flow. For some reason the fluids are just not coupling. I used a similar approach to rayleighTaylor3D (in showCases folder).
Running this code does give answers for permeability. However, the program is just running two separate single phase/component simulations (1 for each fluid). Could someone please suggest how I could modify my code in order to make the lattices couple? Thank you.

[code=“cpp”]
/* This file is part of the Palabos library.
*

• Route d’Oron 2
• 1010 Lausanne, Switzerland
• E-mail contact: contact@flowkit.com
• http://www.palabos.org/
• The library Palabos is free software: you can redistribute it and/or
• modify it under the terms of the GNU Affero General Public License as
• published by the Free Software Foundation, either version 3 of the
• The library is distributed in the hope that it will be useful,
• but WITHOUT ANY WARRANTY; without even the implied warranty of
• MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
• GNU Affero General Public License for more details.
• You should have received a copy of the GNU Affero General Public License
• along with this program. If not, see http://www.gnu.org/licenses/.
*/

/* Main author: Wim Degruyter */

#include “palabos3D.h”
#include “palabos3D.hh”
#include
#include

using namespace plb;
using namespace std;

typedef double T;
// EDITED
#define DESCRIPTOR descriptors::ShanChenD3Q19Descriptor

// This function object returns a zero velocity, and a pressure which decreases
// linearly in x-direction. It is used to initialize the particle populations.
public:
PressureGradient(T deltaP_, plint nx_) : deltaP(deltaP_), nx(nx_)
{ }
void operator() (plint iX, plint iY, plint iZ, T& density, Array<T,3>& velocity) const
{
velocity.resetToZero();
density = 1. - deltaP*DESCRIPTOR::invCs2 / (T)(nx-1) * (T)iX;

``````    }
``````

private:
T deltaP;
plint nx;
};

void porousMediaSetup( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition,
MultiScalarField3D& geometry, T deltaP)
{
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);
setBoundaryDensity(lattice, inlet, 1.);

Box3D outlet(nx-1,nx-1, 1,ny-2, 1,nz-2);
setBoundaryDensity(lattice, outlet, 1. - deltaP*DESCRIPTOR<T>::invCs2);

pcout << "Definition of the geometry." << endl;
// Where "geometry" evaluates to 1, use bounce-back.
defineDynamics(lattice, geometry, new BounceBack<T,DESCRIPTOR>(), 1);
// Where "geometry" evaluates to 2, use no-dynamics (which does nothing).
defineDynamics(lattice, geometry, new NoDynamics<T,DESCRIPTOR>(), 2);

pcout << "Initilization of rho and u." << endl;
initializeAtEquilibrium( lattice, lattice.getBoundingBox(), PressureGradient(deltaP, nx) );

lattice.initialize();
delete boundaryCondition;
``````

}

void writeGifs(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint nx = lattice.getNx();
const plint ny = lattice.getNy();
const plint nz = lattice.getNz();

``````    const plint imSize = 600;
ImageWriter<T> imageWriter("leeloo");

// Write velocity-norm at x=0.
imageWriter.writeScaledGif( createFileName("ux_inlet", iter, 6),
*computeVelocityNorm(lattice, Box3D(0,0, 0,ny-1, 0,nz-1)),
imSize, imSize );

// Write velocity-norm at x=nx/2.
imageWriter.writeScaledGif( createFileName("ux_half", iter, 6),
*computeVelocityNorm(lattice, Box3D(nx/2,nx/2, 0,ny-1, 0,nz-1)),
imSize, imSize );
``````

}

void writeVTK(MultiBlockLattice3D<T,DESCRIPTOR>& lattice, plint iter)
{
VtkImageOutput3D vtkOut(createFileName(“vtk”, iter, 6), 1.);
vtkOut.writeData(*computeVelocityNorm(lattice), “velocityNorm”, 1.);
vtkOut.writeData<3,float>(*computeVelocity(lattice), “velocity”, 1.);
}

T computePermeability (
MultiBlockLattice3D<T,DESCRIPTOR>& lattice, T nu, T deltaP, Box3D domain )
{
pcout << “Computing the permeability.” << endl;

``````    // Compute only the x-direction of the velocity (direction of the flow).
plint xComponent = 0;
plint nx = lattice.getNx();

T meanU = computeAverage (
*computeVelocityComponent (lattice, domain, xComponent ) );

pcout << "Average velocity     = " << meanU             << endl;
pcout << "Lattice viscosity nu = " << nu                << endl;
pcout << "Grad P               = " << deltaP/(T)(nx-1)            << endl;
pcout << "Permeability         = " << nu*meanU / (deltaP/(T)(nx-1)) << endl;

return meanU;
``````

}

int main(int argc, char **argv)
{
plbInit(&argc, &argv);

``````    if (argc!=7)
{
pcout << "Error missing some input parameter\n";
pcout << "The structure is :\n";
pcout << "1. Input file name.\n";
pcout << "2. Output directory name.\n";
pcout << "3. number of cells in X direction.\n";
pcout << "4. number of cells in Y direction.\n";
pcout << "5. number of cells in Z direction.\n";
pcout << "6. Delta P .\n";
pcout << "Example: " << argv << " twoSpheres.dat tmp/ 48 64 64 0.00005\n";
exit (EXIT_FAILURE);
}
std::string fNameIn  = argv;
std::string fNameOut = argv;

const plint nx = atoi(argv);
const plint ny = atoi(argv);
const plint nz = atoi(argv);
const T   deltaP = atof(argv);
const T   G = 1.2;

global::directories().setOutputDir(fNameOut+"/");

const T omega1 = 1.0;
const T omega2 = 1.9;
const T nu1    = ((T)1/omega1-0.5)/DESCRIPTOR<T>::invCs2;
const T nu2    = ((T)1/omega2-0.5)/DESCRIPTOR<T>::invCs2;

pcout << "Creation of the lattices." << endl;
MultiBlockLattice3D<T,DESCRIPTOR> lattice1(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega1));
MultiBlockLattice3D<T,DESCRIPTOR> lattice2(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega2));
// Switch off periodicity.
lattice1.periodicity().toggleAll(false);
lattice2.periodicity().toggleAll(false);

///////

vector<MultiBlockLattice3D<T, DESCRIPTOR>* > blockLattices;
blockLattices.push_back(&lattice1);
blockLattices.push_back(&lattice2);

///////

std::vector<T> constOmegaValues;
constOmegaValues.push_back(omega1);
constOmegaValues.push_back(omega2);
plint processorLevel = 1;
integrateProcessingFunctional (
new ShanChenMultiComponentProcessor3D<T,DESCRIPTOR>(G,constOmegaValues),
Box3D(0,nx-1,1,ny-2,0,nz-1),
blockLattices, processorLevel );

///////

MultiScalarField3D<int> geometry(nx,ny,nz);
MultiScalarField3D<int> slice(1,ny,nz);
for (plint iX=0; iX<nx-1; ++iX) {
string fname = createFileName("slice_", iX, 4)+"_truc.dat";
pcout << "Reading slice " << fname;
plb_ifstream geometryFile(fNameIn.c_str());
if (!geometryFile.is_open()) {
pcout << "Error: could not open geometry file " << fNameIn << endl;
return -1;
}
geometryFile >> slice;
copy(slice, slice.getBoundingBox(), geometry, Box3D(iX,iX,0,ny-1,0,nz-1));
}

////////

pcout << "Reading the geometry file." << endl;
//MultiScalarField3D<int> geometry(nx,ny,nz);
plb_ifstream geometryFile(fNameIn.c_str());
if (!geometryFile.is_open()) {
pcout << "Error: could not open geometry file " << fNameIn << endl;
return -1;
}
geometryFile >> geometry;

pcout << "nu = " << nu1 << " " << nu2 << endl;
pcout << "deltaP = " << deltaP << endl;
pcout << "omega = " << omega1 << " " << omega2 << endl;
pcout << "nx = " << lattice1.getNx() << " " << lattice2.getNx() << endl;
pcout << "ny = " << lattice1.getNy() << " " << lattice2.getNy() << endl;
pcout << "nz = " << lattice1.getNz() << " " << lattice2.getNz() << endl;

porousMediaSetup(lattice2, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);
porousMediaSetup(lattice1, createLocalBoundaryCondition3D<T,DESCRIPTOR>(), geometry, deltaP);

// The value-tracer is used to stop the simulation once is has converged.
// 1st parameter:velocity
// 2nd parameter:size
// 3rd parameters:threshold
// 1st and second parameters are used for the length of the time average (size/velocity)
util::ValueTracer<T> converge1(1.0,1000.0,1.0e-4);
util::ValueTracer<T> converge2(1.0,1000.0,1.0e-4);

pcout << "Simulation begins" << endl;
plint iT=0;

const plint maxT = 20;
for (;iT<maxT; ++iT)
{
if (iT % 1 == 0) {
pcout << "Iteration " << iT << endl;
}
if (iT % 50 == 0 && iT>0) {
writeGifs(lattice1,iT);
}

lattice1.collideAndStream();
lattice2.collideAndStream();

converge1.takeValue(getStoredAverageEnergy(lattice1),true);
converge2.takeValue(getStoredAverageEnergy(lattice2),true);

if ((converge1.hasConverged())&&(converge2.hasConverged()))
{
break;
}
}

pcout << "End of simulation at iteration " << iT << endl;

pcout << "Permeability of First fluid:" << endl << endl;
computePermeability(lattice1, nu1, deltaP, lattice1.getBoundingBox());
pcout << endl;

pcout << "Permeability of Second fluid:" << endl << endl;
computePermeability(lattice2, nu2, deltaP, lattice2.getBoundingBox());
pcout << endl;

pcout << "Writing VTK file ..." << endl << endl;
writeVTK(lattice1,iT);
pcout << "Finished!" << endl << endl;
``````

}

Hi,

Hopefully you’ve already solved the problem.
But if not, I had a quick look at your code, and I think the key thing here is that in you code

[code=“cpp”]
MultiBlockLattice3D<T,DESCRIPTOR> lattice1(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega1));
MultiBlockLattice3D<T,DESCRIPTOR> lattice2(nx,ny,nz, new BGKdynamics<T,DESCRIPTOR>(omega2));

``````

Instead of defining "BGKdynamics", you should have "ExternalMomentBGKdynamics" or "ExternalMomentRegularizedBGKdynamics" for better stability.

Yeah I was using that "BGKdynamics" before and couldn't get it right. I guess that "ExternalMoment" means some extra coupling of shan-chen force to calculate equilibrium velocity (if this reminds you how the shan-chen model works :-P)

Hope it helps !
rex``````

Hello, mlak. I’m trying to do a simulation to calculate the permeability of a 3D multicomponent flow, have you found a solution? may i know how to combine the multicomponent3d code into the permeability code on palabos?