MRT vs BGK

[kolotinsky]

Hello everybody!
I have tested my openlb program using BGK and MRT dynamics. In my program I simulate the following equation with periodical boundary conditions:

df/dt +( az – b*Vflow) * df/dv = collision operator

az – constant acceleration directed along z axis

b – coefficient of viscosity

The analytical result for flow velocity is: Vflow(t) = a/b * [1 – exp(-b*t)]

Both MRT and BGK dynamics give this result correctly, however, when I use BGK operator density remains constant and when I use MRT operator density increases. I thought MRT must give more precise results than BGK. What is wrong?

My code you can see bellow

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#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 ForcedMRTD3Q19Descriptor

// Parameters for the simulation setup
const int maxIter = 500;
const int nx = 30;
const int ny = 30;
const int nz = 30;

// Stores geometry information in form of material numbers
void prepareGeometry( SuperGeometry3D<T>& superGeometry ) {

OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

// Sets material number for fluid
superGeometry.rename( 0,1 );

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();

superGeometry.print();

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

// Set up the geometry of the simulation
void prepareLattice( SuperLattice3D<T, DESCRIPTOR>& sLattice,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
SuperGeometry3D<T>& superGeometry ) {

// Material=1 –>bulk dynamics
sLattice.defineDynamics( superGeometry, 1, &bulkDynamics1 );

// Initial conditions

std::vector<T> v( 3,T() );
AnalyticalConst3D<T,T> zeroVelocity( v );
AnalyticalConst3D<T,T> Rho( 1. );

// Initialize force

std::vector<T> force(3, T());
force[0] = 0;
force[1] = 0;
force[2] = 0.00075;

//olb::BlockVTIreader3D< T,BaseType > datareader(“forcefield.vti”,”ForceField”);
//BlockDataF3D<T, BaseType> blockfield(datareader.getBlockData());
//AnalyticalF3D<T,T>* field;
//field = new SpecialAnalyticalFfromBlockF3D<T,T>(blockfield, datareader.getCuboid(),spacing);
AnalyticalConst3D<T,T> constant(force);
sLattice.defineExternalField( superGeometry, 1,
DESCRIPTOR<T>::ExternalField::forceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfForce, constant );

// Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( superGeometry, 1, Rho, zeroVelocity );
sLattice.iniEquilibrium( superGeometry, 1, Rho, zeroVelocity );

// Make the lattice ready for simulation
sLattice.initialize();

//delete field;
//field = 0;
}

// Output to console and files
void getResults( SuperLattice3D<T, DESCRIPTOR>& sLattice, int iT,
SuperGeometry3D<T>& superGeometry, Timer<T>& timer ) {

OstreamManager clout( std::cout,”getResults” );

SuperVTMwriter3D<T> vtmWriter( “ionwake3d” );
SuperLatticeDensity3D<T, DESCRIPTOR> density( sLattice );
vtmWriter.addFunctor( density );

const int statIter = 1;
vtmWriter.createMasterFile();
vtmWriter.write();

// Writes output on the console
if ( iT%statIter==0 ) {
// Timer console output
timer.update( iT );
timer.printStep();

// Lattice statistics console output
sLattice.getStatistics().print( iT,iT );
}
}

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

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

const T omega1 = 0.1;

// === 2rd Step: Prepare Geometry ===

const int noOfCuboids = 1;

CuboidGeometry3D<T> cuboidGeometry( 0, 0, 0, 1, nx, ny, nz, noOfCuboids );

// Periodic boundaries in x- and y- and z-direction
cuboidGeometry.setPeriodicity( true, true, true );

// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
SuperGeometry3D<T> superGeometry( cuboidGeometry,loadBalancer,2 );

prepareGeometry( superGeometry );

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

ForcedMRTdynamics<T, DESCRIPTOR> bulkDynamics1 ( omega1, instances::getBulkMomenta<T, DESCRIPTOR>() );

prepareLattice( sLattice, bulkDynamics1, superGeometry );

// === 4th Step: Main Loop ===
int iT = 0;
clout << “starting simulation…” << endl;
Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();

//std::vector<T> spacing(3, T());
//spacing[0] = 1;
//spacing[1] = 1;
//spacing[2] = 1;
ofstream fout(“velocity.txt”);
for ( iT = 0; iT < maxIter; ++iT ) {

// === 5th Step: Collide and Stream Execution ===
//olb::BlockVTIreader3D< T,BaseType > datareader(“forcefield.vti”,”ForceField”);
//BlockDataF3D<T, BaseType> blockfield(datareader.getBlockData());
//AnalyticalF3D<T,T>* field;
//field = new SpecialAnalyticalFfromBlockF3D<T,T>(blockfield, datareader.getCuboid(),spacing);
//sLattice.defineExternalField( superGeometry, 1,
// DESCRIPTOR<T>::ExternalField::forceBeginsAt,
// DESCRIPTOR<T>::ExternalField::sizeOfForce, *field );
std::vector<T> force(3, T());
force[0] = 0;
force[1] = 0;
force[2] = 0.001-0.011*sLattice.getStatistics().getMaxU();
AnalyticalConst3D<T,T> f(force);
sLattice.defineExternalField( superGeometry, 1,
DESCRIPTOR<T>::ExternalField::forceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfForce, f );
sLattice.collideAndStream();
fout<<sLattice.getStatistics().getMaxU()<<endl;
//delete field;
//field = 0;

// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, iT, superGeometry, timer );
}

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

[mathias]

Actually in this case with periodic BC, pressure=density is obtained up to a constant (as also in the NSE systems). You can use a pressure filter/projection to get a unique solution (see https://www.openlb.net/wp-content/uploads/2011/12/olb-tr3.pdf). Maybe the MRT and BGK version uses a different force modell. You can adapt them by writing a new dynamics if you like.

Best
Mathias

[Author]

Posts