Parallel Execution of WellBalancedCahnHilliardBGKdynamics with solid boundaries

[tmhys]

Dear OpenLB Community,

I am trying to run a simulation of a liquid bridge between two (sand) particles. Eventually I would like to scale this up to some hundred particles so I would like to run simulations in parallel. I am using Ver 1.8.1.
Basically my starting point is examples/multiComponent/contactAngle2d.cpp. I added two circles using indicators and put a liquid droplet between the circles. (using Nx=320, Ny= 400,diameter=140)
// —-
IndicatorCircle2D<T> particle( {Nx/2., 240.}, diameter/2. );
superGeometry.rename( 1, 2, particle );
IndicatorCircle2D<T> particle2( {Nx/2., 80.}, diameter/2. );
superGeometry.rename( 1, 2, particle2 );
// —
This worked as I expected for the simulation with a single core – it formed a reasonable liquid bridge between the two circles but if I run in parallel, there is a sharp increase in the pressure and velocity at the intersection between “boundary of the surface of circle” and “MPI boundary”. Please find the screenshots from Paraview.
https://drive.google.com/drive/folders/1xAtUWIiaHUuKOnAxFOCadkV_IEOoC0O8?usp=share_link

I confirmed that the original contactAngle2d.cpp worked without problems in parallel but this might be because it does not have an oblique solid boundary crossing the MPI boundaries. I thought this is a problem related to WellBalancedCahnHilliardPostProcessor – if the BOUNDARY field is not communicated, there would be a problem in the gradient calculation in mu and phi. However, this problem actually persists even though I explicitly requested the BOUNDARY field to be communicated.
//–
{
auto& communicator = sLatticeCH.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(2);
communicator.requestField<STATISTIC>();
communicator.requestField<PHIWETTING>();
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestField<BOUNDARY>(); // Added
communicator.exchangeRequests();
}
//–

Please find the full script below. Any input is really appreciated.

//—
#include <olb.h>

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

using T = FLOATING_POINT_TYPE;
using NSDESCRIPTOR = D2Q9<RHO,NABLARHO,FORCE,EXTERNAL_FORCE,TAU_EFF,STATISTIC>;
using CHDESCRIPTOR = D2Q9<FORCE,SOURCE,SOURCE_OLD,PHIWETTING,VELOCITY,OLD_PHIU,STATISTIC,CHEM_POTENTIAL,BOUNDARY>;
using NSBulkDynamics = MultiPhaseIncompressbileBGKdynamics<T,NSDESCRIPTOR>;
using CHBulkDynamics = WellBalancedCahnHilliardBGKdynamics<T,CHDESCRIPTOR>;
using Coupling = WellBalancedCahnHilliardPostProcessor;

// Parameters for the simulation domain
const int Nx = 320; // domain resolution x [lattice units]
const int Ny = 400; // domain resolution y [lattice units]
int diameter = 140; // droplet diameter [lattice units]
const int maxIter = 5000000; // number of iterations to perform
const int vtkIter = 1000; // interval to save vtk output
const int statIter = 1000; // interval to print statistics

// Characteristic physical parameters
const T L_char = 70e-6*4000; // charPhysLength [physical units]
const T DeltaRho = 1000.; // conversion factor density [physical units]
const T viscosityH2O = 9e-7; // physViscosity H2O liquid [physical units]

// Lattice parameters for fluid properties
const T tau_mobil = 5.5; // relaxation time for interface mobility in Cahn-Hilliard equation (mobility=(tau_mobil-0.5)/3) [lattice units]
const T tau_l = 1.0; // relaxation time H2O liquid lattice (kinematic viscosity=(tau_l-0.5)/3) [lattice units]
const T tau_g = 0.52; // relaxation time gas lattice (kinematic viscosity=(tau_g-0.5)/3) [lattice units]
const T sigma = 0.01; // liquid-gas surface tension [lattice units]
const T w = 4.; // diffuse interface width [lattice units]
const T g = 0*9.81; // gravitational force magnitude [lattice units]
const std::vector<T> rhos = {1., 1.}; // densities of the {gas, liquid} [lattice units]
const T theta = 10.; // equilibrium contact angle [degrees]

// Labels for boundary and fluid locations
void prepareGeometry( SuperGeometry<T,2>& superGeometry )
{
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;
// Fluid nodes labelled 2
superGeometry.rename( 0,2 );
// Label edges as 1
superGeometry.rename( 2, 1, {0, 1} );

// Insert circle
IndicatorCircle2D<T> particle( {Nx/2., 240.}, diameter/2. );
superGeometry.rename( 1, 2, particle );
IndicatorCircle2D<T> particle2( {Nx/2., 80.}, diameter/2. );
superGeometry.rename( 1, 2, particle2 );

//IndicatorCircle2D<T> sourceCircle( {Nx/2., 160.}, diameter/2. );
//superGeometry.rename( 1, 3, sourceCircle );

superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
clout << “Prepare Geometry … OK” << std::endl;
}

template <typename SuperLatticeCoupling>
void prepareLattice( SuperLattice<T,NSDESCRIPTOR>& sLatticeNS,
SuperLattice<T,CHDESCRIPTOR>& sLatticeCH,
SuperLatticeCoupling& coupling,
UnitConverter<T,NSDESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry, int diameter )
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

// lattice dynamics for Navier-Stokes equation, only needed in fluid domain with mn 1
sLatticeNS.defineDynamics<NoDynamics>(superGeometry, 0);
sLatticeNS.defineDynamics<NSBulkDynamics>(superGeometry, 1);
sLatticeNS.defineDynamics<NoDynamics>(superGeometry, 2);

// lattice dynamics for Cahn-Hilliard equation, only needed in fluid domain with mn 1
sLatticeCH.defineDynamics<NoDynamics>(superGeometry, 0);
sLatticeCH.defineDynamics<CHBulkDynamics>(superGeometry, 1);
sLatticeCH.defineDynamics<NoDynamics>(superGeometry, 2);

// set velocity to zero everywhere initially
Vector<T,2> u(0., 0.);
AnalyticalConst2D<T,T> zeroVelocity( u );

// analytical constants to use in initialisation
AnalyticalConst2D<T,T> one ( 1. );
AnalyticalConst2D<T,T> two ( 2. );
AnalyticalConst2D<T,T> zero ( 0. );
AnalyticalConst2D<T,T> rhov ( rhos[0] );
AnalyticalConst2D<T,T> rhol ( rhos[1] );
AnalyticalConst2D<T,T> tauv ( tau_g );
AnalyticalConst2D<T,T> taul ( tau_l );

// regions with different material ids
auto bulk = superGeometry.getMaterialIndicator( 1 );
//auto bulk = superGeometry.getMaterialIndicator( {1,3} );
auto walls = superGeometry.getMaterialIndicator( 2 );
auto all = superGeometry.getMaterialIndicator( {0,1,2} );

// initialise phi with a circle with centre at (Nx/2,1) and radius (diameter/2), smoothed by w/2
SmoothIndicatorFactoredCircle2D<T,T> circle( {Nx/2., 160.}, diameter/4., w/2, 0, {0,0}, 0, 1. );
AnalyticalIdentity2D<T,T> phi( one – circle );

// initial values for interpolated rho and tau across interfaces
AnalyticalIdentity2D<T,T> rho( rhov + (rhol-rhov)*phi );
AnalyticalIdentity2D<T,T> tau( tauv + (taul-tauv)*phi );

// initial (hydrodynamic) pressure
AnalyticalIdentity2D<T,T> pressure( zero );

// set the initial values for the fields
sLatticeNS.defineField<descriptors::RHO>( all, rho ); // density
sLatticeNS.defineField<descriptors::TAU_EFF>( bulk, tau ); // relaxation time for navier-stokes lattice
sLatticeCH.defineField<descriptors::SOURCE>( bulk, zero ); // source term for Cahn-Hilliard lattice
sLatticeCH.defineField<descriptors::PHIWETTING>( bulk, phi ); // fluid concentration used by wetting boundaries
sLatticeCH.defineField<descriptors::CHEM_POTENTIAL>( bulk, zero ); // chemical potential
sLatticeCH.defineField<descriptors::BOUNDARY>( walls, two ); // boundary field used by wetting boundaries
sLatticeCH.defineField<descriptors::BOUNDARY>( bulk, zero ); // boundary field used by wetting boundaries

