Reply To: Multiphase flow in 3D using FreeEnergy LBM

[mehrdadvasheghanifarahani]

Dear Sam,

Many thanks for your response. Sure, I’d be happy to share the code:

#include “olb3D.h”
#include “olb3D.hh” // use only generic version!

#include <cstdlib>
#include <iostream>
#include <fstream>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace std;

typedef double T;
#define DESCRIPTOR D3Q19<CHEM_POTENTIAL,FORCE>

// Parameters for the simulation setup
const int N = 100;
const T nx = 20.;
const T ny = 100.;
const T nz = 100.;
const T dx = ny / N;

const T inlet1Velocity = 0.00056; // [lattice units]
const T outletDensity = 1.; // [lattice units]
const T alpha = 1.; // Interfacial width [lattice units]

const T kappa1 = 0.001;
const T kappa2 = 0.001;
const T kappa3 = 0.001;

const T gama = 1.; // For mobility of interfaces [lattice units]
const T h1 = 0.0;
const T h2 = 0.0;
const T h3 = 0.0;

const int maxIter = 1000000;
const int vtkIter = 1000;
const int statIter = 2000;

void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry3D<T>& superGeometry ) {
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

superGeometry.rename(0, 2);
superGeometry.rename(2, 1, 1, 1, 1);

// Inlets and outlet

Vector<T,3> origin = superGeometry.getStatistics().getMinPhysR( 2 );
origin[1] += converter.getConversionFactorLength()/2.;
origin[2] += converter.getConversionFactorLength()/2.;

Vector<T,3> extend = superGeometry.getStatistics().getMaxPhysR( 2 );
extend[1] = extend[1]-origin[1]-converter.getConversionFactorLength()/2.;
extend[2] = extend[2]-origin[2]-converter.getConversionFactorLength()/2.;

// Set material number for inflow
origin[0] = superGeometry.getStatistics().getMinPhysR( 2 )[0]-converter.getConversionFactorLength();
extend[0] = 2*converter.getConversionFactorLength();
IndicatorCuboid3D<T> inflow( extend,origin );
superGeometry.rename( 2,3,1,inflow );

// Set material number for outflow
origin[0] = superGeometry.getStatistics().getMaxPhysR( 2 )[0]-converter.getConversionFactorLength();
extend[0] = 2*converter.getConversionFactorLength();
IndicatorCuboid3D<T> outflow( extend,origin );
superGeometry.rename( 2,8,1,outflow );

superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();

clout << “Prepare Geometry … OK” << std::endl;
}

void prepareLattice( SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
Dynamics<T, DESCRIPTOR>& bulkDynamics2,
Dynamics<T, DESCRIPTOR>& bulkDynamics3,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC1,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC2,
sOnLatticeBoundaryCondition3D<T,DESCRIPTOR>& sOnBC3,
UnitConverter<T, DESCRIPTOR>& converter,
SuperGeometry3D<T>& superGeometry ) {

OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

// define lattice dynamics
sLattice1.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 1, &bulkDynamics1 );
sLattice2.defineDynamics( superGeometry, 1, &bulkDynamics2 );
sLattice3.defineDynamics( superGeometry, 1, &bulkDynamics3 );

sLattice1.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 2, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 3, &instances::getNoDynamics<T, DESCRIPTOR>() );

sLattice1.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice2.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );
sLattice3.defineDynamics( superGeometry, 8, &instances::getNoDynamics<T, DESCRIPTOR>() );

// add wall boundary
sOnBC1.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 1 );
sOnBC2.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 2 );
sOnBC3.addFreeEnergyWallBoundary( superGeometry, 2, alpha, kappa1, kappa2, kappa3, h1, h2, h3, 3 );

// add inlet boundaries
T omega = converter.getLatticeRelaxationFrequency();
auto inlet1Indicator = superGeometry.getMaterialIndicator(3);
sOnBC1.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 1 );
sOnBC2.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 2 );
sOnBC3.addFreeEnergyInletBoundary( inlet1Indicator, omega, “velocity”, 3 );

// add outlet boundary
auto outletIndicator = superGeometry.getMaterialIndicator(8);
sOnBC1.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 1 );
sOnBC2.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 2 );
sOnBC3.addFreeEnergyOutletBoundary( outletIndicator, omega, “density”, 3 );

// bulk initial conditions
std::vector<T> v( 3,T() );
AnalyticalConst3D<T,T> zeroVelocity( v );

