Simulation giving out: absoluteErrorL2(line)=0.0981649 and relativeErrorL2(line)
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Simulation giving out: absoluteErrorL2(line)=0.0981649 and relativeErrorL2(line)
- This topic has 7 replies, 2 voices, and was last updated 10 months, 2 weeks ago by Adrian.
-
AuthorPosts
-
April 2, 2024 at 9:32 am #8438atanuchaudhury1994@gmail.comParticipant
Dear OpenLB Team,
I actually tried to trun the Cavity2D code in Laminar Example into a Turbulent code. But during the simulation its showing an Error: “[getResults] absoluteErrorL2(line)=0.0981633; relativeErrorL2(line)=0.989684”. I ran it using single core CPU. Please help me to solve this. I am providing the Cavity2D code that I changed below:
/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2006 – 2012 Mathias J. Krause, Jonas Fietz,
* Jonas Latt, Jonas Kratzke
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/>
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*//* cavity2d.cpp:
* This example illustrates a flow in a cuboid, lid-driven cavity.
* It also shows how to use the XML parameter files and has an
* example description file for OpenGPI. This version is for parallel
* use. A version for sequential use is also available.
*/#include “olb2D.h”
#include “olb2D.hh”using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<>;
using BulkDynamics = ConstRhoBGKdynamics<T,DESCRIPTOR>;#define USE_SMAGORINSKY //default
void prepareGeometry( UnitConverter<T,DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry )
{
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();
Vector<T,2> extend( T( 1 ) + 2*eps, 2*eps );
Vector<T,2> origin( T() – eps, T( 1 ) – eps );
IndicatorCuboid2D<T> lid( extend, origin );
// Set material number for lid
superGeometry.rename( 2,3,1,lid );// 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.getStatistics().print();
clout << “Prepare Geometry … OK” << std::endl;
}void prepareLattice( UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice<T, DESCRIPTOR>& sLattice,
SuperGeometry<T,2>& superGeometry )
{OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;const T omega = converter.getLatticeRelaxationFrequency();
// Material=1 –>bulk dynamics
//sLattice.defineDynamics<BulkDynamics>(superGeometry, 1);
#ifdef USE_SMAGORINSKY
using BulkDynamics = SmagorinskyBGKdynamics<T,DESCRIPTOR>;
sLattice.defineDynamics<BulkDynamics>(superGeometry, 1);
//sLattice.setParameter<collision::LES::Smagorinsky>(T(0.15));
sLattice.setParameter<collision::LES::Smagorinsky>(0.15);
#endif// Material=2,3 –>bulk dynamics, velocity boundary
setInterpolatedVelocityBoundary<T,DESCRIPTOR,BulkDynamics>(sLattice, omega, superGeometry, 2);
setInterpolatedVelocityBoundary<T,DESCRIPTOR,BulkDynamics>(sLattice, omega, superGeometry, 3);sLattice.setParameter<descriptors::OMEGA>(omega);
clout << “Prepare Lattice … OK” << std::endl;
}void setBoundaryValues( UnitConverter<T,DESCRIPTOR> const& converter,
SuperLattice<T, DESCRIPTOR>& sLattice,
int iT, SuperGeometry<T,2>& superGeometry )
{if ( iT==0 ) {
// set initial values: v = [0,0]
AnalyticalConst2D<T,T> rhoF( 1 );
std::vector<T> velocity( 2,T() );
AnalyticalConst2D<T,T> uF( velocity );auto bulkIndicator = superGeometry.getMaterialIndicator({1, 2, 3});
sLattice.iniEquilibrium( bulkIndicator, rhoF, uF );
sLattice.defineRhoU( bulkIndicator, rhoF, uF );// set non-zero velocity for upper boundary cells
velocity[0] = converter.getCharLatticeVelocity();
AnalyticalConst2D<T,T> u( velocity );
sLattice.defineU( superGeometry, 3, u );// Make the lattice ready for simulation
sLattice.initialize();
}
}void getResults( SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T,DESCRIPTOR> const& converter, std::size_t iT, util::Timer<T>* timer,
const T logT, const T maxPhysT, const T imSave, const T vtkSave, const T gnuplotSave,
std::string filenameGif, std::string filenameVtk, std::string filenameGnuplot,
const int timerPrintMode,
SuperGeometry<T,2>& superGeometry, bool converged )
{OstreamManager clout( std::cout,”getResults” );
SuperVTMwriter2D<T> vtmWriter( filenameVtk );
if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice );
SuperLatticeRank2D<T, DESCRIPTOR> rank( sLattice );
SuperLatticeDiscreteNormal2D<T, DESCRIPTOR> discreteNormal( sLattice, superGeometry, superGeometry.getMaterialIndicator({2, 3}) );
SuperLatticeDiscreteNormalType2D<T, DESCRIPTOR> discreteNormalType( sLattice, superGeometry, superGeometry.getMaterialIndicator({2, 3}) );vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
vtmWriter.write( discreteNormal );
vtmWriter.write( discreteNormalType );
vtmWriter.createMasterFile();
}// Get statistics
if ( iT%converter.getLatticeTime( logT )==0 || converged ) {
sLattice.getStatistics().print( iT, converter.getPhysTime( iT ) );
timer->print( iT,timerPrintMode );
}// Writes the VTK files
if ( ( iT%converter.getLatticeTime( vtkSave )==0 && iT>0 ) || converged ) {
sLattice.setProcessingContext(ProcessingContext::Evaluation);SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity( sLattice, converter );
SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure( sLattice, converter );vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( pressure );
vtmWriter.write( iT );
}// Writes the Gif files
if ( ( iT%converter.getLatticeTime( imSave )==0 && iT>0 ) || converged ) {
SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity( sLattice, converter );
SuperEuklidNorm2D<T,DESCRIPTOR> normVel( velocity );
BlockReduction2D2D<T> planeReduction( normVel, 600, BlockDataSyncMode::ReduceOnly );
// write output of velocity as JPEG
heatmap::write(planeReduction, iT);
}// Output for x-velocity along y-position at the last time step
//if ( iT == converter.getLatticeTime( maxPhysT ) || converged ) {
if ( (iT%converter.getLatticeTime( gnuplotSave ) == 0 && iT>0 ) || converged ) {
// Gives access to velocity information on lattice
SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocityField( sLattice, converter );
// Interpolation functor with velocityField information
AnalyticalFfromSuperF2D<T> interpolation( velocityField, true, 1 );Vector<T,17> y_coord( {128, 125, 124, 123, 122, 109, 94, 79, 64, 58, 36, 22, 13, 9, 8, 7, 0} );
// Ghia, Ghia and Shin, 1982: “High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method”; Table 1
Vector<T,17> vel_ghia_RE1000( { 1.0, 0.65928, 0.57492, 0.51117, 0.46604,
0.33304, 0.18719, 0.05702,-0.06080,-0.10648,
-0.27805,-0.38289,-0.29730,-0.22220,-0.20196,
-0.18109, 0.0
} );
Vector<T,17> vel_ghia_RE100( {1.0, 0.84123, 0.78871, 0.73722, 0.68717,
0.23151, 0.00332,-0.13641,-0.20581,-0.21090,
-0.15662,-0.10150,-0.06434,-0.04775,-0.04192,
-0.03717, 0.0
} );
Vector<T,17> vel_simulation;// Gnuplot interface to create plots
Gnuplot<T> gplot( filenameGnuplot + “_iT” + std::to_string(iT) );
// Define comparison values
Vector<T,17> comparison = vel_ghia_RE1000;for ( int nY = 0; nY < 17; ++nY ) {
// 17 data points evenly distributed between 0 and 1 (height)
T position[2] = {0.5, y_coord[nY]/T(128)};
T velocity[2] = {T(), T()};
// Interpolate velocityField at “position” and save it in “velocity”
interpolation( velocity, position );
// Save value of velocity (in x-direction) in “vel_simulation” for every position “nY”
vel_simulation[nY] = velocity[0];
// Set data for plot output
gplot.setData( position[1], {vel_simulation[nY],comparison[nY]}, {“simulated”,”Ghia”} );
}
// Create PNG file
gplot.writePNG();
// Console output with results
clout << “absoluteErrorL2(line)=” << norm(vel_simulation – comparison) / 17. << “; relativeErrorL2(line)=” << norm(vel_simulation – comparison) / norm(comparison) << std::endl;
}
}int main( int argc, char* argv[] )
{// === 1st Step: Initialization ===
olbInit( &argc, &argv );
OstreamManager clout( std::cout,”main” );std::string fName( “cavity2d.xml” );
XMLreader config( fName );std::string olbdir, outputdir;
config[“Application”][“OlbDir”].read( olbdir );
config[“Output”][“OutputDir”].read( outputdir );
singleton::directories().setOlbDir( olbdir );
singleton::directories().setOutputDir( outputdir );UnitConverter<T,DESCRIPTOR>* converter = createUnitConverter<T,DESCRIPTOR>( config );
// Prints the converter log as console output
converter->print();
// Writes the converter log in a file
converter->write(“cavity2d”);int N = converter->getLatticeLength(1) + 1; // number of voxels in x,y,z direction
util::Timer<T>* timer = util::createTimer<T>( config, *converter, N*N );T logT = config[“Output”][“Log”][“SaveTime”].get<T>();
T imSave = config[“Output”][“VisualizationImages”][“SaveTime”].get<T>();
T vtkSave = config[“Output”][“VisualizationVTK”][“SaveTime”].get<T>();
T gnuplotSave = config[“Output”][“VisualizationGnuplot”][“SaveTime”].get<T>();
T maxPhysT = config[“Application”][“PhysParameters”][“PhysMaxTime”].get<T>();
int timerPrintMode = config[“Output”][“Timer”][“PrintMode”].get<int>();std::string filenameGif = config[“Output”][“VisualizationImages”][“Filename”].get<std::string>();
std::string filenameVtk = config[“Output”][“VisualizationVTK”][“Filename”].get<std::string>();
std::string filenameGnuplot = config[“Output”][“VisualizationGnuplot”][“Filename”].get<std::string>();// === 2nd Step: Prepare Geometry ===
Vector<T,2> extend( 1,1 );
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(), 1 );
#endifcuboidGeometry.print();
HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry );
SuperGeometry<T,2> superGeometry( cuboidGeometry, loadBalancer );
prepareGeometry( *converter, superGeometry );// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLattice( superGeometry );
prepareLattice( *converter, sLattice, superGeometry );
// === 4th Step: Main Loop with Timer ===
int interval = converter->getLatticeTime( 1 /*config[“Application”][“ConvergenceCheck”][“interval”].get<T>()*/ );
T epsilon = 1e-3;
util::ValueTracer<T> converge( interval, epsilon );timer->start();
for ( std::size_t iT=0; iT <= converter->getLatticeTime( maxPhysT ); ++iT ) {
if ( converge.hasConverged() ) {
clout << “Simulation converged.” << std::endl;
getResults( sLattice, *converter, iT, timer, logT, maxPhysT, imSave, vtkSave, gnuplotSave, filenameGif, filenameVtk,
filenameGnuplot, timerPrintMode, superGeometry, converge.hasConverged() );
break;
}
// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( *converter, sLattice, iT, superGeometry );
// === 6th Step: Collide and Stream Execution ===
sLattice.collideAndStream();
// === 7th Step: Computation and Output of the Results ===
getResults( sLattice, *converter, iT, timer, logT, maxPhysT, imSave, vtkSave, gnuplotSave, filenameGif, filenameVtk,
filenameGnuplot, timerPrintMode, superGeometry, converge.hasConverged() );
converge.takeValue( sLattice.getStatistics().getAverageEnergy(), true );
}
timer->stop();
timer->printSummary();
delete converter;
delete timer;
}April 2, 2024 at 9:36 am #8439AdrianKeymasterThis is not an error – the same message (with different results) is also printed for the unmodified lid driven cavity case. As the description suggests it is not a program error but the numerical error between simulation and reference solution.
April 2, 2024 at 12:23 pm #8441atanuchaudhury1994@gmail.comParticipantSo can you suggest what changes I have to make to get the required output. As I tired with different Reynolds number, its giving the same result. Please tell me how to come out of this situation.
Thank youApril 2, 2024 at 1:04 pm #8442AdrianKeymasterWhat do you mean by “required output”? Going by your first message I guess you just want to change the physical parameters which you can achieve by changing the parameters given to the unit converter.
April 2, 2024 at 1:27 pm #8443atanuchaudhury1994@gmail.comParticipantI mean by the message that, the Cavity2D example is in Laminar Form and I want to change it in Turbulent Form and so I made some changes in the code:
1st Change:
using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<>;
using BulkDynamics = ConstRhoBGKdynamics<T,DESCRIPTOR>;#define USE_SMAGORINSKY //default
2nd Change:
// Material=1 –>bulk dynamics
//sLattice.defineDynamics<BulkDynamics>(superGeometry, 1);
#ifdef USE_SMAGORINSKY
using BulkDynamics = SmagorinskyBGKdynamics<T,DESCRIPTOR>;
sLattice.defineDynamics<BulkDynamics>(superGeometry, 1);
//sLattice.setParameter<collision::LES::Smagorinsky>(T(0.15));
sLattice.setParameter<collision::LES::Smagorinsky>(0.15);
#endifand so I am an error message:
[ValueTracer] average=0.000966059; stdDev/average=0.00105275
[main] Simulation converged.
[LatticeStatistics] step=133153; t=53.2612; uMax=0.16; avEnergy=0.000968583; avRho=0.988066
[Timer] step=133153; percent=26.6306; passedTime=3119.43; remTime=8594.28; MLUPs=6.68937
We could not find a gnuplot distribution at your system.
We still write the data files s.t. you can plot the data yourself.
We could not find a gnuplot distribution at your system.
We still write the data files s.t. you can plot the data yourself.
[getResults] absoluteErrorL2(line)=0.0278967; relativeErrorL2(line)=0.281255
[Timer]
[Timer] —————-Summary:Timer—————-
[Timer] measured time (rt) : 3120.689s
[Timer] measured time (cpu): 3104.719s
[Timer] average MLUPs : 6.861
[Timer] average MLUPps: 6.861
[Timer] ———————————————So is there any solution for that? Or there is there any other way to change Cavity2D in Turbulent Form. If there please suggest.
April 2, 2024 at 1:50 pm #8444AdrianKeymasterSolution for what? Did you actually change the physical parameters or just the collision model?
I repeat: The error here is the difference between simulation and expected reference solution, not a program error. If you mean the gnuplot warning: just install gnuplot if you want to directly generate result plots for your simulation.
April 2, 2024 at 2:02 pm #8445atanuchaudhury1994@gmail.comParticipantYes Sir I tried to change the Physical Parameter of the model, from Laminar to Turbulent Form. And so the simulation is showing converged after completing 28% only rather completing 100% as given in the previous reply.
April 2, 2024 at 3:42 pm #8446AdrianKeymasterThe convergence criterion of the laminar version quite likely is not what you want for a turbulent setup. You need to decide yourself which criteria is suitable for your specific situation. I.e. you need to actually decide on your exact simulation model before trying to adapt any code.
-
AuthorPosts
- You must be logged in to reply to this topic.