From 681a99ab3466a5199b6f5a62541f0016a0a9c9a4 Mon Sep 17 00:00:00 2001
From: Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
Date: Wed, 1 Feb 2023 14:19:42 +0100
Subject: [PATCH] Fix incorrect treatment of body force density in FSLBM module
---
.../FreeSurface/BubblyPoiseuille.cpp | 14 +-
.../FreeSurface/BubblyPoiseuille.prm | 11 +-
apps/showcases/FreeSurface/CapillaryWave.cpp | 14 +-
apps/showcases/FreeSurface/CapillaryWave.prm | 3 +-
.../FreeSurface/DamBreakCylindrical.cpp | 14 +-
.../FreeSurface/DamBreakCylindrical.prm | 11 +-
.../FreeSurface/DamBreakRectangular.cpp | 14 +-
.../FreeSurface/DamBreakRectangular.prm | 11 +-
apps/showcases/FreeSurface/DropImpact.cpp | 14 +-
apps/showcases/FreeSurface/DropImpact.prm | 11 +-
apps/showcases/FreeSurface/DropWetting.cpp | 14 +-
apps/showcases/FreeSurface/DropWetting.prm | 11 +-
apps/showcases/FreeSurface/GravityWave.cpp | 69 ++++++--
apps/showcases/FreeSurface/GravityWave.prm | 9 +-
.../FreeSurface/GravityWaveCodegen.cpp | 43 ++---
.../GravityWaveLatticeModelGeneration.py | 49 +++---
apps/showcases/FreeSurface/MovingDrop.cpp | 14 +-
apps/showcases/FreeSurface/MovingDrop.prm | 11 +-
apps/showcases/FreeSurface/RisingBubble.cpp | 14 +-
apps/showcases/FreeSurface/RisingBubble.prm | 3 +-
apps/showcases/FreeSurface/TaylorBubble.cpp | 14 +-
apps/showcases/FreeSurface/TaylorBubble.prm | 3 +-
apps/showcases/FreeSurface/ZalesakDisk.cpp | 14 +-
apps/showcases/FreeSurface/ZalesakDisk.prm | 3 +-
.../communication/SimpleCommunication.h | 18 +-
src/lbm/free_surface/InitFunctions.h | 92 ++++++++++
src/lbm/free_surface/VtkWriter.h | 29 +++-
.../FreeSurfaceBoundaryHandling.impl.h | 24 +--
src/lbm/free_surface/dynamics/CMakeLists.txt | 2 +-
.../free_surface/dynamics/ForceDensitySweep.h | 159 ++++++++++++++++++
.../dynamics/ForceWeightingSweep.h | 96 -----------
.../dynamics/SurfaceDynamicsHandler.h | 50 ++++--
.../dynamics/functionality/AdvectMass.h | 4 +-
.../free_surface/surface_geometry/Utility.cpp | 6 +-
.../dynamics/CellConversionTest.cpp | 8 +-
.../lbm/free_surface/dynamics/CodegenTest.cpp | 81 +++++----
.../lbm/free_surface/dynamics/InflowTest.cpp | 8 +-
.../LatticeModelGenerationFreeSurface.py | 44 ++---
.../dynamics/WettingConversionTest.cpp | 10 +-
39 files changed, 631 insertions(+), 388 deletions(-)
create mode 100644 src/lbm/free_surface/dynamics/ForceDensitySweep.h
delete mode 100644 src/lbm/free_surface/dynamics/ForceWeightingSweep.h
diff --git a/apps/showcases/FreeSurface/BubblyPoiseuille.cpp b/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
index f3726c245..e27edf893 100644
--- a/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
+++ b/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
@@ -176,7 +176,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const bool enableBubbleModel = modelParameters.getParameter< bool >("enableBubbleModel");
const bool enableBubbleSplits = modelParameters.getParameter< bool >("enableBubbleSplits");
@@ -187,7 +186,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
@@ -205,11 +203,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel,
@@ -276,7 +274,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -316,9 +314,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -331,7 +329,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/BubblyPoiseuille.prm b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
index cfd43e279..e4fb01b3d 100644
--- a/apps/showcases/FreeSurface/BubblyPoiseuille.prm
+++ b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
@@ -30,7 +30,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
enableBubbleModel true;
enableBubbleSplits false; // only used if enableBubbleModel=true
@@ -75,16 +74,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/CapillaryWave.cpp b/apps/showcases/FreeSurface/CapillaryWave.cpp
index 2a6ebb277..64148c5b3 100644
--- a/apps/showcases/FreeSurface/CapillaryWave.cpp
+++ b/apps/showcases/FreeSurface/CapillaryWave.cpp
@@ -240,7 +240,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -251,7 +250,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -273,11 +271,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -346,7 +344,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -386,9 +384,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -408,7 +406,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/CapillaryWave.prm b/apps/showcases/FreeSurface/CapillaryWave.prm
index b73208639..a2d60bd59 100644
--- a/apps/showcases/FreeSurface/CapillaryWave.prm
+++ b/apps/showcases/FreeSurface/CapillaryWave.prm
@@ -28,7 +28,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -85,7 +84,7 @@ VTK
//obstacle_normal;
//pdf;
//flag;
- //force;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/DamBreakCylindrical.cpp b/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
index d082b5fcc..f48e23232 100644
--- a/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
+++ b/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
@@ -379,7 +379,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -391,7 +390,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -420,12 +418,12 @@ int main(int argc, char** argv)
const CollisionModel_T collisionModel = lbm::collision_model::SRTField< ScalarField_T >(relaxationRateFieldID);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
const LatticeModel_T latticeModel =
- LatticeModel_T(collisionModel, lbm::force_model::GuoField< VectorField_T >(forceFieldID));
+ LatticeModel_T(collisionModel, lbm::force_model::GuoField< VectorField_T >(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -455,7 +453,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -500,9 +498,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold, relaxationRateFieldID, smagorinskyConstant);
dynamicsHandler.addSweeps(timeloop);
@@ -531,7 +529,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/DamBreakCylindrical.prm b/apps/showcases/FreeSurface/DamBreakCylindrical.prm
index aee842982..ccf8a13ff 100644
--- a/apps/showcases/FreeSurface/DamBreakCylindrical.prm
+++ b/apps/showcases/FreeSurface/DamBreakCylindrical.prm
@@ -28,7 +28,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -78,16 +77,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/DamBreakRectangular.cpp b/apps/showcases/FreeSurface/DamBreakRectangular.cpp
index 6c7109b16..e7abd62e6 100644
--- a/apps/showcases/FreeSurface/DamBreakRectangular.cpp
+++ b/apps/showcases/FreeSurface/DamBreakRectangular.cpp
@@ -336,7 +336,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -348,7 +347,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -377,11 +375,11 @@ int main(int argc, char** argv)
const CollisionModel_T collisionModel = lbm::collision_model::SRTField< ScalarField_T >(relaxationRateFieldID);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -444,7 +442,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -489,9 +487,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold, relaxationRateFieldID, smagorinskyConstant);
dynamicsHandler.addSweeps(timeloop);
@@ -516,7 +514,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/DamBreakRectangular.prm b/apps/showcases/FreeSurface/DamBreakRectangular.prm
index 8c4bc96e7..8e424f919 100644
--- a/apps/showcases/FreeSurface/DamBreakRectangular.prm
+++ b/apps/showcases/FreeSurface/DamBreakRectangular.prm
@@ -28,7 +28,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -78,16 +77,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/DropImpact.cpp b/apps/showcases/FreeSurface/DropImpact.cpp
index e5df0d56d..896c5b098 100644
--- a/apps/showcases/FreeSurface/DropImpact.cpp
+++ b/apps/showcases/FreeSurface/DropImpact.cpp
@@ -183,7 +183,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const bool enableBubbleModel = modelParameters.getParameter< bool >("enableBubbleModel");
const bool enableBubbleSplits = modelParameters.getParameter< bool >("enableBubbleSplits");
@@ -194,7 +193,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
@@ -212,11 +210,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID =
@@ -269,7 +267,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -314,9 +312,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -329,7 +327,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/DropImpact.prm b/apps/showcases/FreeSurface/DropImpact.prm
index 76580d98f..8437db98e 100644
--- a/apps/showcases/FreeSurface/DropImpact.prm
+++ b/apps/showcases/FreeSurface/DropImpact.prm
@@ -32,7 +32,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
enableBubbleModel false;
enableBubbleSplits false; // only used if enableBubbleModel=true
@@ -77,16 +76,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/DropWetting.cpp b/apps/showcases/FreeSurface/DropWetting.cpp
index c2aca47f5..bdba59086 100644
--- a/apps/showcases/FreeSurface/DropWetting.cpp
+++ b/apps/showcases/FreeSurface/DropWetting.cpp
@@ -238,7 +238,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const bool enableBubbleModel = modelParameters.getParameter< bool >("enableBubbleModel");
const bool enableBubbleSplits = modelParameters.getParameter< bool >("enableBubbleSplits");
@@ -249,7 +248,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
@@ -271,11 +269,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -306,7 +304,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -351,9 +349,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -372,7 +370,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/DropWetting.prm b/apps/showcases/FreeSurface/DropWetting.prm
index e10ae9273..5e93daa21 100644
--- a/apps/showcases/FreeSurface/DropWetting.prm
+++ b/apps/showcases/FreeSurface/DropWetting.prm
@@ -28,7 +28,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
enableBubbleModel false;
enableBubbleSplits false; // only used if enableBubbleModel=true
@@ -75,16 +74,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/GravityWave.cpp b/apps/showcases/FreeSurface/GravityWave.cpp
index 10b52e266..acca70e45 100644
--- a/apps/showcases/FreeSurface/GravityWave.cpp
+++ b/apps/showcases/FreeSurface/GravityWave.cpp
@@ -308,8 +308,9 @@ int main(int argc, char** argv)
const real_t reynoldsNumber = physicsParameters.getParameter< real_t >("reynoldsNumber");
const real_t waveNumber = real_c(2) * math::pi / real_c(domainSize[0]);
const real_t waveFrequency = reynoldsNumber * viscosity / real_c(domainSize[0]) / initialAmplitude;
- const real_t forceY = -(waveFrequency * waveFrequency) / waveNumber / std::tanh(waveNumber * liquidDepth);
- const Vector3< real_t > force(real_c(0), forceY, real_c(0));
+ const real_t gravitationalAccelerationY =
+ -(waveFrequency * waveFrequency) / waveNumber / std::tanh(waveNumber * liquidDepth);
+ const Vector3< real_t > acceleration(real_c(0), gravitationalAccelerationY, real_c(0));
const bool enableWetting = physicsParameters.getParameter< bool >("enableWetting");
const real_t contactAngle = physicsParameters.getParameter< real_t >("contactAngle");
@@ -320,7 +321,7 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(contactAngle);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(acceleration);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
@@ -329,7 +330,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -340,7 +340,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -362,11 +361,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
- field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
+ const BlockDataID forceDensityFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Force density field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -435,8 +434,48 @@ int main(int argc, char** argv)
// might not detect solid flags correctly
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
+ // initialize sine profile such that there is exactly one period in the domain, i.e., with wavelength=domainSize[0];
+ // every length is normalized with domainSize[0]
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ ScalarField_T* const fillField = blockIt->getData< ScalarField_T >(fillFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(fillFieldIt, fillField, {
+ // cell in block-local coordinates
+ const Cell localCell = fillFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = fillFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // Monte-Carlo like estimation of the fill level:
+ // create uniformly-distributed sample points in each cell and count the number of points below the sine
+ // profile; this fraction of points is used as the fill level to initialize the profile
+ uint_t numPointsBelow = uint_c(0);
+
+ for (uint_t xSample = uint_c(0); xSample <= fillLevelInitSamples; ++xSample)
+ {
+ // value of the sine-function
+ const real_t functionValue = initializationProfile(real_c(globalCell[0]) + real_c(xSample) * stepsize,
+ initialAmplitude, liquidDepth, real_c(domainSize[0]));
+
+ for (uint_t ySample = uint_c(0); ySample <= fillLevelInitSamples; ++ySample)
+ {
+ const real_t yPoint = real_c(globalCell[1]) + real_c(ySample) * stepsize;
+ // with operator <, a fill level of 1 can not be reached when the line is equal to the cell's top border;
+ // with operator <=, a fill level of 0 can not be reached when the line is equal to the cell's bottom
+ // border
+ if (yPoint < functionValue) { ++numPointsBelow; }
+ }
+ }
+
+ // fill level is fraction of points below sine profile
+ fillField->get(localCell) = real_c(numPointsBelow) / real_c(numTotalPoints);
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -457,7 +496,11 @@ int main(int argc, char** argv)
else { bubbleModel = std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1)); }
// initialize hydrostatic pressure
- initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, force, liquidDepth);
+ initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, acceleration, liquidDepth);
+
+ // initialize force density field
+ initForceDensityField< PdfField_T, FlagField_T, VectorField_T, ScalarField_T >(
+ blockForest, forceDensityFieldID, fillFieldID, pdfFieldID, flagFieldID, flagInfo, acceleration);
// set density in non-liquid or non-interface cells to one (after initializing with hydrostatic pressure)
setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
@@ -485,9 +528,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, acceleration, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -513,7 +556,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/GravityWave.prm b/apps/showcases/FreeSurface/GravityWave.prm
index 0400d62da..ffc5a4c49 100644
--- a/apps/showcases/FreeSurface/GravityWave.prm
+++ b/apps/showcases/FreeSurface/GravityWave.prm
@@ -8,9 +8,9 @@ BlockForestParameters
DomainParameters
{
- domainWidth 200; // equivalent to wavelength
- liquidDepth 100;
- initialAmplitude 2;
+ domainWidth 50; // equivalent to wavelength
+ liquidDepth 25;
+ initialAmplitude 0.5;
}
PhysicsParameters
@@ -28,7 +28,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -85,7 +84,7 @@ VTK
//obstacle_normal;
//pdf;
//flag;
- //force;
+ force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/GravityWaveCodegen.cpp b/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
index 1658f9290..dd603d62c 100644
--- a/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
+++ b/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
@@ -46,8 +46,9 @@ namespace free_surface
{
namespace GravityWaveCodegen
{
-using ScalarField_T = GhostLayerField< real_t, 1 >;
-using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+using VectorFieldFlattened_T = GhostLayerField< real_t, 3 >;
using LatticeModel_T = lbm::GravityWaveLatticeModel;
using LatticeModelStencil_T = LatticeModel_T::Stencil;
@@ -309,8 +310,8 @@ int main(int argc, char** argv)
const real_t reynoldsNumber = physicsParameters.getParameter< real_t >("reynoldsNumber");
const real_t waveNumber = real_c(2) * math::pi / real_c(domainSize[0]);
const real_t waveFrequency = reynoldsNumber * viscosity / real_c(domainSize[0]) / initialAmplitude;
- const real_t forceY = -(waveFrequency * waveFrequency) / waveNumber / std::tanh(waveNumber * liquidDepth);
- const Vector3< real_t > force(real_c(0), forceY, real_c(0));
+ const real_t accelerationY = -(waveFrequency * waveFrequency) / waveNumber / std::tanh(waveNumber * liquidDepth);
+ const Vector3< real_t > acceleration(real_c(0), accelerationY, real_c(0));
const bool enableWetting = physicsParameters.getParameter< bool >("enableWetting");
const real_t contactAngle = physicsParameters.getParameter< real_t >("contactAngle");
@@ -321,7 +322,7 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(contactAngle);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(acceleration);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
@@ -330,7 +331,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -341,7 +341,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -363,11 +362,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
- field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
+ const BlockDataID forceDensityFieldID =
+ field::addToStorage< VectorFieldFlattened_T >(blockForest, "Force field", real_c(0), field::fzyx, uint_c(1));
// create lattice model
- LatticeModel_T latticeModel = LatticeModel_T(force[0], force[1], force[2], relaxationRate);
+ LatticeModel_T latticeModel = LatticeModel_T(forceDensityFieldID, relaxationRate);
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -437,7 +436,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -458,7 +457,11 @@ int main(int argc, char** argv)
else { bubbleModel = std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1)); }
// initialize hydrostatic pressure
- initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, force, liquidDepth);
+ initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, acceleration, liquidDepth);
+
+ // initialize force density field
+ initForceDensityFieldCodegen< PdfField_T, FlagField_T, VectorFieldFlattened_T, ScalarField_T >(
+ blockForest, forceDensityFieldID, fillFieldID, pdfFieldID, flagFieldID, flagInfo, acceleration);
// set density in non-liquid or non-interface cells to 1 (after initializing with hydrostatic pressure)
setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
@@ -485,11 +488,12 @@ int main(int argc, char** argv)
const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
// add boundary handling for standard boundaries and free surface boundaries
- const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T, true > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
- freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
- cellConversionForceThreshold);
+ const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T, true,
+ VectorFieldFlattened_T >
+ dynamicsHandler(blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID,
+ curvatureFieldID, freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel,
+ pdfRefillingModel, excessMassDistributionModel, relaxationRate, acceleration, surfaceTension,
+ useSimpleMassExchange, cellConversionThreshold, cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -513,8 +517,9 @@ int main(int argc, char** argv)
timeloop.addFuncAfterTimeStep(symmetryEvaluator, "Evaluator: symmetry norm");
// add VTK output
- addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T,
+ true, VectorFieldFlattened_T >(
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py b/apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py
index 22f474fcf..aa61ee797 100644
--- a/apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py
+++ b/apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py
@@ -1,5 +1,5 @@
import sympy as sp
-
+import pystencils as ps
from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
from lbmpy.enums import ForceModel, Method, Stencil
from lbmpy.stencils import LBStencil
@@ -7,28 +7,31 @@ from lbmpy.stencils import LBStencil
from pystencils_walberla import CodeGeneration
from lbmpy_walberla import generate_lattice_model
-# general parameters
-stencil = LBStencil(Stencil.D3Q19)
-omega = sp.Symbol('omega')
-force = sp.symbols('force_:3')
-layout = 'fzyx'
-
-# method definition
-lbm_config = LBMConfig(stencil=stencil,
- method=Method.SRT,
- relaxation_rate=omega,
- compressible=True,
- force=force,
- force_model=ForceModel.GUO,
- zero_centered=False,
- streaming_pattern='pull') # free surface implementation only works with pull pattern
-
-# optimizations to be used by the code generator
-lbm_opt = LBMOptimisation(cse_global=True,
- field_layout=layout)
-
-collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
- lbm_optimisation=lbm_opt)
with CodeGeneration() as ctx:
+ # general parameters
+ layout = 'fzyx'
+ data_type = "float64" if ctx.double_accuracy else "float32"
+
+ stencil = LBStencil(Stencil.D3Q19)
+ omega = sp.Symbol('omega')
+ force_field = ps.fields(f"force(3): {data_type}[3D]", layout='fzyx')
+
+ # method definition
+ lbm_config = LBMConfig(stencil=stencil,
+ method=Method.SRT,
+ relaxation_rate=omega,
+ compressible=True,
+ force=force_field,
+ force_model=ForceModel.GUO,
+ zero_centered=False,
+ streaming_pattern='pull') # free surface implementation only works with pull pattern
+
+ # optimizations to be used by the code generator
+ lbm_opt = LBMOptimisation(cse_global=True,
+ field_layout=layout)
+
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
+ lbm_optimisation=lbm_opt)
+
generate_lattice_model(ctx, "GravityWaveLatticeModel", collision_rule, field_layout=layout)
diff --git a/apps/showcases/FreeSurface/MovingDrop.cpp b/apps/showcases/FreeSurface/MovingDrop.cpp
index 9ddba4ae8..352a85c72 100644
--- a/apps/showcases/FreeSurface/MovingDrop.cpp
+++ b/apps/showcases/FreeSurface/MovingDrop.cpp
@@ -170,7 +170,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const bool enableBubbleModel = modelParameters.getParameter< bool >("enableBubbleModel");
const bool enableBubbleSplits = modelParameters.getParameter< bool >("enableBubbleSplits");
@@ -181,7 +180,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
@@ -199,11 +197,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID =
@@ -233,7 +231,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -275,9 +273,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -290,7 +288,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/MovingDrop.prm b/apps/showcases/FreeSurface/MovingDrop.prm
index dff95569f..dceed84e7 100644
--- a/apps/showcases/FreeSurface/MovingDrop.prm
+++ b/apps/showcases/FreeSurface/MovingDrop.prm
@@ -30,7 +30,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
enableBubbleModel false;
enableBubbleSplits false; // only used if enableBubbleModel=true
@@ -75,16 +74,16 @@ VTK
writers
{
+ fill_level;
+ mapped_flag;
velocity;
density;
- //pdf;
- flag;
- fill_level;
- force;
curvature;
normal;
obstacle_normal;
- mapped_flag;
+ //pdf;
+ //flag;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/RisingBubble.cpp b/apps/showcases/FreeSurface/RisingBubble.cpp
index ef48f46cc..09bef1b3f 100644
--- a/apps/showcases/FreeSurface/RisingBubble.cpp
+++ b/apps/showcases/FreeSurface/RisingBubble.cpp
@@ -247,7 +247,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -258,7 +257,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -280,11 +278,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add various fields
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -315,7 +313,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -358,9 +356,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -386,7 +384,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/RisingBubble.prm b/apps/showcases/FreeSurface/RisingBubble.prm
index 62f1053d2..1acf65000 100644
--- a/apps/showcases/FreeSurface/RisingBubble.prm
+++ b/apps/showcases/FreeSurface/RisingBubble.prm
@@ -29,7 +29,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -86,7 +85,7 @@ VTK
//obstacle_normal;
//pdf;
//flag;
- //force;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/TaylorBubble.cpp b/apps/showcases/FreeSurface/TaylorBubble.cpp
index 57436886d..9672337ed 100644
--- a/apps/showcases/FreeSurface/TaylorBubble.cpp
+++ b/apps/showcases/FreeSurface/TaylorBubble.cpp
@@ -253,7 +253,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -264,7 +263,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -286,11 +284,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add various fields
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -345,7 +343,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -388,9 +386,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -414,7 +412,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/TaylorBubble.prm b/apps/showcases/FreeSurface/TaylorBubble.prm
index e413afc7f..c4096891f 100644
--- a/apps/showcases/FreeSurface/TaylorBubble.prm
+++ b/apps/showcases/FreeSurface/TaylorBubble.prm
@@ -31,7 +31,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -88,7 +87,7 @@ VTK
obstacle_normal;
//pdf;
//flag;
- //force;
+ //force_density;
}
inclusion_filters
diff --git a/apps/showcases/FreeSurface/ZalesakDisk.cpp b/apps/showcases/FreeSurface/ZalesakDisk.cpp
index 02beb5706..edbaab294 100644
--- a/apps/showcases/FreeSurface/ZalesakDisk.cpp
+++ b/apps/showcases/FreeSurface/ZalesakDisk.cpp
@@ -163,7 +163,6 @@ int main(int argc, char** argv)
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
- const bool enableForceWeighting = modelParameters.getParameter< bool >("enableForceWeighting");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
@@ -174,7 +173,6 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
@@ -192,11 +190,11 @@ int main(int argc, char** argv)
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
// add force field
- const BlockDataID forceFieldID =
+ const BlockDataID forceDensityFieldID =
field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
@@ -254,7 +252,7 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
@@ -298,9 +296,9 @@ int main(int argc, char** argv)
// add boundary handling for standard boundaries and free surface boundaries
const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+ relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -313,7 +311,7 @@ int main(int argc, char** argv)
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
diff --git a/apps/showcases/FreeSurface/ZalesakDisk.prm b/apps/showcases/FreeSurface/ZalesakDisk.prm
index 259774efd..b9af1280c 100644
--- a/apps/showcases/FreeSurface/ZalesakDisk.prm
+++ b/apps/showcases/FreeSurface/ZalesakDisk.prm
@@ -26,7 +26,6 @@ ModelParameters
pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
- enableForceWeighting false;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
@@ -81,7 +80,7 @@ VTK
//obstacle_normal;
//pdf;
//flag;
- //force;
+ //force_density;
}
inclusion_filters
diff --git a/src/lbm/blockforest/communication/SimpleCommunication.h b/src/lbm/blockforest/communication/SimpleCommunication.h
index 42bffe263..764152da6 100644
--- a/src/lbm/blockforest/communication/SimpleCommunication.h
+++ b/src/lbm/blockforest/communication/SimpleCommunication.h
@@ -42,10 +42,11 @@ using communication::UniformBufferedScheme;
template< typename Stencil_T >
class SimpleCommunication : public communication::UniformBufferedScheme< Stencil_T >
{
- using RealScalarField_T = GhostLayerField< real_t, 1 >;
- using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
- using PdfField_T = GhostLayerField< real_t, Stencil_T::Size >;
- using UintScalarField_T = GhostLayerField< uint_t, 1 >;
+ using RealScalarField_T = GhostLayerField< real_t, 1 >;
+ using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+ using VectorFieldFlattened_T = GhostLayerField< real_t, 3 >;
+ using PdfField_T = GhostLayerField< real_t, Stencil_T::Size >;
+ using UintScalarField_T = GhostLayerField< uint_t, 1 >;
using FlagField16_T = FlagField< uint16_t >;
using FlagField32_T = FlagField< uint32_t >;
@@ -151,7 +152,14 @@ class SimpleCommunication : public communication::UniformBufferedScheme< Stencil
{
this->addPackInfo(make_shared< PackInfo< UintScalarField_T > >(fieldId));
}
- else { WALBERLA_ABORT("Problem with UID"); }
+ else
+ {
+ if (firstBlock.isDataClassOrSubclassOf< VectorFieldFlattened_T >(fieldId))
+ {
+ this->addPackInfo(make_shared< PackInfo< VectorFieldFlattened_T > >(fieldId));
+ }
+ else { WALBERLA_ABORT("Problem with UID"); }
+ }
}
}
}
diff --git a/src/lbm/free_surface/InitFunctions.h b/src/lbm/free_surface/InitFunctions.h
index b98f55d54..22cd09c62 100644
--- a/src/lbm/free_surface/InitFunctions.h
+++ b/src/lbm/free_surface/InitFunctions.h
@@ -200,6 +200,98 @@ void initHydrostaticPressure(const std::weak_ptr< StructuredBlockForest >& block
}
}
+/***********************************************************************************************************************
+ * Initialize the force density field with a given acceleration and density of each cell.
+ **********************************************************************************************************************/
+template< typename PdfField_T, typename FlagField_T, typename VectorField_T, typename ScalarField_T >
+void initForceDensityField(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+ const BlockDataID& forceDensityFieldID, const ConstBlockDataID& fillFieldID,
+ const ConstBlockDataID& pdfFieldID, const ConstBlockDataID& flagFieldID,
+ const FlagInfo< FlagField_T >& flagInfo, const Vector3< real_t >& acceleration)
+{
+ const auto blockForest = blockForestPtr.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ VectorField_T* const forceDensityField = blockIt->getData< VectorField_T >(forceDensityFieldID);
+ const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID);
+ const PdfField_T* const pdfField = blockIt->getData< const PdfField_T >(pdfFieldID);
+ const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(forceDensityFieldIt, forceDensityField, fillFieldIt, fillField, pdfFieldIt, pdfField,
+ flagFieldIt, flagField, {
+ // set force density in cells to acceleration * density * fillLevel (see equation 15
+ // in Koerner et al., 2005);
+
+ *forceDensityFieldIt = Vector3< real_t >(real_c(0));
+
+ if (flagInfo.isInterface(*flagFieldIt))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ *forceDensityFieldIt = acceleration * *fillFieldIt * density;
+ }
+
+ else
+ {
+ if (flagInfo.isLiquid(*flagFieldIt))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ *forceDensityFieldIt = acceleration * density;
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+}
+
+/***********************************************************************************************************************
+ * Initialize the force density field with a given acceleration and density of each cell.
+ * Differs from the version above by using a flattened vector field (GhostLayerField< real_t, 3 >). This is necessary
+ * because Pystencils does not support VectorField_T (GhostLayerField< Vector3<real_t>, 1 >).
+ **********************************************************************************************************************/
+template< typename PdfField_T, typename FlagField_T, typename VectorFieldFlattened_T, typename ScalarField_T >
+void initForceDensityFieldCodegen(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+ const BlockDataID& forceDensityFieldID, const ConstBlockDataID& fillFieldID,
+ const ConstBlockDataID& pdfFieldID, const ConstBlockDataID& flagFieldID,
+ const FlagInfo< FlagField_T >& flagInfo, const Vector3< real_t >& acceleration)
+{
+ const auto blockForest = blockForestPtr.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ VectorFieldFlattened_T* const forceDensityField = blockIt->getData< VectorFieldFlattened_T >(forceDensityFieldID);
+ const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID);
+ const PdfField_T* const pdfField = blockIt->getData< const PdfField_T >(pdfFieldID);
+ const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(forceDensityFieldIt, forceDensityField, fillFieldIt, fillField, pdfFieldIt, pdfField,
+ flagFieldIt, flagField, {
+ // set force density in cells to acceleration * density * fillLevel (see equation 15
+ // in Koerner et al., 2005);
+
+ if (flagInfo.isInterface(*flagFieldIt))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ forceDensityFieldIt[0] = acceleration[0] * *fillFieldIt * density;
+ forceDensityFieldIt[1] = acceleration[1] * *fillFieldIt * density;
+ forceDensityFieldIt[2] = acceleration[2] * *fillFieldIt * density;
+ }
+
+ else
+ {
+ if (flagInfo.isLiquid(*flagFieldIt))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ forceDensityFieldIt[0] = acceleration[0] * density;
+ forceDensityFieldIt[1] = acceleration[1] * density;
+ forceDensityFieldIt[2] = acceleration[2] * density;
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+}
+
/***********************************************************************************************************************
* Set density in non-liquid and non-interface cells to 1.
**********************************************************************************************************************/
diff --git a/src/lbm/free_surface/VtkWriter.h b/src/lbm/free_surface/VtkWriter.h
index f98549a8f..0683479e1 100644
--- a/src/lbm/free_surface/VtkWriter.h
+++ b/src/lbm/free_surface/VtkWriter.h
@@ -47,13 +47,14 @@ namespace free_surface
* config-file.
**********************************************************************************************************************/
template< typename LatticeModel_T, typename FreeSurfaceBoundaryHandling_T, typename PdfField_T, typename FlagField_T,
- typename ScalarField_T, typename VectorField_T >
+ typename ScalarField_T, typename VectorField_T, bool useCodegen = false,
+ typename VectorFieldFlattened_T = GhostLayerField< real_t, 3 > >
void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr, SweepTimeloop& timeloop,
const std::weak_ptr< Config >& configPtr,
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo, const BlockDataID& pdfFieldID,
- const BlockDataID& flagFieldID, const BlockDataID& fillFieldID, const BlockDataID& forceFieldID,
- const BlockDataID& curvatureFieldID, const BlockDataID& normalFieldID,
- const BlockDataID& obstacleNormalFieldID)
+ const BlockDataID& flagFieldID, const BlockDataID& fillFieldID,
+ const BlockDataID& forceDensityFieldID, const BlockDataID& curvatureFieldID,
+ const BlockDataID& normalFieldID, const BlockDataID& obstacleNormalFieldID)
{
const auto blockForest = blockForestPtr.lock();
WALBERLA_CHECK_NOT_NULLPTR(blockForest);
@@ -85,7 +86,15 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(normalFieldID, "normal"));
writers.push_back(
std::make_shared< VTKWriter< VectorField_T, float > >(obstacleNormalFieldID, "obstacle_normal"));
- writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(forceFieldID, "force"));
+ if constexpr (useCodegen)
+ {
+ writers.push_back(
+ std::make_shared< VTKWriter< VectorFieldFlattened_T, float > >(forceDensityFieldID, "force_density"));
+ }
+ else
+ {
+ writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(forceDensityFieldID, "force_density"));
+ }
// map flagIDs to integer values
const auto flagMapper =
@@ -118,6 +127,16 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< VectorField_T > >(normalFieldID));
preVTKComm.addPackInfo(
std::make_shared< field::communication::PackInfo< VectorField_T > >(obstacleNormalFieldID));
+ if constexpr (useCodegen)
+ {
+ preVTKComm.addPackInfo(
+ std::make_shared< field::communication::PackInfo< VectorFieldFlattened_T > >(forceDensityFieldID));
+ }
+ else
+ {
+ preVTKComm.addPackInfo(
+ std::make_shared< field::communication::PackInfo< VectorField_T > >(forceDensityFieldID));
+ }
beforeFuncs["ghost_layer_synchronization"] = preVTKComm;
diff --git a/src/lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.impl.h b/src/lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.impl.h
index 254293cff..1fdb69f98 100644
--- a/src/lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.impl.h
+++ b/src/lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.impl.h
@@ -180,12 +180,12 @@ FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::FreeS
// define IDs (static const variables)
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const field::FlagUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::noSlipFlagID = field::FlagUID("NoSlip");
+const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::noSlipFlagID =
+ field::FlagUID("NoSlip");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const field::FlagUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbFlagID = field::FlagUID("UBB");
+const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbFlagID =
+ field::FlagUID("UBB");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbInflowFlagID =
@@ -200,20 +200,20 @@ const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, S
field::FlagUID("PressureOutflow");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const field::FlagUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::outletFlagID = field::FlagUID("Outlet");
+const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::outletFlagID =
+ field::FlagUID("Outlet");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
const field::FlagUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::freeSlipFlagID =
field::FlagUID("FreeSlip");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const BoundaryUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::noSlipBoundaryID = BoundaryUID("NoSlip");
+const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::noSlipBoundaryID =
+ BoundaryUID("NoSlip");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const BoundaryUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbBoundaryID = BoundaryUID("UBB");
+const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbBoundaryID =
+ BoundaryUID("UBB");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::ubbInflowBoundaryID =
@@ -228,8 +228,8 @@ const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, Scal
BoundaryUID("PressureOutflow");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
-const BoundaryUID
- FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::outletBoundaryID = BoundaryUID("Outlet");
+const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::outletBoundaryID =
+ BoundaryUID("Outlet");
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T >
const BoundaryUID FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >::freeSlipBoundaryID =
diff --git a/src/lbm/free_surface/dynamics/CMakeLists.txt b/src/lbm/free_surface/dynamics/CMakeLists.txt
index 458fabdfe..2d7145db3 100644
--- a/src/lbm/free_surface/dynamics/CMakeLists.txt
+++ b/src/lbm/free_surface/dynamics/CMakeLists.txt
@@ -5,7 +5,7 @@ target_sources( lbm
ExcessMassDistributionModel.h
ExcessMassDistributionSweep.h
ExcessMassDistributionSweep.impl.h
- ForceWeightingSweep.h
+ ForceDensitySweep.h
PdfReconstructionModel.h
PdfRefillingModel.h
PdfRefillingSweep.h
diff --git a/src/lbm/free_surface/dynamics/ForceDensitySweep.h b/src/lbm/free_surface/dynamics/ForceDensitySweep.h
new file mode 100644
index 000000000..afe418e9c
--- /dev/null
+++ b/src/lbm/free_surface/dynamics/ForceDensitySweep.h
@@ -0,0 +1,159 @@
+//======================================================================================================================
+//
+// This file is part of waLBerla. waLBerla is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// waLBerla is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file ForceDensitySweep.h
+//! \ingroup dynamics
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Weight force in interface cells with fill level and density (equation 15 in Koerner et al., 2005).
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/logging/Logging.h"
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include "field/FlagField.h"
+#include "field/GhostLayerField.h"
+
+#include "lbm/field/PdfField.h"
+#include "lbm/free_surface/FlagInfo.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+/***********************************************************************************************************************
+ * Update the force density field as in equation 15 in Koerner et al., 2005. with "acceleration * density * fill level".
+ **********************************************************************************************************************/
+template< typename LatticeModel_T, typename FlagField_T, typename VectorField_T, typename ScalarField_T >
+class ForceDensitySweep
+{
+ public:
+ ForceDensitySweep(BlockDataID forceDensityFieldID, ConstBlockDataID pdfFieldID, ConstBlockDataID flagFieldID,
+ ConstBlockDataID fillFieldID, const FlagInfo< FlagField_T >& flagInfo,
+ const Vector3< real_t >& globalAcceleration)
+ : forceDensityFieldID_(forceDensityFieldID), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID),
+ fillFieldID_(fillFieldID), flagInfo_(flagInfo), globalAcceleration_(globalAcceleration)
+ {}
+
+ void operator()(IBlock* const block)
+ {
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+ VectorField_T* const forceDensityField = block->getData< VectorField_T >(forceDensityFieldID_);
+ const PdfField_T* const pdfField = block->getData< const PdfField_T >(pdfFieldID_);
+ const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID_);
+ const ScalarField_T* const fillField = block->getData< const ScalarField_T >(fillFieldID_);
+
+ WALBERLA_FOR_ALL_CELLS(forceDensityFieldIt, forceDensityField, pdfFieldIt, pdfField, flagFieldIt, flagField,
+ fillFieldIt, fillField, {
+ flag_t flag = *flagFieldIt;
+
+ // set force density in cells to acceleration * density * fillLevel (see equation 15
+ // in Koerner et al., 2005);
+ if (flagInfo_.isInterface(flag))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ *forceDensityFieldIt = globalAcceleration_ * *fillFieldIt * density;
+ }
+
+ else
+ {
+ if (flagInfo_.isLiquid(flag))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ *forceDensityFieldIt = globalAcceleration_ * density;
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ private:
+ using flag_t = typename FlagField_T::flag_t;
+
+ BlockDataID forceDensityFieldID_;
+ ConstBlockDataID pdfFieldID_;
+ ConstBlockDataID flagFieldID_;
+ ConstBlockDataID fillFieldID_;
+ FlagInfo< FlagField_T > flagInfo_;
+ Vector3< real_t > globalAcceleration_;
+}; // class ForceDensitySweep
+
+/***********************************************************************************************************************
+ * Update the force density field as in equation 15 in Koerner et al., 2005. with "acceleration * density * fill level".
+ * Differs from the version above by using a flattened vector field (GhostLayerField< real_t, 3 >). This is necessary
+ * because Pystencils does not support VectorField_T (GhostLayerField< Vector3<real_t>, 1 >).
+ **********************************************************************************************************************/
+template< typename LatticeModel_T, typename FlagField_T, typename VectorFieldFlattened_T, typename ScalarField_T >
+class ForceDensityCodegenSweep
+{
+ public:
+ ForceDensityCodegenSweep(BlockDataID forceDensityFieldID, ConstBlockDataID pdfFieldID, ConstBlockDataID flagFieldID,
+ ConstBlockDataID fillFieldID, const FlagInfo< FlagField_T >& flagInfo,
+ const Vector3< real_t >& globalAcceleration)
+ : forceDensityFieldID_(forceDensityFieldID), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID),
+ fillFieldID_(fillFieldID), flagInfo_(flagInfo), globalAcceleration_(globalAcceleration)
+ {}
+
+ void operator()(IBlock* const block)
+ {
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+ VectorFieldFlattened_T* const forceDensityField = block->getData< VectorFieldFlattened_T >(forceDensityFieldID_);
+ const PdfField_T* const pdfField = block->getData< const PdfField_T >(pdfFieldID_);
+ const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID_);
+ const ScalarField_T* const fillField = block->getData< const ScalarField_T >(fillFieldID_);
+
+ WALBERLA_FOR_ALL_CELLS(forceDensityFieldIt, forceDensityField, pdfFieldIt, pdfField, flagFieldIt, flagField,
+ fillFieldIt, fillField, {
+ flag_t flag = *flagFieldIt;
+
+ // set force density in cells to acceleration * density * fillLevel (see equation 15
+ // in Koerner et al., 2005);
+ if (flagInfo_.isInterface(flag))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ forceDensityFieldIt[0] = globalAcceleration_[0] * *fillFieldIt * density;
+ forceDensityFieldIt[1] = globalAcceleration_[1] * *fillFieldIt * density;
+ forceDensityFieldIt[2] = globalAcceleration_[2] * *fillFieldIt * density;
+ }
+ else
+ {
+ if (flagInfo_.isLiquid(flag))
+ {
+ const real_t density = pdfField->getDensity(pdfFieldIt.cell());
+ forceDensityFieldIt[0] = globalAcceleration_[0] * density;
+ forceDensityFieldIt[1] = globalAcceleration_[1] * density;
+ forceDensityFieldIt[2] = globalAcceleration_[2] * density;
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ private:
+ using flag_t = typename FlagField_T::flag_t;
+
+ BlockDataID forceDensityFieldID_;
+ ConstBlockDataID pdfFieldID_;
+ ConstBlockDataID flagFieldID_;
+ ConstBlockDataID fillFieldID_;
+ FlagInfo< FlagField_T > flagInfo_;
+ Vector3< real_t > globalAcceleration_;
+}; // class ForceDensitySweep
+
+} // namespace free_surface
+} // namespace walberla
diff --git a/src/lbm/free_surface/dynamics/ForceWeightingSweep.h b/src/lbm/free_surface/dynamics/ForceWeightingSweep.h
deleted file mode 100644
index 4a13fb3b4..000000000
--- a/src/lbm/free_surface/dynamics/ForceWeightingSweep.h
+++ /dev/null
@@ -1,96 +0,0 @@
-//======================================================================================================================
-//
-// This file is part of waLBerla. waLBerla is free software: you can
-// redistribute it and/or modify it under the terms of the GNU General Public
-// License as published by the Free Software Foundation, either version 3 of
-// the License, or (at your option) any later version.
-//
-// waLBerla is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
-// for more details.
-//
-// You should have received a copy of the GNU General Public License along
-// with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
-//
-//! \file ForceWeightingSweep.h
-//! \ingroup dynamics
-//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
-//! \brief Weight force in interface cells with fill level and density (equation 15 in Koerner et al., 2005).
-//
-//======================================================================================================================
-
-#pragma once
-
-#include "core/logging/Logging.h"
-
-#include "domain_decomposition/BlockDataID.h"
-
-#include "field/FlagField.h"
-#include "field/GhostLayerField.h"
-
-#include "lbm/field/PdfField.h"
-#include "lbm/free_surface/FlagInfo.h"
-
-namespace walberla
-{
-namespace free_surface
-{
-/***********************************************************************************************************************
- * Weight the specified force in interface cells according to their density and fill level, as in equation 15 in Koerner
- * et al., 2005.
- * In liquid cells, the force is set to the specified constant global force.
- **********************************************************************************************************************/
-template< typename LatticeModel_T, typename FlagField_T, typename VectorField_T, typename ScalarField_T >
-class ForceWeightingSweep
-{
- public:
- ForceWeightingSweep(BlockDataID forceFieldID, ConstBlockDataID pdfFieldID, ConstBlockDataID flagFieldID,
- ConstBlockDataID fillFieldID, const FlagInfo< FlagField_T >& flagInfo,
- const Vector3< real_t >& globalForce)
- : forceFieldID_(forceFieldID), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), fillFieldID_(fillFieldID),
- flagInfo_(flagInfo), globalForce_(globalForce)
- {}
-
- void operator()(IBlock* const block)
- {
- using PdfField_T = lbm::PdfField< LatticeModel_T >;
-
- VectorField_T* const forceField = block->getData< VectorField_T >(forceFieldID_);
- const PdfField_T* const pdfField = block->getData< const PdfField_T >(pdfFieldID_);
- const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID_);
- const ScalarField_T* const fillField = block->getData< const ScalarField_T >(fillFieldID_);
-
- WALBERLA_FOR_ALL_CELLS(forceFieldIt, forceField, pdfFieldIt, pdfField, flagFieldIt, flagField, fillFieldIt,
- fillField, {
- flag_t flag = *flagFieldIt;
-
- // set force in interface cells to globalForce_ * density * fillLevel (see equation 15
- // in Koerner et al., 2005)
- if (flagInfo_.isInterface(flag))
- {
- const real_t density = pdfField->getDensity(pdfFieldIt.cell());
- const real_t fillLevel = *fillFieldIt;
- *forceFieldIt = globalForce_ * fillLevel * density;
- }
- else
- {
- // set force to globalForce_ in all non-interface cells
- *forceFieldIt = globalForce_;
- }
- }) // WALBERLA_FOR_ALL_CELLS
- }
-
- private:
- using flag_t = typename FlagField_T::flag_t;
-
- BlockDataID forceFieldID_;
- ConstBlockDataID pdfFieldID_;
- ConstBlockDataID flagFieldID_;
- ConstBlockDataID fillFieldID_;
- FlagInfo< FlagField_T > flagInfo_;
- Vector3< real_t > globalForce_;
-}; // class ForceWeightingSweep
-
-} // namespace free_surface
-} // namespace walberla
diff --git a/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h b/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
index 93f204d30..d65469940 100644
--- a/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
+++ b/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
@@ -48,7 +48,7 @@
#include "ConversionFlagsResetSweep.h"
#include "ExcessMassDistributionModel.h"
#include "ExcessMassDistributionSweep.h"
-#include "ForceWeightingSweep.h"
+#include "ForceDensitySweep.h"
#include "PdfReconstructionModel.h"
#include "PdfRefillingModel.h"
#include "PdfRefillingSweep.h"
@@ -59,7 +59,7 @@ namespace walberla
namespace free_surface
{
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename VectorField_T,
- bool useCodegen = false >
+ bool useCodegen = false, typename VectorFieldFlattened_T = GhostLayerField< real_t, 3 > >
class SurfaceDynamicsHandler
{
protected:
@@ -74,22 +74,22 @@ class SurfaceDynamicsHandler
public:
SurfaceDynamicsHandler(const std::shared_ptr< StructuredBlockForest >& blockForest, BlockDataID pdfFieldID,
- BlockDataID flagFieldID, BlockDataID fillFieldID, BlockDataID forceFieldID,
+ BlockDataID flagFieldID, BlockDataID fillFieldID, BlockDataID forceDensityFieldID,
ConstBlockDataID normalFieldID, ConstBlockDataID curvatureFieldID,
const std::shared_ptr< FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
const std::shared_ptr< BubbleModelBase >& bubbleModel,
const std::string& pdfReconstructionModel, const std::string& pdfRefillingModel,
const std::string& excessMassDistributionModel, real_t relaxationRate,
- const Vector3< real_t >& globalForce, real_t surfaceTension, bool enableForceWeighting,
+ const Vector3< real_t >& globalAcceleration, real_t surfaceTension,
bool useSimpleMassExchange, real_t cellConversionThreshold,
real_t cellConversionForceThreshold, BlockDataID relaxationRateFieldID = BlockDataID(),
real_t smagorinskyConstant = real_c(0))
: blockForest_(blockForest), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), fillFieldID_(fillFieldID),
- forceFieldID_(forceFieldID), normalFieldID_(normalFieldID), curvatureFieldID_(curvatureFieldID),
+ forceDensityFieldID_(forceDensityFieldID), normalFieldID_(normalFieldID), curvatureFieldID_(curvatureFieldID),
bubbleModel_(bubbleModel), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling),
pdfReconstructionModel_(pdfReconstructionModel), pdfRefillingModel_({ pdfRefillingModel }),
excessMassDistributionModel_({ excessMassDistributionModel }), relaxationRate_(relaxationRate),
- globalForce_(globalForce), surfaceTension_(surfaceTension), enableForceWeighting_(enableForceWeighting),
+ globalAcceleration_(globalAcceleration), surfaceTension_(surfaceTension),
useSimpleMassExchange_(useSimpleMassExchange), cellConversionThreshold_(cellConversionThreshold),
cellConversionForceThreshold_(cellConversionForceThreshold), relaxationRateFieldID_(relaxationRateFieldID),
smagorinskyConstant_(smagorinskyConstant)
@@ -141,16 +141,33 @@ class SurfaceDynamicsHandler
Set< SUID >::emptySet(), StateSweep::onlyGasAndBoundary)
<< Sweep(emptySweep, "Empty sweep: boundary handling", StateSweep::onlyGasAndBoundary);
- if (enableForceWeighting_)
+ if (!(floatIsEqual(globalAcceleration_[0], real_c(0), real_c(1e-14)) &&
+ floatIsEqual(globalAcceleration_[1], real_c(0), real_c(1e-14)) &&
+ floatIsEqual(globalAcceleration_[2], real_c(0), real_c(1e-14))))
{
// add sweep for weighting force in interface cells with fill level and density
- const ForceWeightingSweep< LatticeModel_T, FlagField_T, VectorField_T, ScalarField_T > forceWeightingSweep(
- forceFieldID_, pdfFieldID_, flagFieldID_, fillFieldID_, flagInfo, globalForce_);
- timeloop.add() << Sweep(forceWeightingSweep, "Sweep: force weighting", Set< SUID >::emptySet(),
- StateSweep::onlyGasAndBoundary)
- << Sweep(emptySweep, "Empty sweep: force weighting", StateSweep::onlyGasAndBoundary)
- << AfterFunction(CommunicationCorner_T(blockForest_, forceFieldID_),
- "Communication: after force weighting sweep");
+ if constexpr (useCodegen)
+ {
+ // different versions for codegen because pystencils does not support 'Ghostlayerfield<Vector3(), 1>'
+ const ForceDensityCodegenSweep< LatticeModel_T, FlagField_T, VectorFieldFlattened_T, ScalarField_T >
+ forceDensityCodegenSweep(forceDensityFieldID_, pdfFieldID_, flagFieldID_, fillFieldID_, flagInfo,
+ globalAcceleration_);
+ timeloop.add() << Sweep(forceDensityCodegenSweep, "Sweep: force weighting", Set< SUID >::emptySet(),
+ StateSweep::onlyGasAndBoundary)
+ << Sweep(emptySweep, "Empty sweep: force weighting", StateSweep::onlyGasAndBoundary)
+ << AfterFunction(CommunicationCorner_T(blockForest_, forceDensityFieldID_),
+ "Communication: after force weighting sweep");
+ }
+ else
+ {
+ const ForceDensitySweep< LatticeModel_T, FlagField_T, VectorField_T, ScalarField_T > forceDensitySweep(
+ forceDensityFieldID_, pdfFieldID_, flagFieldID_, fillFieldID_, flagInfo, globalAcceleration_);
+ timeloop.add() << Sweep(forceDensitySweep, "Sweep: force weighting", Set< SUID >::emptySet(),
+ StateSweep::onlyGasAndBoundary)
+ << Sweep(emptySweep, "Empty sweep: force weighting", StateSweep::onlyGasAndBoundary)
+ << AfterFunction(CommunicationCorner_T(blockForest_, forceDensityFieldID_),
+ "Communication: after force weighting sweep");
+ }
}
// sweep for
@@ -412,7 +429,7 @@ class SurfaceDynamicsHandler
BlockDataID pdfFieldID_;
BlockDataID flagFieldID_;
BlockDataID fillFieldID_;
- BlockDataID forceFieldID_;
+ BlockDataID forceDensityFieldID_;
ConstBlockDataID normalFieldID_;
ConstBlockDataID curvatureFieldID_;
@@ -424,9 +441,8 @@ class SurfaceDynamicsHandler
PdfRefillingModel pdfRefillingModel_;
ExcessMassDistributionModel excessMassDistributionModel_;
real_t relaxationRate_;
- Vector3< real_t > globalForce_;
+ Vector3< real_t > globalAcceleration_;
real_t surfaceTension_;
- bool enableForceWeighting_;
bool useSimpleMassExchange_;
real_t cellConversionThreshold_;
real_t cellConversionForceThreshold_;
diff --git a/src/lbm/free_surface/dynamics/functionality/AdvectMass.h b/src/lbm/free_surface/dynamics/functionality/AdvectMass.h
index 5e198e2db..9ec9e655c 100644
--- a/src/lbm/free_surface/dynamics/functionality/AdvectMass.h
+++ b/src/lbm/free_surface/dynamics/functionality/AdvectMass.h
@@ -207,8 +207,8 @@ real_t advectMass(const FlagField_T* flagField, const ConstScalarIt_T& fillSrc,
bool neighborNoFluidNeig = !flagInfo.isLiquid(neighIt.neighbor(relevantDir[0], relevantDir[1], relevantDir[2]));
bool neighborStandardCell = !(neighborNoGasNeig || neighborNoFluidNeig);
bool neighborWettingCell = flagInfo.isKeepInterfaceForWetting(flagFieldIt.neighbor(
- relevantDir[0], relevantDir[1],
- relevantDir[2])); // evaluate flag of this cell (flagFieldIt) and not the neighborhood flags (neighIt)
+ relevantDir[0], relevantDir[1],
+ relevantDir[2])); // evaluate flag of this cell (flagFieldIt) and not the neighborhood flags (neighIt)
if (neighborNoGasNeig && neighborNoFluidNeig &&
!neighborWettingCell) // neighboring cell has only interface neighbors
diff --git a/src/lbm/free_surface/surface_geometry/Utility.cpp b/src/lbm/free_surface/surface_geometry/Utility.cpp
index 9b1582c66..ef9f2520a 100644
--- a/src/lbm/free_surface/surface_geometry/Utility.cpp
+++ b/src/lbm/free_surface/surface_geometry/Utility.cpp
@@ -391,9 +391,9 @@ Vector3< real_t > getInterfacePoint(const Vector3< real_t >& normal, real_t fill
real_t getCellEdgeIntersection(const Vector3< real_t >& edgePoint, const Vector3< real_t >& edgeDirection,
const Vector3< real_t >& normal, const Vector3< real_t >& surfacePoint)
{
- //#ifndef BELOW_CELL
- //# define BELOW_CELL (-10)
- //#endif
+ // #ifndef BELOW_CELL
+ // # define BELOW_CELL (-10)
+ // #endif
#ifndef ABOVE_CELL
# define ABOVE_CELL (-20)
diff --git a/tests/lbm/free_surface/dynamics/CellConversionTest.cpp b/tests/lbm/free_surface/dynamics/CellConversionTest.cpp
index 6f7645c84..b11c64421 100644
--- a/tests/lbm/free_surface/dynamics/CellConversionTest.cpp
+++ b/tests/lbm/free_surface/dynamics/CellConversionTest.cpp
@@ -109,8 +109,8 @@ void testCellConversion()
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill levels", real_c(1.0), field::fzyx, uint_c(1));
- BlockDataID forceFieldID = field::addToStorage< VectorField_T >(
- blockForest, "Force field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
+ BlockDataID forceDensityFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Force density field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
BlockDataID curvatureFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Curvature", real_c(0.0), field::fzyx, uint_c(1));
BlockDataID normalFieldID = field::addToStorage< VectorField_T >(
@@ -162,9 +162,9 @@ void testCellConversion()
// add various sweeps for surface dynamics
SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, "NormalBasedKeepCenter", "EquilibriumRefilling", "EvenlyNewInterface",
- relaxRate, Vector3< real_t >(real_c(0)), real_c(0), false, false, real_c(1e-3), real_c(1e-1));
+ relaxRate, Vector3< real_t >(real_c(0)), real_c(0), false, real_c(1e-3), real_c(1e-1));
dynamicsHandler.addSweeps(timeloop);
real_t initialMass =
diff --git a/tests/lbm/free_surface/dynamics/CodegenTest.cpp b/tests/lbm/free_surface/dynamics/CodegenTest.cpp
index a141a0a3c..9a00890aa 100644
--- a/tests/lbm/free_surface/dynamics/CodegenTest.cpp
+++ b/tests/lbm/free_surface/dynamics/CodegenTest.cpp
@@ -27,6 +27,7 @@
#include "core/debug/TestSubsystem.h"
#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/InitFunctions.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
@@ -44,8 +45,9 @@ namespace free_surface
namespace CodegenTest
{
-using ScalarField_T = GhostLayerField< real_t, 1 >;
-using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+using VectorFieldFlattened_T = GhostLayerField< real_t, 3 >;
template< bool useCodegen >
void runSimulation()
@@ -54,6 +56,7 @@ void runSimulation()
using ForceModel_T = lbm::force_model::GuoField< VectorField_T >;
using LatticeModel_T = typename std::conditional< useCodegen, lbm::GeneratedLatticeModel_FreeSurface,
lbm::D3Q19< CollisionModel_T, true, ForceModel_T, 2 > >::type;
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
using Communication_T = blockforest::SimpleCommunication< typename LatticeModel_T::CommunicationStencil >;
@@ -76,39 +79,45 @@ void runSimulation()
true, true, true); // periodicity
// physics parameters
- const real_t dropDiameter = real_c(2);
- const real_t relaxationRate = real_c(1.8);
- const real_t surfaceTension = real_c(1e-5);
- const bool enableWetting = false;
- const real_t contactAngle = real_c(0);
- const Vector3< real_t > force = Vector3< real_t >(real_c(1e-5), real_c(0), real_c(0));
+ const real_t dropDiameter = real_c(2);
+ const real_t relaxationRate = real_c(1.8);
+ const real_t surfaceTension = real_c(1e-5);
+ const bool enableWetting = false;
+ const real_t contactAngle = real_c(0);
+ const Vector3< real_t > acceleration = Vector3< real_t >(real_c(1e-5), real_c(0), real_c(0));
// model parameters
const std::string pdfReconstructionModel = "NormalBasedKeepCenter";
const std::string pdfRefillingModel = "EquilibriumRefilling";
const std::string excessMassDistributionModel = "EvenlyAllInterface";
const std::string curvatureModel = "FiniteDifferenceMethod";
- const bool enableForceWeighting = false;
const bool useSimpleMassExchange = false;
const real_t cellConversionThreshold = real_c(1e-2);
const real_t cellConversionForceThreshold = real_c(1e-1);
// add force field
- const BlockDataID forceFieldID =
- field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
+ std::shared_ptr< BlockDataID > forceDensityFieldIDPtr = nullptr;
std::shared_ptr< LatticeModel_T > latticeModel;
// create lattice model
if constexpr (useCodegen)
{
- latticeModel = std::make_shared< LatticeModel_T >(force[0], force[1], force[2], relaxationRate);
+ // add force field of type 'GhostLayerField< real_t, 3 >' as required by pystencils
+ forceDensityFieldIDPtr = std::make_shared< BlockDataID >(field::addToStorage< VectorFieldFlattened_T >(
+ blockForest, "Force density field (codegen)", real_c(0), field::fzyx, uint_c(1)));
+ latticeModel = std::make_shared< LatticeModel_T >(*forceDensityFieldIDPtr, relaxationRate);
}
else
{
- latticeModel = std::make_shared< LatticeModel_T >(CollisionModel_T(relaxationRate), ForceModel_T(forceFieldID));
+ forceDensityFieldIDPtr = std::make_shared< BlockDataID >(field::addToStorage< VectorField_T >(
+ blockForest, "Force density field", Vector3< real_t >(0), field::fzyx, uint_c(1)));
+ latticeModel =
+ std::make_shared< LatticeModel_T >(CollisionModel_T(relaxationRate), ForceModel_T(*forceDensityFieldIDPtr));
}
+ WALBERLA_ASSERT_NOT_NULLPTR(forceDensityFieldIDPtr);
+
// add various fields
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", *latticeModel, field::fzyx);
const BlockDataID fillFieldID =
@@ -117,7 +126,8 @@ void runSimulation()
// add boundary handling
const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
- const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
+ const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
+ const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
// add drop to fill level field
const geometry::Sphere sphereDrop(Vector3< real_t >(real_c(0.5) * real_c(domainSize[0]),
@@ -129,8 +139,20 @@ void runSimulation()
// initialize flag field from fill level field
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
- // initial Communication_Tunication
- Communication_T(blockForest, pdfFieldID, fillFieldID, flagFieldID, forceFieldID)();
+ if constexpr (useCodegen)
+ {
+ initForceDensityFieldCodegen< PdfField_T, FlagField_T, VectorFieldFlattened_T, ScalarField_T >(
+ blockForest, *forceDensityFieldIDPtr, fillFieldID, pdfFieldID, flagFieldID, flagInfo, acceleration);
+ }
+ else
+ {
+ // initialize force density field
+ initForceDensityField< PdfField_T, FlagField_T, VectorField_T, ScalarField_T >(
+ blockForest, *forceDensityFieldIDPtr, fillFieldID, pdfFieldID, flagFieldID, flagInfo, acceleration);
+ }
+
+ // initial communication
+ Communication_T(blockForest, pdfFieldID, fillFieldID, flagFieldID, *forceDensityFieldIDPtr)();
// add bubble model
const std::shared_ptr< bubble_model::BubbleModelConstantPressure > bubbleModel =
@@ -154,11 +176,12 @@ void runSimulation()
geometryHandler.addSweeps(timeloop);
// add boundary handling for standard boundaries and free surface boundaries
- SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T, useCodegen > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
- freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
- cellConversionForceThreshold);
+ const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T, useCodegen,
+ VectorFieldFlattened_T >
+ dynamicsHandler(blockForest, pdfFieldID, flagFieldID, fillFieldID, *forceDensityFieldIDPtr, normalFieldID,
+ curvatureFieldID, freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel,
+ pdfRefillingModel, excessMassDistributionModel, relaxationRate, acceleration, surfaceTension,
+ useSimpleMassExchange, cellConversionThreshold, cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
@@ -172,35 +195,35 @@ void runSimulation()
WALBERLA_FOR_ALL_CELLS(fillFieldIt, fillField, {
if (fillFieldIt.cell() == Cell(cell_idx_c(1), cell_idx_c(1), cell_idx_c(1)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.504035), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.512081), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(2), cell_idx_c(1), cell_idx_c(1)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.538017), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.529967), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(1), cell_idx_c(2), cell_idx_c(1)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.504035), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.512081), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(2), cell_idx_c(2), cell_idx_c(1)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.538017), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.529967), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(1), cell_idx_c(1), cell_idx_c(2)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.504035), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.512081), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(2), cell_idx_c(1), cell_idx_c(2)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.538017), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.529967), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(1), cell_idx_c(2), cell_idx_c(2)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.504035), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.512081), real_c(1e-4));
}
if (fillFieldIt.cell() == Cell(cell_idx_c(2), cell_idx_c(2), cell_idx_c(2)))
{
- WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.538017), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.529967), real_c(1e-4));
}
}); // WALBERLA_FOR_ALL_CELLS
}
diff --git a/tests/lbm/free_surface/dynamics/InflowTest.cpp b/tests/lbm/free_surface/dynamics/InflowTest.cpp
index da8b5df18..49423ada6 100644
--- a/tests/lbm/free_surface/dynamics/InflowTest.cpp
+++ b/tests/lbm/free_surface/dynamics/InflowTest.cpp
@@ -90,8 +90,8 @@ void testInflow()
BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill levels", real_c(0.0), field::fzyx, uint_c(2));
- BlockDataID forceFieldID = field::addToStorage< VectorField_T >(
- blockForest, "Force field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
+ BlockDataID forceDensityFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Force density field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
BlockDataID densityAdaptor = field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::Density >(
blockForest, pdfFieldID, "DensityAdaptor");
BlockDataID velocityAdaptor = field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::VelocityVector >(
@@ -155,9 +155,9 @@ void testInflow()
// add boundary handling for standard boundaries and free surface boundaries
SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, "NormalBasedKeepCenter", "EquilibriumRefilling", "EvenlyNewInterface",
- relaxRate, Vector3< real_t >(real_c(0)), real_c(0), false, false, real_c(1e-3), real_c(1e-1));
+ relaxRate, Vector3< real_t >(real_c(0)), real_c(0), false, real_c(1e-3), real_c(1e-1));
dynamicsHandler.addSweeps(timeloop);
diff --git a/tests/lbm/free_surface/dynamics/LatticeModelGenerationFreeSurface.py b/tests/lbm/free_surface/dynamics/LatticeModelGenerationFreeSurface.py
index 0ef55dce4..a1b0b251a 100644
--- a/tests/lbm/free_surface/dynamics/LatticeModelGenerationFreeSurface.py
+++ b/tests/lbm/free_surface/dynamics/LatticeModelGenerationFreeSurface.py
@@ -1,5 +1,5 @@
import sympy as sp
-
+import pystencils as ps
from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
from lbmpy.enums import ForceModel, Method, Stencil
from lbmpy.stencils import LBStencil
@@ -7,28 +7,30 @@ from lbmpy.stencils import LBStencil
from pystencils_walberla import CodeGeneration
from lbmpy_walberla import generate_lattice_model
-# general parameters
-stencil = LBStencil(Stencil.D3Q19)
-omega = sp.Symbol('omega')
-force = sp.symbols('force_:3')
-layout = 'fzyx'
+with CodeGeneration() as ctx:
+ # general parameters
+ layout = 'fzyx'
+ data_type = "float64" if ctx.double_accuracy else "float32"
+
+ stencil = LBStencil(Stencil.D3Q19)
+ omega = sp.Symbol('omega')
+ force_field = ps.fields(f"force(3): {data_type}[3D]", layout=layout)
-# method definition
-lbm_config = LBMConfig(stencil=stencil,
- method=Method.SRT,
- relaxation_rate=omega,
- compressible=True,
- force=force,
- force_model=ForceModel.GUO,
- zero_centered=False,
- streaming_pattern='pull') # free surface implementation only works with pull pattern
+ # method definition
+ lbm_config = LBMConfig(stencil=stencil,
+ method=Method.SRT,
+ relaxation_rate=omega,
+ compressible=True,
+ force=force_field,
+ force_model=ForceModel.GUO,
+ zero_centered=False,
+ streaming_pattern='pull') # free surface implementation only works with pull pattern
-# optimizations to be used by the code generator
-lbm_opt = LBMOptimisation(cse_global=True,
- field_layout=layout)
+ # optimizations to be used by the code generator
+ lbm_opt = LBMOptimisation(cse_global=True,
+ field_layout=layout)
-collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
- lbm_optimisation=lbm_opt)
+ collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
+ lbm_optimisation=lbm_opt)
-with CodeGeneration() as ctx:
generate_lattice_model(ctx, "GeneratedLatticeModel_FreeSurface", collision_rule, field_layout=layout)
diff --git a/tests/lbm/free_surface/dynamics/WettingConversionTest.cpp b/tests/lbm/free_surface/dynamics/WettingConversionTest.cpp
index 8a2590c69..b6e03ede0 100644
--- a/tests/lbm/free_surface/dynamics/WettingConversionTest.cpp
+++ b/tests/lbm/free_surface/dynamics/WettingConversionTest.cpp
@@ -83,8 +83,8 @@ std::vector< Cell > runSimulation(uint_t timesteps, real_t contactAngle)
field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(0.0), field::fzyx, uint_c(1));
// add dummy force field
- BlockDataID forceFieldID = field::addToStorage< VectorField_T >(
- blockForest, "Force field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
+ BlockDataID forceDensityFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Force density field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
// add boundary handling
const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
@@ -103,7 +103,7 @@ std::vector< Cell > runSimulation(uint_t timesteps, real_t contactAngle)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// initial communication
- Communication_T(blockForest, pdfFieldID, fillFieldID, flagFieldID, forceFieldID)();
+ Communication_T(blockForest, pdfFieldID, fillFieldID, flagFieldID, forceDensityFieldID)();
// add (dummy) bubble model
const bool disableSplits = true; // necessary if a gas bubble could split
@@ -124,9 +124,9 @@ std::vector< Cell > runSimulation(uint_t timesteps, real_t contactAngle)
// add surface dynamics handler
SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
freeSurfaceBoundaryHandling, bubbleModel, "NormalBasedKeepCenter", "EquilibriumRefilling", "EvenlyNewInterface",
- relaxRate, Vector3< real_t >(real_c(0)), real_c(1e-2), false, false, real_c(1e-3), real_c(1e-1));
+ relaxRate, Vector3< real_t >(real_c(0)), real_c(1e-2), false, real_c(1e-3), real_c(1e-1));
dynamicsHandler.addSweeps(timeloop);
timeloop.run();
--
GitLab