AnalyticalConst3D<T,T> zero ( 0. );
AnalyticalConst3D<T,T> one ( 1. );
SmoothIndicatorCuboid3D<T,T> section1( {nx * 0.25, ny/2., nz/2.}, nx * 0.5, ny, nz, 0.);
SmoothIndicatorCuboid3D<T,T> section2( {nx * 0.65, ny/2., nz/2.}, nx * 0.3, ny, nz, 0.);

AnalyticalIdentity3D<T,T> c1( section1 );
AnalyticalIdentity3D<T,T> c2( section2 );
AnalyticalIdentity3D<T,T> rho( one );
AnalyticalIdentity3D<T,T> phi( c1 – c2 );
AnalyticalIdentity3D<T,T> psi( rho – c1 – c2 );

auto allIndicator = superGeometry.getMaterialIndicator({1, 2, 3});
sLattice1.iniEquilibrium( allIndicator, rho, zeroVelocity );
sLattice2.iniEquilibrium( allIndicator, phi, zeroVelocity );
sLattice3.iniEquilibrium( allIndicator, psi, zeroVelocity );

// inlet boundary conditions

std::vector<T> maxVelocity( 3,0 );
maxVelocity[0] = 1.5*inlet1Velocity;

T distance2Wall = converter.getConversionFactorLength()/2.;
RectanglePoiseuille3D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall, distance2Wall, distance2Wall );

sLattice1.defineU( inlet1Indicator, poiseuilleU);
sLattice2.defineRho( inlet1Indicator, phi);
sLattice3.defineRho( inlet1Indicator, psi);

// outlet initial / boundary conditions
AnalyticalConst3D<T,T> rhoOutlet( outletDensity );
AnalyticalIdentity3D<T,T> phiOutlet( zero );
AnalyticalIdentity3D<T,T> psiOutlet( rhoOutlet );

sLattice1.defineRho( outletIndicator, rhoOutlet );
sLattice2.defineRho( outletIndicator, phiOutlet );
sLattice3.defineRho( outletIndicator, psiOutlet );

// initialise lattices
sLattice1.initialize();
sLattice2.initialize();
sLattice3.initialize();

sLattice1.communicate();
sLattice2.communicate();
sLattice3.communicate();

clout << “Prepare Lattice … OK” << std::endl;
}

void prepareCoupling(SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3,
SuperGeometry3D<T>& superGeometry) {
OstreamManager clout( std::cout,”prepareCoupling” );
clout << “Add lattice coupling” << endl;

// Bulk couplings
FreeEnergyChemicalPotentialGenerator3D<T,DESCRIPTOR> coupling2( alpha, kappa1, kappa2, kappa3 );
FreeEnergyForceGenerator3D<T,DESCRIPTOR> coupling3;

// Inlet / outlet couplings
FreeEnergyDensityOutletGenerator3D<T,DESCRIPTOR> coupling1( outletDensity );
FreeEnergyInletOutletGenerator3D<T,DESCRIPTOR> coupling4;

// Suppress compiler warnings
coupling2.shift(0, 0, 0);
coupling3.shift(0, 0, 0);
coupling1.shift(0, 0, 0);
coupling4.shift(0, 0, 0);

// The DensityOutlet coupling must be added to the first lattice and come before the ChemicalPotential coupling
// The InletOutlet couplings must be added to the second lattice and come after the Force coupling.
sLattice1.addLatticeCoupling<DESCRIPTOR>( superGeometry, 8, coupling1, {&sLattice2, &sLattice3} );

sLattice1.addLatticeCoupling<DESCRIPTOR>( superGeometry, 1, coupling2, {&sLattice2, &sLattice3} );
sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 1, coupling3, {&sLattice1, &sLattice3} );

sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 3, coupling4, {&sLattice1, &sLattice3} );

sLattice2.addLatticeCoupling<DESCRIPTOR>( superGeometry, 8, coupling4, {&sLattice1, &sLattice3} );

clout << “Add lattice coupling … OK!” << endl;
}

void getResults( SuperLattice3D<T, DESCRIPTOR>& sLattice1,
SuperLattice3D<T, DESCRIPTOR>& sLattice2,
SuperLattice3D<T, DESCRIPTOR>& sLattice3, int iT,
SuperGeometry3D<T>& superGeometry, Timer<T>& timer,
UnitConverter<T, DESCRIPTOR> converter) {

OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter3D<T> vtmWriter( “microFluidics2d” );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry( sLattice1, superGeometry );
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid( sLattice1 );
SuperLatticeRank3D<T, DESCRIPTOR> rank( sLattice1 );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}

// Get statistics
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();
sLattice1.getStatistics().print( iT, converter.getPhysTime(iT) );
sLattice2.getStatistics().print( iT, converter.getPhysTime(iT) );
sLattice3.getStatistics().print( iT, converter.getPhysTime(iT) );
}

