Skip to content

Avoiding spurious peak in residual convergence plot

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Avoiding spurious peak in residual convergence plot

Viewing 7 posts - 1 through 7 (of 7 total)
  • Author
  • #4904

    Dear all

    I am currently extending the existing parallel code available in the openlb examples folder for simulating flow in a 2D lid driven cavity with oscillating boundaries. At Reynolds numbers equal to and exceeding 1000, I observe a spurious peak in the residual convergence plot at a particular time instant (image available at The flow structure immediately shifts to a different pattern at this time instant, which later attains steady state (with converged results). There is a huge discrepancy between these steady state streamline/velocity profiles with the literature. But the profiles at the time instant preceding the spurious peak closely match with the literature, but are not actually converged. For lower Reynolds numbers, there is no such spurious residual convergence peak and the streamline/velocity profiles at the steady state thus match well with the literature.

    So, I was wondering how to avoid this problem at high Reynolds numbers? Both ConstRhoBGKdynamics and MRTdynamics give me similar results. Should I try adopting alternative collision operators implemented in the current version of OpenLB? Or should I try changing the input physical parameters to the code (specially the PhysViscosity or CharPhysVelocity) to avoid this spurious peak. I use the parameters mentioned below in my XML file, which is input to my code (have also tried using a Resolution of 256).


    Below is the code for those interested. The visualization VTK files (which could be loaded in paraview) are shared in this good drive link link

    // Reference : extended from /examples/laminar/cavity2d/ (parallel) code of OLB v1.3-1
    #include "olb2D.h"
    #ifndef OLB_PRECOMPILED // Unless precompiled version is used,
    #include "olb2D.hh"     // include full template code
    #include <vector>
    #include <algorithm>
    #include <cmath>
    #include <iostream>
    using namespace olb;
    using namespace olb::descriptors;
    using namespace olb::graphics;
    using namespace olb::util;
    using namespace std;
    #ifndef M_PI
    #define M_PI 3.14159265358979323846
    typedef double T;
    #define MRT                                                                             // Comment/uncomment this line for BGK/MRT respectively
    #ifdef MRT
    #define DESCRIPTOR MRTD2Q9Descriptor 
    #define DESCRIPTOR D2Q9<>
    T OTime[14] = {0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.7, 0.75, 0.8,                    // time-steps in an oscillating 
      0.85, 0.9, 0.95, 1.0};                                                                // lid-motion to be printed
    void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
                          SuperGeometry2D<T>& superGeometry, int LidType,
                          T ARatio) 
      OstreamManager clout(std::cout,"prepareGeometry");
      clout << "Prepare Geometry ..." << std::endl;
      T eps = converter.getConversionFactorLength();                                        // = deltaX 
      Vector<T,2> extend(T(converter.getCharPhysLength()) + 2*eps, 2*eps);                  // extension length
      Vector<T,2> origin(T() - eps, 
        T(ARatio * converter.getCharPhysLength()) - eps);                                   // origin
      IndicatorCuboid2D<T> lidTop(extend, origin);                                          // center = origin + 0.5*(extension length)
      superGeometry.rename(2, 3, 1, lidTop);                                                // Set material number for lid
      if (LidType == 2)                                                                     // 2 sided lid-motion
        origin[0] = T() - eps;
        origin[1] = T() - eps;
        IndicatorCuboid2D<T> lidBot(extend, origin);                                        // center = origin + 0.5*(extension length)
        superGeometry.rename(2, 4, 1, lidBot);                                              // Set material number for lid
      superGeometry.clean();                                                                // Removes all not needed boundary voxels outside the surface
      superGeometry.innerClean();                                                           // Removes all not needed boundary voxels inside the surface
      clout << "Prepare Geometry ... OK" << std::endl;
    void prepareLattice(UnitConverter<T,DESCRIPTOR> const& converter,
                        SuperLattice2D<T, DESCRIPTOR>& sLattice,
                        Dynamics<T, DESCRIPTOR>& bulkDynamics,
                        sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition,
                        SuperGeometry2D<T>& superGeometry, int LidType) 
      OstreamManager clout(std::cout,"prepareLattice");
      clout << "Prepare Lattice ..." << std::endl;
      const T omega = converter.getLatticeRelaxationFrequency();
      // link lattice with dynamics for collision step
      // Material=0 -->do nothing
      sLattice.defineDynamics(superGeometry, 0, 
        &instances::getNoDynamics<T, DESCRIPTOR>());
      // Material=1 -->bulk dynamics
      sLattice.defineDynamics(superGeometry, 1, &bulkDynamics);
      // Material=2,3 -->bulk dynamics, velocity boundary
      sLattice.defineDynamics(superGeometry, 2, &bulkDynamics);
      sLattice.defineDynamics(superGeometry, 3, &bulkDynamics);
      if (LidType == 2)
        sLattice.defineDynamics(superGeometry, 4, &bulkDynamics);  
      sBoundaryCondition.addVelocityBoundary(superGeometry, 2, omega);
      sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);
      if (LidType == 2)
        sBoundaryCondition.addVelocityBoundary(superGeometry, 4, omega);
      clout << "Prepare Lattice ... OK" << std::endl;
    void setBoundaryValues(UnitConverter<T,DESCRIPTOR> const& converter,
                            SuperLattice2D<T, DESCRIPTOR>& sLattice,
                            int iT, SuperGeometry2D<T>& superGeometry,
                            int LidType, T velRatio, T omegaTop,
                            T omegaBot, char flowtypeTop, 
                            char flowtypeBot) 
      if (iT==0)                                                                                      // initial time step
        AnalyticalConst2D<T,T> rhoF(1);                                                               // set initial values: v = [0,0]
        std::vector<T> velocity(2,T());
        AnalyticalConst2D<T,T> uF(velocity);
        if (LidType == 1)
          auto bulkIndicator = superGeometry.getMaterialIndicator({1, 2, 3});
          sLattice.iniEquilibrium(bulkIndicator, rhoF, uF);
          sLattice.defineRhoU(bulkIndicator, rhoF, uF);
        else if (LidType == 2)
          auto bulkIndicator = superGeometry.getMaterialIndicator({1, 2, 3, 4});
          sLattice.iniEquilibrium(bulkIndicator, rhoF, uF);
          sLattice.defineRhoU(bulkIndicator, rhoF, uF);      
        velocity[0] = converter.getCharLatticeVelocity();                                            // set lid velocity for upper boundary cells
        AnalyticalConst2D<T,T> uTop(velocity);
        sLattice.defineU(superGeometry, 3, uTop);                                                    //define u on material 3
        if (LidType == 2)
          velocity[0] *= velRatio;
          AnalyticalConst2D<T,T> uBot(velocity);
          sLattice.defineU(superGeometry, 4, uBot);                                                  //define u on material 4  
        sLattice.initialize();                                                                       // Make the lattice ready for simulation
        if (flowtypeTop == 'O')
          std::vector<T> velocityTop(2,T());        
          T thetaTop = omegaTop * (T) iT * converter.getPhysDeltaT();                                // oscillating frequency -> angle
          velocityTop[0] = converter.getCharLatticeVelocity() * cos(thetaTop);
          AnalyticalConst2D<T,T> uTop(velocityTop);
          sLattice.defineU(superGeometry, 3, uTop);    
        if (flowtypeBot == 'O' && LidType == 2)
          std::vector<T> velocityBot(2,T());
          T thetaBot = omegaBot * (T) iT * converter.getPhysDeltaT();                                // converter.getPhysTime()
          velocityBot[0] = velRatio *                                                                         
            converter.getCharLatticeVelocity() * cos(thetaBot);
          AnalyticalConst2D<T,T> uBot(velocityBot);
          sLattice.defineU(superGeometry, 4, uBot);   
    void getResults(SuperLattice2D<T, DESCRIPTOR>& sLattice,
                     UnitConverter<T,DESCRIPTOR> const& converter, int iT, 
                     Timer<T>* timer, const T logT, const T maxPhysT, 
                     const T imSave, const T vtkSave, std::string filenameGif,
                     std::string filenameVtk, std::string filenameTec,
                     std::string filenamecVel, std::string filenamexcVel,
                     std::string filenameycVel, const int timerPrintMode,
                     const int timerTimeSteps, const int resolx,              
                     const int resoly, SuperGeometry2D<T>& superGeometry,
                     bool converge, char flowtypeTop, char flowtypeBot,
                     T ARatio, int index = 0) 
      OstreamManager clout(std::cout,"getResults");
      SuperVTMwriter2D<T> vtmWriter(filenameVtk);
      if (iT==0) 
        SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, superGeometry);                            // Write the geometry, cuboid no. and rank no. 
        SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLattice);                                               // as vti file for visualization
        SuperLatticeRank2D<T, DESCRIPTOR> rank(sLattice);
      if (iT%converter.getLatticeTime(logT)==0 || converge) 
        sLattice.getStatistics().print(iT, converter.getPhysTime(iT));                                      // Prints statistics
      if (iT%timerTimeSteps==0 || (converge && flowtypeTop == 'O')) 
        timer->print(iT, timerPrintMode);                                                                   // Prints the passed/remaining time, MLUPS after timerTimesteps
      if ((iT%converter.getLatticeTime(vtkSave)==0 && iT>(int)(0.8 * maxPhysT)) 
        || converge)                                                                                        // Write the VTK files readable in ParaView
        SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity(sLattice, converter);
        SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(sLattice, converter);
    int main(int argc, char* argv[]) 
      // === 1st Step: Initialization ===
      olbInit(&argc, &argv);
      OstreamManager clout(std::cout,"main");
      string fName("cavity2d.xml");                                                                   // XML filename to be read
      XMLreader config(fName);
      std::string olbdir, outputdir;                                                                  // Fix the olb and o/p directories            
      UnitConverter<T,DESCRIPTOR>* converter =                                                        // Read the unit converter related input from the XML file
      converter->print();                                                                              // Prints the converter log as console output
      converter->write("cavity2d");                                                                    // Writes the converter log in a file (.dat format)
      T ARatio = config["Application"]["AspectRatio"].get<T>();                                        // Aspect ratio of the cavity (deep/shallow)
      int N = converter->getResolution() + 1;                                                          // number of voxels in x direction
      int M = (int)ARatio*converter->getResolution() + 1;                                              // number of voxels in y direction
      Timer<T>* timer = createTimer<T>(config, *converter, N*M);
      int LidType = config["Application"]["FlowConfiguration"]["LidType"].get<int>();                   // Flow configuration
      char flowtypeTop = config["Application"]["FlowConfiguration"]["TopLid"]                           // Top lid
      T OscillTimePeriod = config["Application"]["FlowConfiguration"]["TopLid"]
      T omegaTop = 2.0 * M_PI / OscillTimePeriod;
      char flowtypeBot = config["Application"]["FlowConfiguration"]["BottomLid"]                       // Bottom lid
      OscillTimePeriod = config["Application"]["FlowConfiguration"]["BottomLid"]
      T omegaBot = 2.0 * M_PI / OscillTimePeriod;
      T velRatio = config["Application"]["FlowConfiguration"]["BottomLid"]
      T logT = config["Output"]["Log"]["SaveTime"].get<T>();                                            // o/p statistics after (s)
      T imSave = config["Output"]["VisualizationImages"]["SaveTime"].get<T>();                          // save the images after (s)  
      T vtkSave = config["Output"]["VisualizationVTK"]["SaveTime"].get<T>();                            // o/p data to vtk files after (s)
      T maxPhysT = config["Application"]["PhysParameters"]["PhysMaxTime"].get<T>();                     // entire simulation time(s)
      std::string filenameVtk = config["Output"]["VisualizationVTK"]                                    // Filenames for exporting output data
        ["Filename"].get<std::string>();                                                                // name of Vtk files
      std::string filenameGif = config["Output"]["VisualizationImages"]
        ["Filename"].get<std::string>();                                                                // name of Gif images
      std::string filenameTec = config["Output"]
      std::string filenamexcVel = config["Output"]["VisualizationCenterVel"]
      std::string filenameycVel = config["Output"]["VisualizationCenterVel"]
      std::string filenamecVel = config["Output"]["VisualizationCenterVel"]
      std::string errorfilename = "error.txt";
      int timerSkipType = config["Output"]["Timer"]["SkipType"].get<T>();
      int timerPrintMode = config["Output"]["Timer"]["PrintMode"].get<int>();                           // mode for printing time related statistics
      int timerTimeSteps = 1;                                                                           // timerTimeSteps = function of 1/deltaT
      if (timerSkipType == 0) 
        timerTimeSteps = converter->getLatticeTime(config["Output"]["Timer"]
      // === 2rd Step: Prepare Geometry ===
      Vector<T,2> extend(1.0, ARatio);
      extend *= converter->getCharPhysLength();
      Vector<T,2> origin(0,0);
      IndicatorCuboid2D<T> cuboid(extend, origin); 
      CuboidGeometry2D<T> cuboidGeometry(cuboid, 
        converter->getConversionFactorLength(), singleton::mpi().getSize());
      CuboidGeometry2D<T> cuboidGeometry(cuboid, 
        converter->getConversionFactorLength(), 7);
      cuboidGeometry.print();                                                                             // prints cuboid structure statistics
      HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
      SuperGeometry2D<T> superGeometry(cuboidGeometry, loadBalancer, 2);
      prepareGeometry(*converter, superGeometry, LidType, ARatio);
      // === 3rd Step: Prepare Lattice ===
      SuperLattice2D<T, DESCRIPTOR> sLattice(superGeometry);
      Dynamics<T, DESCRIPTOR>* bulkDynamics;
    #ifndef MRT
        bulkDynamics = new ConstRhoBGKdynamics<T, DESCRIPTOR> (
        bulkDynamics = new MRTdynamics<T, DESCRIPTOR> (
      sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryCondition(sLattice);
    #ifndef MRT
          ConstRhoBGKdynamics<T,DESCRIPTOR> > (sBoundaryCondition);
        createInterpBoundaryCondition2D<T,DESCRIPTOR> (sBoundaryCondition);
      prepareLattice(*converter, sLattice, *bulkDynamics, sBoundaryCondition, 
        superGeometry, LidType);
      // === 4th Step: Main Loop with Timer ===
      int interval = converter->getLatticeTime(config["Application"]
        ["ConvergenceCheck"]["interval"].get<T>());                                               // check convergence after interval no. of time-steps
      T epsilon = config["Application"]["ConvergenceCheck"]["residuum"].get<T>();                 // convergence criterion
      bool converge = false;                                                                      // Track convergence
      int convergeLatTime = 0;
      T energy0 = sLattice.getStatistics().getAverageEnergy(); 
      T energy1 = 0.;
      int diffLatTime, incTime, index;
      incTime = converter->getLatticeTime(OscillTimePeriod);                                   // increment time = 1 osc. cycle time-period
      std::vector<int> LatOTime(14,int());   
      for (int k = 0; k < 14; ++k) 
        LatOTime[k] = converter->getLatticeTime(OTime[k]*OscillTimePeriod);
      for (int iT=0; iT <= converter->getLatticeTime(maxPhysT); ++iT) 
        if (converge)                                                                           // attained convergence
          if (iT == convergeLatTime + 1)
            clout << "Simulation converged." << endl;
          if (flowtypeTop == 'O')                                                               // oscillating lid motion
            diffLatTime = iT - convergeLatTime  - 1;                                            // -1 to account for the next iteration step
            auto itr = std::find(LatOTime.begin(), LatOTime.end(), diffLatTime);                // find if diffLatTime exists in LatOTime
            if (itr != LatOTime.cend())                                                          
              index = std::distance(LatOTime.begin(), itr);                                     // find the position of the LatOTime element
              clout << "time-step in the oscillating cycle = " << OTime[index] 
                << "T"<< endl;
              getResults(sLattice, *converter, iT, timer, logT, maxPhysT, imSave, 
                          vtkSave, filenameGif, filenameVtk, filenameTec,
                          filenamecVel, filenamexcVel, filenameycVel,
                          timerPrintMode, timerTimeSteps, N, M, superGeometry, 
                          converge, flowtypeTop, flowtypeBot, ARatio, index);
            if (diffLatTime == incTime) break;                                                  // break the loop once the full cycle is over
          else if (flowtypeTop == 'S')                                                          // steady lid motion
            getResults(sLattice, *converter, iT, timer, logT, maxPhysT, imSave,  
                        vtkSave, filenameGif, filenameVtk, filenameTec,  
                        filenamecVel, filenamexcVel, filenameycVel, timerPrintMode, 
                        timerTimeSteps, N, M, superGeometry, converge, flowtypeTop, 
                        flowtypeBot, ARatio);  
        // === 5th Step: Definition of Initial and Boundary Conditions ===
        setBoundaryValues(*converter, sLattice, iT, superGeometry, LidType, velRatio, 
          omegaTop, omegaBot, flowtypeTop, flowtypeBot);
        // === 6th Step: Collide and Stream Execution ===
        // === 7th Step: Computation and Output of the Results ===
        if (!converge)
          getResults(sLattice, *converter, iT, timer, logT, maxPhysT, imSave, vtkSave, 
                      filenameGif, filenameVtk, filenameTec, filenamecVel, 
                      filenamexcVel, filenameycVel, timerPrintMode, timerTimeSteps,
                      N, M, superGeometry, converge, flowtypeTop, flowtypeBot, ARatio);
        if (iT%interval == 0)                                                                   //check for convergence after interval time-steps
          energy1 = sLattice.getStatistics().getAverageEnergy();
          clout << "step=" << iT << "; Error=" << fabs(energy1 - energy0)/energy1
            << std::endl;     
          if (singleton::mpi().isMainProcessor())
            ofstream myfile;
   (errorfilename.c_str(), ios::app);
            myfile << iT << "\t" << std::setprecision(6) 
              << fabs(energy1 - energy0)/energy1 << endl;
          if (fabs((energy1 - energy0)/energy1) < epsilon && !converge) 
            converge = true;
            convergeLatTime = iT;                                                               // converged time-step
          energy0 = energy1;
      delete converter;
      delete timer;

    Any advice or comments would be deeply appreciated.



    Which spurious peak are you talking about? I dont understand your plot. Best Mathias


    Dear Dr. Mathias

    Thanks for your fast reply.

    I was actually referring to the peak in the residual convergence plot at around 100,000 time steps (plot of residual error vs number of time steps) at the link. I compute the residual error as energy1 - energy0)/energy1 where energy1 and energy0 are sLattice.getStatistics().getAverageEnergy() computed at the current time step and previous time step respectively.

    Though sudden peaks and troughs are possible in such residual convergence plots, the reason I call this as a “Spurious peak” is because I notice a typical change in the flow structure at a time instant around this peak only at high Reynolds numbers. This was not noticed at low Reynolds numbers.



    Maybe that peak is qhysical. For higher Re one can not even expect a steady state solution. Whould be interesting to see the very prot added a case with the same Re but half grid size.


    Dear Mathias

    Thanks for the suggestion. Here is the updated plot for the same Re and different grid sizes. The streamline profiles from 64 x 64 don’t exactly match with that at higher grid sizes due to grid dependence. However, the profiles at higher grid sizes (512 x 512) are still off from the literature.

    I do agree that one cannot expect a steady state solution for higher Re. But Re = 1000 also is not very high enough to expect an unsteady solution state.

    I am now planning to try this by adopting MRTdynamics2, as you state in the code that this uses relaxation times for higher stability. So, is replacing MRTdynamics with MRTdynamics2 in the code and using the same descriptor MRTD2Q9Descriptor the right approach to invoke this?



    The plot confirmes that this is a physical effect. If you plot in physical time instead the time steps you will see that the curves are very similar. So I guess, everything is right. The simulations seem stable, so no need for an MRT but also no reason not to take one.


    Dear Dr. Mathias

    Thanks for your comments and sorry for the delayed reply. I tried both BGK and MRT, along with other boundary conditions – LocalBC, InterpBC and Zou-HeBC. All of them yield similar converged streamline/velocity profiles. So, I believe that this is how the physical solution should be.


Viewing 7 posts - 1 through 7 (of 7 total)
  • You must be logged in to reply to this topic.