Reply To: Particle periodicity running with MPI.
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Particle periodicity running with MPI. › Reply To: Particle periodicity running with MPI.
May 15, 2024 at 10:02 am
#8697
Rookie
Participant
Dear Jan,
Thank you for your patience with this issue. Below is my code. To speed up the simulation for you to test, I’ve shortened the simulation time. Please note that you must change the periodic boundaries to the code I posted earlier in order for it to compile successfully. If you enable parallel mode, you can join me in discovering issues together. Actually, regarding this code, I found that my “geometry2” should also contain particles, and particle collisions should occur in a layer of cells outside “geometry2”. However, it seems that the bidirectional coupling setting cannot be applied to two regions simultaneously.
#include "olb3D.h"
#include "olb3D.hh" // include full template code
using namespace olb;
using namespace olb::descriptors;
using namespace olb::util;
using namespace olb::graphics;
typedef double T;
typedef WallFunctionForcedD3Q19Descriptor DESCRIPTOR;
#define PARTICLE Particle3D
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
// Mathmatical constants
const T pi = util::acos(-1);
// Parameters for the simulation setup
const int N = 59;
const T physRefL = 0.02; // half channel height in meters
const T lx = 2. * pi * physRefL; // streamwise length in meters
const T ly = 2. * physRefL; // spanwise length in meters
const T lz = 2. * physRefL; // wall-normal length in meters
const T radius = 2.5e-5; // particles radius
const T partRho = 2500.; // particles density
// Choose friction reynolds number ReTau
// #define Case_ReTau_1000
// #define Case_ReTau_2000
#define Case_ReTau_180
// Wallfunction parameters
const T latticeWallDistance = 0.5; // lattice distance to boundary
const int rhoMethod = 2; // method for density reconstruction
// 0: Zou-He
// 1: extrapolation
// 2: constant
const int fneqMethod = 0; // method for fneq reconstruction
// 0: regularized NEBB (Latt)
// 1: extrapolation NEQ (Guo Zhaoli)
// 2: regularized second order finite Differnce
// 3: equilibrium scheme
const int wallProfile = 0; // wallfunction profile
// 0: Musker profile
// 1: power law profile
// Reynolds number based on the friction velocity
#if defined(Case_ReTau_1000)
T ReTau = 1000.512;
#elif defined(Case_ReTau_2000)
T ReTau = 1999.756;
#elif defined(Case_ReTau_180)
T ReTau = 180;
#endif
// Characteristic physical kinematic viscosity from DNS Data
// http://turbulence.ices.utexas.edu/channel2015/data/LM_Channel_1000_mean_prof.dat
#if defined(Case_ReTau_1000)
T charPhysNu = 5. / 100000.;
#elif defined(Case_ReTau_2000)
T charPhysNu = 2.3 / 100000.;
#elif defined(Case_ReTau_180)
T charPhysNu = 1.5 / 100000.;
#endif
// number of forcing updates over simulation time
#if defined(Case_ReTau_1000)
T fluxUpdates = 2000;
#elif defined(Case_ReTau_2000)
T fluxUpdates = 4000;
#elif defined(Case_ReTau_180)
T fluxUpdates = 2000;
#endif
// physical simulated length adapted for lattice distance to boundary in meters
const T adaptedPhysSimulatedLength = 2 * physRefL - 2 * ((0.04 / T(N + 2 * latticeWallDistance)) * latticeWallDistance);
const T Re = 6594;
// Characteristic physical mean bulk velocity from Dean correlations in meters - Malaspinas and Sagaut (2014)
const T charPhysU = Re * charPhysNu / (2. * physRefL);
// Time of the simulation in seconds
const T charPhysT = physRefL / (ReTau * charPhysNu / physRefL);
const T physConvergeTime = 12. * charPhysT; // time until until statistics sampling in seconds
const T physStatisticsTime = 4. * charPhysT; // statistics sampling time in seconds
const T particleMaxPhysT = 1. * charPhysT;
const T singleMaxPhysT = physConvergeTime + physStatisticsTime;
const T fluidMaxPhysT = physConvergeTime + physStatisticsTime + particleMaxPhysT; // max. simulation time in seconds
const T statisticsSave = 1. / 25.; // time between statistics samples in seconds
const int noOfParticles = 30000;
const T checkstatistics = (T)fluidMaxPhysT / 200.;
// seed the rng with time if SEED_WITH_TIME is set, otherwise just use a fixed seed.
#if defined SEED_WITH_TIME
#include <chrono>
auto seed = std::chrono::system_clock().now().time_since_epoch().count();
std::default_random_engine generator(seed);
#else
std::default_random_engine generator(0x1337533DAAAAAAAA);
#endif
// Compute mean lattice velocity from musker wallfunction
T computeLatticeVelocity()
{
T Ma_max = 0.1;
T c_s = 1 / util::sqrt(3.0);
T latticeUMax = Ma_max * c_s;
Musker<T, T> musker_tmp(charPhysNu, physRefL, 1.205);
T charPhysU_tau = ReTau * charPhysNu / physRefL;
T tau_w[1];
tau_w[0] = util::pow(charPhysU_tau, 2.);
T charPhysUMax[1];
musker_tmp(charPhysUMax, tau_w);
T latticeU = charPhysU * latticeUMax / charPhysUMax[0];
return latticeU;
}
template <typename T, typename S>
class Channel3D : public AnalyticalF3D<T, S>
{
protected:
T turbulenceIntensity;
T maxVelocity;
T distanceToWall;
T obst_z;
T obst_r;
T a;
T b;
public:
Channel3D(UnitConverter<T, DESCRIPTOR> const &converter, T frac) : AnalyticalF3D<T, S>(3)
{
turbulenceIntensity = 0.05;
distanceToWall = -converter.getPhysDeltaX() / 2.;
maxVelocity = converter.getLatticeVelocity(converter.getCharPhysVelocity() * (8. / 7.)); // Centerline Velocity
obst_z = physRefL + distanceToWall;
obst_r = physRefL;
a = -1.;
b = 1.;
};
bool operator()(T output[], const S input[])
{
std::uniform_real_distribution<T> distribution(a, b);
T nRandom1 = distribution(generator);
T nRandom2 = distribution(generator);
T nRandom3 = distribution(generator);
T u_calc = maxVelocity * util::pow(((obst_r - util::abs(input[2] - obst_z)) / obst_r), 1. / 7.);
output[0] = turbulenceIntensity * nRandom1 * maxVelocity + u_calc;
output[1] = turbulenceIntensity * nRandom2 * maxVelocity;
output[2] = turbulenceIntensity * nRandom3 * maxVelocity;
return true;
};
};
template <typename T, typename S>
class TrackedForcing3D : public AnalyticalF3D<T, S>
{
protected:
T um;
T utau;
T h2;
T aveVelocity;
public:
TrackedForcing3D(UnitConverter<T, DESCRIPTOR> const &converter, int ReTau) : AnalyticalF3D<T, S>(3)
{
um = converter.getCharPhysVelocity();
utau = ReTau * converter.getPhysViscosity() / (converter.getCharPhysLength() / 2.);
h2 = converter.getCharPhysLength() / 2.;
aveVelocity = um;
};
void updateAveVelocity(T newVel)
{
aveVelocity = newVel;
}
bool operator()(T output[], const S input[])
{
output[0] = util::pow(utau, 2) / h2 + (um - aveVelocity) * um / h2;
output[1] = 0;
output[2] = 0;
return true;
};
};
void prepareGeometry(SuperGeometry<T, 3> &superGeometry, IndicatorF3D<T> &indicator,
UnitConverter<T, DESCRIPTOR> const &converter)
{
OstreamManager clout(std::cout, "prepareGeometry");
clout << "Prepare Geometry ..." << std::endl;
superGeometry.rename(0, 2, indicator);
superGeometry.rename(2, 1, {0, 0, 1});
superGeometry.clean();
superGeometry.innerClean();
superGeometry.checkForErrors();
superGeometry.print();
olb::Vector<T, 3> PhyMax = superGeometry.getStatistics().getMaxPhysR(2);
olb::Vector<T, 3> PhyMin = superGeometry.getStatistics().getMinPhysR(2);
clout << "Dimension of the channel in meters: x = " << PhyMax[0] - PhyMin[0];
clout << " ; y = " << PhyMax[1] - PhyMin[1];
clout << " ; z = " << PhyMax[2] - PhyMin[2] << std::endl;
clout << "Prepare Geometry ... OK" << std::endl;
}
// set up initial conditions
void setInitialConditions(SuperLattice<T, DESCRIPTOR> &sLattice,
UnitConverter<T, DESCRIPTOR> const &converter,
SuperGeometry<T, 3> &superGeometry,
AnalyticalScaled3D<T, T> &forceSolScaled,
TrackedForcing3D<T, T> &forceSol)
{
OstreamManager clout(std::cout, "setInitialConditions");
clout << "Set initial conditions ..." << std::endl;
AnalyticalConst3D<T, T> rho(1.205);
Channel3D<T, T> uSol(converter, 1.);
sLattice.defineRhoU(superGeometry, 1, rho, uSol);
sLattice.iniEquilibrium(superGeometry, 1, rho, uSol);
sLattice.defineRhoU(superGeometry, 2, rho, uSol);
sLattice.iniEquilibrium(superGeometry, 2, rho, uSol);
AnalyticalConst3D<T, T> TauEff(1. / converter.getLatticeRelaxationFrequency());
sLattice.defineField<TAU_EFF>(superGeometry, 1, TauEff);
sLattice.defineField<TAU_EFF>(superGeometry, 2, TauEff);
// Force Initialization
forceSol.updateAveVelocity(converter.getCharPhysVelocity()); // New average velocity
// Initialize force
sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
// Tau_w Initialization
T tau_w_guess = 0.0; // Wall shear stress in phys units
AnalyticalConst3D<T, T> tau_w_ini(tau_w_guess);
AnalyticalScaled3D<T, T> tau_w_ini_scaled(tau_w_ini, 1. / (converter.getConversionFactorForce() * util::pow(converter.getConversionFactorLength(), 2.)));
sLattice.defineField<TAU_W>(superGeometry, 1, tau_w_ini_scaled);
sLattice.defineField<TAU_W>(superGeometry, 2, tau_w_ini_scaled);
clout << "Set initial conditions ... OK" << std::endl;
}
// Set up the geometry of the simulation
void prepareLattice(SuperLattice<T, DESCRIPTOR> &sLattice,
UnitConverter<T, DESCRIPTOR> const &converter,
SuperGeometry<T, 3> &superGeometry,
AnalyticalScaled3D<T, T> &forceSolScaled,
TrackedForcing3D<T, T> &forceSol,
wallFunctionParam<T> const &wallFunctionParam)
{
OstreamManager clout(std::cout, "prepareLattice");
clout << "Prepare Lattice ..." << std::endl;
/// Material=1 -->bulk dynamics
sLattice.defineDynamics<SmagorinskyForcedBGKdynamics>(superGeometry, 1);
/// Material = 2 --> boundary node + wallfunction
sLattice.defineDynamics<ExternalTauEffLESForcedBGKdynamics>(superGeometry, 2);
setWallFunctionBoundary<T, DESCRIPTOR>(sLattice, superGeometry, 2, converter, wallFunctionParam);
/// === Set Initial Conditions == ///
setInitialConditions(sLattice, converter, superGeometry, forceSolScaled, forceSol);
sLattice.setParameter<descriptors::OMEGA>(converter.getLatticeRelaxationFrequency());
sLattice.setParameter<collision::LES::Smagorinsky>(0.1);
// Make the lattice ready for simulation
sLattice.initialize();
clout << "Prepare Lattice ... OK" << std::endl;
}
/// Computes the pressure drop between the voxels before and after the cylinder
bool getResults(SuperLattice<T, DESCRIPTOR> &sLattice,
UnitConverter<T, DESCRIPTOR> const &converter, size_t iT, int iTperiod,
SuperGeometry<T, 3> &superGeometry, Timer<double> &fluidTimer,
SuperParticleSystem3D<T, PARTICLE> &supParticleSystem,
T radii, T partRho, Timer<double> &particleTimer,
SuperParticleSysVtuWriter<T, PARTICLE> &supParticleWriter,
bool fluidExists, SuperLatticeTimeAveragedF3D<T> &sAveragedVel)
{
OstreamManager clout(std::cout, "getResults");
std::list<int> materialslist;
materialslist.push_back(1);
materialslist.push_back(2);
using BulkDynamics = ShearSmagorinskyBGKdynamics<T,DESCRIPTOR>;
ParametersOfOperatorD<T,DESCRIPTOR,BulkDynamics> bulkDynamicsParams{};
SuperVTMwriter3D<T> vtmWriter("channel3d");
SuperVTMwriter3D<T> vtmWriterStartTime("startingTimechannel3d");
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
SuperLatticePhysPressure3D<T, DESCRIPTOR> pressure(sLattice, converter);
vtmWriter.addFunctor(geometry);
vtmWriter.addFunctor(velocity);
vtmWriter.addFunctor(pressure);
std::size_t singleMaxT = converter.getLatticeTime(singleMaxPhysT);
if (iT == 0)
{
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLattice, superGeometry);
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLattice);
SuperLatticeRank3D<T, DESCRIPTOR> rank(sLattice);
vtmWriter.write(geometry);
vtmWriter.write(cuboid);
vtmWriter.write(rank);
vtmWriter.createMasterFile();
vtmWriterStartTime.createMasterFile();
// Print some output of the chosen simulation setup
clout << "N=" << N << "; maxTimeSteps(fluid)="
<< converter.getLatticeTime(fluidMaxPhysT) << "; noOfCuboid="
<< superGeometry.getCuboidGeometry().getNc() << "; Re=" << Re
<< "; noOfParticles=" << noOfParticles << "; maxTimeSteps(particle)="
<< converter.getLatticeTime(particleMaxPhysT)
<< "; St=" << (2. * partRho * radius * radius * converter.getCharPhysVelocity()) / (9. * converter.getPhysViscosity() * converter.getPhysDensity() * converter.getCharPhysLength()) << std::endl;
}
// Writes output on the console for the fluid phase
if (iT < converter.getLatticeTime(fluidMaxPhysT) && iT % iTperiod == 0)
{
// Timer console output
fluidTimer.update(iT);
fluidTimer.printStep(2);
// Lattice statistics console output
sLattice.getStatistics().print(iT, iT * converter.getPhysDeltaT());
clout << "Max. physical velocity(m/s): " << converter.getPhysVelocity(sLattice.getStatistics().getMaxU()) << std::endl;
clout << "Max u+:" << converter.getPhysVelocity(sLattice.getStatistics().getMaxU()) / (ReTau * charPhysNu / (converter.getCharPhysLength() / 2.)) << std::endl;
}
if (iT < converter.getLatticeTime(fluidMaxPhysT) && iT % converter.getLatticeTime(statisticsSave) == 0 && iT > converter.getLatticeTime(physConvergeTime))
{
// Add ensemble to temporal averaged velocity
sLattice.communicate();
sAveragedVel.addEnsemble();
}
if (iT < converter.getLatticeTime(singleMaxPhysT) && (iT % converter.getLatticeTime(singleMaxPhysT / 16) == 0 || iT == converter.getLatticeTime(fluidMaxPhysT) - 1))
{
// Writes the vtk files
vtmWriter.write(iT);
}
// Writes output on the console for the fluid phase
if (iT >= converter.getLatticeTime(singleMaxPhysT) &&
(iT % (iTperiod) == 0 || iT == converter.getLatticeTime(singleMaxPhysT) || iT == converter.getLatticeTime(fluidMaxPhysT) - 1))
{
vtmWriter.write(iT);
particleTimer.print(iT - singleMaxT);
// console output number of particles at different material numbers mat
supParticleSystem.print({1, 2});
// only write .vtk-files after the fluid calculation is finished
supParticleWriter.write(iT - singleMaxT);
// true as long as certain amount of active particles
if (supParticleSystem.globalNumOfActiveParticles() < 0.0001 * noOfParticles && iT > 0.9 * converter.getLatticeTime(fluidMaxPhysT))
{
return false;
}
}
return true;
}
int main(int argc, char *argv[])
{
/// === 1st Step: Initialization ===
olbInit(&argc, &argv);
singleton::directories().setOutputDir("./tmp/");
OstreamManager clout(std::cout, "main");
// display messages from every single mpi process
// clout.setMultiOutput(true);
UnitConverterFromResolutionAndLatticeVelocity<T, DESCRIPTOR> converter(
int{N}, // resolution: number of voxels per charPhysL
(T)computeLatticeVelocity(), // latticeU : mean lattice velocity
(T)adaptedPhysSimulatedLength, // charPhysLength: reference length of simulation geometry
(T)2.19, // charPhysVelocity: mean bulk velocity in __m / s__
(T)charPhysNu, // physViscosity: physical kinematic viscosity in __m^2 / s__
(T)1.205 // physDensity: physical density in __kg / m^3__
);
converter.print();
converter.write("channelpraticle");
clout << "----------------------------------------------------------------------" << std::endl;
clout << "Converge time(s): " << physConvergeTime << std::endl;
clout << "Lattice converge time: " << converter.getLatticeTime(physConvergeTime) << std::endl;
clout << "Max. Phys. simulation time(s): " << fluidMaxPhysT << std::endl;
clout << "Max. Lattice simulation time: " << converter.getLatticeTime(fluidMaxPhysT) << std::endl;
clout << "Frequency Statistics Save(Hz): " << 1. / statisticsSave << std::endl;
clout << "Statistics save period(s): " << statisticsSave << std::endl;
clout << "Lattice statistics save period: " << converter.getLatticeTime(statisticsSave) << std::endl;
clout << "----------------------------------------------------------------------" << std::endl;
clout << "Channel height(m): " << adaptedPhysSimulatedLength << std::endl;
clout << "y+ value: " << (ReTau * converter.getPhysViscosity() / (physRefL)) * ((0.04 / T(N + 2 * latticeWallDistance)) * latticeWallDistance) / converter.getPhysViscosity() << std::endl;
clout << "y+ value spacing: " << (ReTau * converter.getPhysViscosity() / (physRefL)) * (converter.getPhysDeltaX()) / converter.getPhysViscosity() << std::endl;
clout << "----------------------------------------------------------------------" << std::endl;
clout << "N=" << N << "; maxTimeSteps(fluid)="
<< converter.getLatticeTime(fluidMaxPhysT) << "; Re=" << Re
<< "; noOfParticles=" << noOfParticles << "; maxTimeSteps(particle)="
<< converter.getLatticeTime(particleMaxPhysT) << "; singleMaxPhysT="
<< converter.getLatticeTime(singleMaxPhysT)
<< "; St=" << (2. * partRho * radius * radius * converter.getCharPhysVelocity()) / (9. * converter.getPhysViscosity() * converter.getPhysDensity() * converter.getCharPhysLength()) << std::endl;
clout << "----------------------------------------------------------------------" << std::endl;
clout << "iTperiod=" << converter.getLatticeTime(checkstatistics) << std::endl;
clout << "utau=" << ReTau * converter.getPhysViscosity() / (converter.getCharPhysLength() / 2.) << std::endl;
clout << "----------------------------------------------------------------------" << std::endl;
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = singleton::mpi().getSize();
#else
const int noOfCuboids = 1;
#endif
Vector<T, 3> extend(lx, ly, adaptedPhysSimulatedLength);
extend[2] += (1. / 8.) * converter.getPhysDeltaX();
Vector<T, 3> origin(0., 0., 0.);
IndicatorCuboid3D<T> cuboid(extend, origin);
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids);
cuboidGeometry.setPeriodicity(true, true, false);
HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry<T, 3> superGeometry(cuboidGeometry, loadBalancer, 3);
prepareGeometry(superGeometry, cuboid, converter);
clout << "noOfCuboid=" << superGeometry.getCuboidGeometry().getNc() << std::endl;
/// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLattice(superGeometry);
// forcing of the channel
TrackedForcing3D<T, T> forceSol(converter, ReTau);
AnalyticalScaled3D<T, T> forceSolScaled(forceSol, 1. / (converter.getConversionFactorForce() / converter.getConversionFactorMass()));
int input[3];
T output[5];
Vector<T, 3> normal(1, 0, 0);
std::vector<int> normalvec{1, 0, 0};
Vector<T, 3> center;
for (int i = 0; i < 3; ++i)
{
center[i] = origin[i] + extend[i] / 2;
}
SuperLatticePhysVelocity3D<T, DESCRIPTOR> velocity(sLattice, converter);
std::vector<int> materials;
materials.push_back(1);
materials.push_back(2);
std::list<int> materialslist;
materialslist.push_back(1);
materialslist.push_back(2);
wallFunctionParam<T> wallFunctionParam;
wallFunctionParam.bodyForce = true;
wallFunctionParam.wallProfile = wallProfile;
wallFunctionParam.rhoMethod = rhoMethod;
wallFunctionParam.fneqMethod = fneqMethod;
wallFunctionParam.latticeWalldistance = latticeWallDistance;
wallFunctionParam.vonKarman = 0.4;
wallFunctionParam.curved = false;
prepareLattice(sLattice, converter, superGeometry,
forceSolScaled, forceSol, wallFunctionParam);
SuperPlaneIntegralFluxVelocity3D<T> velFlux(sLattice, converter, superGeometry, center, normal, materials,
BlockDataReductionMode::Discrete);
// === 3.1 Step: Particles ===
clout << "Prepare Particles ..." << std::endl;
// SuperParticleSystems3D
SuperParticleSystem3D<T, PARTICLE> supParticleSystem(superGeometry);
// define which properties are to be written in output data
SuperParticleSysVtuWriter<T, PARTICLE> supParticleWriter(supParticleSystem,
"particles", SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::velocity | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::mass | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::radius
| SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::active | SuperParticleSysVtuWriter<T, PARTICLE>::particleProperties::force);
SuperLatticeInterpPhysVelocity3D<T, DESCRIPTOR> getVel(sLattice, converter);
T dynVisc = converter.getPhysViscosity() * converter.getPhysDensity();
T physDensity = converter.getPhysDensity();
auto schillerNaumannDragForce = std::make_shared<SchillerNaumannDragForce3D<T, PARTICLE, DESCRIPTOR>>(getVel, dynVisc, physDensity);
supParticleSystem.addForce(schillerNaumannDragForce);
std::vector<T> direction{1, 0, 0};
const T g = 9.81;
auto weightForce = std::make_shared<WeightForce3D<T, PARTICLE>>(direction, g);
supParticleSystem.addForce(weightForce);
T dT = converter.getConversionFactorTime();
std::set<int> reflBMat = {2};
auto materialreflectBoundary = std::make_shared<SimpleReflectBoundary3D<T, PARTICLE>>(dT, superGeometry, reflBMat);
supParticleSystem.addBoundary(materialreflectBoundary);
auto materialperiodicBoundary = std::make_shared
< PeriodicBoundary3D<T, PARTICLE>
> (superGeometry, true, true, false);
supParticleSystem.addBoundary(materialperiodicBoundary);
supParticleSystem.setOverlap(2. * converter.getConversionFactorLength());
auto dragModel = std::make_shared<SchillerNaumannDragModel<T, DESCRIPTOR, PARTICLE>>(converter);
NaiveForwardCouplingModel<T, DESCRIPTOR, PARTICLE> forwardCoupling(converter, sLattice, superGeometry, dragModel);
int overlap = 2;
LocalBackCouplingModel<T, DESCRIPTOR, PARTICLE> backCoupling(converter, sLattice, superGeometry, overlap);
int subSteps = 10;
// particles generation at inlet3
Vector<T, 3> extendP = {lx, ly, adaptedPhysSimulatedLength};
Vector<T, 3> originP = {0, 0, 0};
IndicatorCuboid3D<T> inletCuboid(extendP, originP);
supParticleSystem.addParticle(inletCuboid, 4. / 3. * M_PI * util::pow(radius, 3) * partRho, radius, noOfParticles);
clout << "Prepare Particles ... OK" << std::endl;
/// === 5th Step: Definition of turbulent Statistics Objects ===
SuperLatticePhysVelocity3D<T, DESCRIPTOR> sVel(sLattice, converter);
SuperLatticeTimeAveragedF3D<T> sAveragedVel(sVel);
SuperLatticePhysPressure3D<T, DESCRIPTOR> sPre(sLattice, converter);
/// === 4th Step: Main Loop with Timer ===
Timer<double> fluidTimer(converter.getLatticeTime(fluidMaxPhysT), superGeometry.getStatistics().getNvoxel());
Timer<double> particleTimer(converter.getLatticeTime(particleMaxPhysT),
noOfParticles);
fluidTimer.start();
std::size_t iT = 0;
int iTperiod = converter.getLatticeTime(checkstatistics);
bool fluidExists = true;
// checks whether there is already data of the fluid from an earlier calculation
if (!(sLattice.load("fluidSolution")))
{
fluidExists = false;
for (; iT <= converter.getLatticeTime(singleMaxPhysT); ++iT)
{
if (iT % converter.getLatticeTime(fluidMaxPhysT / fluxUpdates) == 0 || iT == 0)
{
velFlux(output, input);
T flux = output[0];
T area = output[1];
forceSol.updateAveVelocity(flux / area);
sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
}
/// === 6th Step: Computation and Output of the Results ===
getResults(sLattice, converter, iT, iTperiod, superGeometry, fluidTimer,
supParticleSystem, radius, partRho, particleTimer,
supParticleWriter, fluidExists, sAveragedVel);
/// === 7th Step: Collide and Stream Execution ===
sLattice.collideAndStream();
}
sLattice.save("fluidSolution");
}
else
{
iT = converter.getLatticeTime(singleMaxPhysT);
getResults(sLattice, converter, iT,
iTperiod, superGeometry, fluidTimer,
supParticleSystem, radius, partRho, particleTimer,
supParticleWriter, fluidExists, sAveragedVel);
}
// when the fluid simulation time reaches singleMaxPhysT, the particle simulation begins
supParticleSystem.setVelToFluidVel( getVel );
particleTimer.start();
for ( ; iT <= converter.getLatticeTime( fluidMaxPhysT ); ++iT ) {
if (iT % converter.getLatticeTime(fluidMaxPhysT / fluxUpdates) == 0)
{
velFlux(output, input);
T flux = output[0];
T area = output[1];
forceSol.updateAveVelocity(flux / area);
sLattice.defineField<FORCE>(superGeometry, 1, forceSolScaled);
sLattice.defineField<FORCE>(superGeometry, 2, forceSolScaled);
}
// particles simulation starts after run up time is over
// supParticleSystem.simulate( converter.getConversionFactorTime());
supParticleSystem.simulateWithTwoWayCoupling_Mathias ( dT, forwardCoupling, backCoupling, 1, subSteps, true);
if ( !getResults( sLattice, converter, iT,
iTperiod, superGeometry, fluidTimer,
supParticleSystem, radius, partRho, particleTimer,
supParticleWriter, fluidExists, sAveragedVel) )
{
break;
}
/// === 7th Step: Collide and Stream Execution ===
sLattice.collideAndStream();
}
fluidTimer.stop();
fluidTimer.printSummary();
particleTimer.stop();
particleTimer.printSummary();
}
Best regards,
Rookie
- This reply was modified 2 months, 1 week ago by Rookie.