// use interpolated bounce back conditions for the distribution function
std::vector<T> origin = { -1.5, 0.5 };
std::vector<T> extend = { Nx+2, Ny-2 };
IndicatorCuboid2D<T> rectangle( extend, origin ); // rectangle spans from (0,0.5) to (Nx-1,Ny-1.5), in the x direction we add padding for the periodic boundaries
setBouzidiBoundary( sLatticeNS, superGeometry, 2, rectangle );
setBouzidiWellBalanced( sLatticeCH, superGeometry, 2, rectangle );

// apply gravitational force
std::vector<T> F( 2,T() );
F[1] = -g/(converter.getPhysDeltaX()/converter.getPhysDeltaT()/converter.getPhysDeltaT());
AnalyticalConst2D<T,T> f( F );
sLatticeNS.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 1, f );

// initialise distribution functions in equilibrium with initial values of phi, velocity and pressure
sLatticeCH.defineRhoU( all, phi, zeroVelocity );
sLatticeCH.iniEquilibrium( all, phi, zeroVelocity );
sLatticeNS.defineRhoU( all, pressure, zeroVelocity );
sLatticeNS.iniEquilibrium( all, pressure, zeroVelocity );

// postprocessor to update rhowetting
sLatticeCH.addPostProcessor<stage::PostStream>(bulk,meta::id<RhoWettingStatistics>());

// give coupling the values of tau_l, tau_g, rho_l and rho_g so that it can update the viscosity and density
coupling.template setParameter<Coupling::TAUS>({tau_l,tau_g});
coupling.template setParameter<Coupling::RHOS>(rhos);

// postprocessor to calculate the chemical potential
sLatticeCH.addPostProcessor<stage::ChemPotCalc>(meta::id<ChemPotentialPhaseFieldProcessor>());

// parameters needed by models
sLatticeNS.setParameter<descriptors::OMEGA>( 1./tau_l ); // collision rate for Navier-Stokes lattice
sLatticeCH.setParameter<descriptors::OMEGA>( 1./tau_mobil ); // collision rate for Cahn-Hilliard lattice
sLatticeCH.setParameter<descriptors::THETA>( (M_PI-theta*M_PI/180.) ); // contact angle
sLatticeCH.setParameter<descriptors::INTERFACE_WIDTH>( w ); // diffuse interface width
sLatticeCH.setParameter<descriptors::SCALAR>( sigma ); // surface tension

// define fields that must be communicated by mpi (needed for finite difference gradients)
{
auto& communicator = sLatticeCH.getCommunicator(stage::PreCoupling());
communicator.requestOverlap(2);
communicator.requestField<STATISTIC>();
communicator.requestField<PHIWETTING>();
communicator.requestField<CHEM_POTENTIAL>();
communicator.requestField<BOUNDARY>(); // Added
communicator.exchangeRequests();
}

// some initial steps
sLatticeCH.executePostProcessors(stage::PreCoupling());
sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();
sLatticeCH.executePostProcessors(stage::ChemPotCalc());
sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();
sLatticeNS.initialize();
sLatticeCH.initialize();
sLatticeCH.iniEquilibrium( all, phi, zeroVelocity );
sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();

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

// for statistic output and vtk saving
T getResults( SuperLattice<T,NSDESCRIPTOR>& sLatticeNS,
SuperLattice<T,CHDESCRIPTOR>& sLatticeCH,
int iT, SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer,
UnitConverter<T,NSDESCRIPTOR> converter )
{
OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter2D<T> vtmWriter( “contactAngle2d” );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid( sLatticeNS );
SuperLatticeRank2D<T, NSDESCRIPTOR> rank( sLatticeNS );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.createMasterFile();
}

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

