Problem with MRT compiling allModels2d

Hello,

I am having problems compiling the example codesByTopic/navierStokesModels/allModels2d.cpp when I comment and uncomment lines to active MRT. If I leave the code as-is, it compiles and runs flawlessly, but when I change lines 107 to 116, I get problems. I’ll post the full code in the bottom of the thread. The lines in question are:

[code=“cpp”]
107 // Comment these 3 lines for MRT
108
109 MultiBlockLattice2D<T, DESCRIPTOR> lattice (
110 parameters.getNx(), parameters.getNy(),
111 new BackgroundDynamics(parameters.getOmega()) );
112
113 //Uncomment these 4 lines for MRT
114
115 //MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
116 //MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), new BackgroundDynamics(mrt_param) );



This way, it compiles straight away. But when I change it to

[code="cpp"]
107 // Comment these 3 lines for MRT
108 
109 //   MultiBlockLattice2D<T, DESCRIPTOR> lattice (
110 //           parameters.getNx(), parameters.getNy(),
111 //          new BackgroundDynamics(parameters.getOmega()) );
112
113 //Uncomment these 4 lines for MRT
114   
115    MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
116    MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), new    BackgroundDynamics(mrt_param) );

I receive the following errors:


eduardo@Tardis:~/Documents/UFSC/PG/Doutorado/Palabos/palabos-v1.5r1/examples/codesByTopic/navierStokesModels$ make
python ../../../scons/scons.py -j 6 -f ../../../SConstruct palabosRoot=../../.. projectFiles="allModels2d_lixo.cpp" optimize=true debug=true profile=false MPIparallel=true SMPparallel=false usePOSIX=true serialCXX=g++ parallelCXX=mpicxx compileFlags="-Wall -Wnon-virtual-dtor" linkFlags="" optimFlags="-O3" debugFlags="-g" profileFlags="-pg" libraryPaths="" includePaths="../include" libraries=""
scons: Reading SConscript files ...
scons: done reading SConscript files.
scons: Building targets ...
mpicxx -o allModels2d_lixo.o -c -Wall -Wnon-virtual-dtor -O3 -g -DPLB_DEBUG -DPLB_MPI_PARALLEL -DPLB_USE_POSIX -I/home/eduardo/Documents/UFSC/PG/Doutorado/Palabos/palabos-v1.5r1/src -I/home/eduardo/Documents/UFSC/PG/Doutorado/Palabos/palabos-v1.5r1/externalLibraries -I/home/eduardo/Documents/UFSC/PG/Doutorado/Palabos/palabos-v1.5r1/examples/codesByTopic/include allModels2d_lixo.cpp
allModels2d_lixo.cpp: In function ‘int main(int, char**)’:
allModels2d_lixo.cpp:115:5: error: ‘MRTparam’ was not declared in this scope
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
     ^
allModels2d_lixo.cpp:115:15: error: expected primary-expression before ‘,’ token
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
               ^
allModels2d_lixo.cpp:115:26: error: missing template arguments before ‘>’ token
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
                          ^
allModels2d_lixo.cpp:115:30: error: ‘mrt_param’ was not declared in this scope
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
                              ^
allModels2d_lixo.cpp:115:46: error: ‘MRTparam’ does not name a type
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
                                              ^
allModels2d_lixo.cpp:115:56: error: expected primary-expression before ‘,’ token
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
                                                        ^
allModels2d_lixo.cpp:115:67: error: missing template arguments before ‘>’ token
     MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
                                                                   ^

I am trying my way around this using the documentation but I can’t solve this problem. I must use MRT because I have to tune different relaxation parameters for viscosity and diffusivity.

In a related question, can I define multiple sets of MRTparam just like I can define multiple sets of IncomprFlowParam, to model several components?

The full code is:

[code=“cpp”]
/* This file is part of the Palabos library.
*

  • Copyright © 2011-2015 FlowKit Sarl
  • Route d’Oron 2
  • 1010 Lausanne, Switzerland
  • E-mail contact: contact@flowkit.com
  • The most recent release of Palabos can be downloaded at
  • http://www.palabos.org/
  • The library Palabos is free software: you can redistribute it and/or
  • modify it under the terms of the GNU Affero General Public License as
  • published by the Free Software Foundation, either version 3 of the
  • License, or (at your option) any later version.
  • The library is distributed in the hope that it will be useful,
  • but WITHOUT ANY WARRANTY; without even the implied warranty of
  • MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  • GNU Affero General Public License for more details.
  • You should have received a copy of the GNU Affero General Public License
  • along with this program. If not, see http://www.gnu.org/licenses/.
    */

#include “palabos2D.h”
#include “palabos2D.hh”
#include “complexDynamics/mrtDynamics.h”
#include
#include
#include
#include
#include
#include “poiseuille.h”
#include “poiseuille.hh”
#include “cylinder.h”
#include “cylinder.hh”

