diff --git a/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp b/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
index 7d7c35b673396e082fd46e31be505a5837a55bac..90c6326aa36e061a00a28f778d9e6d6d4b75fb39 100644
--- a/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
+++ b/apps/showcases/PorousMediaGPU/PorousMediaGPU.cpp
@@ -29,6 +29,7 @@
#include "core/logging/Logging.h"
#include "core/math/Vector3.h"
#include "core/math/Constants.h"
+#include "core/math/Random.h"
#include "core/Environment.h"
#include "core/mpi/MPIManager.h"
#include "core/MemoryUsage.h"
@@ -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 ) );
}
+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){
Environment env( argc, argv );
mpi::MPIManager::instance()->useWorldComm();
@@ -123,12 +260,10 @@ int main(int argc, char** argv){
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
auto parameters = config->getOneBlock("Parameters");
- auto tpmsParameters = config->getOneBlock("TPMS");
- auto domainParameters = config->getOneBlock("DomainSetup");
auto loggingParameters = config->getOneBlock("Logging");
auto boundariesConfig = config->getBlock("Boundaries");
- Setup setup(parameters, tpmsParameters, domainParameters, loggingParameters, infoMap);
+ Setup setup(config, infoMap);
TPMS tpms(setup);
WALBERLA_LOG_INFO_ON_ROOT( "Simulation run data:"
"\n- simulation parameters:" << std::setprecision(16) <<
@@ -194,14 +329,7 @@ int main(int argc, char** argv){
}
sweepCollection.initialiseBlockPointer();
- for (auto& block : *blocks){
- 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);
-
+ Inflow<VelocityField_T> inflow(blocks, setup, ids);
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// LB SWEEPS AND BOUNDARY HANDLING ///
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -250,6 +378,8 @@ int main(int argc, char** argv){
auto defaultStream = gpu::StreamRAII::newPriorityStream(streamLowPriority);
SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
+ timeloop.addFuncBeforeTimeStep(inflow, "Update inflow velocity");
+
lbm_generated::BasicRecursiveTimeStepGPU< GPUPdfField_T, SweepCollection_T, BoundaryCollection_T >
LBMMeshRefinement(blocks, ids.pdfFieldGPU, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
diff --git a/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm b/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
index 60899bae570ffab5329b640d2d43db90d1880639..877f46b3db8d69208e127a2fbe5d37a5d3bf059c 100644
--- a/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
+++ b/apps/showcases/PorousMediaGPU/PorousMediaGPU.prm
@@ -9,8 +9,6 @@ Parameters
timesteps 11;
- useGrid false;
-
processMemoryLimit 512; // MiB
innerOuterSplit <1, 1, 1>;
@@ -18,6 +16,15 @@ Parameters
gpuEnabledMPI false;
}
+Turbulence
+{
+ useGrid false;
+ U_rms 0.1;
+ characteristicLength 0.3333333333333333;
+ numberOfFourrierModes 5;
+ U_c < 1.0, 0.0, 0.0 >;
+}
+
TPMS
{
waveNumber 6.2831853072;
diff --git a/apps/showcases/PorousMediaGPU/Setup.h b/apps/showcases/PorousMediaGPU/Setup.h
index ac994ad3f063ed6423e2f52a3400c68a8455a5e7..d45dfa2dc9d63fb0c1643c8d1006a0c5e143fdde 100644
--- a/apps/showcases/PorousMediaGPU/Setup.h
+++ b/apps/showcases/PorousMediaGPU/Setup.h
@@ -32,9 +32,16 @@ namespace walberla
struct Setup
{
- Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & tpms,
- const Config::BlockHandle & domainParameters, const Config::BlockHandle & logging, std::map<std::string, std::string>& infoMap)
+ Setup(std::shared_ptr<Config>& config, 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");
numProcesses = domainParameters.getParameter< uint_t >("numberProcesses");
domainSize = domainParameters.getParameter< Vector3< real_t > >("domainSize");
@@ -52,7 +59,12 @@ struct Setup
reynoldsNumber = parameters.getParameter< real_t >("reynoldsNumber");
referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
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");
scaling = tpms.getParameter< real_t >("scaling");
@@ -112,7 +124,12 @@ struct Setup
real_t referenceVelocity;
real_t referenceLength;
real_t latticeVelocity;
+
bool useGrid;
+ real_t urms;
+ real_t characteristicLength;
+ uint_t numberOfFourrierModes;
+ Vector3<real_t> uc;
real_t waveNumber;
real_t scaling;