T contact_angle = 0;
// Writes the VTK files
if ( iT%vtkIter==0 ) {

SuperLatticeDensity2D<T, NSDESCRIPTOR> p_hydro( sLatticeNS );
p_hydro.getName() = “p_hydro”;

SuperLatticeField2D<T, CHDESCRIPTOR, STATISTIC> phi( sLatticeCH );
phi.getName() = “phi”;

SuperLatticeExternalScalarField2D<T, NSDESCRIPTOR, RHO> rho_L( sLatticeNS );
AnalyticalConst2D<T,T> ConversionDensity_( converter.getConversionFactorDensity() );
SuperLatticeFfromAnalyticalF2D<T, NSDESCRIPTOR> ConversionDensity(ConversionDensity_, sLatticeNS);
SuperIdentity2D<T,T> rho ( rho_L * ConversionDensity );
rho.getName() = “rho”;

SuperLatticeVelocity2D<T, NSDESCRIPTOR> velocity( sLatticeNS );
velocity.getName() = “u”;

vtmWriter.addFunctor( p_hydro );
vtmWriter.addFunctor( phi );
vtmWriter.addFunctor( rho );
vtmWriter.addFunctor( velocity );

// Source term:
SuperLatticeField2D<T, CHDESCRIPTOR, SOURCE> sourceField( sLatticeCH );
sourceField.getName() = “source”;
vtmWriter.addFunctor( sourceField );

// Existing phi, chem pot, etc. already created above
// Add PHIWETTING field (CH lattice, descriptor field PHIWETTING maps to tag PHIWETTING)
SuperLatticeField2D<T, CHDESCRIPTOR, PHIWETTING> phiWettingField( sLatticeCH );
phiWettingField.getName() = “phiWetting”;
vtmWriter.addFunctor( phiWettingField );

// Add BOUNDARY field (CH lattice, BOUNDARY tag)
SuperLatticeField2D<T, CHDESCRIPTOR, BOUNDARY> boundaryField( sLatticeCH );
boundaryField.getName() = “boundaryFlag”;
vtmWriter.addFunctor( boundaryField );

// If CHEM_POTENTIAL exists but not yet added, add it similarly:
SuperLatticeField2D<T, CHDESCRIPTOR, CHEM_POTENTIAL> chemPotField( sLatticeCH );
chemPotField.getName() = “chemPot”;
vtmWriter.addFunctor( chemPotField );
vtmWriter.write( iT );

// function to interpolate phi at a given position
AnalyticalFfromSuperF2D<T,T> interpolPhi( phi, true, 1 );
T point1 = 0, point2 = 0, point3 = 0;

// contact angle fitting, will be done by finding three points on the circlular droplet

// point 1, x = ? y = 2
T pos[2] = {0., 2.};
for (int ix=0; ix<0.5*Nx; ix++) {
T phi1, phi2;
pos[0] = ix;
interpolPhi( &phi1, pos );
// if phi at ix is greater than 0.5, we have passed the interface of the liquid phase
if (phi1 < 0.5) {
pos[0] = (ix-1); // check the value at the previous point
interpolPhi( &phi2, pos );
point1 = ix – 1 + (0.5-phi2)/(phi1-phi2); // interpolate between these values
break;
}
}

// point 2, x = ? y = 2
for (int ix=0.5*Nx; ix<Nx; ix++) {
T phi1, phi2;
pos[0] = ix;
interpolPhi( &phi1, pos );
if (phi1 > 0.5) {
pos[0] = (ix-1);
interpolPhi( &phi2, pos );
point2 = ix – 1 + (0.5-phi2)/(phi1-phi2);
break;
}
}

// point 3, x = Nx/2 y = ?
pos[0] = 0.5*(Nx);
for (int iy=3; iy<Ny; iy++) {
T phi1, phi2;
pos[1] = iy;
interpolPhi( &phi1, pos );
if (phi1 > 0.5) {
pos[0] = (iy-1);
interpolPhi( &phi2, pos );
point3 = iy – 1 + (0.5-phi2)/(phi1-phi2);
break;
}
}

T x1=point1;
T y1=2;
T x2=point2;
T y2=2;
T x3=0.5*(Nx);
T y3=point3;

// estimate centre and radius from three points, we must solve three simulatneous equations
T s1 = x1*x1 + y1*y1;
T s2 = x2*x2 + y1*y1;
T s3 = x3*x3 + y3*y3;
T M11 = x1*y2 + x2*y3 + x3*y1 – (x2*y1 + x3*y2 + x1*y3);
T M12 = s1*y2 + s2*y3 + s3*y1 – (s2*y1 + s3*y2 + s1*y3);
T M13 = s1*x2 + s2*x3 + s3*x1 – (s2*x1 + s3*x2 + s1*x3);
T xc = 0.5*M12/M11;
T yc = -0.5*M13/M11;
T r = sqrt(pow(x2 – xc,2) + pow(y2 – yc,2));

contact_angle = util::acos(-(yc-0.5) / r) * 180 / M_PI;

clout << “Numerical Contact angle: ” << contact_angle << std::endl;
clout << “Analytical contact angle: ” << theta << std::endl;
}
return contact_angle;
}

