diff options
Diffstat (limited to 'examples/laminar/powerLaw2d/powerLaw2d.cpp')
| -rw-r--r-- | examples/laminar/powerLaw2d/powerLaw2d.cpp | 456 |
1 files changed, 456 insertions, 0 deletions
diff --git a/examples/laminar/powerLaw2d/powerLaw2d.cpp b/examples/laminar/powerLaw2d/powerLaw2d.cpp new file mode 100644 index 0000000..fd5e833 --- /dev/null +++ b/examples/laminar/powerLaw2d/powerLaw2d.cpp @@ -0,0 +1,456 @@ +/* Lattice Boltzmann sample, written in C++, using the OpenLB + * library + * + * Copyright (C) 2016 Vojtech Cvrcek, Mathias J. Krause + * 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. + */ + +/* powerLaw2d.cpp: + * This example examines a steady flow of a non-newtonian fluid in a channel. + * At the inlet, a profile for non-newtonian fluid is imposed on the velocity, + * where as the outlet implements an outflow condition grad_x p = 0. + * The power law model is for n=1 and m=charNu in fact the classical poiseuille flow. + * One can validate the error with using functors in void error. + * + * + */ + +#include "olb2D.h" +#include "olb2D.hh" // include full template code + +#include <vector> +#include <cmath> +#include <iostream> +#include <fstream> + +using namespace olb; +using namespace olb::descriptors; +using namespace olb::graphics; +using namespace olb::util; +using namespace std; + +typedef double T; +#define DESCRIPTOR DynOmegaD2Q9Descriptor + +// Parameters for the simulation setup +int N = 40; // resolution of the model +T Re = 10.; // Reynolds number +T tau = 0.8; +T lx = 2.; // channel lenght +T ly = 1.; // channel width +T maxU = 1; // Max velocity +T Tmax = 20; // max. phys. time in s +T Tprint = 1; // Phys time at which the status of the system is print +// set the changes for n and m in powerLawBGKdynamics.h +T n = .2; // parameter in power law model (n=1 Newtonian fluid) +bool bcTypePeriodic = false; //true works only with one core + +const T residuum = 1e-5; // residuum for the convergence check + +void prepareGeometry( PowerLawUnitConverter<T,DESCRIPTOR> const& converter, + SuperGeometry2D<T>& superGeometry ) { + OstreamManager clout( std::cout,"prepareGeometry" ); + clout << "Prepare Geometry ..." << std::endl; + + Vector<T,2> extend( lx, ly ); + Vector<T,2> origin; + + superGeometry.rename( 0,2 ); + superGeometry.rename( 2,1,1,1 ); + + // Set material number for inflow + extend[0] = 1.2*converter.getConversionFactorLength(); + origin[0] = -converter.getConversionFactorLength(); + IndicatorCuboid2D<T> inflow( extend, origin ); + if (bcTypePeriodic) + superGeometry.rename( 1,3,inflow ); + else + superGeometry.rename( 2,3,1,inflow ); + // Set material number for outflow + origin[0] = lx-.5*converter.getConversionFactorLength(); + IndicatorCuboid2D<T> outflow( extend, origin ); + if (bcTypePeriodic) + superGeometry.rename( 1,4,outflow ); + else + superGeometry.rename( 2,4,1,outflow ); + // 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; + return; +} + +// Set up the geometry of the simulation +void prepareLattice( SuperLattice2D<T,DESCRIPTOR>& sLattice, + PowerLawUnitConverter<T,DESCRIPTOR> const& converter, + Dynamics<T, DESCRIPTOR>& bulkDynamics, + Dynamics<T, DESCRIPTOR>& inDynamics, + Dynamics<T, DESCRIPTOR>& outDynamics, + sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition, + SuperGeometry2D<T>& superGeometry ) { + + OstreamManager clout( std::cout,"prepareLattice" ); + clout << "Prepare Lattice ..." << std::endl; + + const T omega = converter.getLatticeRelaxationFrequency(); + + // Material=0 -->do nothing + sLattice.defineDynamics( superGeometry.getMaterialIndicator(0), &instances::getNoDynamics<T, DESCRIPTOR>() ); + + // Material=1 -->bulk dynamics + sLattice.defineDynamics( superGeometry.getMaterialIndicator(1), &bulkDynamics ); + + // Material=2 -->bounce back + sLattice.defineDynamics( superGeometry.getMaterialIndicator(2), &instances::getBounceBack<T, DESCRIPTOR>() ); + + // Material=3 -->bulk dynamics (inflow) + if (bcTypePeriodic) + sLattice.defineDynamics( superGeometry.getMaterialIndicator(3), &inDynamics ); + else { + sLattice.defineDynamics( superGeometry.getMaterialIndicator(3), &bulkDynamics ); + // Setting of the boundary conditions + sBoundaryCondition.addVelocityBoundary( superGeometry, 3, omega ); + } + + // Material=4 -->bulk dynamics (outflow) + if (bcTypePeriodic) + sLattice.defineDynamics( superGeometry.getMaterialIndicator(4), &outDynamics ); + else { + sLattice.defineDynamics( superGeometry.getMaterialIndicator(4), &bulkDynamics ); + // Setting of the boundary conditions + sBoundaryCondition.addPressureBoundary( superGeometry, 4, omega ); + } + clout << "Prepare Lattice ... OK" << std::endl; +} + + +void setBoundaryValues( SuperLattice2D<T, DESCRIPTOR>& sLattice, + PowerLawUnitConverter<T,DESCRIPTOR> const& converter, + int iT, SuperGeometry2D<T>& superGeometry ) { + + OstreamManager clout( std::cout,"setBoundaryValues" ); + + // Set initial and steady boundary conditions + if ( iT==0 ) { + + // Define the analytical solutions for pressure and velocity + T maxVelocity = converter.getCharLatticeVelocity(); + T distance2Wall = converter.getConversionFactorLength()/2.; + + T p0 = converter.getPhysConsistencyCoeff()*pow( converter.getCharPhysVelocity(),n )*pow( ( n + 1. )/n,n )*pow( 2./( ly-distance2Wall*2 ),n + 1. ); + + AnalyticalLinear2D<T,T> rho( converter.getLatticeDensityFromPhysPressure( -p0 ) - 1., 0, converter.getLatticeDensityFromPhysPressure( p0*(lx + distance2Wall*2.)/2. ) ); + + PowerLaw2D<T> u( superGeometry, 3, maxVelocity, distance2Wall, ( n + 1. )/n ); + + // Set the analytical solutions for pressure and velocity + AnalyticalConst2D<T,T> omega0( converter.getLatticeRelaxationFrequency() ); + sLattice.defineField<descriptors::OMEGA>( superGeometry, 1, omega0 ); + sLattice.defineField<descriptors::OMEGA>( superGeometry, 3, omega0 ); + sLattice.defineField<descriptors::OMEGA>( superGeometry, 4, omega0 ); + + // Set the analytical solutions for pressure and velocity + // Initialize all values of distribution functions to their local equilibrium + + sLattice.defineRhoU( superGeometry, 1, rho, u ); + sLattice.iniEquilibrium( superGeometry, 1, rho, u ); + + sLattice.iniEquilibrium( superGeometry, 3, rho, u ); + sLattice.defineRhoU( superGeometry, 3, rho, u ); + + sLattice.iniEquilibrium( superGeometry, 4, rho, u ); + sLattice.defineRhoU( superGeometry, 4, rho, u ); + + // Make the lattice ready for simulation + sLattice.initialize(); + } +} + +// Compute error norms +void error( SuperGeometry2D<T>& superGeometry, + SuperLattice2D<T, DESCRIPTOR>& sLattice, + PowerLawUnitConverter<T,DESCRIPTOR> const& converter, + Dynamics<T, DESCRIPTOR>& bulkDynamics ) { + OstreamManager clout( std::cout,"error" ); + + int input[1] = { }; + T result[1] = { }; + + T distance2Wall = converter.getConversionFactorLength()/2.; + + PowerLaw2D<T> uSol( superGeometry,3,converter.getCharPhysVelocity(),distance2Wall,( n + 1. )/n ); + SuperLatticePhysVelocity2D<T,DESCRIPTOR> u( sLattice,converter ); + + T p0 = converter.getPhysConsistencyCoeff()*pow( converter.getCharPhysVelocity(),n )*pow( ( n + 1. )/n,n )*pow( 2./( ly-distance2Wall*2 ),n + 1. ); + AnalyticalLinear2D<T,T> pressureSol( -p0, 0, p0*(lx + distance2Wall*2.)/2. ); + SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure( sLattice,converter ); + + auto indicatorF = superGeometry.getMaterialIndicator(1); + + // velocity error + SuperAbsoluteErrorL1Norm2D<T> absVelocityErrorNormL1(u, uSol, indicatorF); + absVelocityErrorNormL1(result, input); + clout << "velocity-L1-error(abs)=" << result[0]; + SuperRelativeErrorL1Norm2D<T> relVelocityErrorNormL1(u, uSol, indicatorF); + relVelocityErrorNormL1(result, input); + clout << "; velocity-L1-error(rel)=" << result[0] << std::endl; + + SuperAbsoluteErrorL2Norm2D<T> absVelocityErrorNormL2(u, uSol, indicatorF); + absVelocityErrorNormL2(result, input); + clout << "velocity-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm2D<T> relVelocityErrorNormL2(u, uSol, indicatorF); + relVelocityErrorNormL2(result, input); + clout << "; velocity-L2-error(rel)=" << result[0] << std::endl; + + SuperAbsoluteErrorLinfNorm2D<T> absVelocityErrorNormLinf(u, uSol, indicatorF); + absVelocityErrorNormLinf(result, input); + clout << "velocity-Linf-error(abs)=" << result[0]; + SuperRelativeErrorLinfNorm2D<T> relVelocityErrorNormLinf(u, uSol, indicatorF); + relVelocityErrorNormLinf(result, input); + clout << "; velocity-Linf-error(rel)=" << result[0] << std::endl; + + // pressure error + SuperAbsoluteErrorL1Norm2D<T> absPressureErrorNormL1(pressure, pressureSol, indicatorF); + absPressureErrorNormL1(result, input); + clout << "pressure-L1-error(abs)=" << result[0]; + SuperRelativeErrorL1Norm2D<T> relPressureErrorNormL1(pressure, pressureSol, indicatorF); + relPressureErrorNormL1(result, input); + clout << "; pressure-L1-error(rel)=" << result[0] << std::endl; + + SuperAbsoluteErrorL2Norm2D<T> absPressureErrorNormL2(pressure, pressureSol, indicatorF); + absPressureErrorNormL2(result, input); + clout << "pressure-L2-error(abs)=" << result[0]; + SuperRelativeErrorL2Norm2D<T> relPressureErrorNormL2(pressure, pressureSol, indicatorF); + relPressureErrorNormL2(result, input); + clout << "; pressure-L2-error(rel)=" << result[0] << std::endl; + + SuperAbsoluteErrorLinfNorm2D<T> absPressureErrorNormLinf(pressure, pressureSol, indicatorF); + absPressureErrorNormLinf(result, input); + clout << "pressure-Linf-error(abs)=" << result[0]; + SuperRelativeErrorLinfNorm2D<T> relPressureErrorNormLinf(pressure, pressureSol, indicatorF); + relPressureErrorNormLinf(result, input); + clout << "; pressure-Linf-error(rel)=" << result[0] << std::endl; +} + +// Output to console and files +void getResults( SuperLattice2D<T, DESCRIPTOR>& sLattice, + Dynamics<T, DESCRIPTOR>& bulkDynamics, + PowerLawUnitConverter<T,DESCRIPTOR> const& converter, int iT, + SuperGeometry2D<T>& superGeometry, Timer<double>& timer ) { + OstreamManager clout( std::cout,"getResults" ); + + SuperVTMwriter2D<T> vtmWriter( "powerLaw2d" ); + SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity( sLattice, converter ); + SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure( sLattice, converter ); + + vtmWriter.addFunctor( velocity ); + vtmWriter.addFunctor( pressure ); + + if ( iT==0 ) { + SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice ); + SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice,superGeometry ); + SuperLatticeRank2D<T, DESCRIPTOR> rank( sLattice ); + vtmWriter.write( geometry ); + vtmWriter.write( cuboid ); + vtmWriter.write( rank ); + vtmWriter.createMasterFile(); + } + + if ( iT%converter.getLatticeTime( Tprint )==0 ) { + vtmWriter.write( iT ); + + SuperEuklidNorm2D<T, DESCRIPTOR> normVel( velocity ); + BlockReduction2D2D<T> planeReduction( normVel, 600, BlockDataSyncMode::ReduceOnly ); + // write output of velocity as JPEG + heatmap::write(planeReduction, iT); + } + + // Writes output on the console + if ( iT%converter.getLatticeTime( Tprint )==0 ) { + timer.update( iT ); + timer.printStep(); + sLattice.getStatistics().print( iT,converter.getPhysTime( iT ) ); + error( superGeometry, sLattice, converter, bulkDynamics ); + } + return; +} + + +int main( int argc, char* argv[] ) { + + // === 1st Step: Initialization === + olbInit( &argc, &argv ); + singleton::directories().setOutputDir( "./tmp/" ); + OstreamManager clout( std::cout,"main" ); + + /* + if ( argc > 1 ) { + N = atoi( argv[1] ); + } + if ( argc > 2 ) { + n = atof( argv[2] ); + } + */ + + singleton::directories().setOutputDir( "./tmp/" ); + + XMLreader config("input.xml"); + config["geometry"]["N"].read(N); + config["geometry"]["lx"].read(lx); + config["geometry"]["ly"].read(ly); + config["dynamics"]["maxU"].read(maxU); + config["dynamics"]["Re"].read(Re); + config["dynamics"]["tau"].read(tau); + config["dynamics"]["n"].read(n); + config["time"]["Tmax"].read(Tmax); + config["time"]["Tprint"].read(Tprint); + + /* + T physDeltaX = (ly/N); + T m = ly*maxU*pow( maxU/(2*ly),1-n )/Re; + T physViscosity = m*pow( maxU/(2*ly),n-1 ); + T physDeltaT = (tau - 0.5) / DESCRIPTOR::invCs2 * pow(physDeltaX,2) / physViscosity; + PowerLawUnitConverter<T, DESCRIPTOR> const converter( + (T) physDeltaX, + (T) physDeltaT, + (T) ly, + (T) maxU, + (T) m, + (T) n, // power-law index + (T) 1.0 // physDensity: physical density in __kg / m^3__ + ); + */ + PowerLawUnitConverterFrom_Resolution_RelaxationTime_Reynolds_PLindex<T, DESCRIPTOR> const converter( + int {N}, // resolution: number of voxels per charPhysL + (T) tau, // latticeRelaxationTime: relaxation time, have to be greater than 0.5! + (T) ly, // charPhysLength: reference length of simulation geometry + (T) maxU, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__ + (T) Re, // Reynolds number + (T) n, // power-law index + (T) 1.0 // physDensity: physical density in __kg / m^3__ + ); + // Prints the converter log as console output + converter.print(); + // Writes the converter log in a file + converter.write("powerLaw2d"); + + + // === 2rd Step: Prepare Geometry === + // Instantiation of a cuboidGeometry with weights + + Vector<T,2> extend( lx, ly ); + Vector<T,2> origin; + 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 + + // Periodic boundaries in x-direction + if (bcTypePeriodic) + cuboidGeometry.setPeriodicity( true, false ); + + //cuboidGeometry.printExtended(); + + HeuristicLoadBalancer<T> loadBalancer( cuboidGeometry ); + SuperGeometry2D<T> superGeometry( cuboidGeometry, loadBalancer, 2 ); + prepareGeometry( converter, superGeometry ); + + // === 3rd Step: Prepare Lattice === + SuperLattice2D<T, DESCRIPTOR> sLattice( superGeometry ); + + T distance2Wall = converter.getConversionFactorLength()/2.; + T p0 = converter.getPhysConsistencyCoeff()*pow( converter.getCharPhysVelocity(),n )*pow( ( n + 1. )/n,n )*pow( 2./( ly-distance2Wall*2 ),n + 1. ); + + clout << "Dimensionalized version-1." << std::endl; + PowerLawBGKdynamics<T, DESCRIPTOR> bulkDynamics( converter.getLatticeRelaxationFrequency(), instances::getBulkMomenta<T, DESCRIPTOR>(), converter.getLatticeConsistencyCoeff(), n ); + + PeriodicPressureDynamics<T, DESCRIPTOR, PowerLawBGKdynamics<T,DESCRIPTOR>> outDynamics( bulkDynamics,converter.getLatticeDensityFromPhysPressure( p0*(lx + distance2Wall*2.))-1,1,0); + PeriodicPressureDynamics<T, DESCRIPTOR, PowerLawBGKdynamics<T,DESCRIPTOR>> inDynamics( bulkDynamics,-converter.getLatticeDensityFromPhysPressure( p0*(lx + distance2Wall*2. ))+1,-1,0); + std::cout << -converter.getLatticeDensityFromPhysPressure( p0 )+1 << std::endl; + + sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> sBoundaryCondition( sLattice ); + createLocalBoundaryCondition2D<T, DESCRIPTOR, PowerLawBGKdynamics<T,DESCRIPTOR> > ( sBoundaryCondition ); + + prepareLattice( sLattice, converter, bulkDynamics, inDynamics, outDynamics, sBoundaryCondition, superGeometry ); + + // === 4th Step: Main Loop with Timer === + Timer<double> timer( converter.getLatticeTime( Tmax ), superGeometry.getStatistics().getNvoxel() ); + util::ValueTracer<T> converge( converter.getLatticeTime( 0.1*Tprint ), residuum ); + timer.start(); + + for ( int iT=0; iT<converter.getLatticeTime( Tmax ); ++iT ) { + if ( converge.hasConverged() ) { + clout << "Simulation converged." << endl; + getResults( sLattice, bulkDynamics, converter, iT, superGeometry, timer ); + + break; + } + + // === 5th Step: Definition of Initial and Boundary Conditions === + setBoundaryValues( sLattice, converter, iT, superGeometry ); + // === 7th Step: Computation and Output of the Results === + getResults( sLattice, bulkDynamics, converter, iT, superGeometry, timer ); + converge.takeValue( sLattice.getStatistics().getAverageEnergy(), true ); + // === 6th Step: Collide and Stream Execution === + sLattice.collideAndStream(); + + if (bcTypePeriodic) + sLattice.stripeOffDensityOffset ( sLattice.getStatistics().getAverageRho()-(T)1 ); + } + timer.stop(); + timer.printSummary(); + + // === 7th Step: Gnuplot === + Gnuplot<T> gplot( "centerVelocity" ); + T Ly = converter.getLatticeLength( ly ); + for ( int iY=0; iY<=Ly; ++iY ) { + T dx = 1. / T(converter.getResolution()); + // const T maxVelocity = converter.getPhysVelocity( converter.getCharLatticeVelocity() ); + T point[2]= {T(),T()}; + point[0] = .9*lx; + point[1] = ( T )iY/Ly; + // const T radius = ly/2.; + std::vector<T> axisPoint( 2,T() ); + axisPoint[0] = lx/2.; + axisPoint[1] = ly/2.; + std::vector<T> axisDirection( 2,T() ); + axisDirection[0] = 1; + axisDirection[1] = 0; + T distance2Wall = converter.getConversionFactorLength()/2.; + PowerLaw2D<T> uSol( superGeometry,3,converter.getCharPhysVelocity(),distance2Wall,( n + 1. )/n ); + //Poiseuille2D<T> uSol( axisPoint, axisDirection, maxVelocity, radius ); // <--- + T analytical[2] = {T(),T()}; + uSol( analytical,point ); + SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity( sLattice, converter ); + AnalyticalFfromSuperF2D<T> intpolateVelocity( velocity, true ); + T numerical[2] = {T(),T()}; + intpolateVelocity( numerical,point ); + gplot.setData( iY*dx, {analytical[0],numerical[0]}, {"analytical","numerical"} ); + } + // Create PNG file + gplot.writePNG(); +} |