Skip to content
Snippets Groups Projects
Commit 3282d7ea authored by Markus Holzer's avatar Markus Holzer
Browse files

[skip ci] First Draft Turbulent inflow

parent 721943d3
No related branches found
No related tags found
No related merge requests found
Pipeline #72627 skipped
Hide whitespace changes
Inline Side-by-side
apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
+ 141
− 11
View file @ 3282d7ea
...@@ -29,6 +29,7 @@ ...@@ -29,6 +29,7 @@
#include "core/logging/Logging.h" #include "core/logging/Logging.h"
#include "core/math/Vector3.h" #include "core/math/Vector3.h"
#include "core/math/Constants.h" #include "core/math/Constants.h"
#include "core/math/Random.h"
#include "core/Environment.h" #include "core/Environment.h"
#include "core/mpi/MPIManager.h" #include "core/mpi/MPIManager.h"
#include "core/MemoryUsage.h" #include "core/MemoryUsage.h"
...@@ -109,6 +110,142 @@ inline Vector3< uint_t > inflowSliceSize( const shared_ptr< StructuredBlockStora ...@@ -109,6 +110,142 @@ inline Vector3< uint_t > inflowSliceSize( const shared_ptr< StructuredBlockStora
return Vector3<uint_t>( uint_c(1), blocks->getNumberOfYCells( *block ), blocks->getNumberOfZCells( *block ) ); return Vector3<uint_t>( uint_c(1), blocks->getNumberOfYCells( *block ), blocks->getNumberOfZCells( *block ) );
} }
template< typename Field_T >
class Inflow
{
public:
Inflow(const std::shared_ptr< StructuredBlockForest >& blocks, const Setup& setup, const IDs& ids)
: blocks_(blocks), setup_(setup), ids_(ids), executionCounter_(uint_c(0))
{
nu_ = setup_.kinViscosity;
urms_ = setup_.urms;
le_ = setup_.characteristicLength;
n_ = setup_.numberOfFourrierModes;
dt_ = setup_.dt;
uc_ = setup_.uc;
alpha_ = (real_c(55.0) * std::tgamma(real_c(5.0) / real_c(6.0))) / (real_c(9.0) * std::sqrt(math::pi) * std::tgamma(real_c(1.0) / real_c(3.0)));
kappaVKP_ = (real_c(9.0) * math::pi * alpha_) / (real_c(55.0) * le_);
epsilonVKP_ = (std::pow(alpha_ * std::tgamma(real_c(2.0)/real_c(3.0)), (real_c(3.0)/real_c(2.0)))) * (real_c(0.5) * kappaVKP_ * std::pow(urms_, real_c(3.0)));
kappaKolmogorov_ = std::pow(epsilonVKP_ / std::pow(nu_, real_c(3.0)), real_c(1.0)/real_c(4.0));
s_ = std::log10(kappaVKP_ / real_c(2.0));
e_ = std::log10(kappaKolmogorov_ * real_c(2.0));
firstMode_ = std::pow(real_c(10.0), s_);
lastMode_ = std::pow(real_c(10.0), e_);
deltaKappa_ = (lastMode_ - firstMode_) / (real_c(n_) - real_c(1.0));
logging();
}
void operator()(){
++executionCounter_;
// only evaluate in given intervals
// if (executionCounter_ % interval_ != uint_c(0)) { return; }
// update the inflow velocity plane
update();
}
void logging(){
WALBERLA_LOG_INFO_ON_ROOT( "Using turbulent injection with:" << std::setprecision(16) <<
"\n + kin. viscosity: " << nu_ << " [m^2/s]" <<
"\n + uRMS: " << urms_ <<
"\n + characteristicLength: " << le_ <<
"\n + numberOfFourrierModes: " << n_ <<
"\n + alpha: " << alpha_ <<
"\n + kappa VKP: " << kappaVKP_ <<
"\n + epsilon VKP: " << epsilonVKP_ <<
"\n + kappa Kolmogorov: " << kappaKolmogorov_ <<
"\n + firstMode: " << firstMode_ <<
"\n + lastMode: " << lastMode_)
}
private:
const std::shared_ptr< StructuredBlockForest > blocks_;
const Setup setup_;
const IDs ids_;
real_t nu_;
real_t urms_;
real_t le_;
uint_t n_;
real_t alpha_;
real_t kappaVKP_;
real_t epsilonVKP_;
real_t kappaKolmogorov_;
real_t dt_;
Vector3<real_t> uc_;
real_t s_;
real_t e_;
real_t firstMode_;
real_t lastMode_;
real_t deltaKappa_;
uint_t executionCounter_; // number of times operator() has been called
void update(){
WALBERLA_LOG_INFO_ON_ROOT("Update inflow")
for (auto& block : *blocks_){
auto velocityFieldInflowSlice = block.template getData< Field_T >(ids_.velocityFieldInflowSlice);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(velocityFieldInflowSlice, Cell globalCell;
velocityFieldInflowSlice->get(x, y, z, 0) = 0.05;
)
}
gpu::fieldCpy< gpu::GPUField< typename Field_T::value_type >, Field_T >(blocks_, ids_.velocityFieldInflowSliceGPU, ids_.velocityFieldInflowSlice);
}
real_t getKarmanPaoEnergy(const real_t kappa){
const real_t a1 = alpha_ * std::pow(urms_, real_c(2.0)) / kappaVKP_;
const real_t a2 = std::pow(kappa / kappaVKP_, real_c(4.0));
const real_t a3 = std::pow(real_c(1.0) + std::pow(kappa / kappaVKP_, real_c(2.0)), real_c(17.0) / real_c(6.0));
const real_t a4 = std::exp(real_c(-2.0) * std::pow(kappa / kappaKolmogorov_, real_c(2.0)));
return a1 * a2 / a3 * a4;
}
Vector3<real_t> summOfFourrierModes(Vector3<real_t> p){
Vector3<real_t> result(real_c(0.0), real_c(0.0), real_c(0.0));
for(uint_t i = 0; i < n_; ++i){
const real_t a = math::realRandom(real_c(0.0), real_c(1.0));
const real_t theta = std::acos(real_c(2.0) * a - real_c(1.0));
const real_t phi = real_c(2.0) * math::pi * math::realRandom(real_c(0.0), real_c(1.0));
const real_t psi = real_c(2.0) * math::pi * math::realRandom(real_c(0.0), real_c(1.0));
const real_t kappaN = firstMode_ + real_c(0.5) * deltaKappa_ + real_c(i) * deltaKappa_;
Vector3<real_t> kappaVector(std::cos(phi) * std::sin(theta) * kappaN,
std::sin(theta) * std::sin(phi) * kappaN,
std::cos(theta) * kappaN);
Vector3<real_t> zetta(math::realRandom(real_c(-1.0), real_c(1.0)),
math::realRandom(real_c(-1.0), real_c(1.0)),
math::realRandom(real_c(-1.0), real_c(1.0)));
zetta = zetta.getNormalized();
Vector3<real_t> crossProduct = kappaVector%=zetta;
Vector3<real_t> sigma = crossProduct / crossProduct.getNormalized();
const real_t ut = std::sqrt(getKarmanPaoEnergy(kappaN) * deltaKappa_);
const real_t t = executionCounter_ * dt_;
const Vector3<real_t> x = p - t*uc_;
const real_t arg = std::cos(kappaVector * x + psi);
result += real_c(2.0) * ut * arg * sigma;
}
return result;
}
}; // class Inflow
int main(int argc, char** argv){ int main(int argc, char** argv){
Environment env( argc, argv ); Environment env( argc, argv );
mpi::MPIManager::instance()->useWorldComm(); mpi::MPIManager::instance()->useWorldComm();
...@@ -123,12 +260,10 @@ int main(int argc, char** argv){ ...@@ -123,12 +260,10 @@ int main(int argc, char** argv){
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
auto parameters = config->getOneBlock("Parameters"); auto parameters = config->getOneBlock("Parameters");
auto tpmsParameters = config->getOneBlock("TPMS");
auto domainParameters = config->getOneBlock("DomainSetup");
auto loggingParameters = config->getOneBlock("Logging"); auto loggingParameters = config->getOneBlock("Logging");
auto boundariesConfig = config->getBlock("Boundaries"); auto boundariesConfig = config->getBlock("Boundaries");
Setup setup(parameters, tpmsParameters, domainParameters, loggingParameters, infoMap); Setup setup(config, infoMap);
TPMS tpms(setup); TPMS tpms(setup);
WALBERLA_LOG_INFO_ON_ROOT( "Simulation run data:" WALBERLA_LOG_INFO_ON_ROOT( "Simulation run data:"
"\n- simulation parameters:" << std::setprecision(16) << "\n- simulation parameters:" << std::setprecision(16) <<
...@@ -194,14 +329,7 @@ int main(int argc, char** argv){ ...@@ -194,14 +329,7 @@ int main(int argc, char** argv){
} }
sweepCollection.initialiseBlockPointer(); sweepCollection.initialiseBlockPointer();
for (auto& block : *blocks){ Inflow<VelocityField_T> inflow(blocks, setup, ids);
auto velocityFieldInflowSlice = block.getData< VelocityField_T >(ids.velocityFieldInflowSlice);
WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(velocityFieldInflowSlice, Cell globalCell;
velocityFieldInflowSlice->get(x, y, z, 0) = setup.latticeVelocity;
)
}
gpu::fieldCpy< gpu::GPUField< VelocityField_T::value_type >, VelocityField_T >(blocks, ids.velocityFieldInflowSliceGPU, ids.velocityFieldInflowSlice);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// LB SWEEPS AND BOUNDARY HANDLING /// /// LB SWEEPS AND BOUNDARY HANDLING ///
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
...@@ -250,6 +378,8 @@ int main(int argc, char** argv){ ...@@ -250,6 +378,8 @@ int main(int argc, char** argv){
auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority); auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority);
SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps); SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
timeloop.addFuncBeforeTimeStep(inflow, "Update inflow velocity");
lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T > lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T >
LBMMeshRefinement(blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo); LBMMeshRefinement(blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
......
apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
+ 9
− 2
View file @ 3282d7ea
...@@ -9,8 +9,6 @@ Parameters ...@@ -9,8 +9,6 @@ Parameters
timesteps 11; timesteps 11;
useGrid false;
processMemoryLimit 512; // MiB processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>; innerOuterSplit <1, 1, 1>;
...@@ -18,6 +16,15 @@ Parameters ...@@ -18,6 +16,15 @@ Parameters
gpuEnabledMPI false; gpuEnabledMPI false;
} }
Turbulence
{
useGrid false;
U_rms 0.1;
characteristicLength 0.3333333333333333;
numberOfFourrierModes 5;
U_c < 1.0, 0.0, 0.0 >;
}
TPMS TPMS
{ {
waveNumber 6.2831853072; waveNumber 6.2831853072;
......
apps/showcases/PorousMediaGPU/Setup.h
+ 20
− 3
View file @ 3282d7ea
...@@ -32,9 +32,16 @@ namespace walberla ...@@ -32,9 +32,16 @@ namespace walberla
struct Setup struct Setup
{ {
Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & tpms, Setup(std::shared_ptr<Config>& config, std::map<std::string, std::string>& infoMap)
const Config::BlockHandle & domainParameters, const Config::BlockHandle & logging, std::map<std::string, std::string>& infoMap)
{ {
auto parameters = config->getOneBlock("Parameters");
auto domainParameters = config->getOneBlock("DomainSetup");
auto logging = config->getOneBlock("Logging");
auto tpms = config->getOneBlock("TPMS");
auto turbulence = config->getOneBlock("Turbulence");
blockForestFilestem = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest"); blockForestFilestem = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest");
numProcesses = domainParameters.getParameter< uint_t >("numberProcesses"); numProcesses = domainParameters.getParameter< uint_t >("numberProcesses");
domainSize = domainParameters.getParameter< Vector3< real_t > >("domainSize"); domainSize = domainParameters.getParameter< Vector3< real_t > >("domainSize");
...@@ -52,7 +59,12 @@ struct Setup ...@@ -52,7 +59,12 @@ struct Setup
reynoldsNumber = parameters.getParameter< real_t >("reynoldsNumber"); reynoldsNumber = parameters.getParameter< real_t >("reynoldsNumber");
referenceVelocity = parameters.getParameter< real_t >("referenceVelocity"); referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
latticeVelocity = parameters.getParameter< real_t >("latticeVelocity"); latticeVelocity = parameters.getParameter< real_t >("latticeVelocity");
useGrid = parameters.getParameter< bool >("useGrid");
useGrid = turbulence.getParameter< bool >("useGrid");
urms = turbulence.getParameter< real_t >("U_rms");
characteristicLength = turbulence.getParameter< real_t >("characteristicLength");
numberOfFourrierModes = turbulence.getParameter< uint_t >("numberOfFourrierModes");
uc = turbulence.getParameter< Vector3< real_t > >("U_c");
waveNumber = tpms.getParameter< real_t >("waveNumber"); waveNumber = tpms.getParameter< real_t >("waveNumber");
scaling = tpms.getParameter< real_t >("scaling"); scaling = tpms.getParameter< real_t >("scaling");
...@@ -112,7 +124,12 @@ struct Setup ...@@ -112,7 +124,12 @@ struct Setup
real_t referenceVelocity; real_t referenceVelocity;
real_t referenceLength; real_t referenceLength;
real_t latticeVelocity; real_t latticeVelocity;
bool useGrid; bool useGrid;
real_t urms;
real_t characteristicLength;
uint_t numberOfFourrierModes;
Vector3<real_t> uc;
real_t waveNumber; real_t waveNumber;
real_t scaling; real_t scaling;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment