After I started with openLB, unit conversion gave me quite some headache for a while, so I decided to grab all the things I need together and form a small library out of it, and therefore never have to take a look at it again - as long as it works hopefully 
Although the class is very small at the moment and only usefull for an applied pressure gradient, I post it here in case anybody else might find it usefull.
It’s basically very simple:
the class Units inherits from IncomprFlowParam and therefore it replaces it in the code.
So that’s what the code can look like - instead of the IncomprFlowParam:
// Define physical and numerical units
Units<T> units (
(T) (4.E-3), // lx_p [m]
(T) 4.E-3, // ly_p [m]
(T) (4.E-3), // lref_p [m]
(T) (5.E-6)/7., // visc_p [m²s?¹]
(T) 7000., // rho_p [kg m?³]
(T) 0.01, // Pressure gradient [Pa m?¹]
(T) 0.9, // force tau to 0.9
200 // Lattice Cells
);
With
units.WriteUnitsFile("Units.dat");
a more detailed units file is written. Also I intend to add some post-processing calculations to the library, but at the moment there exists only one.
There are 3 Constructors available, whit 2 very simple ones. The third one is the one I use - because with it you can directly define the value of tau.
As you can see in the code, there are a few values which are always 1, and don’t do anything. I just added them since I followed J. Latts paper on unit conversion when I wrote it, and added them just in case
So far, it’s only for 2D, but I plan to add 3D as well.
Also, I really don’t like the way I coded the call of the IncomprFlowParam - Constructor - but I can’t think of another 
so if anybody knows a better way to solve this please let me know
Here is the code:
/**
** Units conversion for the Palabos library.
** Author: Christian Hoelzl, christian.hoelzl@unileoben.ac.at
**/
#include <iostream>
#include "palabos2D.h"
#include <cmath>
using namespace plb;
using namespace plb::util;
using namespace std;
#define DESCRIPTOR descriptors::ForcedD2Q9Descriptor
typedef double T;
/****************************************************************************************************************************************
* Class: Units *
* Offers different Constructors to simplify unit conversion. *
* At the moment, only works for an applied pressure gradient. *
* *
*****************************************************************************************************************************************/
template <typename T>
class Units : public plb::IncomprFlowParam<T> {
private:
plint tSave;
//Variables in physical units
T lx_p; // lenght in x-direction [m]
T ly_p; // lenght in y-direction [m]
T lref_p; // reference lenght [m]
T tref_p; // reference time [s]
T uref_p; // reference velocity [ms?¹]
T visc_p; // kinematic viscosity of the liquid [m²s?¹]
T rho_p; // Density [kg m?³]
T pgrad_p; // pressure gradient [Pa/m]
T dp_p; // pressure difference along x-axis [Pa] bzw. [kgm?¹s?²] = lx_p * pgrad_p
T Re; // Reynolds-number []
//dimensionless:
T lx_d; // dimensionless lenght in x-direction = lx_p / lref_p
T ly_d; // dimensionless lenght in y-direction = ly_p / lref_p
T lref_d; // dimensionless reference lenght (is always 1) = 1
T tref_d; // dimensionless reference time (is always 1) = 1
T uref_d; // dimensionless reference velocity(is always 1) = lref_d / t_d = 1 !!
T visc_d; // dimensionless kinematic viscosity = 1 / Re
T dp_d; // dimensionless pressure difference = 1 / rho_p * tref_p² / lref_p² * dp_p
//Lattice units:
plint N; // Number of lattice cells per reference lenght
T xstep_lb; // discrete space intervall = 1 / N
T tstep_lb; // discrete timestep = 1 / N_iter
T uref_lb; // velocity in LB units = tstep_lb / xstep_lb or = N / N_iter
T dp_lb; // pressure difference in LB units = tstep_lb² / xstep_lb² * dp_d
T K; // applied Force per LB lattice cell = dp_lb / N
T rho_lb; // Density in LB units = 1
T visc_lb; // Viscosity in LB units = tstep_lb / xstep_lb² * visc_d
T visc_plb; // Viscosity calculated by Palabos = plb::IncomprFlowParam<T>::getLatticeNu()
T tau_lb; // Tau of the LB system = 3 * visc_lb + 0.5
T tau_plb; // Tau calculated by Palabos = plb::IncomprFlowParam<T>::getTau()
T omega_lb; // Omega of the LB system = 1 / tau_lb
T omega_plb; // Omega calculated by Palabos = plb::IncomprFlowParam<T>::getOmega()
public:
Units(T lx_p_, T ly_p_, T visc_p_, T rho_p_, T pgrad_p_, plint N_);
Units(T lx_p_, T ly_p_, T lref_p_, T visc_p_, T rho_p_, T pgrad_p_, plint N_);
Units(T lx_p_, T ly_p_, T lref_p, T visc_p_, T rho_p_, T pgrad_p_, T tau_lb_, plint N_);
void WriteUnitsFile(string fName);
void WriteInletVelocity ( string fName, MultiBlockLattice2D<T, DESCRIPTOR>& lattice, T deltaE);
T get_dp_p() { return dp_p; }
T get_K() { return K; }
//void Units<T>WriteMyLogFile (Units<T> const& parameters, std::string const& title)
};
/****************************************************************************************************************************************
* Constructor: (lx_p_, ly_p_, visc_p_, rho_p_, pgrad_p_, N_) *
* *
* This is the most simple constructor. The lattice velocity will be calculated from the relation between the timestep and the *
* space intervall: tstep_lb = xstep_lb². *
* the reference lenght lref_p is the lenght in x-direction: lref_p = lx_p *
* lx_d = lx_p / lref_p *
* lref_d = lx_d ( =1) *
* *
* With only N_ as input parameter for the lattice, u_lb will be calculated from u_lb = tstep_lb / xstep_lb *
* *
* -> the relation tstep_lb = xstep_lb² directly leads to u_lb = 1 / N *
* *
* IncomprFlowParam (T latticeU_, T Re_, plint resolution_, T lx_, T ly_, T lz_=T()) *
*****************************************************************************************************************************************/
template <typename T>
Units<T>::Units(T lx_p_, T ly_p_, T visc_p_, T rho_p_, T pgrad_p_, plint N_) :
plb::IncomprFlowParam<T> ((T) 1. / N_, // Lattice velocity
sqrt(2*lx_p_*pgrad_p_ / rho_p_) * lx_p_ / visc_p_, // Reynoldsnumber
N_, // number of lattice cells per reference lenght
1, // lx, in this case also the reference lenght
ly_p_ / lx_p_) // ly
{
// Input variables
lx_p = lx_p_; // physical length along x-axis
ly_p = ly_p_; // physical lenght along y-axis
visc_p = visc_p_; // physical viscosity
rho_p = rho_p_; // physical density
pgrad_p = pgrad_p_; // physical pressure gradient
uref_lb = (T) 1. / N_; // velocity in LB units
N = N_; // Number of Lattice Cells per reference lenght
lref_p = lx_p_; // reference lenght
//Variables in physical units
dp_p = lx_p * pgrad_p; // kg/(m s^2) = 1 Pa
uref_p = sqrt(2*dp_p/rho_p); // m/s reference velocity
tref_p = lref_p / uref_p; // s reference time
Re = uref_p * lref_p / visc_p; // Reynolds number
//dimensionless:
lx_d = lx_p/lref_p; // dimensionless length in x-direction, in this case = 1
ly_d = ly_p/lref_p; // dimensionless lenght in y-direction
uref_d = tref_p/lref_p * uref_p; // dimensionless velocity (= 1 !!)
visc_d = 1./Re; // dimensionless viscosity
dp_d = 1./rho_p * (tref_p*tref_p) / (lref_p * lref_p) * dp_p; // dimensionless pressure difference
lref_d = lx_d; // dimensionless reference lenght, is always 1
tref_d = 1; // dimensionless reference time, also always 1
//Lattice units:
tstep_lb = plb::IncomprFlowParam<T>::getDeltaT(); // get the timestep from IncomprFlowParam
xstep_lb = plb::IncomprFlowParam<T>::getDeltaX(); // same for the discrete space intervall
dp_lb = tstep_lb*tstep_lb / (xstep_lb * xstep_lb) * dp_d; // calculate the pressure difference in LB units
K = dp_lb / N; // and finally the value for the applied external force
rho_lb = 1.; // LB density = 1 !!
visc_lb = tstep_lb / (xstep_lb * xstep_lb) * visc_d; // viscosity in LB units, to calculate tau
visc_plb = plb::IncomprFlowParam<T>::getLatticeNu(); // compare the value with palabos result
tau_lb = 3.*visc_lb + 0.5; // LB constant for the simulation, where 3 = 1/(soundspeed)² in the LB system
tau_plb = plb::IncomprFlowParam<T>::getTau(); // compare the value with palabos result
omega_lb = 1./tau_lb; // omega from tau
omega_plb = plb::IncomprFlowParam<T>::getOmega(); // compare the value with palabos result
}
/****************************************************************************************************************************************
* Reference Lenght <lref_p> is an required input value for this constructor: *
* *
* lref_p = input value lref_p_ *
* lx_d = lx_p / lref_p *
* lref_d = lx_d *
* *
*****************************************************************************************************************************************/
template <typename T>
Units<T>::Units(T lx_p_, T ly_p_, T lref_p_, T visc_p_, T rho_p_, T pgrad_p_, plint N_) :
plb::IncomprFlowParam<T> ((T) 1. / N_, // Lattice velocity
sqrt(2*lx_p_*pgrad_p_ / rho_p_) * lref_p_ / visc_p_, // Reynoldsnumber
N_, // number of lattice cells per reference lenght
lx_p_ / lref_p_, // lx_d = dimensionless lenght in x-direction
ly_p_ / lref_p_) // ly_d = dimensionless lenght in y-direction
{
// Input variables
lx_p = lx_p_; // physical length along x-axis
ly_p = ly_p_; // physical lenght along y-axis
visc_p = visc_p_; // physical viscosity
rho_p = rho_p_; // physical density
pgrad_p = pgrad_p_; // physical pressure gradient
uref_lb = (T) 1. / N_; // velocity in LB units
N = N_; // Number of Lattice Cells per reference lenght
lref_p = lref_p_; // reference lenght
//Variables in physical units
dp_p = lx_p * pgrad_p; // kg/(m s^2) = 1 Pa
uref_p = sqrt(2*dp_p/rho_p); // m/s reference velocity
tref_p = lref_p / uref_p; // s reference time
Re = uref_p * lref_p / visc_p; // Reynolds number
//dimensionless:
lx_d = lx_p/lref_p; // dimensionless length in x-direction, in this case = 1
ly_d = ly_p/lref_p; // dimensionless lenght in y-direction
uref_d = tref_p/lref_p * uref_p; // dimensionless velocity (= 1 !!)
visc_d = 1./Re; // dimensionless viscosity
dp_d = 1./rho_p * (tref_p*tref_p) / (lref_p * lref_p) * dp_p; // dimensionless pressure difference
lref_d = 1; // dimensionless reference lenght, is always 1
tref_d = 1; // dimensionless reference time, also always 1
//Lattice units:
tstep_lb = plb::IncomprFlowParam<T>::getDeltaT(); // get the timestep from IncomprFlowParam
xstep_lb = plb::IncomprFlowParam<T>::getDeltaX(); // same for the discrete space intervall
dp_lb = tstep_lb*tstep_lb / (xstep_lb * xstep_lb) * dp_d; // calculate the pressure difference in LB units
K = dp_lb / N; // and finally the value for the applied external force
rho_lb = 1.; // LB density = 1 !!
visc_lb = tstep_lb / (xstep_lb * xstep_lb) * visc_d; // viscosity in LB units, to calculate tau
visc_plb = plb::IncomprFlowParam<T>::getLatticeNu(); // compare the value with palabos result
tau_lb = 3.*visc_lb + 0.5; // LB constant for the simulation, where 3 = 1/(soundspeed)² in the LB system
tau_plb = plb::IncomprFlowParam<T>::getTau(); // compare the value with palabos result
omega_lb = 1./tau_lb; // omega from tau
omega_plb = plb::IncomprFlowParam<T>::getOmega(); // compare the value with palabos result
}
/****************************************************************************************************************************************
* Constructor: (lx_p_, ly_p_, visc_p_, rho_p_, pgrad_p_, u_lb_, N_) *
* *
* the reference lenght lref_p is the lenght in x-direction: lref_p = lx_p *
* lx_d = lx_p / lref_p *
* lref_d = lx_d ( =1) *
* *
* Tau is an input parameter. Therefore the lattice velocity uref_lb is calculated. *
* with: *
* visc_d = 1 / Re = visc_p / (uref_p * lref_p), uref_p = sqrt(2 * lx_p * pgrad_p / rho_p) *
* visc_d = xstep_lb² / tstep_lb * visc_lb, xstep_lb = 1 / N *
* uref_lb = tstep_lb / xstep_lb *
* *
* -> uref_lb = ( (tau_lb - 0.5) * sqrt(2 * lx_p * pgrad_p / rho_p) * lref_p ) / (3 * N * visc_p) *
* IncomprFlowParam (T latticeU_, T Re_, plint resolution_, T lx_, T ly_, T lz_=T()) *
*****************************************************************************************************************************************/
template <typename T>
Units<T>::Units(T lx_p_, T ly_p_, T lref_p_, T visc_p_, T rho_p_, T pgrad_p_, T tau_lb_, plint N_) :
plb::IncomprFlowParam<T> ( ( (tau_lb_ - 0.5) * sqrt(2.*lx_p_*pgrad_p_ / rho_p_) * lref_p_) / (3. * (T) N_ * visc_p_), // Lattice velocity
sqrt(2*lx_p_*pgrad_p_ / rho_p_) * lref_p_ / visc_p_, // Reynoldsnumber
N_, // number of lattice cells per reference lenght
lx_p_ / lref_p_, // lx, in this case also the reference lenght
ly_p_ / lref_p_) // ly
{
// Input variables
lx_p = lx_p_; // physical length along x-axis
ly_p = ly_p_; // physical lenght along y-axis
lref_p = lref_p_; // reference lenght
visc_p = visc_p_; // physical viscosity
rho_p = rho_p_; // physical density
pgrad_p = pgrad_p_; // physical pressure gradient
tau_lb = tau_lb_; // tau
N = N_; // Number of Lattice Cells per reference lenght
uref_lb = ( (tau_lb_ - 0.5) * sqrt(2.*lx_p_*pgrad_p_ / rho_p_) * lref_p_) / (3. * (T) N_ * visc_p_);
cout << plb::IncomprFlowParam<T>::getLatticeU() << endl;
//Variables in physical units
dp_p = lx_p * pgrad_p; // kg/(m s^2) = 1 Pa
uref_p = sqrt(2.*dp_p/rho_p); // m/s reference velocity
tref_p = lref_p / uref_p; // s reference time
Re = uref_p * lref_p / visc_p; // Reynolds number
//dimensionless:
lx_d = lx_p/lref_p; // dimensionless length in x-direction, in this case = 1
ly_d = ly_p/lref_p; // dimensionless lenght in y-direction
uref_d = tref_p/lref_p * uref_p; // dimensionless velocity (= 1 !!)
visc_d = 1./Re; // dimensionless viscosity
dp_d = 1./rho_p * (tref_p*tref_p) / (lref_p * lref_p) * dp_p; // dimensionless pressure difference
lref_d = 1; // dimensionless reference lenght, is always 1
tref_d = 1; // dimensionless reference time, also always 1
//Lattice units:
tstep_lb = plb::IncomprFlowParam<T>::getDeltaT(); // get the timestep from IncomprFlowParam
xstep_lb = plb::IncomprFlowParam<T>::getDeltaX(); // same for the discrete space intervall
dp_lb = tstep_lb*tstep_lb / (xstep_lb * xstep_lb) * dp_d; // calculate the pressure difference in LB units
K = dp_lb / N; // and finally the value for the applied external force
rho_lb = 1.; // LB density = 1 !!
visc_lb = tstep_lb / (xstep_lb * xstep_lb) * visc_d; // viscosity in LB units, to calculate tau
visc_plb = plb::IncomprFlowParam<T>::getLatticeNu(); // compare the value with palabos result
tau_lb = 3.*visc_lb + 0.5; // LB constant for the simulation, where 3 = 1/(soundspeed)² in the LB system
tau_plb = plb::IncomprFlowParam<T>::getTau(); // compare the value with palabos result
omega_lb = 1./tau_lb; // omega from tau
omega_plb = plb::IncomprFlowParam<T>::getOmega(); // compare the value with palabos result
}
template <typename T>
void Units<T>::WriteUnitsFile(string fName)
{
std::string fullName = global::directories().getLogOutDir() + fName.c_str();
plb_ofstream variablesOut(fullName.c_str());
variablesOut << "Variables used in Simulation: " << endl;
variablesOut << "---------------------------------------" << endl << endl;
variablesOut << "Physical Units: " << endl << endl;
variablesOut << "Lenght in x-direction: lx_p = " << setw(15) << setfill(' ') << lx_p << " [m]" << endl;
variablesOut << "Lenght in y-direction: ly_p = " << setw(15) << setfill(' ') << ly_p << " [m]" << endl;
variablesOut << "Reference lenght lref_p = " << setw(15) << setfill(' ') << lref_p << " [m]" << endl;
variablesOut << "Reference time tref_p = " << setw(15) << setfill(' ') << tref_p << " [s]" << endl;
variablesOut << "Reference velocity uref_p = " << setw(15) << setfill(' ') << uref_p << " [m/s]" << endl;
variablesOut << "Kinematic viscosity visc_p = " << setw(15) << setfill(' ') << visc_p << " [m^2/s]" << endl;
variablesOut << "Density rho_p = " << setw(15) << setfill(' ') << rho_p << " [kg/m^3]" << endl;
variablesOut << "Pressure gradient pgrad_p = " << setw(15) << setfill(' ') << pgrad_p << " [Pa/m]" << endl;
variablesOut << "Pressure difference dp_p = " << setw(15) << setfill(' ') << dp_p << " [Pa]" << endl;
variablesOut << "Reynolds number Re = " << setw(15) << setfill(' ') << Re << " []" << endl << endl << endl;
variablesOut << "Dimensionless Units: " << endl << endl;
variablesOut << "Length in x-direction: lx_d = " << setw(15) << setfill(' ') << lx_d << endl;
variablesOut << "Lenght in y-direction: ly_d = " << setw(15) << setfill(' ') << ly_d << endl;
variablesOut << "Reference lenght lref_d = " << setw(15) << setfill(' ') << lref_d << endl;
variablesOut << "Reference time tref_d = " << setw(15) << setfill(' ') << tref_d << endl;
variablesOut << "Reference velocity uref_d = " << setw(15) << setfill(' ') << uref_d << endl;
variablesOut << "Kinematic viscosity visc_d = " << setw(15) << setfill(' ') << visc_d << endl;
variablesOut << "Pressure difference dp_d = " << setw(15) << setfill(' ') << dp_d << endl << endl << endl;
variablesOut << "Lattice Units: " << endl << endl;
variablesOut << "Lattice Cells in x-direction: Nx = " << setw(15) << setfill(' ') << plb::IncomprFlowParam<T>::getNx() << endl;
variablesOut << "Lattice Cells in y-direction: Ny = " << setw(15) << setfill(' ') << plb::IncomprFlowParam<T>::getNy() << endl;
variablesOut << "Discrete space Intervall xstep_lb = " << setw(15) << setfill(' ') << xstep_lb << endl;
variablesOut << "Discrete time step tstep_lb = " << setw(15) << setfill(' ') << tstep_lb << endl;
variablesOut << "Lattice Velocity u_lb = " << setw(15) << setfill(' ') << uref_lb << endl;
variablesOut << "Viscosity visc_lb = " << setw(15) << setfill(' ') << visc_lb << endl;
variablesOut << "Viscosity from Palabos visc_plb = " << setw(15) << setfill(' ') << visc_plb << endl;
variablesOut << "Pressure difference dp_lb = " << setw(15) << setfill(' ') << dp_lb << endl;
variablesOut << "Applied force per Lattice cell K = " << setw(15) << setfill(' ') << K << endl;
variablesOut << "Tau tau_lb = " << setw(15) << setfill(' ') << tau_lb << endl;
variablesOut << "Tau from Palabos tau_plb = " << setw(15) << setfill(' ') << tau_plb << endl;
variablesOut << "Omega omega_lb = " << setw(15) << setfill(' ') << omega_lb << endl;
variablesOut << "Omega from Palabos omega_plb = " << setw(15) << setfill(' ') << omega_plb << endl;
}
template <typename T>
void Units<T>::WriteInletVelocity ( string fName,
MultiBlockLattice2D<T, DESCRIPTOR>& lattice, T error)
{
std::string fullName = global::directories().getLogOutDir() + fName.c_str();
plb_ofstream ofile(fullName.c_str());
ofile << "Error at last step: " << error << endl;
ofile << "iY[Lattice Cell Nr.]\t iY[Lattice Units]\t iY[m]\t ux_lb[Lattice Units]\t ux_d[dimless Units] uref_p[m/s]\n";
for (int iY=0; iY<plb::IncomprFlowParam<T>::getNy(); ++iY)
{
Array<T, 2> velocity;
lattice.get(0, iY).computeVelocity(velocity);
ofile << iY << " " << iY*xstep_lb << " " << iY*xstep_lb*lref_p << " " << velocity[0] << " " << velocity[0]*xstep_lb/tstep_lb <<
" " << velocity[0]*xstep_lb/tstep_lb*lref_p/tref_p << "\n";
}
}
PS: I just looked at the code for some time and I think the descriptor definition at the begining is not required here.
Edith found some typos