Avoiding spurious peak in residual convergence plot

[sthavishtha]

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 https://ibb.co/cCQt8GF). 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).

    <Discretization>
      <Resolution>128</Resolution>
     <LatticeRelaxationTime>0.5384</LatticeRelaxationTime>
    </Discretization>
    <PhysParameters>
      <PhysMaxTime>250</PhysMaxTime>
      <CharPhysLength>1</CharPhysLength>
      <CharPhysVelocity>1</CharPhysVelocity>
      <PhysViscosity>0.001</PhysViscosity>
      <PhysDensity>1</PhysDensity>
    </PhysParameters>

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

typedef double T;

#define MRT                                                                             // Comment/uncomment this line for BGK/MRT respectively
#ifdef MRT
#define DESCRIPTOR MRTD2Q9Descriptor 
#else
#define DESCRIPTOR D2Q9<>
#endif

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;

  superGeometry.rename(0,2);
  superGeometry.rename(2,1,1,1);
  superGeometry.clean();

  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
  superGeometry.checkForErrors();   
  superGeometry.getStatistics().print();
  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
  }

  else
  {
    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);
    vtmWriter.write(geometry);
    vtmWriter.write(cuboid);
    vtmWriter.write(rank);
    vtmWriter.createMasterFile();
  }
  
  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);
    vtmWriter.addFunctor(velocity);
    vtmWriter.addFunctor(pressure);
    vtmWriter.write(iT);
  }
}

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            
  config["Application"]["OlbDir"].read(olbdir);
  config["Output"]["OutputDir"].read(outputdir);
  singleton::directories().setOlbDir(olbdir);
  singleton::directories().setOutputDir(outputdir);

  UnitConverter<T,DESCRIPTOR>* converter =                                                        // Read the unit converter related input from the XML file
    createUnitConverter<T,DESCRIPTOR>(config);  
  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
    ["FlowType"].get<char>();
  T OscillTimePeriod = config["Application"]["FlowConfiguration"]["TopLid"]
    ["OscillTimePeriod"].get<T>();
  T omegaTop = 2.0 * M_PI / OscillTimePeriod;
  
  char flowtypeBot = config["Application"]["FlowConfiguration"]["BottomLid"]                       // Bottom lid
    ["FlowType"].get<char>();
  OscillTimePeriod = config["Application"]["FlowConfiguration"]["BottomLid"]
    ["OscillTimePeriod"].get<T>();
  T omegaBot = 2.0 * M_PI / OscillTimePeriod;
  T velRatio = config["Application"]["FlowConfiguration"]["BottomLid"]
    ["VelocityRatio"].get<T>();

  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"]
    ["VisualizationTecplot"].get<std::string>();  
  std::string filenamexcVel = config["Output"]["VisualizationCenterVel"]
    ["xCenterlineVelocity"].get<std::string>();
  std::string filenameycVel = config["Output"]["VisualizationCenterVel"]
    ["yCenterlineVelocity"].get<std::string>();
  std::string filenamecVel = config["Output"]["VisualizationCenterVel"]
    ["CenterVelocity"].get<std::string>();
  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"]
      ["PhysTime"].get<T>());
  }
  else
    config["Output"]["Timer"]["TimeSteps"].read(timerTimeSteps);

  // === 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); 

#ifdef PARALLEL_MODE_MPI
  CuboidGeometry2D<T> cuboidGeometry(cuboid, 
    converter->getConversionFactorLength(), singleton::mpi().getSize());
#else
  CuboidGeometry2D<T> cuboidGeometry(cuboid, 
    converter->getConversionFactorLength(), 7);
#endif

  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> (
      converter->getLatticeRelaxationFrequency(), 
      instances::getBulkMomenta<T,DESCRIPTOR>());  
#else
    bulkDynamics = new MRTdynamics<T, DESCRIPTOR> (
      converter->getLatticeRelaxationFrequency(), 
      instances::getBulkMomenta<T,DESCRIPTOR>());      
#endif

  sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryCondition(sLattice);

#ifndef MRT
    createInterpBoundaryCondition2D<T,DESCRIPTOR,
      ConstRhoBGKdynamics<T,DESCRIPTOR> > (sBoundaryCondition);
#else
    createInterpBoundaryCondition2D<T,DESCRIPTOR> (sBoundaryCondition);
#endif

  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;

  timer->start();
  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);  
        break;
      }      
    }

    // === 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 ===
    sLattice.collideAndStream();
    
    // === 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;
        myfile.open (errorfilename.c_str(), ios::app);
        myfile << iT << "\t" << std::setprecision(6) 
          << fabs(energy1 - energy0)/energy1 << endl;
        myfile.close();          
      }       

      if (fabs((energy1 - energy0)/energy1) < epsilon && !converge) 
      {
        converge = true;
        convergeLatTime = iT;                                                               // converged time-step
      }

      energy0 = energy1;
    }
  }

  timer->stop();
  timer->printSummary();
  delete converter;
  delete timer;
}

Any advice or comments would be deeply appreciated.

Regards
Sthavishtha

[mathias]

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

[sthavishtha]

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.

Regards
Sthavishtha

[mathias]

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.

[sthavishtha]

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?

Regards
Sthavishtha

[mathias]

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.

[sthavishtha]

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.

Regards
Sthavishtha

[Author]

Posts