// run the simulation
void simulate( int diameter )
{
OstreamManager clout( std::cout,”main” );
// === 1st Step: Initialization ===
UnitConverterFromResolutionAndRelaxationTime<T,NSDESCRIPTOR> converter(
int {diameter}, // resolution
(T) tau_l, // lattice relaxation time
(T) L_char, // charPhysLength: reference length of simulation geometry
(T) 0, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) viscosityH2O, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T) DeltaRho // 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 };
std::vector<T> origin = { 0., 0. };
IndicatorCuboid2D<T> cuboid(extend,origin);
#ifdef PARALLEL_MODE_MPI
CuboidDecomposition2D<T> cuboidDecomposition( cuboid, 1., singleton::mpi().getSize() );
#else
CuboidDecomposition2D<T> cuboidDecomposition( cuboid, 1. );
#endif

// Set periodic boundaries to the domain
cuboidDecomposition.setPeriodicity({ true, false });

// Instantiation of loadbalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidDecomposition );
loadBalancer.print();
// Instantiation of superGeometry
SuperGeometry<T,2> superGeometry( cuboidDecomposition,loadBalancer );
prepareGeometry( superGeometry );

// === 3rd Step: Prepare Lattice ===
SuperLattice<T,NSDESCRIPTOR> sLatticeNS( superGeometry );
SuperLattice<T,CHDESCRIPTOR> sLatticeCH( superGeometry );
SuperLatticeCoupling coupling(
WellBalancedCahnHilliardPostProcessor{},
names::NavierStokes{}, sLatticeNS,
names::Component1{}, sLatticeCH);
coupling.restrictTo(superGeometry.getMaterialIndicator({1}));
prepareLattice( sLatticeNS, sLatticeCH, coupling, converter, superGeometry, diameter );

// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();
std::ofstream outfile;
outfile.open (“contactAngleVsTime.dat”);
T old_angle = 1.;

for ( iT=0; iT<=maxIter; ++iT ) {
// Collide and stream (and coupling) execution
sLatticeNS.collideAndStream();
sLatticeCH.collideAndStream();

sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();
sLatticeCH.executePostProcessors(stage::PreCoupling());
sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();
sLatticeCH.executePostProcessors(stage::ChemPotCalc());
sLatticeCH.getCommunicator(stage::PreCoupling()).communicate();

coupling.execute();

AnalyticalConst2D<T,T> one ( -1e-5 );
AnalyticalConst2D<T,T> zero ( 0. );
IndicatorCircle2D<T> sourceCircle( {Nx/2., 160.}, 5. );

if (iT > 10000) {
sLatticeCH.defineField<descriptors::SOURCE>( superGeometry, sourceCircle, one );
}

T angle = getResults( sLatticeNS, sLatticeCH, iT, superGeometry, timer, converter );
if ( iT%vtkIter == 0 ) {
outfile << iT*converter.getPhysDeltaT() << “,” ;
outfile << angle << “\n”;
/*
if ( fabs(angle-old_angle)/fabs(old_angle) < 5e-5 ) {
clout << “contact angle converged…” << std::endl;
break;
}
*/
old_angle = angle;
}

if ( std::isnan( sLatticeNS.getStatistics().getAverageEnergy() ) ) {
clout << “Simulation terminated: ” << iT << std::endl;
break;
}
}
outfile.close();
timer.stop();
timer.printSummary();
}

int main( int argc, char *argv[] )
{
initialize( &argc, &argv );
if (argc > 1) {
diameter = atof(argv[1]);
}

singleton::directories().setOutputDir( “./tmp/” );
simulate( diameter );
}

[mathias]

It seems like one of the flields is not updated correctly at the overlab, there might one field which has not sufficient overlab or is missing a communication step..

[Author]

Posts