using namespace plb;
using namespace plb::descriptors;
using namespace std;

typedef double T;

// Uncomment the two following lines for BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef BGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for so-called incompressible BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef IncBGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for Regularized BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef RegularizedBGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for MRT dynamics
#define DESCRIPTOR MRTD2Q9Descriptor
typedef MRTdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for Entropic dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef EntropicDynamics<T,DESCRIPTOR> BackgroundDynamics;

void defineCylinderGeometry( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();

plint cx      = nx/4;
plint cy      = ny/2+ny/10;
plint radius  = cy/4;

createCylinder(lattice, cx, cy, radius);

}

void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;

ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif(createFileName("u", iter, 6),
                           *computeVelocityNorm(lattice),
                           imSize, imSize );

}

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

global::directories().setOutputDir("./tmp/");

IncomprFlowParam<T> parameters(
        (T) 1e-2,  // uMax
        (T) 300.,  // Re
        100,       // N
        5.,        // lx
        1.         // ly 
);
const T logT     = (T)0.02;
const T imSave   = (T)0.1;
const T maxT     = (T)10.1;

writeLogFile(parameters, "Poiseuille flow");

// Comment these 3 lines for MRT

// MultiBlockLattice2D<T, DESCRIPTOR> lattice (
// parameters.getNx(), parameters.getNy(),
// new BackgroundDynamics(parameters.getOmega()) );

//Uncomment these 4 lines for MRT

MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
MultiBlockLattice2D<T, DESCRIPTOR> lattice (parameters.getNx(), parameters.getNy(), new BackgroundDynamics(mrt_param) );

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
    boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

defineCylinderGeometry(lattice, parameters);

createPoiseuilleBoundaries(lattice, parameters, *boundaryCondition);
createPoiseuilleInitialValues(lattice, parameters);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
    if (iT%parameters.nStep(logT)==0) {
        pcout << "step " << iT
              << "; lattice time=" << lattice.getTimeCounter().getTime()
              << "; t=" << iT*parameters.getDeltaT()
              << "; av energy="
              << setprecision(10) << getStoredAverageEnergy<T>(lattice)
              << "; av rho="
              << getStoredAverageDensity<T>(lattice) << endl;
    }

    if (iT%parameters.nStep(imSave)==0) {
        pcout << "Saving Gif ..." << endl;
        writeGifs(lattice, iT);
    }

    // Lattice Boltzmann iteration step.
    lattice.collideAndStream();
}

delete boundaryCondition;

}

Dear eduardowoj,

searching the code revealed that the class MRTparam simply does not exist. Also, if you look at the implementation of the MRTdynamics ( ./src/complexDynamics/mrtDynamics.h ), you will realize that it only needs omega in the constructor and NOT an object of the MRTparam class as done in the allModels2d.cpp.

To get this to work, you can change the line

[code=“cpp”]
MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega()));
MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameter.getNx(), parameters.getNy(),
new BackgroundDynamics(mrt_param) );



to

[code="cpp"]
// MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega()));
MultiBlockLattice2D<T, DESCRIPTOR> lattice (
     parameter.getNx(), parameters.getNy(),
     new BackgroundDynamics(parameters.getOmega()) );

and thus simply avoid using the (non-existing) MRTparam class. When I apply these changes, the code compiles and leads to stable simulations.
I am assuming that at some point the MRTparam class got obsolete and was removed from the code, but they forgot to update the example. Also, since this part of the code is obviously not maintained properly and no developer/admin is answering to this problem, I would highly recommend that you dig into the implementation and make some tests before using it for an actual problem.

Regards,

kk

Dear kk,

First of all, thanks for your prompt response. The example code had both implementations, with the one using the MRTparam<T,DESCRIPTOR> function already commented. The other implementation, the one you suggested to uncomment, has no explicit mention to MRT, so I thought it wasn’t the most appropriated way to implement a MRT dynamics to the lattice.

Let me see if I understood it correctly. The lines

[code=“cpp”]
#define DESCRIPTOR MRTD2Q9Descriptor
typedef MRTdynamics<T,DESCRIPTOR> BackgroundDynamics;



make "BackgroundDynamics" be an alias to MRTdynamics, so when I declare

MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(),
  new BackgroundDynamics(parameters.getOmega()) );

I am actually indicating this lattice [b]will[/b] use MRTdynamics as parsed by BackgroundDynamics?

The MRT-related functions for setting and getting the epsilon, lambda and other relaxation parameters are all written for the MRTparam. How can I, then, set and get lambda, epsilon and the other relaxation times from the relaxation matrix with this implementation of the MRT?

Thanks again!!

