From a4c9e800ce2d59bbd93931251fc9106f800ba35a Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Fri, 22 Jul 2022 13:06:53 +0200
Subject: [PATCH] Add drop impact and Taylor bubble showcase for phase-field
LBM
---
.../PhaseFieldAllenCahn/CPU/CMakeLists.txt | 2 +-
.../CPU/InitializerFunctions.cpp | 278 +++++++++++++-----
.../CPU/InitializerFunctions.h | 16 +-
.../PhaseFieldAllenCahn/CPU/multiphase.cpp | 97 +++++-
.../CPU/multiphase_RTI_3D.py | 2 +-
.../CPU/multiphase_codegen.py | 5 +-
.../CPU/multiphase_drop.py | 108 +++++++
.../CPU/multiphase_rising_bubble.py | 2 +-
.../CPU/multiphase_taylor_bubble.py | 184 ++++++++++++
.../PhaseFieldAllenCahn/GPU/CMakeLists.txt | 2 +-
10 files changed, 605 insertions(+), 91 deletions(-)
create mode 100755 apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py
create mode 100755 apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
index abce2c997..383c9d31a 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
@@ -18,4 +18,4 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenCPU
waLBerla_add_executable(NAME multiphaseCPU
FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp multiphase_codegen.py
- DEPENDS blockforest core field postprocessing python_coupling lbm geometry timeloop gui PhaseFieldCodeGenCPU)
+ DEPENDS blockforest core field postprocessing python_coupling lbm geometry timeloop PhaseFieldCodeGenCPU)
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
index d7cb146a0..da5d24f34 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
@@ -23,38 +23,82 @@
#include "field/FlagField.h"
#include "field/communication/PackInfo.h"
-#include "field/vtk/VTKWriter.h"
-
-#include "python_coupling/DictWrapper.h"
namespace walberla
{
-using PhaseField_T = GhostLayerField< real_t, 1 >;
-using FlagField_T = FlagField< uint8_t >;
+using PhaseField_T = GhostLayerField< real_t, 1 >;
+using FlagField_T = FlagField< uint8_t >;
+using VelocityField_T = GhostLayerField< real_t, 3 >;
void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t R = real_c(10.0),
+ BlockDataID velocityFieldID, const real_t R = real_c(10.0),
const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0),
- const bool bubble = true, const real_t W = real_c(5.0))
+ const bool bubble = true, const real_t W = real_c(5.0),
+ const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0)))
+{
+ for (auto& block : *blocks)
+ {
+ auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+ auto velocityField = block.getData< VelocityField_T >(velocityFieldID);
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ real_t Ri =
+ real_c(sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])));
+ if (bubble) { phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); }
+ else
+ {
+ phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W));
+
+ if (Ri < R + W)
+ {
+ velocityField->get(x, y, z, 0) = initialVelocity[0];
+ velocityField->get(x, y, z, 1) = initialVelocity[1];
+ velocityField->get(x, y, z, 2) = initialVelocity[2];
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+ }
+}
+
+void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W,
+ real_t poolDepth)
{
for (auto& block : *blocks)
{
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- // clang-format off
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell;
+
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+
+ phaseField->get(x, y, z) +=
+ real_c(0.5) - real_c(0.5) * real_c(tanh(real_c(2.0) * (real_c(globalCell[1]) - poolDepth) / W));
+ if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+ }
+}
+
+void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID densityFieldID,
+ real_t GravitationalAcceleration, real_t poolDepth)
+{
+ for (auto& block : *blocks)
+ {
+ auto densityField = block.getData< PhaseField_T >(densityFieldID);
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(densityField, {
+ Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
- (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
- (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
- if (bubble) phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W));
- else phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W));
- )
- // clang-format on
+ real_t pressure = real_c(3.0) * GravitationalAcceleration * (real_c(globalCell[1]) - poolDepth);
+ densityField->get(x, y, z) += pressure;
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0), const real_t Donut_x0 = real_c(40.0))
+ const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0),
+ const real_t Donut_x0 = real_c(40.0))
{
auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
@@ -62,16 +106,66 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
for (auto& block : *blocks)
{
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
- phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t Ri = D * real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2)));
+ real_t Ri =
+ D *
+ real_c(sqrt(pow(H, 2) -
+ pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2)));
- real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx));
+ real_t shifter = real_c(atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx)));
if (shifter < 0) shifter = shifter + 2 * math::pi;
- if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) {
+ if ((real_c(globalCell[1]) < Donut_x0 + Ri * real_c(sin(shifter / 2.0))) && (real_c(globalCell[1]) > Donut_x0 - Ri))
+ {
phaseField->get(x, y, z) = 0.0;
- } else { phaseField->get(x, y, z) = real_c(1.0); })
+ }
+ else { phaseField->get(x, y, z) = real_c(1.0); }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+ }
+}
+
+void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& phaseFieldID,
+ real_t R, real_t H, real_t L, real_t W)
+{
+ auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
+ auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
+
+ for (auto& block : *blocks)
+ {
+ PhaseField_T* phaseField = block.getData< PhaseField_T >(phaseFieldID);
+
+ // initialize top and bottom walls of cylinder
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+
+ const real_t Ri = real_c(sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)));
+
+ const real_t cylinderMidY = real_c(H) + real_c(0.5) * real_c(L);
+
+ const real_t distance = std::abs(real_c(globalCell[1]) - cylinderMidY);
+ const real_t lHalf = real_c(0.5) * real_c(L);
+
+ // cylinder side walls
+ if (real_c(globalCell[1]) >= H - W && real_c(globalCell[1]) < H + L + W)
+ {
+ phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (Ri - R) / W));
+
+ // cylinder top and bottom walls (must be added to not overwrite the side wall initialization)
+ if (Ri < R + W)
+ {
+ phaseField->get(x, y, z) +=
+ real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (distance - lHalf) / W));
+
+ // limit maximum values of phase field
+ if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); }
+ if (phaseField->get(x, y, z) < real_c(0)) { phaseField->get(x, y, z) = real_c(0); }
+ }
+ }
+ else { phaseField->get(x, y, z) = real_c(1); }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
@@ -92,8 +186,8 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block
// distance in between the bubbles
real_t dist = real_c(4.0);
- auto nx = uint_c(floor(X / (D + dist * W)));
- auto nz = uint_c(floor(Z / (D + dist * W)));
+ auto nx = uint_c(floor(X / (D + dist * W)));
+ auto nz = uint_c(floor(Z / (D + dist * W)));
// fluctuation of the bubble radii
std::vector< std::vector< real_t > > fluctuation_radius(nx, std::vector< real_t >(nz, 0.0));
@@ -113,39 +207,50 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block
for (auto& block : *blocks)
{
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
- phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- for (unsigned int i = 0; i < nx; ++i) {
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ for (unsigned int i = 0; i < nx; ++i)
+ {
for (unsigned int j = 0; j < nz; ++j)
{
- bubbleMidPoint[0] = real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
+ bubbleMidPoint[0] =
+ real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
bubbleMidPoint[1] = R + W + 4;
- bubbleMidPoint[2] = real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
+ bubbleMidPoint[2] =
+ real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j];
- real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
- (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
- (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]));
- if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) &&
- real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W))
- phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W);
+ real_t Ri = real_c(
+ sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) +
+ (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) +
+ (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])));
+ if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) &&
+ real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) &&
+ real_c(globalCell[2]) >= real_c(j) * (D + dist * W) &&
+ real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W))
+ phaseField->get(x, y, z) =
+ real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W));
if (real_c(globalCell[0]) > real_c(nx) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0);
if (real_c(globalCell[2]) > real_c(nz) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0);
}
}
- if (real_c(globalCell[1]) > gas_top) {
- phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W);
- })
+ if (real_c(globalCell[1]) > gas_top)
+ {
+ phaseField->get(x, y, z) =
+ real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W));
+ }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
const real_t W = real_c(5.0), const bool pipe = true)
{
- auto X = real_c(blocks->getDomainCellBB().xMax());
- auto Z = real_c(blocks->getDomainCellBB().zMax());
- auto halfY = real_c(blocks->getDomainCellBB().yMax()) / real_c(2.0);
+ auto X = real_c(blocks->getDomainCellBB().xMax());
+ auto Z = real_c(blocks->getDomainCellBB().zMax());
+ auto halfY = real_c(blocks->getDomainCellBB().yMax()) / real_c(2.0);
auto perturbation = real_c(0.05);
if (pipe)
@@ -153,12 +258,16 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
for (auto& block : *blocks)
{
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
- phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) +
- (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2));
- if (R > X) R = X; real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X);
- phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)));)
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ real_t R = real_c(sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) +
+ (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)));
+ if (R > X) R = X;
+ real_t tmp = perturbation * X * real_c(cos((real_c(2.0) * math::pi * R) / X));
+ phaseField->get(x, y, z) =
+ real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0))));
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
else
@@ -166,11 +275,15 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc
for (auto& block : *blocks)
{
auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
- phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t tmp = perturbation * X *
- (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X));
- phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)));)
+ (real_c(cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X)) +
+ real_c(cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X)));
+ phaseField->get(x, y, z) =
+ real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0))));
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
}
@@ -193,31 +306,41 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
{
auto flagField = block.template getData< FlagField_T >(flagFieldID);
auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell;
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+ Cell globalCell;
blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
real_t R1;
- if (real_c(globalCell[1]) <= start_transition) {
- R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
- } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ if (real_c(globalCell[1]) <= start_transition)
+ {
+ R1 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
+ }
+ else if (real_c(globalCell[1]) > start_transition &&
+ real_c(globalCell[1]) < start_transition + length_transition)
+ {
real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
- real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp));
- R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) +
- (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
- } else {
- R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) +
- (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
+ real_t shift_tmp = shift * real_c(0.5) * (1 - real_c(cos(tmp)));
+ R1 = real_c(sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
+ }
+ else
+ {
+ R1 = real_c(sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
}
- real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
+ real_t R2 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag);
- if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
+ if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
else
{
auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0);
auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0);
-
+
auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0);
auto const shift = real_c(eccentricity * R_in);
@@ -227,21 +350,24 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl
{
auto flagField = block.template getData< FlagField_T >(flagFieldID);
auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID);
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
- flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- if (real_c(globalCell[1]) <= start_transition) {
- R_tmp = R_in;
- } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) {
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+ Cell globalCell;
+ blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
+ if (real_c(globalCell[1]) <= start_transition) { R_tmp = R_in; }
+ else if (real_c(globalCell[1]) > start_transition &&
+ real_c(globalCell[1]) < start_transition + length_transition)
+ {
real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition);
- real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp));
- R_tmp = R_in + shift_tmp;
- } else {
- R_tmp = R_in + shift;
+ real_t shift_tmp = shift * real_c(0.5) * (1 - real_c(cos(tmp)));
+ R_tmp = R_in + shift_tmp;
}
+ else { R_tmp = R_in + shift; }
- real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz));
+ real_t R2 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) +
+ (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)));
if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag);
- if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);)
+ if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
}
}
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
index 585639ac9..ae99d3ab0 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
@@ -31,12 +31,24 @@
namespace walberla
{
-void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
- Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5);
+void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks,
+ BlockDataID phaseFieldID, BlockDataID velocityFieldID,
+ real_t R,
+ Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5,
+ const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0)));
+
+void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W,
+ real_t poolDepth);
+
+void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID densityFieldID,
+ real_t GravitationalAcceleration, real_t poolDepth);
void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t D = 5,
real_t H = 2, real_t DT = 20, real_t Donut_x0 = 40);
+void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& phaseFieldID,
+ real_t R, real_t H, real_t L, real_t W);
+
void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
real_t W = 5);
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
index 8e0f2cb62..2ed8706f3 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
@@ -30,9 +30,12 @@
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
+#include "geometry/mesh/TriangleMeshIO.h"
#include "lbm/PerformanceEvaluation.h"
+#include "postprocessing/FieldToSurfaceMesh.h"
+
#include "python_coupling/CreateConfig.h"
#include "python_coupling/PythonCallback.h"
@@ -68,8 +71,9 @@ int main(int argc, char** argv)
// ADD DOMAIN PARAMETERS //
//////////////////////////
- auto domainSetup = config->getOneBlock("DomainSetup");
- const bool tube = domainSetup.getParameter< bool >("tube", true);
+ auto domainSetup = config->getOneBlock("DomainSetup");
+ Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize");
+ const bool tube = domainSetup.getParameter< bool >("tube", false);
////////////////////////////////////////
// ADD GENERAL SIMULATION PARAMETERS //
@@ -92,7 +96,8 @@ int main(int argc, char** argv)
field::addToStorage< PdfField_hydro_T >(blocks, "LB velocity field", real_c(0.0), field::fzyx);
BlockDataID velocity_field_ID =
field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
- BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
+ BlockDataID density_field_ID = field::addToStorage< PhaseField_T >(blocks, "density", real_c(1.0), field::fzyx);
+ BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
@@ -108,15 +113,45 @@ int main(int argc, char** argv)
const real_t interface_thickness = physical_parameters.getParameter< real_t >("interface_thickness");
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field")
- if (scenario == 1)
+ switch (scenario)
{
+ case 1: {
auto bubbleParameters = config->getOneBlock("Bubble");
const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
const bool bubble = bubbleParameters.getParameter< bool >("bubble", true);
- initPhaseField_sphere(blocks, phase_field_ID, bubbleRadius, bubbleMidPoint, bubble);
+ initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble);
+ break;
+ }
+ case 2: {
+ initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube);
+ break;
+ }
+ case 3: {
+ auto dropParameters = config->getOneBlock("Drop");
+ const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("drop_mid_point");
+ const real_t dropRadius = dropParameters.getParameter< real_t >("drop_radius");
+ const real_t poolDepth = dropParameters.getParameter< real_t >("pool_depth");
+ const Vector3< real_t > impact_velocity = dropParameters.getParameter< Vector3< real_t > >("impact_velocity");
+ init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth);
+
+ initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false,
+ interface_thickness, impact_velocity);
+ initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth);
+ break;
+ }
+ case 4: {
+ auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble");
+ const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("bubble_radius");
+ const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("initial_height");
+ const real_t Length = TaylorBubbleParameters.getParameter< real_t >("length");
+ init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length,
+ real_c(interface_thickness));
+ break;
+ }
+ default:
+ WALBERLA_ABORT("Scenario is not defined.")
}
- else if (scenario == 2) { initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); }
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done")
/////////////////
@@ -125,7 +160,8 @@ int main(int argc, char** argv)
pystencils::initialize_phase_field_distributions init_h(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID,
interface_thickness);
- pystencils::initialize_velocity_based_distributions init_g(hydrodynamic_PDFs_ID, velocity_field_ID);
+ pystencils::initialize_velocity_based_distributions init_g(density_field_ID, hydrodynamic_PDFs_ID,
+ velocity_field_ID);
pystencils::phase_field_LB_step phase_field_LB_step(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID,
mobility, interface_thickness);
@@ -180,6 +216,29 @@ int main(int argc, char** argv)
if (boundariesConfig)
{
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
+ if (tube)
+ {
+ // initialize cylindrical domain walls
+ const Vector3< real_t > domainCylinderBottomEnd =
+ Vector3< real_t >(real_c(domainSize[0]) * real_c(0.5), real_c(0), real_c(domainSize[2]) * real_c(0.5));
+ const Vector3< real_t > domainCylinderTopEnd = Vector3< real_t >(
+ real_c(domainSize[0]) * real_c(0.5), real_c(domainSize[1]), real_c(domainSize[2]) * real_c(0.5));
+ const real_t radius = real_c(domainSize[0]) * real_c(0.5);
+ const geometry::Cylinder cylinderTube(domainCylinderBottomEnd, domainCylinderTopEnd, radius);
+
+ for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+ {
+ FlagField_T* flagField = blockIt->template getData< FlagField_T >(flagFieldID);
+ auto wallFlag = flagField->getOrRegisterFlag(wallFlagUID);
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+ Cell globalCell = Cell(x, y, z);
+ blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
+ const Vector3< real_t > globalPoint = blocks->getCellCenter(globalCell);
+
+ if (!geometry::contains(cylinderTube, globalPoint)) { addFlag(flagField->get(x, y, z), wallFlag); }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+ }
+ }
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID);
}
@@ -248,6 +307,30 @@ int main(int argc, char** argv)
},
"Python callback");
+ int meshWriteFrequency = parameters.getParameter< int >("meshWriteFrequency", 0);
+ int counter = 0;
+ int targetRank = 0;
+ if (meshWriteFrequency > 0)
+ {
+ timeloop.addFuncAfterTimeStep(
+ [&]() {
+ if (timeloop.getCurrentTimeStep() % uint_t(meshWriteFrequency) == 0)
+ {
+ auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field_ID, 0.5, 0,
+ true, targetRank, MPI_COMM_WORLD);
+ WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank)
+ {
+ std::string path = "";
+ std::ostringstream out;
+ out << std::internal << std::setfill('0') << std::setw(6) << counter;
+ geometry::writeMesh(path + out.str() + ".obj", *mesh);
+ counter++;
+ }
+ }
+ },
+ "Mesh writer");
+ }
+
// remaining time logger
timeloop.addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
index d7e2f1959..58ac2fbac 100755
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
@@ -61,7 +61,7 @@ class Scenario:
# everything else
self.dbFile = "RTI.csv"
- self.scenario = 2 # 1 rising bubble or droplet, 2 RTI
+ self.scenario = 2 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
self.config_dict = self.config()
@wlb.member_callback
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
index 1b5113454..0b117fc23 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
@@ -28,7 +28,7 @@ with CodeGeneration() as ctx:
stencil_hydro = LBStencil(Stencil.D3Q27)
assert (stencil_phase.D == stencil_hydro.D)
- # In the codegneration skript we only need the symbolic parameters
+ # In the code-generation skript, we only need the symbolic parameters
parameters = AllenCahnParameters(0, 0, 0, 0, 0)
########################
@@ -37,6 +37,7 @@ with CodeGeneration() as ctx:
# velocity field
u = fields(f"vel_field({stencil_hydro.D}): {field_type}[{stencil_hydro.D}D]", layout='fzyx')
+ density_field = fields(f"density_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
# phase-field
C = fields(f"phase_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
C_tmp = fields(f"phase_field_tmp: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
@@ -87,7 +88,7 @@ with CodeGeneration() as ctx:
# create the kernels for the initialization of the g and h field
h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
- g_updates = initializer_kernel_hydro_lb(method_hydro, 1.0, u, g)
+ g_updates = initializer_kernel_hydro_lb(method_hydro, density_field, u, g)
force_h = interface_tracking_force(C, stencil_phase, parameters)
hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py
new file mode 100755
index 000000000..2845605b1
--- /dev/null
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py
@@ -0,0 +1,108 @@
+import waLBerla as wlb
+import math
+
+
+class Scenario:
+ def __init__(self):
+ # physical parameters
+ self.drop_diameter = 20
+ self.pool_depth = self.drop_diameter * 0.5
+
+ self.bond_number = 3.18
+ self.weber_number = 2010
+ self.ohnesorge_number = 0.0384
+
+ self.relaxation_rate_heavy = 1.988
+ self.density_heavy = 1
+ self.density_ratio = 1000
+ self.dynamic_viscosity_ratio = 100
+
+ self.impact_angle_degree = 0 # drop impact angle in degree
+
+ self.interface_thickness = 4
+ self.mobility = 0.03
+
+ self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy
+ self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5)
+ self.dynamic_viscosity_heavy = self.kinematic_viscosity_heavy * self.density_heavy
+
+ self.density_light = self.density_heavy / self.density_ratio
+ self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio
+ self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light
+ self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5
+
+ self.surface_tension = ((self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter
+ / self.density_heavy)
+
+ self.gravitational_acceleration = (- self.bond_number * self.surface_tension / self.drop_diameter**2 /
+ self.density_heavy)
+
+ self.impact_velocity_magnitude = (self.weber_number * self.surface_tension / self.drop_diameter /
+ self.density_heavy)**0.5
+
+ self.impact_velocity = (self.impact_velocity_magnitude * math.sin(self.impact_angle_degree * math.pi / 180),
+ -self.impact_velocity_magnitude * math.cos(self.impact_angle_degree * math.pi / 180),
+ 0)
+
+ self.reference_time = self.drop_diameter / self.impact_velocity_magnitude
+ self.timesteps = 15 * self.reference_time
+ self.vtkWriteFrequency = self.reference_time
+ self.meshWriteFrequency = self.reference_time
+
+ # domain parameters
+ self.size = (self.drop_diameter * 10,
+ self.drop_diameter * 5,
+ self.drop_diameter * 10)
+ self.blocks = (1, 1, 1)
+ self.periodic = (1, 0, 1)
+ self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2])
+ self.drop_mid_point = (self.size[0] / 2, self.pool_depth + self.drop_diameter / 2, self.size[2] / 2)
+
+ self.scenario = 3 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
+
+ self.counter = 0
+ self.yPositions = []
+
+ @wlb.member_callback
+ def config(self):
+ return {
+ 'DomainSetup': {
+ 'blocks': self.blocks,
+ 'domainSize': self.size,
+ 'cellsPerBlock': self.cells,
+ 'periodic': self.periodic,
+ },
+ 'Parameters': {
+ 'timesteps': self.timesteps,
+ 'vtkWriteFrequency': self.vtkWriteFrequency,
+ 'meshWriteFrequency': self.meshWriteFrequency,
+ 'remainingTimeLoggerFrequency': 10.0,
+ 'scenario': self.scenario,
+ },
+ 'PhysicalParameters': {
+ 'density_liquid': self.density_heavy,
+ 'density_gas': self.density_light,
+ 'surface_tension': self.surface_tension,
+ 'mobility': self.mobility,
+ 'gravitational_acceleration': self.gravitational_acceleration,
+ 'relaxation_time_liquid': self.relaxation_time_heavy - 0.5,
+ 'relaxation_time_gas': self.relaxation_time_light - 0.5,
+ 'interface_thickness': self.interface_thickness,
+ },
+ 'Boundaries': {
+ 'Border': [
+ {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
+ ]
+ },
+ 'Drop': {
+ 'drop_mid_point': self.drop_mid_point,
+ 'drop_radius': self.drop_diameter / 2,
+ 'pool_depth': self.pool_depth,
+ 'impact_velocity': self.impact_velocity,
+ },
+ }
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
index 5caf4b343..213dd9b1b 100755
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
@@ -50,7 +50,7 @@ class Scenario:
# everything else
self.dbFile = "risingBubble3D.db"
- self.scenario = 1 # 1 rising bubble, 2 RTI
+ self.scenario = 1 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
self.counter = 0
self.yPositions = []
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py
new file mode 100755
index 000000000..f07adfe57
--- /dev/null
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py
@@ -0,0 +1,184 @@
+import os
+import waLBerla as wlb
+import numpy as np
+import pandas as pd
+
+from waLBerla.core_extension import makeSlice
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+
+class Scenario:
+ def __init__(self):
+ # physical parameters
+ self.tube_diameter = 64
+
+ self.bond_number = 100
+ self.morton_number = 0.015
+
+ self.relaxation_rate_heavy = 1.76
+ self.density_heavy = 1
+ self.density_ratio = 744
+ self.dynamic_viscosity_ratio = 4236
+
+ self.interface_width = 3
+ self.mobility = 0.08
+
+ self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy
+ self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5)
+ self.dynamic_viscosity_heavy = self.density_heavy * self.kinematic_viscosity_heavy
+
+ self.density_light = self.density_heavy / self.density_ratio
+ self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio
+ self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light
+ self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5
+
+ self.surface_tension = (self.bond_number * self.dynamic_viscosity_heavy**4 / self.morton_number /
+ self.density_heavy**2 / self.tube_diameter**2)**0.5
+
+ self.gravitational_acceleration = - (self.morton_number * self.density_heavy * self.surface_tension**3 /
+ self.dynamic_viscosity_heavy**4)
+
+ self.reference_time = (self.tube_diameter / abs(self.gravitational_acceleration))**0.5
+ self.timesteps = 20 * self.reference_time
+
+ self.vtkWriteFrequency = self.reference_time
+ self.dbWriteFrequency = self.reference_time
+ self.meshWriteFrequency = self.reference_time
+
+ # simulation parameters
+ self.size = (self.tube_diameter, self.tube_diameter * 10, self.tube_diameter)
+ self.blocks = (1, 1, 1)
+ self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2])
+ self.periodic = (0, 0, 0)
+ self.inner_radius = self.tube_diameter
+
+ self.center_x = self.size[0] / 2
+ self.center_y = self.size[1] / 2
+ self.center_z = self.size[2] / 2
+
+ self.scenario = 4 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
+
+ self.counter = 0
+ self.yPositions = []
+
+ self.eccentricity_or_pipe_ratio = False # if True eccentricity is conducted otherwise pipe ratio
+ self.ratio = 0
+
+ self.start_transition = (self.size[1] // 2) - 2 * self.tube_diameter
+ self.length_transition = 4 * self.tube_diameter
+
+ setup = "eccentricity" if self.eccentricity_or_pipe_ratio else "ratio"
+
+ self.csv_file = f"Taylor_bubble_D_{self.tube_diameter}_C_{setup}_{self.ratio}_W_" \
+ f"{self.interface_width}_M_{self.mobility}.csv"
+
+ d = self.tube_diameter / 2
+ dh = self.tube_diameter - d
+
+ resolution = self.tube_diameter / 128
+
+ self.Donut_D = 0.1 * self.tube_diameter / resolution
+ self.Donut_h = dh / 6
+ self.DonutTime = 0.5 * (self.tube_diameter + d) / 2
+
+ self.config_dict = self.config()
+
+ @wlb.member_callback
+ def config(self):
+ return {
+ 'DomainSetup': {
+ 'blocks': self.blocks,
+ 'domainSize': self.size,
+ 'cellsPerBlock': self.cells,
+ 'diameter': self.tube_diameter,
+ 'periodic': self.periodic,
+ 'inner_radius': self.inner_radius,
+ 'ratio': self.ratio,
+ 'start_transition': self.start_transition,
+ 'length_transition': self.length_transition,
+ 'eccentricity_or_pipe_ration': self.eccentricity_or_pipe_ratio,
+ 'tube': True
+ },
+ 'Parameters': {
+ 'timesteps': self.timesteps,
+ 'vtkWriteFrequency': self.vtkWriteFrequency,
+ 'dbWriteFrequency': self.dbWriteFrequency,
+ 'meshWriteFrequency': self.meshWriteFrequency,
+ 'remainingTimeLoggerFrequency': 60.0,
+ 'scenario': self.scenario,
+ },
+ 'PhysicalParameters': {
+ 'density_liquid': self.density_heavy,
+ 'density_gas': self.density_light,
+ 'surface_tension': self.surface_tension,
+ 'mobility': self.mobility,
+ 'gravitational_acceleration': self.gravitational_acceleration,
+ 'relaxation_time_liquid': self.relaxation_time_heavy - 0.5,
+ 'relaxation_time_gas': self.relaxation_time_light - 0.5,
+ 'interface_thickness': self.interface_width
+ },
+ 'Boundaries': {
+ 'Border': [
+ {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'W', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'E', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'},
+ {'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'},
+ ],
+ },
+ 'TaylorBubble': {
+ 'bubble_radius': 0.75 * 0.5 * self.tube_diameter,
+ 'initial_height': 1 * self.tube_diameter,
+ 'length': 3 * self.tube_diameter,
+ }
+ }
+
+ @wlb.member_callback
+ def at_end_of_time_step(self, blocks, **kwargs):
+ t = kwargs["timeStep"]
+ if t % self.dbWriteFrequency == 0:
+ wlb_field = wlb.field.gather(blocks, 'phase', makeSlice[:, :, self.size[2] // 2])
+ if wlb_field:
+ data = {'timestep': t}
+ data.update(self.config_dict['Parameters'])
+ data.update(self.config_dict['DomainSetup'])
+ data.update(self.config_dict['PhysicalParameters'])
+ data.update(self.config_dict['TaylorBubble'])
+ data.update(kwargs)
+
+ phase_field = np.asarray(wlb_field).squeeze()
+ location_of_gas = np.where(phase_field < 0.5)
+ center_of_mass = np.mean(location_of_gas, axis=1)
+
+ assert np.isfinite(np.sum(phase_field)), "NaN detected in the phasefield"
+
+ self.yPositions.append(center_of_mass[1])
+ if len(self.yPositions) > 1:
+ speed = self.yPositions[-1] - self.yPositions[-2]
+ else:
+ speed = 0
+
+ data['center_of_mass_x'] = center_of_mass[0]
+ data['center_of_mass_y'] = center_of_mass[1]
+ # data['center_of_mass_z'] = center_of_mass[2]
+
+ data['xCells'] = self.size[0]
+ data['yCells'] = self.size[1]
+ data['zCells'] = self.size[2]
+
+ data['rising_velocity'] = speed
+
+ data['StencilHydroLB'] = kwargs["stencil_hydro"]
+ data['StencilPhaseLB'] = kwargs["stencil_phase"]
+ sequenceValuesToScalars(data)
+
+ df = pd.DataFrame.from_records([data])
+ if not os.path.isfile(self.csv_file):
+ df.to_csv(self.csv_file, index=False)
+ else:
+ df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
index ac4feda33..9116f0b19 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
@@ -18,4 +18,4 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenGPU
waLBerla_add_executable(NAME multiphaseGPU
FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp util.cpp multiphase_codegen.py
- DEPENDS blockforest core cuda field postprocessing python_coupling lbm geometry timeloop gui PhaseFieldCodeGenGPU)
+ DEPENDS blockforest core cuda field postprocessing python_coupling lbm geometry timeloop PhaseFieldCodeGenGPU)
--
GitLab