EDIT I posted the solution below my original text, and removed the original, bugged code.
Hi guys,
I really don’t know how to express my problem better than how’s written on the subject. I got sdea’s sample advection-diffusion code for the advective transport of a component (http://www.palabos.org/forum/read.php?4,8390) and modified it to consider two transported components - my ultimate goal with this is to replicate Duncan & Toor’s experiment with LBM.
Since heat and mass transfer are analogous, where it says “temperature” I’m considering “concentration”. My full code is below the text. Sorry for the huge amount of comments, I’ve been testing a lot of stuff trying to figure out how to solve this problem.
I was able to pinpoint the source of the problem. On the function concentrationSetup (again, concentration is modelled as temperature, due to the mathematical analogy between heat and mass transfer), there’s the definition onf the boundary conditions:
[code="cpp"]
Box2D leftwall(0, 0 , 0,ny-1);
Box2D rightwall(nx ,nx, 0,ny-1);
// Set the boundary-conditions for lattice1
BC.addTemperatureBoundary1N(leftwall, lattice1);
BC.addTemperatureBoundary1P(rightwall, lattice1);
// Set the boundary-conditions for lattice2
BC.addTemperatureBoundary1N(leftwall, lattice2);
BC.addTemperatureBoundary1P(rightwall, lattice2);
The problem lies within the leftwall boundary condition. With this definition, and initializing the concentration (temperature) of 0.5 (green) for each components at the first and last sixth of the channel and 0 (blue) for the rest (reffer to the full code below), this leftwall column of nodes initializes to 1.0 (dark red):
https://preview.ibb.co/fxZq6Q/ADV_conc1_000000.gif
A bit hard to see, but it's the very first column of nodes, on the left.
When I change leftwall a bit to
[code="cpp"]
Box2D leftwall(10, 10 , 0,ny-1);
the red column shifts its position to the tenth column:
https://preview.ibb.co/gKtTt5/ADV_conc1_000000.gif
And when I increase its width, the red column gets a bit fatter too:
[code=“cpp”]
Box2D leftwall(10, 15, 0, ny-1)
https://preview.ibb.co/dNa56Q/ADV_conc1_000000.gif
And this column is present also on the snapshots for the concentration of the second component, represented by lattice2:
https://preview.ibb.co/hkpXmQ/ADV_conc2_000000.gif
This is obviously a problem, because this spurious concentration/temperature gets advected with the rest of the fluid, and I don't want this behaviour! For example, this is the following iteration of the program:
https://preview.ibb.co/kh8Q6Q/ADV_conc1_000100.gif
This is driving me nuts!!
##############################
[b]SOLUTION[/b]
Here's what I did. I starter from "scratch", from sdea's code. Configured two components as I wanted and I still had the fixed-values column problem. Then I decided to try something a little bit different: I commented the boundary condition definition lines, and replaced these definitions with explicit definitions over the region I wanted. And, voilà, problem solved!
However, I am not very comfortable with simply commenting boundary conditions definitions. So this is something I must correct in the future. I set up the advective field with a null velocity, i.e., it's a still fluid, and configured three diffusive components. I then initialized random velocities on each lattice site for these components, to try to mimic the thermodynamic definition of temperature as thermal kinetic energy of the molecules and put this to run, writing the concentration profile of each component (at constant y = ny/2). Below is an animation of the evolution of the concentration profiles, showing that the components do diffuse on the system.
I just need to find how to couple them now.
https://media.giphy.com/media/3o7btQfmlDsSHP1Epa/giphy.gif
The code is:
[code="cpp"]
#include "palabos2D.h"
#include "palabos2D.hh"
#include <iostream>
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>
#include <cstdio>
#include <cstdlib>
using namespace std;
using namespace plb;
typedef double T;
#define NSDES descriptors::D2Q9Descriptor
#define AVDES descriptors::AdvectionDiffusionD2Q5Descriptor
// ************************************************************************
// Part for the NS part of the code. A classic poiseuille flow is simulated.
// Most of the code is taken from Palabos tutorial 1.5
// ************************************************************************
// Poiselle velocity analitical
// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) {
T y = (T)iY / parameters.getResolution();
return 4.*parameters.getLatticeU() * (y-y*y);
}
// A functional, used to initialize the velocity for the boundary conditions
template<typename T>
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(IncomprFlowParam<T> parameters_)
: parameters(parameters_)
{ }
// This version of the operator returns the velocity only,
// to instantiate the boundary condition.
void operator()(plint iX, plint iY, Array<T,2>& u) const {
u[0] = 0.*poiseuilleVelocity(iY, parameters);
u[1] = T();
}
// This version of the operator returns also a constant value for
// the density, to create the initial condition.
void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const {
u[0] = 0.*poiseuilleVelocity(iY, parameters);
u[1] = T();
rho = (T)1;
}
private:
IncomprFlowParam<T> parameters;
};
// Function to set up the channel
void channelSetup (
MultiBlockLattice2D<T,NSDES>& lattice,
IncomprFlowParam<T> const& parameters,
OnLatticeBoundaryCondition2D<T,NSDES>& boundaryCondition )
{
// Create Velocity boundary conditions.
boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);
// Specify the boundary velocity.
setBoundaryVelocity (
lattice, lattice.getBoundingBox(),
PoiseuilleVelocity<T>(parameters) );
// Create the initial condition.
initializeAtEquilibrium (
lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );
lattice.initialize();
}
// ************************************************************************
// Advection-diffusion setup *********************************************
// ************************************************************************
void concentrationSetup(MultiBlockLattice2D<T,AVDES> &lattice1,
MultiBlockLattice2D<T,AVDES> &lattice2,
MultiBlockLattice2D<T,AVDES> &lattice3,
OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> &BC,
IncomprFlowParam<T> parameters1,
IncomprFlowParam<T> parameters2,
IncomprFlowParam<T> parameters3) {
// Get lattice dimensions
const plint nx = parameters1.getNx();
const plint ny = parameters1.getNy();
plint xx, yy;
T ran1, ran2;
// Set-up the geometry
//Box2D topwall(0, nx-1, ny-1, ny-1);
//Box2D bottomwall(0, nx-1, 0, 0);
Box2D leftwall(0, nx/6, 0,ny-1);
Box2D rightwall(nx-nx/6, nx-1, 0,ny-1);
// Set the boundary-conditions for lattice1
//BC.addTemperatureBoundary0P(rightwall, lattice1);
//BC.addTemperatureBoundary0N(leftwall, lattice1);
// Set the boundary-conditions for lattice2
//BC.addTemperatureBoundary0P(rightwall, lattice2);
//BC.addTemperatureBoundary0N(leftwall, lattice2);
// Try to impose a different temperature for lattice1
setBoundaryTemperature(lattice1, rightwall, 0.0);
setBoundaryTemperature(lattice1, leftwall, 0.5);
// Try to impose a different temperature for lattice2
setBoundaryTemperature(lattice2, rightwall, 0.5);
setBoundaryTemperature(lattice2, leftwall, 0.0);
// Try to impose a different temperature for lattice3
setBoundaryTemperature(lattice3, rightwall, 0.0);
setBoundaryTemperature(lattice3, leftwall, 0.25);
//Box2D domain(0, nx-1, 0, ny-1);
//
//Array<T,2> u1(1.0,0.);
//Array<T,2> u2(-1.0,0.);
//
//// Init lattice1
//initializeAtEquilibrium(lattice1, domain, 0.0, u1,0.0);
//initializeAtEquilibrium(lattice1, leftwall, 0.5,u1,0.5);
//// Init lattice2
//initializeAtEquilibrium(lattice2, domain, 0.0, u2,0.0);
//initializeAtEquilibrium(lattice2, rightwall, 0.5,u2,0.5);
pcout << "Preparing random velocities on lattices" << endl;
for(xx = 0; xx < nx; xx++)
{
for(yy = 0; yy < ny; yy++)
{
Box2D domain(xx, xx, yy, yy);
ran1 = 2*((double)random()/RAND_MAX-0.5);
ran2 = 2*((double)random()/RAND_MAX-0.5);
Array<T,2> u1(ran1,ran2);
// Init lattice1
initializeAtEquilibrium(lattice1, domain, 0.0, u1,0.0);
initializeAtEquilibrium(lattice1, leftwall, 0.5,u1,0.5);
// Init lattice2
ran1 = 2*((double)random()/RAND_MAX-0.5);
ran2 = 2*((double)random()/RAND_MAX-0.5);
Array<T,2> u2(-ran1,-ran2);
initializeAtEquilibrium(lattice2, domain, 0.0, u2,0.0);
initializeAtEquilibrium(lattice2, rightwall, 0.5,u2,0.5);
}
}
Box2D domain2(0, nx-1, 0, ny-1);
Array<T,2> u3(0.0,0.0);
initializeAtEquilibrium(lattice3, domain2, 0.0, u3,0.0);
initializeAtEquilibrium(lattice3, leftwall, 0.25,u3,0.25);
// Command to init the lattices
lattice1.periodicity().toggle(1,true);
lattice1.initialize();
lattice2.periodicity().toggle(1,true);
lattice2.initialize();
lattice3.periodicity().toggle(1,true);
lattice3.initialize();
}
// ************************************************************************
// ************************************************************************
// ************************************************************************
int main(int argc, char *argv[])
{
// Init simulation
plbInit(&argc, &argv);
global::directories().setOutputDir("./tmp/");
// Parameters for the dynamics
IncomprFlowParam<T> parameters1(
(T) 1e-2, // Reference velocity (the maximum velocity
// in the Poiseuille profile) in lattice units.
(T) 5000., // Reynolds number
150, // Resolution of the reference length (channel height).
2, // Channel length in dimensionless variables
1. // Channel height in dimensionless variables
);
IncomprFlowParam<T> parameters2(
(T) 1e-2,
(T) 50.,
150,
2,
1.
);
IncomprFlowParam<T> parameters3(
(T) 1e-2,
(T) 50.,
150,
2,
1.
);
stringstream buffer;
string filename1, filename2, filename3;
plint nx = parameters1.getNx();
const plint maxIter = 30000;
Array<T,2> u2(2.0, 0.0);
// Instantiate the NS lattice
MultiBlockLattice2D<T, NSDES> nsLattice(parameters1.getNx(), parameters1.getNy(), new BGKdynamics<T, NSDES>(parameters1.getOmega()));
// Instantiate adv-diff lattice1
MultiBlockLattice2D<T, AVDES> advLattice1(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(1.));
// Instantiate adv-diff lattice2
MultiBlockLattice2D<T, AVDES> advLattice2(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(1.));
// Instantiate adv-diff lattice3
MultiBlockLattice2D<T, AVDES> advLattice3(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(0.7));
// Instantiate the boundary conditions for both
OnLatticeBoundaryCondition2D<T,NSDES> *nsBC = createLocalBoundaryCondition2D<T, NSDES>();
OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> *avBC = createLocalAdvectionDiffusionBoundaryCondition2D<T, AVDES>();
// Set-up channel for NS equations
channelSetup(nsLattice, parameters1, *nsBC);
//Set-up concentration field for both lattices
concentrationSetup(advLattice1, advLattice2, advLattice3, *avBC, parameters1, parameters2,parameters3);
// Main LBM cycle
ImageWriter<T> image("leeloo");
for (plint iT=0; iT <= maxIter; ++iT) {
if(iT%100==0) {
pcout << "Navier-Stokes step " << iT << endl;
}
nsLattice.collideAndStream();
}
// First attempt to couple the two physics
latticeToPassiveAdvDiff(nsLattice, advLattice1, advLattice1.getBoundingBox());
latticeToPassiveAdvDiff(nsLattice, advLattice2, advLattice1.getBoundingBox());
latticeToPassiveAdvDiff(nsLattice, advLattice3, advLattice1.getBoundingBox());
for (plint iT=0; iT <= maxIter; ++iT) {
if(iT%100==0) {
pcout << "Writing GIF for advection - timestep = " << iT << endl;
image.writeScaledGif(createFileName("ADV-conc1_",iT,6), *computeDensity(advLattice1), parameters1.getNx(),parameters1.getNy());
image.writeScaledGif(createFileName("ADV-conc2_",iT,6), *computeDensity(advLattice2), parameters1.getNx(),parameters1.getNy());
image.writeScaledGif(createFileName("ADV-conc3_",iT,6), *computeDensity(advLattice3), parameters1.getNx(),parameters1.getNy());
buffer << "./ascii_output/conc_profile1_";
buffer << setprecision(5) << iT << ".dat";
filename1 = buffer.str();
//pcout << filename << endl; // Debug - Check filename
plb_ofstream concProfile1(filename1.c_str());
buffer.str("");
buffer << "./ascii_output/conc_profile2_";
buffer << setprecision(5) << iT << ".dat";
filename2 = buffer.str();
//pcout << filename << endl; // Debug - Check filename
plb_ofstream concProfile2(filename2.c_str());
buffer.str("");
buffer << "./ascii_output/conc_profile3_";
buffer << setprecision(5) << iT << ".dat";
filename3 = buffer.str();
//pcout << filename << endl; // Debug - Check filename
plb_ofstream concProfile3(filename3.c_str());
buffer.str("");
for(nx = 0; nx < parameters1.getNx(); nx++)
{
Box2D concSection(nx, nx, 0.5*parameters1.getNy(), 0.5*parameters1.getNy());
concProfile1 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice1,concSection) << endl;
concProfile2 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice2,concSection) << endl;
concProfile3 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice3,concSection) << endl;
}
//image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy());
}
advLattice1.collideAndStream();
advLattice2.collideAndStream();
advLattice3.collideAndStream();
}
//pcout << "Hello palabos";
return 0;
}