Dear eduardowoj,

yes, the typedef command will replace “BackgroundDynamics” with the respective definition at the top of the code. Thus, if you have the MRTdynamics uncommented, you will use MRTdynamics. To run the code with MRTdynamics, you have to follow all the comment/uncomment instructions and apply the changes i posted earlier. I will post my working allModels2d.cpp down below.

To your second question: The MRTD2Q9Descriptor contains all the relevant information. The vector set for the d2q9 stencil, the weights and some additional constants used throughout the code (e.g. cs^2). Moreover, the diagonal matrix with the other relaxation times is also defined in there (see: ./src/latticeBoltzmann/mrtLattice.hh). Sadly, these are hardcoded… Getting them to be variable and accessable from your main function would involve some serious changes in the code. The other alternative would be to adjust them in there and recompile all the time…

My (working) allModels2d.cpp with MRTdynamics:

[code=“cpp”]
#include “palabos2D.h”
#include “palabos2D.hh”
#include
#include
#include
#include
#include
#include “poiseuille.h”
#include “poiseuille.hh”
#include “cylinder.h”
#include “cylinder.hh”

using namespace plb;
using namespace plb::descriptors;
using namespace std;

typedef double T;

// Uncomment the two following lines for BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef BGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for so-called incompressible BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef IncBGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for Regularized BGK dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef RegularizedBGKdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for MRT dynamics
#define DESCRIPTOR MRTD2Q9Descriptor
typedef MRTdynamics<T,DESCRIPTOR> BackgroundDynamics;

// Uncomment the two following lines for Entropic dynamics
//#define DESCRIPTOR D2Q9Descriptor
//typedef EntropicDynamics<T,DESCRIPTOR> BackgroundDynamics;

void defineCylinderGeometry( MultiBlockLattice2D<T,DESCRIPTOR>& lattice,
IncomprFlowParam const& parameters )
{
const plint nx = parameters.getNx();
const plint ny = parameters.getNy();

plint cx      = nx/4;
plint cy      = ny/2+ny/10;
plint radius  = cy/4;

createCylinder(lattice, cx, cy, radius);

}

void writeGifs(MultiBlockLattice2D<T,DESCRIPTOR>& lattice, plint iter)
{
const plint imSize = 600;

ImageWriter<T> imageWriter("leeloo");
imageWriter.writeScaledGif(createFileName("u", iter, 6),
                           *computeVelocityNorm(lattice),
                           imSize, imSize );

}

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

global::directories().setOutputDir("./tmp/");

IncomprFlowParam<T> parameters(
        (T) 1e-2,  // uMax
        (T) 300.,  // Re
        100,       // N
        5.,        // lx
        1.         // ly 
);
const T logT     = (T)0.02;
const T imSave   = (T)0.1;
const T maxT     = (T)10.1;

writeLogFile(parameters, "Poiseuille flow");

// Comment these 3 lines for MRT

// MultiBlockLattice2D<T, DESCRIPTOR> lattice (
// parameters.getNx(), parameters.getNy(),
// new BackgroundDynamics(parameters.getOmega()) );

// Uncomment these 4 lines for MRT

// MRTparam<T,DESCRIPTOR> * mrt_param = new MRTparam<T,DESCRIPTOR>(parameters.getOmega());
MultiBlockLattice2D<T, DESCRIPTOR> lattice (
parameters.getNx(), parameters.getNy(),
new BackgroundDynamics(parameters.getOmega()) );
// new BackgroundDynamics(mrt_param) );

OnLatticeBoundaryCondition2D<T,DESCRIPTOR>*
    boundaryCondition = createLocalBoundaryCondition2D<T,DESCRIPTOR>();

defineCylinderGeometry(lattice, parameters);

createPoiseuilleBoundaries(lattice, parameters, *boundaryCondition);
createPoiseuilleInitialValues(lattice, parameters);

// Main loop over time iterations.
for (plint iT=0; iT*parameters.getDeltaT()<maxT; ++iT) {
    if (iT%parameters.nStep(logT)==0) {
        pcout << "step " << iT
              << "; lattice time=" << lattice.getTimeCounter().getTime()
              << "; t=" << iT*parameters.getDeltaT()
              << "; av energy="
              << setprecision(10) << getStoredAverageEnergy<T>(lattice)
              << "; av rho="
              << getStoredAverageDensity<T>(lattice) << endl;
    }

    if (iT%parameters.nStep(imSave)==0) {
        pcout << "Saving Gif ..." << endl;
        writeGifs(lattice, iT);
    }

    // Lattice Boltzmann iteration step.
    lattice.collideAndStream();
}

delete boundaryCondition;

}




Regards,

kk

Many thanks!!! I will take a look at the source, and fiddle with it as necessary for the parameters.