Commit effb6b03 authored by Dominik Schuster's avatar Dominik Schuster
Browse files

Added HCSITS and preliminary timeloop functionality

parent a64c4053
......@@ -158,10 +158,39 @@ namespace fluidizedBed {
CellInterval domainBB = storage->getDomainCellBB();
//domainBB.xMin() -= cell_idx_c(FieldGhostLayers);
//domainBB.xMax() += cell_idx_c(FieldGhostLayers);
//domainBB.zMin() -= cell_idx_c(FieldGhostLayers);
//domainBB.zMax() += cell_idx_c(FieldGhostLayers);
domainBB.xMin() -= cell_idx_c(FieldGhostLayers);
domainBB.xMax() += cell_idx_c(FieldGhostLayers);
if (!xPeriodic_) {
// LEFT
CellInterval left(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMin(), domainBB.yMax(), domainBB.zMax());
storage->transformGlobalToBlockLocalCellInterval(left, *block);
handling->forceBoundary(noSlip, left);
// RIGHT
CellInterval right(domainBB.xMax(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax());
storage->transformGlobalToBlockLocalCellInterval(right, *block);
handling->forceBoundary(noSlip, right);
}
domainBB.zMin() -= cell_idx_c(FieldGhostLayers);
domainBB.zMax() += cell_idx_c(FieldGhostLayers);
if (!zPeriodic_) {
// FRONT
CellInterval front(domainBB.xMin(), domainBB.yMin(), domainBB.zMin(), domainBB.xMax(), domainBB.yMax(), domainBB.zMin());
storage->transformGlobalToBlockLocalCellInterval(front, *block);
handling->forceBoundary(noSlip, front);
// BACK
CellInterval back(domainBB.xMin(), domainBB.yMin(), domainBB.zMax(), domainBB.xMax(), domainBB.yMax(), domainBB.zMax());
storage->transformGlobalToBlockLocalCellInterval(back, *block);
handling->forceBoundary(noSlip, back);
}
domainBB.yMin() -= cell_idx_c(FieldGhostLayers);
domainBB.yMax() += cell_idx_c(FieldGhostLayers);
// BOTTOM
......@@ -376,7 +405,7 @@ namespace fluidizedBed {
const real_t peAccelerator = sim_parameters.getParameter<real_t>("pe_accelerator", real_c(50));
const bool fixedBed = sim_parameters.getParameter<bool>("fixed_bed", false);
const bool runPreliminaryTimeloop = sim_parameters.getParameter<bool>("run_preliminary_timeloop", false);
const bool runPreliminaryTimeloop = sim_parameters.getParameter<bool>("run_preliminary_timeloop", false);
const real_t eps = sim_parameters.getParameter<real_t>("eps", real_c(0.05));
const bool initializeFromCheckPointFile = sim_parameters.getParameter<bool>("initialize_from_checkpoint", false);
......@@ -467,6 +496,11 @@ namespace fluidizedBed {
// MATERIAL PARAMETERS
auto mat_parameters = env.config()->getBlock("Material");
const bool useDEM = mat_parameters.getParameter<bool>("use_DEM", true);
const uint_t HCSITSMaxIterations = mat_parameters.getParameter<uint_t>("HCSITS_max_iterations", uint_c(10));
const real_t HCSITSRelaxationParameter = mat_parameters.getParameter<real_t>("HCSITS_relaxation_parameter", real_c(0.7));
const real_t volumeA = real_c(4.0 / 3.0 * math::PI * radiusA * radiusA * radiusA);
const real_t volumeB = real_c(4.0 / 3.0 * math::PI * radiusB * radiusB * radiusB);
const real_t MijA = real_c(0.5 * config.densityRatioA * volumeA);
......@@ -873,9 +907,20 @@ namespace fluidizedBed {
// set up collision response, here DEM solver
auto fcdID_PE = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
pe::cr::DEM cr_PE(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID_PE);
std::unique_ptr<pe::cr::ICR> cr_PE;
if (useDEM)
{
cr_PE = std::make_unique<pe::cr::DEM>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID_PE);
WALBERLA_LOG_INFO_ON_ROOT("Using DEM!");
} else {
cr_PE = std::make_unique<pe::cr::HCSITS>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID_PE);
pe::cr::HCSITS* hcsits_PE = static_cast<pe::cr::HCSITS*>(cr_PE.get());
hcsits_PE->setMaxIterations(HCSITSMaxIterations);
hcsits_PE->setRelaxationParameter(HCSITSRelaxationParameter);
WALBERLA_LOG_INFO_ON_ROOT("Using HCSITS!");
}
cr_PE.setGlobalLinearAcceleration(pe::Vec3(0.0, - (config.gravity), 0.0));
cr_PE->setGlobalLinearAcceleration(pe::Vec3(0.0, - (config.gravity), 0.0));
real_t avgVel = real_c(100.0);
real_t maxVel = real_c(100.0);
......@@ -885,7 +930,7 @@ namespace fluidizedBed {
while ((currentSettleVelocity > settleVelocity) || (itPE < 4999)) {
cr_PE.timestep( peAccelerator / real_c(config.peSubCycles));
cr_PE->timestep( peAccelerator / real_c(config.peSubCycles));
syncCall_PE();
++itPE;
......@@ -913,7 +958,7 @@ namespace fluidizedBed {
}
cr_PE.setGlobalLinearAcceleration(pe::Vec3(0.0, 0.0, 0.0));
cr_PE->setGlobalLinearAcceleration(pe::Vec3(0.0, 0.0, 0.0));
FBfunc::restSphere(blocks, bodyStorageID);
avgVel = FBfunc::getAvgVel(blocks, bodyStorageID, (numberCreatedParticlesA + numberCreatedParticlesB)) * velocityConversion;
......@@ -937,7 +982,6 @@ namespace fluidizedBed {
std::cout << "\nSetup Information\n" << std::endl;
std::cout << "Physical time step is: " << config.dt_SI << " s" << std::endl;
std::cout << "Physical cell size is: " << config.dx_SI << " m" << std::endl;
std::cout << "Number of particles A is: " << config.nrParticlesA << std::endl;
std::cout << "Particle Reynoldsnumber(LBM) is: " << particleReynoldsNumber_LBM << std::endl;
std::cout << "Froudenumber(LBM) is: " << FroudeNumber_LBM << std::endl;
std::cout << "Minimum Fluidization Re number is: " << REpmf << std::endl;
......@@ -1031,7 +1075,17 @@ namespace fluidizedBed {
// set up collision response, here DEM solver
auto fcdID = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
std::unique_ptr<pe::cr::ICR> cr;
if (useDEM)
{
cr = std::make_unique<pe::cr::DEM>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
} else {
cr = std::make_unique<pe::cr::HCSITS>(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID);
pe::cr::HCSITS* hcsits = static_cast<pe::cr::HCSITS*>(cr.get());
hcsits->setMaxIterations(HCSITSMaxIterations);
hcsits->setRelaxationParameter(HCSITSRelaxationParameter);
}
// add flag field
BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>(blocks, "flag field");
......@@ -1040,12 +1094,7 @@ namespace fluidizedBed {
BlockDataID bodyFieldID = field::addToStorage<BodyField_T>(blocks, "body field", NULL, field::zyxf);
// Object for keeping track of time
shared_ptr<lbm::TimeTracker> preTimeTrack = make_shared<lbm::TimeTracker>();
shared_ptr<lbm::TimeTracker> timeTrack = make_shared<lbm::TimeTracker>();
BlockDataID preBoundaryHandlingID = blocks->addStructuredBlockData<BoundaryHandling_T>(
MyBoundaryHandling_T(flagFieldID, pdfFieldID, bodyFieldID, xPeriodic, zPeriodic, spotInflow, spotDiameter, real_c(config.xlength)*0.5,
real_c(config.zlength)*0.5, config.velocity * initialRampingVelocityFactor, real_t(1), real_t(0), preTimeTrack), "pre boundary handling");
shared_ptr<lbm::TimeTracker> timeTrack = make_shared<lbm::TimeTracker>();
BlockDataID boundaryHandlingID = blocks->addStructuredBlockData<BoundaryHandling_T>(
MyBoundaryHandling_T(flagFieldID, pdfFieldID, bodyFieldID, xPeriodic, zPeriodic, spotInflow, spotDiameter, real_c(config.xlength)*0.5,
......@@ -1053,10 +1102,10 @@ namespace fluidizedBed {
// map pe bodies into the LBM simulation
// moving bodies are handled by the momentum exchange method
pe_coupling::mapMovingBodies<BoundaryHandling_T>(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, FBfunc::MO_Flag);
pe_coupling::mapMovingBodies<BoundaryHandling_T>(*blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, bodyFieldID, FBfunc::MO_Flag, pe_coupling::selectRegularBodies);
// map planes into the LBM simulation -> act as no-slip boundaries
pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, FBfunc::NoSlip_Flag, pe_coupling::selectGlobalBodies );
//pe_coupling::mapBodies< BoundaryHandling_T >( *blocks, boundaryHandlingID, bodyStorageID, *globalBodyStorage, FBfunc::NoSlip_Flag, pe_coupling::selectGlobalBodies );
/////////////////////
// PRE TIME LOOP //
......@@ -1065,22 +1114,21 @@ namespace fluidizedBed {
SweepTimeloop timeloopPRE(blocks->getBlockStorage(), 1000000);
// add LBM communication function and boundary handling sweep
timeloopPRE.add() << BeforeFunction(commFunction, "LBM Communication") << Sweep(BoundaryHandling_T::getBlockSweep(preBoundaryHandlingID), "Boundary Handling");
timeloopPRE.add() << BeforeFunction(commFunction, "LBM Communication") << Sweep(BoundaryHandling_T::getBlockSweep(boundaryHandlingID), "Boundary Handling");
// LBM stream collide sweep
timeloopPRE.add() << Sweep(makeSharedSweep(lbm::makeCellwiseSweep<LatticeModel_T, FlagField_T>(pdfFieldID, flagFieldID, FBfunc::Fluid_Flag)), "LBM DEFAULT");
// Reset Forces
timeloopPRE.addFuncAfterTimeStep(pe_coupling::ForceTorqueOnBodiesResetter(blocks, bodyStorageID), "Reset Body Forces");
// Velocity Check
timeloopPRE.add() << Sweep(FBfunc::VelocityCheckMEM<LatticeModel_T, FlagField_T, BodyField_T>(blocks, bodyStorageID, bodyFieldID, pdfFieldID, flagFieldID, FBfunc::Fluid_Flag, config),
"LBM Velocity Check");
shared_ptr<FBfunc::PressureDropper> PressureDropper = make_shared<FBfunc::PressureDropper>(blocks, bodyStorageID, config);
shared_ptr<FBfunc::PressureDropper> PressureDropper = make_shared<FBfunc::PressureDropper>(blocks, bodyStorageID, config, true);
timeloopPRE.addFuncAfterTimeStep( SharedFunctor<FBfunc::PressureDropper>(PressureDropper), "Calculate Pressure Drop");
WcTimingPool timeloopTimingPRE;
// Reset Forces
timeloopPRE.addFuncAfterTimeStep(pe_coupling::ForceTorqueOnBodiesResetter(blocks, bodyStorageID), "Reset Body Forces");
WcTimingPool timeloopTimingPRE;
////////////////////
// MAIN TIME LOOP //
......@@ -1156,7 +1204,7 @@ namespace fluidizedBed {
}
if (calculatePressureDrop) {
timeloop.addFuncAfterTimeStep(FBfunc::PressureDropper(blocks, bodyStorageID, config), "Calculate Pressure Drop");
timeloop.addFuncAfterTimeStep(FBfunc::PressureDropper(blocks, bodyStorageID, config, false), "Calculate Pressure Drop");
}
if (calculateParticleFlux) {
......@@ -1188,7 +1236,7 @@ namespace fluidizedBed {
// advance pe rigid body simulation
if (!fixedBed) {
timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, cr, syncCall, config.lbmSubCycles, config.peSubCycles ), "pe Time Step" );
timeloop.addFuncAfterTimeStep( pe_coupling::TimeStep( blocks, bodyStorageID, *cr, syncCall, config.lbmSubCycles, config.peSubCycles ), "pe Time Step" );
}
////////////////
......@@ -1219,11 +1267,10 @@ namespace fluidizedBed {
timeloop.addFuncAfterTimeStep(vtk::writeFiles(flagFieldVTK), "VTK (flag field data)");
}
if (writeFluidField) {
// pdf field (ghost layers cannot be written because re-sampling/coarsening is applied)
auto pdfFieldVTK = vtk::createVTKOutput_BlockData(blocks, "fluid_field", frequencyFluidField, 0, false); // 0 => no fieldghostlayer
pdfFieldVTK->setSamplingResolution(samplingResolutionFluidField);
pdfFieldVTK->setSamplingResolution(samplingResolutionFluidField);
blockforest::communication::UniformBufferedScheme<stencil::D3Q27> pdfGhostLayerSync(blocks);
pdfGhostLayerSync.addPackInfo(make_shared<field::communication::PackInfo<PdfField_T> >(pdfFieldID));
......@@ -1281,20 +1328,35 @@ namespace fluidizedBed {
std::cout << std::endl << "||||||||||||| PRELIMINARY SIM ||||||||||||||" << std::endl;
}
real_t pressureDifference = real_t(1000);
real_t relativePressureDifference = real_t(1000);
real_t previousPressure = real_t(0);
// perform a single simulation step
timeloopPRE.singleStep(timeloopTimingPRE);
previousPressure = PressureDropper->getPressureDrop();
while (relativePressureDifference > eps) {
while (pressureDifference > eps) {
// perform a single simulation step
timeloopPRE.singleStep(timeloopTimingPRE);
pressureDifference = fabs(PressureDropper->getPressureDrop() - pressureDifference) / pressureDifference;
if (timeloopPRE.getCurrentTimeStep() % uint_t(10) == uint_t(0)) {
relativePressureDifference = fabs(PressureDropper->getPressureDrop() - previousPressure) / previousPressure;
if (timeloopPRE.getCurrentTimeStep() % 20) {
WALBERLA_LOG_INFO_ON_ROOT(
"Current pressure drop over whole bed: " << PressureDropper->getPressureDrop()
<< " [Pa] after "
<< timeloopPRE.getCurrentTimeStep()
<< " timesteps!");
WALBERLA_LOG_INFO_ON_ROOT(
"Relativ change: " << relativePressureDifference );
previousPressure = PressureDropper->getPressureDrop();
WALBERLA_MPI_SECTION() {
mpi::broadcastObject(relativePressureDifference);
}
}
}
......
......@@ -168,10 +168,11 @@ namespace FBfunc {
PressureDropper(const shared_ptr<StructuredBlockStorage> &blocks,
const BlockDataID &bodyStorageID,
SetupFB &config) :
SetupFB &config, const bool preliminaryTimeloop ) :
blocks_(blocks), bodyStorageID_(bodyStorageID),
checkFrequency_((config.evaluationFreq > 0) ? config.evaluationFreq : uint_c(1)), executionCounter_(uint_c(0)),
preFac_((config.densityFluid_SI * config.dx_SI * config.dx_SI / (config.dt_SI * config.dt_SI)) / (real_c(config.xlength * config.zlength))) {
preFac_((config.densityFluid_SI * config.dx_SI * config.dx_SI / (config.dt_SI * config.dt_SI)) / (real_c(config.xlength * config.zlength))),
preliminaryTimeloop_(preliminaryTimeloop), currentPressureDrop_(real_t(0)){
std::ofstream output;
output.open("PressureDrop.txt", std::ofstream::out | std::ofstream::trunc);
......@@ -180,30 +181,34 @@ namespace FBfunc {
void operator()() {
++executionCounter_;
if ((executionCounter_ - uint_c(1)) % checkFrequency_ != 0)
if (((executionCounter_ - uint_c(1)) % checkFrequency_ != 0) && !preliminaryTimeloop_)
return;
pe::Vec3 forceSum = pe::Vec3(0.0, 0.0, 0.0);
pe::Vec3 force = pe::Vec3(0.0, 0.0, 0.0);
for (auto blockIt = blocks_->begin(); blockIt != blocks_->end(); ++blockIt) {
for (auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt) {
forceSum += bodyIt->getForce();
force += bodyIt->getForce();
}
}
WALBERLA_MPI_SECTION() {
mpi::reduceInplace(forceSum[0], mpi::SUM);
mpi::reduceInplace(forceSum[1], mpi::SUM);
mpi::reduceInplace(forceSum[2], mpi::SUM);
mpi::reduceInplace(force[0], mpi::SUM);
mpi::reduceInplace(force[1], mpi::SUM);
mpi::reduceInplace(force[2], mpi::SUM);
}
WALBERLA_ROOT_SECTION() {
forceSum = forceSum * preFac_;
currentPressureDrop_ = forceSum.length();
std::ofstream output;
output.open("PressureDrop.txt", std::ofstream::out | std::ofstream::app);
//output << executionCounter_ << "\tx: " << forceSum[0]<<"\ty: " << forceSum[1]<<"\tz: "<< forceSum[2]<<"\tAbs: "<< forceSum.length() << "\t[Pa]" << std::endl;
output << executionCounter_ << "\t" << currentPressureDrop_ << "\t[Pa]" << std::endl;
output.close();
currentPressureDrop_ = force.length() * preFac_;
if (!preliminaryTimeloop_) {
std::ofstream output;
output.open("PressureDrop.txt", std::ofstream::out | std::ofstream::app);
//output << executionCounter_ << "\tx: " << forceSum[0]<<"\ty: " << forceSum[1]<<"\tz: "<< forceSum[2]<<"\tAbs: "<< forceSum.length() << "\t[Pa]" << std::endl;
output << executionCounter_ << "\t" << currentPressureDrop_ << "\t[Pa]" << std::endl;
output.close();
}
}
}
......@@ -219,6 +224,7 @@ namespace FBfunc {
const uint_t checkFrequency_;
uint_t executionCounter_;
const real_t preFac_;
const bool preliminaryTimeloop_;
real_t currentPressureDrop_;
......@@ -1155,10 +1161,6 @@ namespace FBfunc {
}; //
////////////////
// MEM //
////////////////
template<typename LatticeModel_T, typename FlagField_T, typename BodyField_T>
class VelocityCheckMEM {
public:
......@@ -1268,115 +1270,6 @@ namespace FBfunc {
}; // VelocityCheckMEM
//TODO unfinished (see MEM)
template< typename LatticeModel_T, typename FlagField_T , typename BodyAndVolumeFractionField_T>
class VelocityCheckPSM
{
public:
typedef lbm::PdfField< LatticeModel_T > PdfField;
VelocityCheckPSM( const shared_ptr< StructuredBlockStorage > & blocks, const BlockDataID & fractionFieldID, const ConstBlockDataID & pdfFieldId, const ConstBlockDataID & flagFieldId,
const Set< FlagUID > & cellsToCheck, FBfunc::SetupFB & config) :
blocks_(blocks), fractionFieldID_( fractionFieldID ), pdfFieldId_( pdfFieldId ),
flagFieldId_( flagFieldId ), cellsToCheck_( cellsToCheck ),
uMax_(config.stopVel * config.stopVel), checkFrequency_((config.velCheckFreq > 0) ? config.velCheckFreq : uint_c(1)), dx_SI_(config.dx_SI), executionCounter_(uint_c(0)){
std::ofstream output;
output.open("MaximalLBMvelocity.txt", std::ofstream::out | std::ofstream::trunc);
output.close();
}
void operator()( const IBlock * const block )
{
++executionCounter_;
if( ( executionCounter_ - uint_c(1) ) % checkFrequency_ != 0 )
return;
const FlagField_T * flagField = block->getData< FlagField_T >( flagFieldId_ );
const PdfField * pdfField = block->getData< PdfField >( pdfFieldId_ );
const BodyAndVolumeFractionField_T * fractionField = block->getData<BodyAndVolumeFractionField_T>( fractionFieldID_ );
//const pe::BodyStorage * bodyStorage = block->getData< pe::BodyStorage >( bodyStorageID_ );
const CellInterval & cellBB = pdfField->xyzSize();
WALBERLA_ASSERT_EQUAL( flagField->xyzSize(), cellBB );
real_t currentMax = real_c(0.0);
typename FlagField_T::flag_t mask = 0;
for( auto flag = cellsToCheck_.begin(); flag != cellsToCheck_.end(); ++flag )
mask = static_cast< typename FlagField_T::flag_t >( mask | flagField->getFlag( *flag ) );
for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z ) {
for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y ) {
for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x ){
if( flagField->isPartOfMaskSet(x,y,z,mask) )
{
if( pdfField->getVelocity( x, y, z ).sqrLength() > uMax_ ){
real_t sumSolidFrac = real_c(0);
Cell globalPos;
blocks_->transformBlockLocalToGlobalCell(globalPos, *block , Cell(x,y,z));
std::ofstream output;
output.open("BrokenArea" + std::to_string(MPIManager::instance()->rank()) + ".txt");
output << "Position of maximum value is x: " << globalPos.x()*dx_SI_*1000 << " y: " << globalPos.y()*dx_SI_*1000 << " z: " << globalPos.z()*dx_SI_*1000 <<" [mm]"<< std::endl;
output << "Broken velocity is " << pdfField->getVelocity( x, y, z ).length() << std::endl << std::endl;
for ( auto dir = stencil::D3Q27::beginNoCenter(); dir!= stencil::D3Q27::end(); ++dir ){
if( flagField->isPartOfMaskSet(x + dir.cx(),y + dir.cy(),z + dir.cz(), mask)){
output << "Cell in direction\t"<< dir.cx() << "\t" << dir.cy() << "\t" << dir.cz() << "\thas velocity " << pdfField->getVelocity( x + dir.cx(), y + dir.cy(), z + dir.cz() ).length() << std::endl;
}
else{
output << "Cell in direction\t"<< dir.cx() << "\t" << dir.cy() << "\t" << dir.cz() << "\tis solid" << std::endl << std::endl;
}
}
for( auto bodyFracIt = fractionField->get(x,y,z).begin(); bodyFracIt != fractionField->get(x,y,z).end(); ++bodyFracIt ){
output << "Body with ID: "<< (*bodyFracIt).first << " has a solid volume fraction of: " << (*bodyFracIt).second << std::endl << std::endl;
sumSolidFrac += (*bodyFracIt).second;
}
output.close();
WALBERLA_CHECK( false , "Position of maximum value is x: "
<< globalPos.x()*dx_SI_*1000 << " y: " << globalPos.y()*dx_SI_*1000 << " z: " << globalPos.z()*dx_SI_*1000 <<" [mm]"
<< std::endl << "Total solid volume fraction:\t" << sumSolidFrac << std::endl );
}
if (currentMax < pdfField->getVelocity( x, y, z ).sqrLength()){
currentMax = pdfField->getVelocity( x, y, z ).sqrLength();
}
}
}
}
}
WALBERLA_MPI_SECTION(){
mpi::reduceInplace( currentMax, mpi::MAX );
}
WALBERLA_ROOT_SECTION(){
std::cout << "Maximal LBM Velocity: " << sqrt(currentMax) << std::endl;
std::ofstream output;
output.open( "MaximalLBMvelocity.txt", std::ofstream::out | std::ofstream::app);
output << executionCounter_ << "\t" << sqrt(currentMax) << std::endl;
output.close();
}
}
private:
shared_ptr< StructuredBlockStorage > blocks_;
const BlockDataID fractionFieldID_;
const ConstBlockDataID pdfFieldId_;
const ConstBlockDataID flagFieldId_;
const Set< FlagUID > cellsToCheck_;
const real_t uMax_;
const uint_t checkFrequency_;
const real_t dx_SI_;
uint_t executionCounter_;
}; // VelocityCheckPSM
// CALCULATE SOLID FRACTION
template<typename FlagField_T>
class FractionCalculatorMEM {
......
......@@ -3,14 +3,14 @@
Model
{
turbulence_model true;
turbulence_model false;
}
Geometry
{
width 0.10;
width 0.04;
depth 0.02;
height 0.40;
height 0.08;
offset_bot 0.00;
offset_top 0.00;
......@@ -25,8 +25,8 @@ Geometry
diameter_standard_distribution_B 0.00006;
cell_size 0.00020;
blocks_width 5;
blocks_height 20;
blocks_width 2;
blocks_height 4;
blocks_depth 1;
initial_sphere_distance 0.00400;
......@@ -39,18 +39,18 @@ Geometry
FluidizedBed
{
#particles_A 5000;
#particles_B 3000;
density_particle_A 22.5;
density_particle_B 22.5;
#particles_A 200;
#particles_B 100;
density_particle_A 5;
density_particle_B 5;
density_fluid 1.2;
dynamic_viscosity 1.84e-5;
gravity 9.81;
volume_flow_rate 2.5;
//volume_flow_rate 2.5;
//velocity 0.4;
//particle_RE 50.0;
particle_RE 10.0;
}
Simulation
......@@ -62,7 +62,7 @@ Simulation
path_to_checkpoint_pdf /simdata/schuster/i10hpc/dsfd/checkpoint_pdf.dat;
enable_checkpointing true;
evaluation_frequency 100;
evaluation_frequency 2;
velocity_check_frequency 10;
checkpoint_frequency 100000;
......@@ -77,18 +77,18 @@ Simulation
smagorinsky_constant 0.1;
simulation_time 8.0;
simulation_time 0.5;
//#time_steps 80000;
ramping_time 2.0;
ramping_time 0.5;
//#ramping_steps 20000;
start_ramping_scaling 0.75;
start_ramping_scaling 0.50;
run_preliminary_timeloop false;
eps 0.05;
run_preliminary_timeloop true;
eps 0.001;
create_packed_bed true;
check_maximum_particle_velocity true; //else check average particle velocity
settle_velocity 0.001;
check_maximum_particle_velocity true; // else check for average particle velocity
settle_velocity 0.02;
settle_check_frequency 200;
pe_accelerator 1;
fixed_bed false;
......@@ -96,6 +96,11 @@ Simulation
Material
{
use_DEM false; // else use HCSITS
HCSITS_max_iterations 10;
HCSITS_relaxation_parameter 0.7;
static_friction_particles_A 0.5;
static_friction_particles_B 0.5;
static_friction_particles_wall 0.5;
......@@ -165,7 +170,7 @@ VTK
frequency_spheres 200;
write_flag_field false;
frequency_flag_field 200000;
frequency_flag_field 10000000;
write_fluid_field false;
frequency_fluid_field 50000;
......@@ -174,4 +179,4 @@ VTK
write_omega false;
frequency_omega 50000;
sampling_resolution omega 1.0;
}
\ No newline at end of file
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment