Currently, as far as i can tell the default Shan-Chen data processor uses the Shan-Chen force implementation. I would like to use the Guo implementation instead.
Does this mean i need to modify the data processor or is there already an implementation for the Guo force which can be used?
Dear nlooije,
I know this thread is old, but have you worked on this problem any further? I wanted to try out the same thing
and adjusted the single component data processor to work with the guo forcing scheme (see code below).
I have initialized a static droplet with the densities obtained from the mechanical stability condition.
The simulation is stable, but the results are wrong. At equilibrium, I end up with the following liquid
and vapor densities:
GuoForce: Vapor: 3.578450e-01 Liquid: 2.467872e+00
ShanChenForce: Vapor: 1.590207e-01 Liquid: 1.953080e+00
Obviously something is wrong, but I can’t figure out what. Have you (or anybody else) encountered similar problems
or see where I went wrong?
Modified Shan-Chen Processor:
[code=“cpp”]
template<typename T, template class Descriptor>
void singleComponentWithGuoScheme2D<T,Descriptor>::process ( Box2D domain, BlockLattice2D<T,Descriptor>& lattice ) {
// Short-hand notation for the lattice descriptor
typedef Descriptor<T> D;
// Handle to external scalars
enum {
densityOffset = D::ExternalField::densityBeginsAt,
forceOffset = D::ExternalField::forceBeginsAt
} ;
// Include Boundary layer
plint nx = domain.getNx() + 2*D::vicinity ;
plint ny = domain.getNy() + 2*D::vicinity ;
plint offsetX = domain.x0 - D::vicinity ;
plint offsetY = domain.y0 - D::vicinity ;
// Scalarfield for values of potential function
ScalarField2D<T> psiField(nx,ny);
// Calculate moments
for (plint iX=domain.x0-D::vicinity; iX<=domain.x1+D::vicinity; ++iX) {
for (plint iY=domain.y0-D::vicinity; iY<=domain.y1+D::vicinity; ++iY) {
// Get cell
Cell<T,Descriptor>& cell = lattice.get(iX,iY);
// Get density
T rho = cell.computeDensity();
// Get psi
psiField.get(iX-offsetX, iY-offsetY) = Psi->compute(rho) ;
// Store externally
*cell.getExternal(densityOffset) = rho ;
}
}
// Calculate force
for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
Array<T,D::d> rhoContribution;
rhoContribution.resetToZero();
// Calculate finite difference
for (plint iPop = 0; iPop < D::q; ++iPop) {
// Get neighbour
plint nextX = iX + D::c[iPop][0] ;
plint nextY = iY + D::c[iPop][1] ;
// Get psi
T psi = psiField.get(nextX-offsetX, nextY-offsetY) ;
// Get finite difference
for (plint iD = 0; iD < D::d; ++iD) {
rhoContribution[iD] += D::t[iPop] * psi * D::c[iPop][iD] ;
}
}
// Access external force
Cell<T,Descriptor>& cell = lattice.get(iX,iY);
T *force = cell.getExternal(forceOffset) ;
force[0] = (T) 0 ;
force[1] = (T) 0 ;
for (plint iD = 0; iD < D::d; ++iD) {
// Add interaction term.
T psi = psiField.get(iX-offsetX, iY-offsetY);
// Store force vector in external field
force[iD] -= G * psi * rhoContribution[iD] ;
}
}
}
}
And my main:
[code="cpp"]
#include "palabos2D.h"
#include "palabos2D.hh"
#include "palabosAuxiliary.h"
#include "math.h"
#include <cstdlib>
#include <iostream>
using namespace plb ;
using namespace std ;
typedef double T ;
#define DESCRIPTOR descriptors::ForcedShanChenD2Q9Descriptor
#define DYNAMICS GuoExternalForceBGKdynamics
//#define DYNAMICS ExternalMomentBGKdynamics
// Initializer
template<typename T, template<typename U> class Descriptor>
class staticDropletInitializer2D : public OneCellIndexedFunctional2D<T,Descriptor> {
public :
// Constructor
staticDropletInitializer2D ( T rhoInside_, T rhoOutside_, plint rad_, plint xCenter_, plint yCenter_)
: rhoInside(rhoInside_), rhoOutside(rhoOutside_), rad(rad_), xCenter(xCenter_), yCenter(yCenter_) { } ;
// Execute
virtual void execute(plint iX, plint iY, Cell<T,Descriptor>& cell) const {
// Zero-velocity condition
Array<T,Descriptor<T>::d> zeroVelocity ((T)0,(T)0) ;
T insideCircle = sqrt(pow(((T)(iX-xCenter)),(T)2) + pow(((T)(iY-yCenter)),(T)2)) ;
if(insideCircle<=rad){
iniCellAtEquilibrium(cell,rhoInside,zeroVelocity) ;
}
else{
iniCellAtEquilibrium(cell,rhoOutside,zeroVelocity) ;
}
} ;
// Clone
virtual staticDropletInitializer2D<T,Descriptor>* clone() const {
return new staticDropletInitializer2D<T,Descriptor>(*this) ;
} ;
private :
T rhoInside ;
T rhoOutside ;
plint rad ;
plint xCenter ;
plint yCenter ;
} ;
int main(int argc, char *argv[]){
plbInit(&argc, &argv);
global::directories().setOutputDir("./results-test/") ;
// Read console parameters
const T tau = (T)1 ;
const T omega = (T)1/tau ;
// Parameters
const T rhoVapor = 1.564130e-01 ;
const T rhoLiquid = 1.932442e+00 ;
// Constant for potential function
const T rho0 = 1. ;
const T G = -5. ;
// set output format
pcout << scientific ;
// Problem description
plint nx = 100 ;
plint ny = 100 ;
plint radius = 25 ;
plint centerX = nx/2 ;
plint centerY = ny/2 ;
// Parameters and variables for the convergence analysis
plint maxIter = 20001 ;
plint statIter = 100 ;
plint saveIter = 100 ;
// Initialize lattice
MultiBlockLattice2D<T,DESCRIPTOR> lattice (nx,ny,new DYNAMICS<T,DESCRIPTOR>(omega)) ;
// Initialize all boundaries to be periodic
lattice.periodicity().toggleAll(true) ;
// Define processor level
const plint processorLevel = 1;
// Apply data processor
// Shan-Chen Processor with Guo force
integrateProcessingFunctional ( new singleComponentWithGuoScheme2D<T,DESCRIPTOR> (G, new interparticlePotential::PsiShanChen93<T>(rho0)) ,
lattice.getBoundingBox() ,
lattice ,
processorLevel ) ;
// Shan-Chen Processor with Shan-Chen force
// integrateProcessingFunctional ( new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> ( G, new interparticlePotential::PsiShanChen93<T>(rho0) ),
// lattice.getBoundingBox(),
// lattice,
// processorLevel ) ;
// Setup problem
applyIndexed(lattice,lattice.getBoundingBox(), new staticDropletInitializer2D<T,DESCRIPTOR>(rhoLiquid,rhoVapor,radius,centerX,centerY) ) ;
// Execute actual initialization
lattice.initialize();
// Simulate
for(plint i = 0 ; i<maxIter ; i++){
// Lattice Boltzmann iteration
lattice.collideAndStream() ;
if (i%statIter==0) {
auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
pcout << i << ":\t Minimum density: " << computeMin(*rho) ;
pcout << "\t Maximum density: " << computeMax(*rho) << endl ;
}
if (i%saveIter == 0) {
ImageWriter<T>("leeloo").writeScaledGif ( createFileName("rho", i, 6), *computeDensity(lattice) );
}
}
}
Regards,
kk
Hi kk,
I was a bit annoyed that no replies came to this post, it seems to me that something like this should already have been implemented for some time.
I went ahead and implemented the Guo scheme along with some other things (Yuan-Schaefer interaction potentials, Lycett-Brown forcing, etc.), you can find them at my github repo. There is an example code included with which you can compare your results. Feel free to fork 
I don’t know why your code is not performing as you expect and unfortunately don’t have time to go into it at great depth but what i noticed was you don’t actually give the Guo code in ‘singleComponentWithGuoScheme2D’. Are you sure it is implemented correctly? Compare it to ‘GuoExternalMomentBGKdynamics’
nlooije
Hi,
thank you for the quick reply. Some more EOS were also the first thing I added.
To my understanding, when using the ‘GuoExternalMomentBGKdynamics’ the force that is stored in the corresponding external field is applied via the guo forcing scheme. Therefore, I only calculated the force in the shan chen data processor, stored it in said field and used the ‘GuoExternalMomentBGKdynamics’ (defined at the top of the main-file). Maybe I am getting this wrong (obviously something bad happens).
Thanks for your time and I will look into your implementation.
Regards,
kk
No problem, let me know if you find the bug (or if you don’t)…
It would be nice if we don’t reinvent the wheel so feel free to send a pull request if you make any edits to the files
especially if you make any improvements to them
nlooije