// Writes the VTK files
if ( iT%vtkIter==0 ) {
SuperLatticeVelocity3D<T, DESCRIPTOR> velocity( sLattice1 );
SuperLatticeDensity3D<T, DESCRIPTOR> density1( sLattice1 );
density1.getName() = “rho”;
SuperLatticeDensity3D<T, DESCRIPTOR> density2( sLattice2 );
density2.getName() = “phi”;
SuperLatticeDensity3D<T, DESCRIPTOR> density3( sLattice3 );
density3.getName() = “density-fluid-3”;

AnalyticalConst3D<T,T> half_( 0.5 );
SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR> half(half_, sLattice1);

SuperIdentity3D<T,T> c1 (half*(density1+density2-density3));
c1.getName() = “density-fluid-1”;
SuperIdentity3D<T,T> c2 (half*(density1-density2-density3));
c2.getName() = “density-fluid-2”;

vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( density1 );
vtmWriter.addFunctor( density2 );
vtmWriter.addFunctor( density3 );
vtmWriter.addFunctor( c1 );
vtmWriter.addFunctor( c2 );
vtmWriter.write( iT );
}
}

int main( int argc, char *argv[] ) {

// === 1st Step: Initialization ===

olbInit( &argc, &argv );
singleton::directories().setOutputDir( “./tmp/” );
OstreamManager clout( std::cout,”main” );

UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR> converter(
(T) N, // resolution
(T) 1., // lattice relaxation time (tau)
(T) nz, // charPhysLength: reference length of simulation geometry
(T) 1.e-6, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 0.1, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) 1. // physDensity: physical density in __kg / m^3__
);

// Prints the converter log as console output
converter.print();

// === 2nd Step: Prepare Geometry ===
std::vector<T> extend = {nx, ny, nz};
std::vector<T> origin = {0, 0, 0};
IndicatorCuboid3D<T> cuboid(extend,origin);
#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cGeometry( cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize() );
#else
CuboidGeometry3D<T> cGeometry( cuboid, converter.getPhysDeltaX() );
#endif

// Instantiation of loadbalancer
HeuristicLoadBalancer<T> loadBalancer( cGeometry );
loadBalancer.print();

// Instantiation of superGeometry
SuperGeometry3D<T> superGeometry( cGeometry,loadBalancer );

prepareGeometry( converter, superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice3D<T, DESCRIPTOR> sLattice1( superGeometry );
SuperLattice3D<T, DESCRIPTOR> sLattice2( superGeometry );
SuperLattice3D<T, DESCRIPTOR> sLattice3( superGeometry );

ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
converter.getLatticeRelaxationFrequency(),
instances::getBulkMomenta<T,DESCRIPTOR>() );

FreeEnergyBGKdynamics<T, DESCRIPTOR> bulkDynamics23 (
converter.getLatticeRelaxationFrequency(), gama,
instances::getBulkMomenta<T,DESCRIPTOR>() );

sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC1( sLattice1 );
sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC2( sLattice2 );
sOnLatticeBoundaryCondition3D<T, DESCRIPTOR> sOnBC3( sLattice3 );
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC1);
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC2);
createLocalBoundaryCondition3D<T, DESCRIPTOR> (sOnBC3);

prepareLattice( sLattice1, sLattice2, sLattice3, bulkDynamics1, bulkDynamics23,
bulkDynamics23, sOnBC1, sOnBC2, sOnBC3, converter, superGeometry );

prepareCoupling( sLattice1, sLattice2, sLattice3, superGeometry );

SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal1 (superGeometry, sLattice1, sLattice1.getOverlap() );
SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal2 (superGeometry, sLattice2, sLattice2.getOverlap() );
SuperExternal3D<T,DESCRIPTOR,CHEM_POTENTIAL> sExternal3 (superGeometry, sLattice3, sLattice3.getOverlap() );

// === 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 ) {
// Computation and output of the results
getResults( sLattice1, sLattice2, sLattice3, iT, superGeometry, timer, converter );

// Collide and stream execution
sLattice1.collideAndStream();
sLattice2.collideAndStream();
sLattice3.collideAndStream();

// MPI communication for lattice data
sLattice1.communicate();
sLattice2.communicate();
sLattice3.communicate();

// Execute coupling between the two lattices
sLattice1.executeCoupling();
sExternal1.communicate();
sExternal2.communicate();
sExternal3.communicate();
sLattice2.executeCoupling();
}

timer.stop();
timer.printSummary();

}

As I mentioned, my results show that “phi” is monotonously increasing in the outlet boundary and I believe this might cause the issue.

Kind Regards,
Mehrdad