Commit 444caaca authored by Christian Godenschwager's avatar Christian Godenschwager
Browse files

Use _r for all real_t literals

parent 4c98649d
......@@ -211,7 +211,7 @@ int main( int argc, char * argv[] )
meshWorkloadMemory.setOutsideCellWorkload(1);
bfc.setWorkloadMemorySUIDAssignmentFunction( meshWorkloadMemory );
bfc.setPeriodicity( Vector3<bool>(true) );
bfc.setRefinementSelectionFunction( makeRefinementSelection( distanceOctree, numLevels - uint_t(1), dx, dx * real_t(1) ) );
bfc.setRefinementSelectionFunction( makeRefinementSelection( distanceOctree, numLevels - uint_t(1), dx, dx * 1_r ) );
auto structuredBlockforest = bfc.createStructuredBlockForest( blockSize );
......@@ -225,7 +225,7 @@ int main( int argc, char * argv[] )
static const uint_t NUM_GHOSTLAYERS = 4;
BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( structuredBlockforest, "pdf field", latticeModel, Vector3<real_t>(0), real_t(1), NUM_GHOSTLAYERS, field::fzyx );
BlockDataID pdfFieldId = lbm::addPdfFieldToStorage( structuredBlockforest, "pdf field", latticeModel, Vector3<real_t>(0), 1_r, NUM_GHOSTLAYERS, field::fzyx );
BlockDataID flagFieldId = field::addFlagFieldToStorage< FlagField_T >( structuredBlockforest, "flag field", NUM_GHOSTLAYERS );
const FlagUID fluidFlagUID( "Fluid" );
......
......@@ -307,7 +307,7 @@ static shared_ptr< SetupBlockForest > createSetupBlockForest( const blockforest:
if( blocksPerProcess != 0 )
numberOfProcesses = uint_c( std::ceil( real_c( forest->getNumberOfBlocks() ) / real_c( blocksPerProcess ) ) );
forest->balanceLoad( blockforest::StaticLevelwiseCurveBalance(true), numberOfProcesses, real_t(0), processMemoryLimit, true );
forest->balanceLoad( blockforest::StaticLevelwiseCurveBalance(true), numberOfProcesses, 0_r, processMemoryLimit, true );
if( outputSetupForest )
{
......@@ -408,7 +408,7 @@ MyBoundaryHandling<LatticeModel_T>::operator()( IBlock * const block ) const
return new BoundaryHandling_T( "boundary handling", flagField, fluid,
boost::tuples::make_tuple( NoSlip_T( "no slip", NoSlip_Flag, pdfField ),
UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, real_t(0), real_t(0) ) ) );
UBB_T( "velocity bounce back", UBB_Flag, pdfField, topVelocity_, 0_r, 0_r ) ) );
}
......@@ -570,7 +570,7 @@ real_t exactFlowRate( const real_t flowRate )
Vector3< real_t > exactVelocity( const Vector3< real_t > & p, const math::AABB & domain, const real_t maxLatticeVelocity )
{
return Vector3< real_t >( maxLatticeVelocity * ( p[1] - domain.yMin() ) / domain.ySize(), real_t(0), real_t(0) );
return Vector3< real_t >( maxLatticeVelocity * ( p[1] - domain.yMin() ) / domain.ySize(), 0_r, 0_r );
}
......@@ -667,8 +667,8 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
setup.viscosity_L = latticeModel.collisionModel().viscosity( uint_t(0) );
setup.meanVelocity_L = ( setup.Re * setup.viscosity_L ) / real_c( setup.yBlocks * setup.yCells );
setup.maxVelocity_L = real_t(2) * setup.meanVelocity_L;
setup.flowRate_L = ( setup.maxVelocity_L * real_c( setup.yBlocks * setup.yCells ) * real_c( setup.zBlocks * setup.zCells ) ) / real_t(2);
setup.maxVelocity_L = 2_r * setup.meanVelocity_L;
setup.flowRate_L = ( setup.maxVelocity_L * real_c( setup.yBlocks * setup.yCells ) * real_c( setup.zBlocks * setup.zCells ) ) / 2_r;
// creating the block structure
......@@ -676,13 +676,13 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
// add pdf field to blocks
const real_t initVelocity = ( configBlock.getParameter< bool >( "initWithMeanVelocity", false ) ) ? setup.meanVelocity_L : real_t(0);
const real_t initVelocity = ( configBlock.getParameter< bool >( "initWithMeanVelocity", false ) ) ? setup.meanVelocity_L : 0_r;
BlockDataID pdfFieldId = fzyx ? lbm::addPdfFieldToStorage( blocks, "pdf field (fzyx)", latticeModel,
Vector3< real_t >( initVelocity, real_c(0), real_c(0) ), real_t(1),
Vector3< real_t >( initVelocity, real_c(0), real_c(0) ), 1_r,
FieldGhostLayers, field::fzyx ) :
lbm::addPdfFieldToStorage( blocks, "pdf field (zyxf)", latticeModel,
Vector3< real_t >( initVelocity, real_c(0), real_c(0) ), real_t(1),
Vector3< real_t >( initVelocity, real_c(0), real_c(0) ), 1_r,
FieldGhostLayers, field::zyxf );
using VelocityAdaptor_T = typename lbm::Adaptor< LatticeModel_T >::VelocityVector;
......@@ -743,18 +743,18 @@ void run( const shared_ptr< Config > & config, const LatticeModel_T & latticeMod
flagFieldId, Fluid_Flag,
std::bind( exactFlowRate, setup.flowRate_L ),
exactSolutionFunction );
volumetricFlowRate->setNormalizationFactor( real_t(1) / setup.maxVelocity_L );
volumetricFlowRate->setDomainNormalization( Vector3<real_t>( real_t(1) ) );
volumetricFlowRate->setNormalizationFactor( 1_r / setup.maxVelocity_L );
volumetricFlowRate->setDomainNormalization( Vector3<real_t>( 1_r ) );
timeloop.addFuncBeforeTimeStep( makeSharedFunctor( volumetricFlowRate ), "volumetric flow rate evaluation" );
auto accuracyEvaluation = field::makeAccuracyEvaluation< VelocityAdaptor_T >( configBlock, blocks, velocityAdaptorId, exactSolutionFunction );
accuracyEvaluation->setNormalizationFactor( real_t(1) / setup.maxVelocity_L );
accuracyEvaluation->setNormalizationFactor( 1_r / setup.maxVelocity_L );
timeloop.addFuncBeforeTimeStep( makeSharedFunctor( accuracyEvaluation ), "accuracy evaluation" );
auto linePlot = field::makeAccuracyEvaluationLinePlot< VelocityAdaptor_T >( configBlock, blocks, velocityAdaptorId, exactSolutionFunction );
linePlot->setNormalizationFactor( real_t(1) / setup.maxVelocity_L );
linePlot->setNormalizationFactor( 1_r / setup.maxVelocity_L );
timeloop.addFuncBeforeTimeStep( makeSharedFunctor( field::makeAccuracyEvaluationLinePlotter( configBlock, linePlot ) ), "accuracy evaluation (line plot)" );
......@@ -1024,7 +1024,7 @@ int main( int argc, char **argv )
setup.yCells = configBlock.getParameter< uint_t >( "yCells", uint_t(50) );
setup.zCells = configBlock.getParameter< uint_t >( "zCells", uint_t(10) );
setup.Re = configBlock.getParameter< real_t >( "Re", real_t(10) );
setup.Re = configBlock.getParameter< real_t >( "Re", 10_r );
// ... in bytes
const memory_t memoryPerCell = configBlock.getParameter< memory_t >( "memoryPerCell", memory_t( 19 * 8 + 1 ) );
......@@ -1047,7 +1047,7 @@ int main( int argc, char **argv )
if( !configBlock.isDefined("borderRefinementLevel") )
WALBERLA_ABORT( "You have to specify \'borderRefinementLevel\' in the \"CouetteFlow\" block of the configuration file (" << argv[1] << ")" );
const real_t borderRefinementBuffer = configBlock.getParameter< real_t >( "borderRefinementBuffer", real_t(0) );
const real_t borderRefinementBuffer = configBlock.getParameter< real_t >( "borderRefinementBuffer", 0_r );
BorderRefinementSelection borderRefinementSelection( setup, configBlock.getParameter< uint_t >( "borderRefinementLevel" ),
borderRefinementBuffer );
......@@ -1144,19 +1144,19 @@ int main( int argc, char **argv )
// executing benchmark
const real_t omega = configBlock.getParameter< real_t >( "omega", real_t(1.4) );
const real_t omega = configBlock.getParameter< real_t >( "omega", 1.4_r );
const real_t magicNumber = configBlock.getParameter< real_t >( "magicNumber", real_t(3) / real_t(16) );
const real_t magicNumber = configBlock.getParameter< real_t >( "magicNumber", 3_r / 16_r );
const real_t lambda_e = configBlock.getParameter< real_t >( "lambda_e", real_t(1.4) );
const real_t lambda_d = configBlock.getParameter< real_t >( "lambda_d", real_t(1.4) );
const real_t lambda_e = configBlock.getParameter< real_t >( "lambda_e", 1.4_r );
const real_t lambda_d = configBlock.getParameter< real_t >( "lambda_d", 1.4_r );
const real_t s1 = configBlock.getParameter< real_t >( "s1", real_t(1.4) );
const real_t s2 = configBlock.getParameter< real_t >( "s2", real_t(1.4) );
const real_t s4 = configBlock.getParameter< real_t >( "s4", real_t(1.4) );
const real_t s9 = configBlock.getParameter< real_t >( "s9", real_t(1.4) );
const real_t s10 = configBlock.getParameter< real_t >( "s10", real_t(1.4) );
const real_t s16 = configBlock.getParameter< real_t >( "s16", real_t(1.4) );
const real_t s1 = configBlock.getParameter< real_t >( "s1", 1.4_r );
const real_t s2 = configBlock.getParameter< real_t >( "s2", 1.4_r );
const real_t s4 = configBlock.getParameter< real_t >( "s4", 1.4_r );
const real_t s9 = configBlock.getParameter< real_t >( "s9", 1.4_r );
const real_t s10 = configBlock.getParameter< real_t >( "s10", 1.4_r );
const real_t s16 = configBlock.getParameter< real_t >( "s16", 1.4_r );
const uint_t relaxationParametersLevel = configBlock.getParameter< uint_t >( "relaxationParametersLevel", uint_t(0) );
......
......@@ -113,7 +113,7 @@ const FlagUID MO_CLI_Flag( "moving obstacle CLI" );
static void refinementSelection( SetupBlockForest& forest, uint_t levels, AABB refinementBox )
{
real_t dx = real_t(1); // dx on finest level
real_t dx = 1_r; // dx on finest level
for( auto block = forest.begin(); block != forest.end(); ++block )
{
uint_t blockLevel = block->getLevel();
......@@ -170,11 +170,11 @@ static shared_ptr< StructuredBlockForest > createBlockStructure( const AABB & do
// refinement area is the complete area along the bottom plane, containing the sphere
// this avoids refinement borders in flow direction
refinementBox = AABB(domainAABB.xMin(), domainAABB.yMin(), domainAABB.zMin(),
domainAABB.xMax(), domainAABB.yMax(), initialPosition[2] + real_t(0.5) * diameter);
domainAABB.xMax(), domainAABB.yMax(), initialPosition[2] + 0.5_r * diameter);
} else{
// refinement area is just around the sphere
refinementBox = AABB(initialPosition[0] - real_t(0.5) * diameter, initialPosition[1] - real_t(0.5) * diameter, domainAABB.zMin(),
initialPosition[0] + real_t(0.5) * diameter, initialPosition[1] + real_t(0.5) * diameter, initialPosition[2] + real_t(0.5) * diameter);
refinementBox = AABB(initialPosition[0] - 0.5_r * diameter, initialPosition[1] - 0.5_r * diameter, domainAABB.zMin(),
initialPosition[0] + 0.5_r * diameter, initialPosition[1] + 0.5_r * diameter, initialPosition[2] + 0.5_r * diameter);
}
WALBERLA_LOG_INFO_ON_ROOT(" - refinement box: " << refinementBox);
......@@ -187,7 +187,7 @@ static shared_ptr< StructuredBlockForest > createBlockStructure( const AABB & do
// calculate process distribution
const memory_t memoryLimit = math::Limits< memory_t >::inf();
sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalance(true), uint_c( MPIManager::instance()->numProcesses() ), real_t(0), memoryLimit, true );
sforest.balanceLoad( blockforest::StaticLevelwiseCurveBalance(true), uint_c( MPIManager::instance()->numProcesses() ), 0_r, memoryLimit, true );
WALBERLA_LOG_INFO_ON_ROOT( sforest );
......@@ -276,8 +276,8 @@ public:
void operator()()
{
Vector3<real_t> force(real_t(0));
Vector3<real_t> torque(real_t(0));
Vector3<real_t> force(0_r);
Vector3<real_t> torque(0_r);
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
......@@ -446,15 +446,15 @@ int main( int argc, char **argv )
std::string baseFolderLogging = ".";
// physical setup
real_t diameter = real_t(20); // cells per diameter -> determines overall resolution
real_t normalizedWallDistance = real_t(1); // distance of the sphere center to the bottom wall, normalized by the diameter
real_t ReynoldsNumberShear = real_t(1); // = shearRate * wallDistance * diameter / viscosity
real_t diameter = 20_r; // cells per diameter -> determines overall resolution
real_t normalizedWallDistance = 1_r; // distance of the sphere center to the bottom wall, normalized by the diameter
real_t ReynoldsNumberShear = 1_r; // = shearRate * wallDistance * diameter / viscosity
//numerical parameters
real_t minimumNonDimTimesteps = real_t(100); // minimum number of non-dimensional time steps before simulation can be terminated by convergence
real_t minimumNonDimTimesteps = 100_r; // minimum number of non-dimensional time steps before simulation can be terminated by convergence
uint_t numberOfLevels = uint_t(4); // number of grid levels for static refinement ( 1 = no refinement)
real_t xOffsetOfSpherePosition = real_t(0); // offset in x-direction of sphere position
real_t yOffsetOfSpherePosition = real_t(0); // offset in y-direction of sphere position
real_t xOffsetOfSpherePosition = 0_r; // offset in x-direction of sphere position
real_t yOffsetOfSpherePosition = 0_r; // offset in y-direction of sphere position
bool useSBB = false; // use simple bounce-back boundary condition for sphere
bool useLargeRefinementRegion = false; // uses the whole area near the bottom plane as the finest grid, else refinement is only applied around the sphere
......@@ -479,41 +479,41 @@ int main( int argc, char **argv )
WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
}
WALBERLA_CHECK_GREATER_EQUAL(normalizedWallDistance, real_t(0.5));
WALBERLA_CHECK_GREATER_EQUAL(ReynoldsNumberShear, real_t(0));
WALBERLA_CHECK_GREATER_EQUAL(diameter, real_t(0));
WALBERLA_CHECK_GREATER_EQUAL(normalizedWallDistance, 0.5_r);
WALBERLA_CHECK_GREATER_EQUAL(ReynoldsNumberShear, 0_r);
WALBERLA_CHECK_GREATER_EQUAL(diameter, 0_r);
//////////////////////////
// NUMERICAL PARAMETERS //
//////////////////////////
const real_t domainLength = real_t(48) * diameter; //x
const real_t domainWidth = real_t(16) * diameter; //y
const real_t domainLength = 48_r * diameter; //x
const real_t domainWidth = 16_r * diameter; //y
const real_t domainHeight = real_t( 8) * diameter; //z
Vector3<uint_t> domainSize( uint_c( std::ceil(domainLength)), uint_c( std::ceil(domainWidth)), uint_c( std::ceil(domainHeight)) );
real_t wallVelocity = real_t(0.01);
if( zeroShearTest ) wallVelocity = real_t(0);
real_t wallVelocity = 0.01_r;
if( zeroShearTest ) wallVelocity = 0_r;
const real_t wallDistance = diameter * normalizedWallDistance;
const real_t shearRate = wallVelocity / domainHeight;
const real_t velAtSpherePosition = shearRate * wallDistance;
const real_t viscosity = ( zeroShearTest ) ? real_t(0.015) : ( velAtSpherePosition * diameter / ReynoldsNumberShear );
const real_t viscosity = ( zeroShearTest ) ? 0.015_r : ( velAtSpherePosition * diameter / ReynoldsNumberShear );
const real_t relaxationTime = real_t(1) / lbm::collision_model::omegaFromViscosity(viscosity);
const real_t relaxationTime = 1_r / lbm::collision_model::omegaFromViscosity(viscosity);
const real_t densityFluid = real_t(1);
const real_t densityFluid = 1_r;
const real_t dx = real_t(1);
const real_t dx = 1_r;
const real_t physicalTimeScale = ( zeroShearTest ) ? real_t(10) : (diameter / velAtSpherePosition);
const real_t physicalTimeScale = ( zeroShearTest ) ? 10_r : (diameter / velAtSpherePosition);
const uint_t minimumLBMtimesteps = uint_c(minimumNonDimTimesteps * physicalTimeScale);
const real_t omega = real_t(1) / relaxationTime;
const real_t omega = 1_r / relaxationTime;
const uint_t finestLevel = numberOfLevels - uint_t(1);
Vector3<real_t> initialPosition( domainLength * real_t(0.5) + xOffsetOfSpherePosition,
domainWidth * real_t(0.5) + yOffsetOfSpherePosition,
Vector3<real_t> initialPosition( domainLength * 0.5_r + xOffsetOfSpherePosition,
domainWidth * 0.5_r + yOffsetOfSpherePosition,
wallDistance );
WALBERLA_LOG_INFO_ON_ROOT("Setup:");
......@@ -556,7 +556,7 @@ int main( int argc, char **argv )
domainSize[1] / ( coarseBlocksPerDirection[1] * levelScalingFactor ),
domainSize[2] / ( coarseBlocksPerDirection[2] * levelScalingFactor ) );
AABB simulationDomain( real_t(0), real_t(0), real_t(0), real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
AABB simulationDomain( 0_r, 0_r, 0_r, real_c(domainSize[0]), real_c(domainSize[1]), real_c(domainSize[2]) );
auto blocks = createBlockStructure( simulationDomain, blockSizeInCells, numberOfLevels, diameter, initialPosition, useLargeRefinementRegion );
//write domain decomposition to file
......@@ -583,13 +583,13 @@ int main( int argc, char **argv )
// create pe bodies
// bounding planes (global)
const auto planeMaterial = pe::createMaterial( "myPlaneMat", real_t(8920), real_t(0), real_t(1), real_t(1), real_t(0), real_t(1), real_t(1), real_t(0), real_t(0) );
const auto planeMaterial = pe::createMaterial( "myPlaneMat", 8920_r, 0_r, 1_r, 1_r, 0_r, 1_r, 1_r, 0_r, 0_r );
pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,1), Vector3<real_t>(0,0,0), planeMaterial );
auto topPlane = pe::createPlane( *globalBodyStorage, 0, Vector3<real_t>(0,0,-1), Vector3<real_t>(0,0,domainHeight), planeMaterial );
topPlane->setLinearVel(wallVelocity, real_t(0), real_t(0));
topPlane->setLinearVel(wallVelocity, 0_r, 0_r);
// add the sphere
pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, initialPosition, real_t(0.5) * diameter, planeMaterial );
pe::createSphere( *globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0, initialPosition, 0.5_r * diameter, planeMaterial );
uint_t minBlockSizeInCells = blockSizeInCells.min();
for( uint_t i = 0; i < uint_c(diameter / real_c(minBlockSizeInCells)) + 1; ++i)
......@@ -604,7 +604,7 @@ int main( int argc, char **argv )
// add PDF field
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
Vector3< real_t >( real_t(0) ), real_t(1),
Vector3< real_t >( 0_r ), 1_r,
FieldGhostLayers, field::zyxf );
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field", FieldGhostLayers );
......@@ -669,12 +669,12 @@ int main( int argc, char **argv )
auto refinementTimestep = lbm::refinement::makeTimeStep< LatticeModel_T, BoundaryHandling_T >( blocks, sweep, pdfFieldID, boundaryHandlingID );
// add force evaluation and logging
real_t normalizationFactor = ( zeroShearTest ) ? real_t(1) : ( math::M_PI / real_t(8) * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter );
real_t normalizationFactor = ( zeroShearTest ) ? 1_r : ( math::M_PI / 8_r * densityFluid * shearRate * shearRate * wallDistance * wallDistance * diameter * diameter );
std::string loggingFileName( baseFolderLogging + "/LoggingForcesNearPlane");
loggingFileName += "_lvl" + std::to_string(numberOfLevels);
loggingFileName += "_D" + std::to_string(uint_c(diameter));
loggingFileName += "_Re" + std::to_string(uint_c(ReynoldsNumberShear));
loggingFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * real_t(1000)));
loggingFileName += "_WD" + std::to_string(uint_c(normalizedWallDistance * 1000_r));
loggingFileName += ".txt";
shared_ptr< SpherePropertyLogger > logger = walberla::make_shared< SpherePropertyLogger >( blocks, bodyStorageID,
loggingFileName, fileIO,
......@@ -693,22 +693,22 @@ int main( int argc, char **argv )
// compute reference values from literature
const real_t normalizedGapSize = normalizedWallDistance - real_t(0.5);
const real_t normalizedGapSize = normalizedWallDistance - 0.5_r;
// drag correlation for the drag coefficient
const real_t standardDragCorrelation = real_t(24) / ReynoldsNumberShear * (real_t(1) + real_t(0.15) * std::pow(ReynoldsNumberShear, real_t(0.687))); // Schiller-Naumann correlation
const real_t dragCorrelationWithGapSizeStokes = real_t(24) / ReynoldsNumberShear * (real_t(1) + real_t(0.138) * std::exp(real_t(-2) * normalizedGapSize) + real_t(9)/( real_t(16) * (real_t(1) + real_t(2) * normalizedGapSize) ) ); // Goldman et al. (1967)
const real_t alphaDragS = real_t(0.15) - real_t(0.046) * ( real_t(1) - real_t(0.16) * normalizedGapSize * normalizedGapSize ) * std::exp( -real_t(0.7) * normalizedGapSize);
const real_t betaDragS = real_t(0.687) + real_t(0.066)*(real_t(1) - real_t(0.76) * normalizedGapSize * normalizedGapSize) * std::exp( - std::pow( normalizedGapSize, real_t(0.9) ) );
const real_t dragCorrelationZeng = dragCorrelationWithGapSizeStokes * ( real_t(1) + alphaDragS * std::pow( ReynoldsNumberShear, betaDragS ) ); // Zeng et al. (2009) - Eqs. (13) and (14)
const real_t standardDragCorrelation = 24_r / ReynoldsNumberShear * (1_r + 0.15_r * std::pow(ReynoldsNumberShear, 0.687_r)); // Schiller-Naumann correlation
const real_t dragCorrelationWithGapSizeStokes = 24_r / ReynoldsNumberShear * (1_r + 0.138_r * std::exp(-2_r * normalizedGapSize) + 9_r/( 16_r * (1_r + 2_r * normalizedGapSize) ) ); // Goldman et al. (1967)
const real_t alphaDragS = 0.15_r - 0.046_r * ( 1_r - 0.16_r * normalizedGapSize * normalizedGapSize ) * std::exp( -0.7_r * normalizedGapSize);
const real_t betaDragS = 0.687_r + 0.066_r*(1_r - 0.76_r * normalizedGapSize * normalizedGapSize) * std::exp( - std::pow( normalizedGapSize, 0.9_r ) );
const real_t dragCorrelationZeng = dragCorrelationWithGapSizeStokes * ( 1_r + alphaDragS * std::pow( ReynoldsNumberShear, betaDragS ) ); // Zeng et al. (2009) - Eqs. (13) and (14)
// lift correlations for the lift coefficient
const real_t liftCorrelationZeroGapStokes = real_t(5.87); // Leighton, Acrivos (1985)
const real_t liftCorrelationZeroGap = real_t(3.663) / std::pow( ReynoldsNumberShear * ReynoldsNumberShear + real_t(0.1173), real_t(0.22) ); // Zeng et al. (2009) - Eq. (19)
const real_t alphaLiftS = - std::exp( -real_t(0.3) + real_t(0.025) * ReynoldsNumberShear);
const real_t betaLiftS = real_t(0.8) + real_t(0.01) * ReynoldsNumberShear;
const real_t lambdaLiftS = ( real_t(1) - std::exp(-normalizedGapSize)) * std::pow( ReynoldsNumberShear / real_t(250), real_t(5) / real_t(2) );
const real_t liftCorrelationZeng = liftCorrelationZeroGap * std::exp( - real_t(0.5) * normalizedGapSize * std::pow( ReynoldsNumberShear / real_t(250), real_t(4)/real_t(3))) *
const real_t liftCorrelationZeroGapStokes = 5.87_r; // Leighton, Acrivos (1985)
const real_t liftCorrelationZeroGap = 3.663_r / std::pow( ReynoldsNumberShear * ReynoldsNumberShear + 0.1173_r, 0.22_r ); // Zeng et al. (2009) - Eq. (19)
const real_t alphaLiftS = - std::exp( -0.3_r + 0.025_r * ReynoldsNumberShear);
const real_t betaLiftS = 0.8_r + 0.01_r * ReynoldsNumberShear;
const real_t lambdaLiftS = ( 1_r - std::exp(-normalizedGapSize)) * std::pow( ReynoldsNumberShear / 250_r, 5_r / 2_r );
const real_t liftCorrelationZeng = liftCorrelationZeroGap * std::exp( - 0.5_r * normalizedGapSize * std::pow( ReynoldsNumberShear / 250_r, 4_r/3_r)) *
( std::exp( alphaLiftS * std::pow( normalizedGapSize, betaLiftS ) ) - lambdaLiftS ); // Zeng et al. (2009) - Eqs. (28) and (29)
////////////////////////
......@@ -717,8 +717,8 @@ int main( int argc, char **argv )
WcTimingPool timeloopTiming;
const real_t relativeChangeConvergenceEps = real_t(1e-3);
const real_t physicalCheckingFrequency = real_t(0.00625);
const real_t relativeChangeConvergenceEps = 1e-3_r;
const real_t physicalCheckingFrequency = 0.00625_r;
const uint_t checkingFrequency = (zeroShearTest) ? uint_t(1) : uint_c( physicalCheckingFrequency * physicalTimeScale );
WALBERLA_LOG_INFO_ON_ROOT("Starting simulation with at least " << timesteps << " (coarse) time steps");
......
......@@ -200,7 +200,7 @@ BoundaryHandling_T * MyBoundaryHandling::operator()( IBlock * const block, const
BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
boost::tuples::make_tuple( UBB_T( "UBB", UBB_Flag, pdfField, velocity_),
Outlet_T( "Outlet", Outlet_Flag, pdfField, real_t(1) ),
Outlet_T( "Outlet", Outlet_Flag, pdfField, 1_r ),
MEM_BB_T ( "MEM_BB", MEM_BB_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
MEM_CLI_T( "MEM_CLI", MEM_CLI_Flag, pdfField, flagField, bodyField, fluid, *storage, *block ),
MEM_MR_T ( "MEM_MR", MEM_MR_Flag, pdfField, flagField, bodyField, fluid, *storage, *block, pdfFieldPreCol ) ) );
......@@ -279,7 +279,7 @@ public:
dragForceOld_ = dragForceNew_;
dragForceNew_ = currentAverage_ / real_c( averageFrequency_ );
currentAverage_ = real_t(0);
currentAverage_ = 0_r;
}
......@@ -319,11 +319,11 @@ private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID bodyStorageID_;
real_t currentAverage_ = real_t(0);
real_t currentAverage_ = 0_r;
uint_t averageFrequency_;
std::string filename_;
real_t dragForceOld_ = real_t(0);
real_t dragForceNew_ = real_t(0);
real_t dragForceOld_ = 0_r;
real_t dragForceNew_ = 0_r;
};
......@@ -358,7 +358,7 @@ public:
fileSetup << "viscosity = " << viscosity << "\n";
fileSetup << "diameter = " << diameter << "\n";
fileSetup << "uInfty = " << u_infty_[2] << "\n";
real_t u_ref = std::sqrt( std::fabs(densityRatio - real_t(1)) * gravity * diameter );
real_t u_ref = std::sqrt( std::fabs(densityRatio - 1_r) * gravity * diameter );
fileSetup << "u_{ref} = " << u_ref << "\n";
fileSetup << "numLBMSubCycles = " << numLBMSubCycles << "\n";
fileSetup.close();
......@@ -392,12 +392,12 @@ public:
{
const uint_t timestep ( timeloop_->getCurrentTimeStep() * numLBMSubCycles_ + 1 );
Vector3<real_t> transVel( real_t(0) );
Vector3<real_t> angularVel( real_t(0) );
Vector3<real_t> pos( real_t(0) );
Vector3<real_t> transVel( 0_r );
Vector3<real_t> angularVel( 0_r );
Vector3<real_t> pos( 0_r );
Vector3<real_t> force( real_t(0) );
Vector3<real_t> torque( real_t(0) );
Vector3<real_t> force( 0_r );
Vector3<real_t> torque( 0_r );
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
......@@ -450,7 +450,7 @@ public:
real_t omega_p_H = std::sqrt( particleAngularVel[0] * particleAngularVel[0] + particleAngularVel[1] * particleAngularVel[1] );
real_t omega_p_V = particleAngularVel[2];
real_t u_ref = std::sqrt( std::fabs(densityRatio_ - real_t(1)) * gravity_ * diameter_ );
real_t u_ref = std::sqrt( std::fabs(densityRatio_ - 1_r) * gravity_ * diameter_ );
real_t Reynolds = u_p_r.length() * diameter_ / viscosity_;
// results
......@@ -528,7 +528,7 @@ public:
void operator()()
{
Vector3<real_t> u_p( real_t(0) );
Vector3<real_t> u_p( 0_r );
for( auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt )
{
for( auto bodyIt = pe::LocalBodyIterator::begin<pe::Sphere>(*blockIt, bodyStorageID_ ); bodyIt != pe::LocalBodyIterator::end<pe::Sphere>(); ++bodyIt )
......@@ -550,9 +550,9 @@ public:
real_t u_p_H = std::sqrt( u_p_r[0] * u_p_r[0] + u_p_r[1] * u_p_r[1]);
Vector3<real_t> e_p_H (u_p_r[0], u_p_r[1], real_t(0));
Vector3<real_t> e_p_H (u_p_r[0], u_p_r[1], 0_r);
e_p_H /= u_p_H;
Vector3<real_t> e_p_Hz_perp (-u_p_r[1], u_p_r[0], real_t(0));
Vector3<real_t> e_p_Hz_perp (-u_p_r[1], u_p_r[0], 0_r);
e_p_Hz_perp /= u_p_H;
Vector3<real_t> e_p_parallel = u_p_r / u_p_r.length();
......@@ -669,8 +669,8 @@ int main( int argc, char **argv )
bool usePSM = false;
PSMVariant psmVariant = PSMVariant::SC3W2;
real_t Galileo = real_t(144);
real_t diameter = real_t(18);
real_t Galileo = 144_r;
real_t diameter = 18_r;
////////////////////////////
// COMMAND LINE ARGUMENTS //
......@@ -695,39 +695,39 @@ int main( int argc, char **argv )
// SIMULATION PROPERTIES //
///////////////////////////
const real_t radius = real_t(0.5) * diameter;
const uint_t xlength = uint_c( diameter * real_t(5.34) );
const real_t radius = 0.5_r * diameter;
const uint_t xlength = uint_c( diameter * 5.34_r );
const uint_t ylength = xlength;
uint_t zlength = uint_c( diameter * real_t(16) );
uint_t zlength = uint_c( diameter * 16_r );
if (longDom)
{
zlength *= uint_t(3);
}
real_t viscosity = real_t(0.01);
real_t Re_target = real_t(1);
real_t timestepsNonDim = real_t(1);
real_t viscosity = 0.01_r;
real_t Re_target = 1_r;
real_t timestepsNonDim = 1_r;
// estimate Reynolds number (i.e. inflow velocity) based on Galileo number
// values are taken from the original simulation of Uhlmann, Dusek
switch( int(Galileo) )
{
case 144:
Re_target = real_t(185.08);
timestepsNonDim = real_t(100);
Re_target = 185.08_r;
timestepsNonDim = 100_r;
break;
case 178:
Re_target = real_t(243.01);
timestepsNonDim = real_t(250);
Re_target = 243.01_r;
timestepsNonDim = 250_r;
break;
case 190:
Re_target = real_t(262.71);
timestepsNonDim = real_t(250);
Re_target = 262.71_r;
timestepsNonDim = 250_r;
break;
case 250:
Re_target = real_t(365.10);
timestepsNonDim = real_t(510);
Re_target = 365.10_r;
timestepsNonDim = 510_r;
break;
default:
WALBERLA_ABORT("Galileo number is different from the usual ones (144, 178, 190, or 250). No estimate of the Reynolds number available. Add this case manually!");
......@@ -738,20 +738,20 @@ int main( int argc, char **argv )
real_t omega = lbm::collision_model::omegaFromViscosity(viscosity);
Vector3<real_t> uInfty = Vector3<real_t>( real_t(0), real_t(0), uIn );
Vector3<real_t> uInfty = Vector3<real_t>( 0_r, 0_r, uIn );
const real_t densityRatio = real_t(1.5);
const real_t densityRatio = 1.5_r;
const uint_t averageFrequency = uint_c( ( ( uint_c(Galileo) >= 200) ? real_t(500) : real_t(2) ) * diameter / uIn ); // for initial simulation
const real_t convergenceLimit = real_t(1e-4);
const real_t convergenceLimitGalileo = real_t(1e-4);
const real_t dx = real_t(1);
const uint_t averageFrequency = uint_c( ( ( uint_c(Galileo) >= 200) ? 500_r : 2_r ) * diameter / uIn ); // for initial simulation
const real_t convergenceLimit = 1e-4_r;
const real_t convergenceLimitGalileo = 1e-4_r;
const real_t dx = 1_r;