diff --git a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
index bbdf8c144bd730dc0ec68c650ae56220f32cd4cf..f9a20e01a0ddcd11a65eade491cf05a35f25cca4 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
+++ b/apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
@@ -86,7 +86,6 @@ using FlagField_T = FlagField< flag_t >;
#if defined(WALBERLA_BUILD_WITH_CUDA)
typedef cuda::GPUField< real_t > GPUField;
#endif
-// using CommScheme_T = cuda::communication::UniformGPUScheme<stencil::D2Q9>;
int main(int argc, char** argv)
{
@@ -185,7 +184,7 @@ int main(int argc, char** argv)
auto Comm_velocity_based_distributions =
make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
auto generatedPackInfo_velocity_based_distributions =
- make_shared< pystencils::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
+ make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
Comm_velocity_based_distributions->addPackInfo(generatedPackInfo_velocity_based_distributions);
auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field_gpu);
Comm_velocity_based_distributions->addPackInfo(generatedPackInfo_phase_field);
@@ -193,7 +192,7 @@ int main(int argc, char** argv)
auto Comm_phase_field_distributions =
make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
auto generatedPackInfo_phase_field_distributions =
- make_shared< pystencils::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
+ make_shared< lbm::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
Comm_phase_field_distributions->addPackInfo(generatedPackInfo_phase_field_distributions);
#else
@@ -202,14 +201,14 @@ int main(int argc, char** argv)
auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field);
auto generatedPackInfo_velocity_based_distributions =
- make_shared< pystencils::PackInfo_velocity_based_distributions >(lb_velocity_field);
+ make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field);
Comm_velocity_based_distributions.addPackInfo(generatedPackInfo_phase_field);
Comm_velocity_based_distributions.addPackInfo(generatedPackInfo_velocity_based_distributions);
blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_phase_field_distributions(blocks);
auto generatedPackInfo_phase_field_distributions =
- make_shared< pystencils::PackInfo_phase_field_distributions >(lb_phase_field);
+ make_shared< lbm::PackInfo_phase_field_distributions >(lb_phase_field);
Comm_phase_field_distributions.addPackInfo(generatedPackInfo_phase_field_distributions);
#endif
diff --git a/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py b/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
index a751f35a43bab045b19e3ec94253f6930aacb839..bbd994de58a0b180293b7e643a2a27e0b6c50068 100644
--- a/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
+++ b/apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
@@ -5,11 +5,12 @@ from pystencils import AssignmentCollection
from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
from lbmpy.stencils import get_stencil
-from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field
+from lbmpy_walberla import generate_lb_pack_info
from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
initializer_kernel_hydro_lb, interface_tracking_force, \
- hydrodynamic_force, get_collision_assignments_hydro
+ hydrodynamic_force, get_collision_assignments_hydro, get_collision_assignments_phase
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
@@ -52,6 +53,7 @@ w_c = 1.0 / (0.5 + (3.0 * M))
u = fields(f"vel_field({dimensions}): [{dimensions}D]", layout='fzyx')
# phase-field
C = fields(f"phase_field: [{dimensions}D]", layout='fzyx')
+C_tmp = fields(f"phase_field_tmp: [{dimensions}D]", layout='fzyx')
# phase-field distribution functions
h = fields(f"lb_phase_field({q_phase}): [{dimensions}D]", layout='fzyx')
@@ -88,32 +90,26 @@ h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
+force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W, fd_stencil=get_stencil("D3Q27"))]
force_model_h = MultiphaseForceModel(force=force_h)
-force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
+force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force,
+ fd_stencil=get_stencil("D3Q27"))
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
+force_model_g = MultiphaseForceModel(force=force_g, rho=density)
####################
# LBM UPDATE RULES #
####################
-method_phase.set_force_model(force_model_h)
+phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
+ velocity_input=u,
+ output={'density': C_tmp},
+ force_model=force_model_h,
+ symbolic_fields={"symbolic_field": h,
+ "symbolic_temporary_field": h_tmp},
+ kernel_type='stream_pull_collide')
-phase_field_LB_step = create_lb_update_rule(lb_method=method_phase,
- velocity_input=u,
- compressible=True,
- optimization={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
-
-
-phase_field_LB_step.set_main_assignments_from_dict({**phase_field_LB_step.main_assignments_dict, **{C.center: sum_h}})
-
-phase_field_LB_step = AssignmentCollection(main_assignments=phase_field_LB_step.main_assignments,
- subexpressions=phase_field_LB_step.subexpressions)
phase_field_LB_step = sympy_cse(phase_field_LB_step)
# ---------------------------------------------------------------------------------------------------------
@@ -121,18 +117,12 @@ phase_field_LB_step = sympy_cse(phase_field_LB_step)
hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
- force=force_g,
- sub_iterations=1,
+ force_model=force_model_g,
+ sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
kernel_type='collide_stream_push')
-# streaming of the hydrodynamic distribution
-stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
- optimization={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='stream_pull_only')
-
###################
# GENERATE SWEEPS #
###################
@@ -161,7 +151,7 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates)
generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
- field_swaps=[(h, h_tmp)],
+ field_swaps=[(h, h_tmp), (C, C_tmp)],
inner_outer_split=True,
cpu_vectorize_info=cpu_vec)
@@ -171,12 +161,13 @@ with CodeGeneration() as ctx:
cpu_vectorize_info=cpu_vec)
# communication
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field_distributions',
- phase_field_LB_step.main_assignments, target='cpu')
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field',
- hydro_LB_step.all_assignments, target='cpu', kind='pull')
- generate_pack_info_from_kernel(ctx, 'PackInfo_velocity_based_distributions',
- hydro_LB_step.all_assignments, target='cpu', kind='push')
+ generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
+ streaming_pattern='pull', target='cpu')
+
+ generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
+ streaming_pattern='push', target='cpu')
+
+ generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='cpu')
ctx.write_file("GenDefines.h", info_header)
@@ -187,7 +178,7 @@ with CodeGeneration() as ctx:
g_updates, target='gpu')
generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
- field_swaps=[(h, h_tmp)],
+ field_swaps=[(h, h_tmp), (C, C_tmp)],
inner_outer_split=True,
target='gpu',
gpu_indexing_params=sweep_params,
@@ -200,12 +191,13 @@ with CodeGeneration() as ctx:
gpu_indexing_params=sweep_params,
varying_parameters=vp)
# communication
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field_distributions',
- phase_field_LB_step.main_assignments, target='gpu')
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field',
- hydro_LB_step.all_assignments, target='gpu', kind='pull')
- generate_pack_info_from_kernel(ctx, 'PackInfo_velocity_based_distributions',
- hydro_LB_step.all_assignments, target='gpu', kind='push')
+ generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
+ streaming_pattern='pull', target='gpu')
+
+ generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
+ streaming_pattern='push', target='gpu')
+
+ generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='gpu')
ctx.write_file("GenDefines.h", info_header)
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
index adacc326397db495352f958afa90fb26ae63e5d4..cf550a8e1f8392dbce80bb709ac71c3857724787 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
@@ -10,11 +10,14 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenCPU
phase_field_LB_NoSlip.cpp phase_field_LB_NoSlip.h
hydro_LB_step.cpp hydro_LB_step.h
hydro_LB_NoSlip.cpp hydro_LB_NoSlip.h
- stream_hydro.cpp stream_hydro.h
+ PackInfo_phase_field_distributions.cpp PackInfo_phase_field_distributions.h
+ PackInfo_velocity_based_distributions.cpp PackInfo_velocity_based_distributions.h
+ PackInfo_phase_field.cpp PackInfo_phase_field.h
+ ContactAngle.cpp ContactAngle.h
GenDefines.h)
waLBerla_add_executable(NAME multiphaseCPU
- FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp contact.cpp CalculateNormals.cpp multiphase_codegen.py
+ FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp multiphase_codegen.py
DEPENDS blockforest core field postprocessing lbm geometry timeloop gui PhaseFieldCodeGenCPU)
set_target_properties(multiphaseCPU PROPERTIES CXX_VISIBILITY_PRESET hidden)
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.cpp
deleted file mode 100644
index e20053e38130262fdf2bb8c556f1fd665888a945..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.cpp
+++ /dev/null
@@ -1,83 +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 CalculateNormals.cpp
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-#include "CalculateNormals.h"
-
-#include "core/Environment.h"
-#include "core/logging/Initialization.h"
-
-#include "field/FlagField.h"
-
-namespace walberla
-{
-using FlagField_T = FlagField< uint8_t >;
-using NormalsField_T = GhostLayerField< int8_t, 3 >;
-
-void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
- ConstBlockDataID flagFieldID, FlagUID domainFlagUID, FlagUID boundaryFlagUID)
-{
- for (auto& block : *blocks)
- {
- CellInterval globalCellBB = blocks->getBlockCellBB(block);
- CellInterval blockLocalCellBB;
- blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, block, globalCellBB);
-
- auto* normalsField = block.getData< NormalsField_T >(normalsFieldID);
- auto* flagField = block.getData< FlagField_T >(flagFieldID);
- auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
- auto domainFlag = flagField->getFlag(domainFlagUID);
-
- // clang-format off
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(normalsField,
-
- if( x < blockLocalCellBB.xMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x + 1, y, z) == domainFlag)
- normalsField->get(x, y, z, 0) = 1;
- }
-
- if( x > blockLocalCellBB.xMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x - 1, y, z) == domainFlag)
- normalsField->get(x, y, z, 0) = - 1;
- }
-
- if( y < blockLocalCellBB.yMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y + 1, z) == domainFlag)
- normalsField->get(x, y, z, 1) = 1;
- }
-
- if( y > blockLocalCellBB.yMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y - 1, z) == domainFlag)
- normalsField->get(x, y, z, 1) = - 1;
- }
-
- if( z < blockLocalCellBB.zMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z + 1) == domainFlag)
- normalsField->get(x, y, z, 2) = 1;
- }
-
- if( z > blockLocalCellBB.zMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z - 1) == domainFlag)
- normalsField->get(x, y, z, 2) = - 1;
- }
-
- )
- // clang-format on
- }
-}
-} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
index 0fd64cb73cefcfefb25a565fee9c80b4c127b155..979c87ff330f669567df87c2e051a75271280b6b 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
@@ -53,24 +53,196 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B
}
}
+void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+ const real_t D = 5, const real_t H = 2, const real_t DT = 20, const real_t Donut_x0 = 40)
+{
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ 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 Ri = D * sqrt(pow(H, 2) - pow(DT - sqrt(pow(globalCell[0] - Mx, 2) + pow(globalCell[2] - Mz, 2)), 2));
+
+ real_t shifter = atan2((globalCell[2] - Mz), (globalCell[0] - Mx));
+ if (shifter < 0) shifter = shifter + 2 * math::pi;
+ if ((globalCell[1] < Donut_x0 + Ri * sin(shifter / 2.0)) && (globalCell[1] > Donut_x0 - Ri)) {
+ phaseField->get(x, y, z) = 0.0;
+ } else { phaseField->get(x, y, z) = 1.0; })
+ }
+}
+
+void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+ real_t W = 5)
+{
+ Vector3< real_t > bubbleMidPoint;
+
+ auto X = blocks->getDomainCellBB().xMax();
+ auto Y = blocks->getDomainCellBB().yMax();
+ auto Z = blocks->getDomainCellBB().zMax();
+
+ // 20 percent from the top are filled with the gas phase
+ real_t gas_top = Y - Y / 5.0;
+
+ // Diameter of the bubble
+ real_t D = R * 2;
+
+ // distance in between the bubbles
+ int dist = 4;
+ auto nx = static_cast< unsigned int >(floor(X / (D + dist * W)));
+ auto nz = static_cast< unsigned int >(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));
+ std::vector< std::vector< real_t > > fluctuation_pos(nx, std::vector< real_t >(nz, 0.0));
+
+ real_t max_fluctuation_radius = R / 5;
+ real_t max_fluctuation_pos = (dist * W) / 3.0;
+ for (unsigned int i = 0; i < nx; ++i)
+ {
+ for (unsigned int j = 0; j < nz; ++j)
+ {
+ fluctuation_radius[i][j] = math::realRandom< real_t >(-max_fluctuation_radius, max_fluctuation_radius);
+ fluctuation_pos[i][j] = math::realRandom< real_t >(-max_fluctuation_pos, max_fluctuation_pos);
+ }
+ }
+
+ 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) {
+ for (unsigned int j = 0; j < nz; ++j)
+ {
+ bubbleMidPoint[0] = (i + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[1] = R + W + 4;
+ bubbleMidPoint[2] = (j + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+
+ real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
+ (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
+ (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
+ if (globalCell[0] >= i * (D + dist * W) && globalCell[0] <= (i + 1) * (D + dist * W) &&
+ globalCell[2] >= j * (D + dist * W) && globalCell[2] <= (j + 1) * (D + dist * W))
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - (R - fluctuation_radius[i][j])) / W);
+
+ if (globalCell[0] > nx * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ if (globalCell[2] > nz * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ }
+ }
+
+ if (globalCell[1] > gas_top) {
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (gas_top + 10 - globalCell[1]) / W);
+ })
+ }
+}
+
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t W = 5)
+ const real_t W = 5, const bool pipe = true)
{
auto X = blocks->getDomainCellBB().xMax();
+ auto Z = blocks->getDomainCellBB().zMax();
auto halfY = (blocks->getDomainCellBB().yMax()) / 2.0;
- double perturbation = 0.05;
+ real_t perturbation = 0.05;
- for (auto& block : *blocks)
+ if (pipe)
{
- auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- // clang-format off
- 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((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));
- )
- // clang-format on
+ 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((globalCell[0] - X / 2) * (globalCell[0] - X / 2) +
+ (globalCell[2] - Z / 2) * (globalCell[2] - Z / 2));
+ if (R > X) R = X; real_t tmp = perturbation * X * cos((2.0 * math::pi * R) / X);
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) + tmp) / (W / 2.0));)
+ }
+ }
+ else
+ {
+ 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 tmp = perturbation * X *
+ (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));)
+ }
+ }
+}
+
+void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID flagFieldID,
+ field::FlagUID boundaryFlagUID, real_t const R_in, real_t const eccentricity,
+ real_t const start_transition, real_t const length_transition,
+ bool const eccentricity_or_pipe_ratio)
+{
+ if (eccentricity_or_pipe_ratio)
+ {
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+
+ real_t const shift = eccentricity * Mx / 2;
+
+ for (auto& block : *blocks)
+ {
+ 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));
+ real_t R1;
+ if (globalCell[1] <= start_transition) {
+ R1 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
+ real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
+ real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ R1 = sqrt((globalCell[0] - Mx - shift_tmp) * (globalCell[0] - Mx - shift_tmp) +
+ (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ } else {
+ R1 = sqrt((globalCell[0] - Mx - shift) * (globalCell[0] - Mx - shift) +
+ (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ }
+
+ real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (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);)
+ }
+ }
+ else
+ {
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+
+ real_t const shift = eccentricity * R_in;
+ real_t R_tmp;
+
+ for (auto& block : *blocks)
+ {
+ 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 (globalCell[1] <= start_transition) {
+ R_tmp = R_in;
+ } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
+ real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
+ real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ R_tmp = R_in + shift_tmp;
+ } else {
+ R_tmp = R_in + shift;
+ }
+
+ real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (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);)
+ }
}
}
} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
index 04bac53f2e75d17ddd393da0e247b63df05d0704..585639ac9078beb048a8ce694ea97552e59b833a 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
@@ -34,6 +34,17 @@ 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_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5);
+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_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+ real_t W = 5);
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5,
+ const bool pipe = true);
+
+void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID flagFieldID,
+ field::FlagUID boundaryFlagUID, real_t R_in, real_t eccentricity, real_t start_transition,
+ real_t length_transition, bool const eccentricity_or_pipe_ratio);
} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/contact.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/contact.cpp
deleted file mode 100644
index 8253e1d4cd26c071e4c1abd775d8aee82e27ce26..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/contact.cpp
+++ /dev/null
@@ -1,116 +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 contact.cpp
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-
-#include "contact.h"
-
-#include "core/DataTypes.h"
-
-#include <cmath>
-
-#define FUNC_PREFIX
-
-namespace walberla
-{
-namespace lbm
-{
-#ifdef __GNUC__
-# pragma GCC diagnostic push
-#endif
-
-namespace internal_boundary_contact
-{
-static FUNC_PREFIX void contact_angle_treatment(uint8_t* WALBERLA_RESTRICT const _data_indexVector, double* WALBERLA_RESTRICT _data_phase,
- int64_t const _stride_phase_0, int64_t const _stride_phase_1,
- int64_t const _stride_phase_2, int64_t indexVectorSize, double alpha)
-{
- for (int ctr_0 = 0; ctr_0 < indexVectorSize; ctr_0 += 1)
- {
- const int32_t x = *((int32_t*) (&_data_indexVector[24 * ctr_0]));
- const int32_t y = *((int32_t*) (&_data_indexVector[24 * ctr_0 + 4]));
- const int32_t z = *((int32_t*) (&_data_indexVector[24 * ctr_0 + 8]));
-
- const int32_t nx = *((int32_t*) (&_data_indexVector[24 * ctr_0 + 12]));
- const int32_t ny = *((int32_t*) (&_data_indexVector[24 * ctr_0 + 16]));
- const int32_t nz = *((int32_t*) (&_data_indexVector[24 * ctr_0 + 20]));
-
- const int32_t x1 = x + nx;
- const int32_t y1 = y + ny;
- const int32_t z1 = z + nz;
-
- const double h = 0.5 * sqrt((float) (nx * nx + ny * ny + nz * nz));
- const double a = cos(alpha);
- const double W = 5;
-
- double* WALBERLA_RESTRICT _phase_wall = _data_phase + _stride_phase_1 * y + _stride_phase_2 * z;
- double* WALBERLA_RESTRICT _phase_interior = _data_phase + _stride_phase_1 * y1 + _stride_phase_2 * z1;
- if (h < 0.001) { _phase_wall[_stride_phase_0 * x] = 1.0; }
- else if (a > 1e-8 || a < -1e-8)
- {
- const double var = -h * (4.0 / W) * a;
- _phase_wall[_stride_phase_0 * x] =
- (1 + var - sqrt((1 + var) * (1 + var) - 4 * var * _phase_interior[_stride_phase_0 * x1])) / (var + 1e-12) -
- _phase_interior[_stride_phase_0 * x1];
- }
- else
- {
- _phase_wall[_stride_phase_0 * x] = _phase_interior[_stride_phase_0 * x1];
- }
- }
-}
-} // namespace internal_boundary_contact
-
-#ifdef __GNUC__
-# pragma GCC diagnostic pop
-#endif
-
-#ifdef __CUDACC__
-# pragma pop
-#endif
-
-void contact::run(IBlock* block, IndexVectors::Type type)
-{
- auto* indexVectors = block->getData< IndexVectors >(indexVectorID);
- int64_t indexVectorSize = int64_c(indexVectors->indexVector(type).size());
- if (indexVectorSize == 0) return;
-
- auto pointer = indexVectors->pointerCpu(type);
-
- auto* _data_indexVector = reinterpret_cast< uint8_t* >(pointer);
-
- auto phaseField = block->getData< field::GhostLayerField< double, 1 > >(phaseFieldID);
- auto& alpha = this->alpha_;
-
- WALBERLA_ASSERT_GREATER_EQUAL(0, -int_c(phaseField->nrOfGhostLayers()))
- double* WALBERLA_RESTRICT _data_phase = phaseField->dataAt(0, 0, 0, 0);
- const auto _stride_pdfs_0 = int64_t(phaseField->xStride());
- const auto _stride_pdfs_1 = int64_t(phaseField->yStride());
- const auto _stride_pdfs_2 = int64_t(phaseField->zStride());
- internal_boundary_contact::contact_angle_treatment(_data_indexVector, _data_phase, _stride_pdfs_0, _stride_pdfs_1,
- _stride_pdfs_2, indexVectorSize, alpha);
-}
-
-void contact::operator()(IBlock* block) { run(block, IndexVectors::ALL); }
-
-void contact::inner(IBlock* block) { run(block, IndexVectors::INNER); }
-
-void contact::outer(IBlock* block) { run(block, IndexVectors::OUTER); }
-
-} // namespace lbm
-} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/contact.h b/apps/showcases/PhaseFieldAllenCahn/CPU/contact.h
deleted file mode 100644
index dbaf4a560b3a717c96b34ae811a58f048776eff8..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/contact.h
+++ /dev/null
@@ -1,132 +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 contact.h
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-
-#include "blockforest/StructuredBlockForest.h"
-
-#include "core/DataTypes.h"
-
-#include "domain_decomposition/BlockDataID.h"
-#include "domain_decomposition/IBlock.h"
-
-#include "field/FlagField.h"
-#include "field/GhostLayerField.h"
-
-#include <set>
-#include <vector>
-
-namespace walberla
-{
-namespace lbm
-{
-class contact
-{
- public:
- struct IndexInfo
- {
- int32_t x1;
- int32_t y1;
- int32_t z1;
- int32_t x2;
- int32_t y2;
- int32_t z2;
- IndexInfo(int32_t x1_, int32_t y1_, int32_t z1_, int32_t x2_, int32_t y2_, int32_t z2_)
- : x1(x1_), y1(y1_), z1(z1_), x2(x2_), y2(y2_), z2(z2_)
- {}
- bool operator==(const IndexInfo& o) const
- {
- return x1 == o.x1 && y1 == o.y1 && z1 == o.z1 && x2 == o.x2 && y2 == o.y2 && z2 == o.z2;
- }
- };
-
- class IndexVectors
- {
- public:
- using CpuIndexVector = std::vector< IndexInfo >;
-
- enum Type { ALL = 0, INNER = 1, OUTER = 2, NUM_TYPES = 3 };
-
- IndexVectors() : cpuVectors_(NUM_TYPES) {}
- bool operator==(IndexVectors& other) { return other.cpuVectors_ == cpuVectors_; }
-
- CpuIndexVector& indexVector(Type t) { return cpuVectors_[t]; }
- IndexInfo* pointerCpu(Type t) { return &(cpuVectors_[t][0]); }
-
- private:
- std::vector< CpuIndexVector > cpuVectors_;
- };
-
- contact(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID phaseFieldID_, double alpha)
- : phaseFieldID(phaseFieldID_), alpha_(alpha)
- {
- auto createIdxVector = [](IBlock* const, StructuredBlockStorage* const) { return new IndexVectors(); };
- indexVectorID = blocks->addStructuredBlockData< IndexVectors >(createIdxVector, "IndexField_contact_angle");
- };
-
- void operator()(IBlock* block);
- void inner(IBlock* block);
- void outer(IBlock* block);
-
- template< typename NormalField_T >
- void fillFromNormalField(const shared_ptr< StructuredBlockForest >& blocks, ConstBlockDataID normalFieldID)
- {
- for (auto& block : *blocks)
- {
- auto* indexVectors = block.getData< IndexVectors >(indexVectorID);
- auto& indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
- auto& indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
- auto& indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
-
- auto* normalField = block.getData< NormalField_T >(normalFieldID);
-
- auto inner = normalField->xyzSize();
- inner.expand(cell_idx_t(-1));
-
- indexVectorAll.clear();
- indexVectorInner.clear();
- indexVectorOuter.clear();
- // clang-format off
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(normalField, Cell globalCell;
- blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- if(normalField->get(x, y, z, 0) != 0 || normalField->get(x, y, z, 1) != 0 || normalField->get(x, y, z, 2) != 0)
- {
- auto element = IndexInfo(x, y, z, normalField->get(x, y, z, 0), normalField->get(x, y, z, 1), normalField->get(x, y, z, 2));
- indexVectorAll.push_back( element );
- if( inner.contains( x, y, z ) )
- indexVectorInner.push_back( element );
- else
- indexVectorOuter.push_back( element );
- }
- )
- // clang-format on
- }
- }
-
- private:
- void run(IBlock* block, IndexVectors::Type type);
-
- BlockDataID indexVectorID;
-
- public:
- BlockDataID phaseFieldID;
- double alpha_;
-};
-
-} // namespace lbm
-} // namespace walberla
\ No newline at end of file
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/droplet_contact_angle.py b/apps/showcases/PhaseFieldAllenCahn/CPU/droplet_contact_angle.py
index 68229ece6fda45f0467a6277548ac23515054a8b..ee297075c1030f87b1d4b3cab4ee505f07dc381b 100755
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/droplet_contact_angle.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/droplet_contact_angle.py
@@ -7,7 +7,7 @@ class Scenario:
self.vtkWriteFrequency = 1000
# simulation parameters
- self.timesteps = 5001
+ self.timesteps = 10001
self.cells = (32, 64, 32)
self.blocks = (4, 1, 4)
self.periodic = (0, 0, 0)
@@ -18,25 +18,24 @@ class Scenario:
self.overlappingWidth = (8, 1, 1)
self.timeStepStrategy = 'normal'
- self.contactAngle = 22
-
# bubble parameters
self.dropletRadius = 24.0
self.dropletMidPoint = (64, 24, 64)
# everything else
- self.scenario = 1 # 1 rising bubble, 2 RTI
+ self.scenario = 1 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up
self.counter = 0
self.yPositions = []
@wlb.member_callback
- def config(self, **kwargs):
+ def config(self):
return {
'DomainSetup': {
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
+ 'tube': False
},
'Parameters': {
'timesteps': self.timesteps,
@@ -45,7 +44,6 @@ class Scenario:
'overlappingWidth': self.overlappingWidth,
'remainingTimeLoggerFrequency': 10.0,
'scenario': self.scenario,
- 'contactAngle': self.contactAngle
},
'PhysicalParameters': {
'density_liquid': 1.0,
@@ -55,6 +53,7 @@ class Scenario:
'gravitational_acceleration': 0.0,
'relaxation_time_liquid': 3 * 0.166,
'relaxation_time_gas': 3 * 0.0166,
+ 'interface_thickness': 5
},
'Boundaries': {
'Border': [
@@ -69,7 +68,7 @@ class Scenario:
'Bubble': {
'bubbleMidPoint': self.dropletMidPoint,
'bubbleRadius': self.dropletRadius,
- 'bubble': False
+ 'bubble': False # this means we are simulating a droplet rather than a bubble
},
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
index 00f5fae4d7caad2a9540aade589e05ea9aa06cc3..07c9542cf3c83de774e3cf39f8cfafed8b9ec597 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp
@@ -27,21 +27,17 @@
#include "field/AddToStorage.h"
#include "field/FlagField.h"
-#include "field/communication/PackInfo.h"
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
#include "python_coupling/CreateConfig.h"
#include "python_coupling/PythonCallback.h"
-#include "python_coupling/export/FieldExports.h"
#include "timeloop/SweepTimeloop.h"
-#include "CalculateNormals.h"
#include "InitializerFunctions.h"
#include "PythonExports.h"
-#include "contact.h"
//////////////////////////////
// INCLUDE GENERATED FILES //
@@ -54,7 +50,10 @@
#include "initialize_velocity_based_distributions.h"
#include "phase_field_LB_NoSlip.h"
#include "phase_field_LB_step.h"
-#include "stream_hydro.h"
+#include "ContactAngle.h"
+#include "PackInfo_phase_field_distributions.h"
+#include "PackInfo_velocity_based_distributions.h"
+#include "PackInfo_phase_field.h"
////////////
// USING //
@@ -62,11 +61,6 @@
using namespace walberla;
-using PdfField_phase_T = GhostLayerField< real_t, Stencil_phase_T::Size >;
-using PdfField_hydro_T = GhostLayerField< real_t, Stencil_hydro_T::Size >;
-using VelocityField_T = GhostLayerField< real_t, Stencil_hydro_T::Dimension >;
-using NormalsField_T = GhostLayerField< int8_t, Stencil_hydro_T::Dimension >;
-using PhaseField_T = GhostLayerField< real_t, 1 >;
using flag_t = walberla::uint8_t;
using FlagField_T = FlagField< flag_t >;
@@ -89,6 +83,7 @@ int main(int argc, char** argv)
auto domainSetup = config->getOneBlock("DomainSetup");
Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const bool tube = domainSetup.getParameter< bool >("tube", true);
////////////////////////////////////////
// ADD GENERAL SIMULATION PARAMETERS //
@@ -100,8 +95,7 @@ int main(int argc, char** argv)
const real_t remainingTimeLoggerFrequency =
parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
- const real_t alpha = parameters.getParameter< real_t >("contactAngle", real_c(90));
- const real_t alpha_rad = alpha * (math::pi / 180);
+
Vector3< int > overlappingWidth =
parameters.getParameter< Vector3< int > >("overlappingWidth", Vector3< int >(1, 1, 1));
@@ -116,9 +110,19 @@ int main(int argc, char** argv)
BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
- BlockDataID normals = field::addToStorage< NormalsField_T >(blocks, "normals", int8_t(0), field::fzyx);
BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+ auto physical_parameters = config->getOneBlock("PhysicalParameters");
+ const real_t density_liquid = physical_parameters.getParameter< real_t >("density_liquid", real_c(1.0));
+ const real_t density_gas = physical_parameters.getParameter< real_t >("density_gas");
+ const real_t surface_tension = physical_parameters.getParameter< real_t >("surface_tension");
+ const real_t mobility = physical_parameters.getParameter< real_t >("mobility");
+ const real_t gravitational_acceleration =
+ physical_parameters.getParameter< real_t >("gravitational_acceleration");
+ const real_t relaxation_time_liquid = physical_parameters.getParameter< real_t >("relaxation_time_liquid");
+ const real_t relaxation_time_gas = physical_parameters.getParameter< real_t >("relaxation_time_gas");
+ 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)
{
@@ -130,7 +134,7 @@ int main(int argc, char** argv)
}
else if (scenario == 2)
{
- initPhaseField_RTI(blocks, phase_field);
+ initPhaseField_RTI(blocks, phase_field, interface_thickness, tube);
}
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done")
@@ -138,43 +142,31 @@ int main(int argc, char** argv)
// ADD SWEEPS //
///////////////
- auto physical_parameters = config->getOneBlock("PhysicalParameters");
- const real_t density_liquid = physical_parameters.getParameter< real_t >("density_liquid", real_c(1.0));
- const real_t density_gas = physical_parameters.getParameter< real_t >("density_gas");
- const real_t surface_tension = physical_parameters.getParameter< real_t >("surface_tension");
- const real_t mobility = physical_parameters.getParameter< real_t >("mobility");
- const real_t gravitational_acceleration =
- physical_parameters.getParameter< real_t >("gravitational_acceleration");
- const real_t relaxation_time_liquid = physical_parameters.getParameter< real_t >("relaxation_time_liquid");
- const real_t relaxation_time_gas = physical_parameters.getParameter< real_t >("relaxation_time_gas");
-
- pystencils::initialize_phase_field_distributions init_h(lb_phase_field, phase_field, vel_field);
+ pystencils::initialize_phase_field_distributions init_h(lb_phase_field, phase_field, vel_field, interface_thickness);
pystencils::initialize_velocity_based_distributions init_g(lb_velocity_field, vel_field);
pystencils::phase_field_LB_step phase_field_LB_step(
- lb_phase_field, phase_field, vel_field, mobility,
+ lb_phase_field, phase_field, vel_field, interface_thickness, mobility,
Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
pystencils::hydro_LB_step hydro_LB_step(lb_velocity_field, phase_field, vel_field, gravitational_acceleration,
- density_liquid, density_gas, surface_tension, relaxation_time_liquid,
+ interface_thickness, density_liquid, density_gas, surface_tension, relaxation_time_liquid,
relaxation_time_gas,
Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
- pystencils::stream_hydro stream_hydro(lb_velocity_field,
- Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
////////////////////////
// ADD COMMUNICATION //
//////////////////////
-
blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_phase_field(blocks);
- Comm_phase_field.addPackInfo(make_shared< field::communication::PackInfo< PhaseField_T > >(phase_field));
+ auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field);
+ Comm_phase_field.addPackInfo(generatedPackInfo_phase_field);
blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_velocity_based_distributions(blocks);
- Comm_velocity_based_distributions.addPackInfo(
- make_shared< field::communication::PackInfo< PdfField_hydro_T > >(lb_velocity_field));
+ auto generatedPackInfo_velocity_based_distributions = make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field);
+ Comm_velocity_based_distributions.addPackInfo(generatedPackInfo_velocity_based_distributions);
blockforest::communication::UniformBufferedScheme< Stencil_hydro_T > Comm_phase_field_distributions(blocks);
- Comm_phase_field_distributions.addPackInfo(
- make_shared< field::communication::PackInfo< PdfField_phase_T > >(lb_phase_field));
+ auto generatedPackInfo_phase_field_distributions = make_shared< lbm::PackInfo_phase_field_distributions >(lb_phase_field);
+ Comm_phase_field_distributions.addPackInfo(generatedPackInfo_phase_field_distributions);
////////////////////////
// BOUNDARY HANDLING //
@@ -190,14 +182,13 @@ int main(int argc, char** argv)
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID);
}
- calculate_normals(blocks, normals, flagFieldID, fluidFlagUID, wallFlagUID);
lbm::phase_field_LB_NoSlip phase_field_LB_NoSlip(blocks, lb_phase_field);
lbm::hydro_LB_NoSlip hydro_LB_NoSlip(blocks, lb_velocity_field);
- lbm::contact contact_angle(blocks, phase_field, alpha_rad);
+ pystencils::ContactAngle contact_angle(blocks, phase_field, interface_thickness);
phase_field_LB_NoSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldID, FlagUID("NoSlip"), fluidFlagUID);
hydro_LB_NoSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldID, FlagUID("NoSlip"), fluidFlagUID);
- contact_angle.fillFromNormalField< NormalsField_T >(blocks, normals);
+ contact_angle.fillFromFlagField< FlagField_T >(blocks, flagFieldID, wallFlagUID, fluidFlagUID);
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the normals-field done")
@@ -209,48 +200,44 @@ int main(int argc, char** argv)
auto normalTimeStep = [&]() {
for (auto& block : *blocks)
{
- phase_field_LB_NoSlip(&block);
Comm_phase_field_distributions();
+ phase_field_LB_NoSlip(&block);
- phase_field_LB_step(&block);
- Comm_phase_field();
+ phase_field_LB_step(&block);
contact_angle(&block);
+ Comm_phase_field();
+
hydro_LB_step(&block);
- hydro_LB_NoSlip(&block);
- Comm_velocity_based_distributions();
- stream_hydro(&block);
+ Comm_velocity_based_distributions();
+ hydro_LB_NoSlip(&block);
}
};
auto simpleOverlapTimeStep = [&]() {
- for (auto& block : *blocks)
- phase_field_LB_NoSlip(&block);
-
- Comm_phase_field_distributions.startCommunication();
- for (auto& block : *blocks)
- phase_field_LB_step.inner(&block);
- Comm_phase_field_distributions.wait();
- for (auto& block : *blocks)
- phase_field_LB_step.outer(&block);
-
- Comm_phase_field.startCommunication();
- for (auto& block : *blocks)
- hydro_LB_step.inner(&block);
- Comm_phase_field.wait();
- for (auto& block : *blocks)
- hydro_LB_step.outer(&block);
-
- for (auto& block : *blocks)
- hydro_LB_NoSlip(&block);
-
- Comm_velocity_based_distributions.startCommunication();
- for (auto& block : *blocks)
- stream_hydro.inner(&block);
- Comm_velocity_based_distributions.wait();
- for (auto& block : *blocks)
- stream_hydro.outer(&block);
+ for (auto& block : *blocks)
+ phase_field_LB_NoSlip(&block);
+
+ Comm_phase_field_distributions.startCommunication();
+ for (auto& block : *blocks)
+ phase_field_LB_step.inner(&block);
+ Comm_phase_field_distributions.wait();
+ for (auto& block : *blocks)
+ phase_field_LB_step.outer(&block);
+
+ for (auto& block : *blocks)
+ contact_angle(&block);
+
+ Comm_phase_field.startCommunication();
+ for (auto& block : *blocks)
+ hydro_LB_step.inner(&block);
+ Comm_phase_field.wait();
+ for (auto& block : *blocks)
+ hydro_LB_step.outer(&block);
+
+ for (auto& block : *blocks)
+ hydro_LB_NoSlip(&block);
};
std::function< void() > timeStep;
if (timeStepStrategy == "overlap")
@@ -294,6 +281,8 @@ int main(int argc, char** argv)
{
callback.data().exposeValue("blocks", blocks);
callback.data().exposeValue( "timeStep", timeLoop->getCurrentTimeStep());
+ callback.data().exposeValue("stencil_phase", stencil_phase_name);
+ callback.data().exposeValue("stencil_hydro", stencil_hydro_name);
callback();
}
}
@@ -313,9 +302,6 @@ int main(int argc, char** argv)
auto phaseWriter = make_shared< field::VTKWriter< PhaseField_T > >(phase_field, "PhaseField");
vtkOutput->addCellDataWriter(phaseWriter);
- // auto normlasWriter = make_shared<field::VTKWriter<NormalsField_T>>(normals, "Normals");
- // vtkOutput->addCellDataWriter(normlasWriter);
-
// auto velWriter = make_shared<field::VTKWriter<VelocityField_T>>(vel_field, "Velocity");
// vtkOutput->addCellDataWriter(velWriter);
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
index 1ef86f585a2ad8e90ae0400321284431eaa71b04..401ce693ac893aba8375c64dd2b56f6853a9c08c 100755
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py
@@ -1,8 +1,11 @@
+import os
+
import waLBerla as wlb
-import waLBerla.tools.sqlitedb as wlbSqlite
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
from waLBerla.core_extension import makeSlice
import numpy as np
+import pandas as pd
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
@@ -10,12 +13,11 @@ class Scenario:
def __init__(self):
# output frequencies
self.vtkWriteFrequency = 1000
- self.dbWriteFrequency = 200
# simulation parameters
- self.timesteps = 27001
+ self.time = 2 # physical time in seconds
- self.cells = (64, 32, 64)
+ self.cells = (128, 64, 128)
self.blocks = (1, 8, 1)
self.periodic = (1, 0, 1)
self.size = (self.cells[0] * self.blocks[0],
@@ -24,25 +26,40 @@ class Scenario:
# physical parameters
self.density_heavy = 1.0
- self.reference_time = 6000
- self.parameters = calculate_parameters_rti(reference_length=128,
+ self.reference_time = 4000
+ self.dbWriteFrequency = self.reference_time // 20
+ self.timesteps = int(self.reference_time * self.time) + 1
+
+ self.capillary_number = 8.7
+ self.reynolds_number = 3000
+ self.atwood_number = 1
+ self.peclet_number = 744
+ self.density_ratio = 1000
+ self.viscosity_ratio = 100
+
+ self.parameters = calculate_parameters_rti(reference_length=self.cells[0],
reference_time=self.reference_time,
density_heavy=self.density_heavy,
- capillary_number=9.1,
- reynolds_number=128,
- atwood_number=1.0,
- peclet_number=140,
- density_ratio=3,
- viscosity_ratio=3)
+ capillary_number=self.capillary_number,
+ reynolds_number=self.reynolds_number,
+ atwood_number=self.atwood_number,
+ peclet_number=self.peclet_number,
+ density_ratio=self.density_ratio,
+ viscosity_ratio=self.viscosity_ratio)
+
+ self.interface_thickness = 5
+ self.tube = False
# everything else
- self.dbFile = "RTI.db"
+ self.dbFile = "RTI.csv"
- self.scenario = 2 # 1 rising bubble, 2 RTI
+ self.scenario = 2 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up
self.counter = 0
self.yPositions = []
+ self.config_dict = self.config()
+
@wlb.member_callback
def config(self):
return {
@@ -50,23 +67,24 @@ class Scenario:
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
+ 'tube': self.tube
},
'Parameters': {
'timesteps': self.timesteps,
'vtkWriteFrequency': self.vtkWriteFrequency,
'dbWriteFrequency': self.dbWriteFrequency,
- 'useGui': 0,
'remainingTimeLoggerFrequency': 10.0,
'scenario': self.scenario,
},
'PhysicalParameters': {
'density_liquid': self.density_heavy,
- 'density_gas': self.parameters["density_light"],
- 'surface_tension': self.parameters["surface_tension"],
- 'mobility': self.parameters.get("mobility", 0.1),
- 'gravitational_acceleration': self.parameters["gravitational_acceleration"],
+ 'density_gas': self.parameters.get("density_light"),
+ 'surface_tension': self.parameters.get("surface_tension"),
+ 'mobility': self.parameters.get("mobility"),
+ 'gravitational_acceleration': self.parameters.get("gravitational_acceleration"),
'relaxation_time_liquid': self.parameters.get("relaxation_time_heavy"),
'relaxation_time_gas': self.parameters.get("relaxation_time_light"),
+ 'interface_thickness': self.interface_thickness
},
'Boundaries': {
'Border': [
@@ -79,55 +97,79 @@ class Scenario:
@wlb.member_callback
def at_end_of_time_step(self, blocks, **kwargs):
t = kwargs['timeStep']
- ny = self.size[1]
- l0 = self.size[0]
if t % self.dbWriteFrequency == 0:
- location_of_spike = -100
- location_of_bubble = -100
- location_of_saddle = -100
- mass = -100
- spike_data = wlb.field.gather(blocks, 'phase', makeSlice[self.size[0] // 2, :, self.size[2] // 2])
- if spike_data:
- spike_field = np.asarray(spike_data).squeeze()
- location_of_spike = (np.argmax(spike_field > 0.5) - ny // 2) / l0
-
- bubble_data = wlb.field.gather(blocks, 'phase', makeSlice[0, :, 0])
- if bubble_data:
- bubble_field = np.asarray(bubble_data).squeeze()
- location_of_bubble = (np.argmax(bubble_field > 0.5) - ny // 2) / l0
-
- saddle_data = wlb.field.gather(blocks, 'phase', makeSlice[0, :, self.size[2] // 2])
- if saddle_data:
- saddle_field = np.asarray(saddle_data).squeeze()
- location_of_saddle = (np.argmax(saddle_field > 0.5) - ny // 2) / l0
-
phase = wlb.field.gather(blocks, 'phase', makeSlice[:, :, :])
if phase:
+ data = {'timestep': t}
+ data.update(self.config_dict['PhysicalParameters'])
+ data.update({'total_timesteps': self.timesteps})
+ data.update({'normalized_time': t / self.reference_time})
+ data.update({'tube_setup': self.tube})
+ data.update({'interface_thickness': self.interface_thickness})
+ data.update({'capillary_number': self.capillary_number})
+ data.update({'reynolds_number': self.reynolds_number})
+ data.update({'atwood_number': self.atwood_number})
+ data.update({'peclet_number': self.peclet_number})
+ data.update({'density_ratio': self.density_ratio})
+ data.update({'viscosity_ratio': self.viscosity_ratio})
+ data.update({'reference_time': self.reference_time})
+ data.update(kwargs)
+
phase_field = np.asarray(phase).squeeze()
+ stable = np.isfinite(np.sum(phase_field))
mass = np.sum(phase_field)
+ rel_max = np.max(phase_field) - 1
+ rel_min = np.min(phase_field)
+ data.update({'mass': mass})
+ data.update({'rel_max': rel_max})
+ data.update({'rel_min': rel_min})
+ data.update({'stable': stable})
+
+ if self.tube:
+ location_of_spike = self.get_interface_location(
+ phase_field[self.size[0] // 2, :, self.size[2] // 2])
+ a = np.where(phase_field < 0.5)
+ value = np.argmax(a[1])
+ location_of_bubble = self.get_interface_location(
+ phase_field[a[0][value], a[1][value] - 10:a[1][value] + 10, a[2][value]], a[1][value] - 10)
+
+ data.update({'location_of_spike': location_of_spike})
+ data.update({'location_of_bubble': location_of_bubble})
+ else:
+ location_of_spike = self.get_interface_location(
+ phase_field[self.size[0] // 2, :, self.size[2] // 2])
+ location_of_bubble = self.get_interface_location(phase_field[0, :, 0])
+ location_of_saddle = self.get_interface_location(phase_field[0, :, self.size[2] // 2])
+
+ data.update({'location_of_spike': location_of_spike})
+ data.update({'location_of_bubble': location_of_bubble})
+ data.update({'location_of_saddle': location_of_saddle})
+
+ sequenceValuesToScalars(data)
+
+ csv_file = f"RTI_{data['stencil_phase']}_{data['stencil_hydro']}_Re_{self.reynolds_number}_tube.csv"
+
+ df = pd.DataFrame.from_records([data])
+ if not os.path.isfile(csv_file):
+ df.to_csv(csv_file, index=False)
+ else:
+ df.to_csv(csv_file, index=False, mode='a', header=False)
+
+ def get_interface_location(self, one_dimensional_array, shift=None):
+ ny = self.size[1]
+ l0 = self.size[0]
- self.write_result_to_database(t, location_of_spike, location_of_bubble, location_of_saddle, mass)
-
- def write_result_to_database(self, t, spike, bubble, saddle, mass):
- """Writes the simulation result stored in the global variables shapeFactors and angles to
- an sqlite3 database, and resets the global variables."""
- result = {'waLBerlaVersion': wlb.build_info.version,
- 'xCells': self.size[0],
- 'yCells': self.size[1],
- 'zCells': self.size[2],
- 'spike': spike,
- 'bubble': bubble,
- 'saddle': saddle,
- 'mass': mass,
- 'timesteps': t,
- 'normalized_time': t / self.reference_time,
- 'processes': wlb.mpi.numProcesses(),
- }
- try:
- wlbSqlite.checkAndUpdateSchema(result, 'interface_location', self.dbFile)
- wlbSqlite.storeSingle(result, 'interface_location', self.dbFile)
- except Exception as e:
- wlb.log_warning("Failed to store run in database " + str(e) + "\n" + str(result))
+ index = np.argmax(one_dimensional_array > 0.5)
+
+ if index > 0:
+ zw1 = one_dimensional_array[index]
+ zw2 = one_dimensional_array[index - 1]
+ absolute_location = (index - 1) + (zw2 - 0.5) / (zw2 - zw1)
+ if shift:
+ absolute_location += shift
+ return (absolute_location - ny // 2) / l0
+ else:
+ return -100
scenarios = wlb.ScenarioManager()
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
index 28f40c5e4a68f477fcd2d8137fa270ca07f0f577..607fc161686c9bec98be1a103d6c70a47e6d0095 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py
@@ -1,26 +1,31 @@
+from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from pystencils import fields
from pystencils.simp import sympy_cse
-from pystencils import AssignmentCollection
from lbmpy.boundaries import NoSlip
-from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.creationfunctions import create_lb_method
from lbmpy.stencils import get_stencil
-from pystencils_walberla import CodeGeneration, generate_sweep
-from lbmpy_walberla import generate_boundary
+import pystencils_walberla
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field
+from lbmpy_walberla import generate_boundary, generate_lb_pack_info
from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb, \
- initializer_kernel_hydro_lb, interface_tracking_force, \
- hydrodynamic_force, get_collision_assignments_hydro
+ initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, get_collision_assignments_hydro, \
+ get_collision_assignments_phase
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
import numpy as np
import sympy as sp
-import pystencils as ps
-stencil_phase = get_stencil("D3Q19")
-stencil_hydro = get_stencil("D3Q27")
+stencil_phase_name = "D3Q27"
+stencil_hydro_name = "D3Q27"
+
+contact_angle_in_degrees = 22
+
+stencil_phase = get_stencil(stencil_phase_name)
+stencil_hydro = get_stencil(stencil_hydro_name)
q_phase = len(stencil_phase)
q_hydro = len(stencil_hydro)
@@ -45,7 +50,7 @@ relaxation_time_gas = sp.Symbol("tau_L")
# phase-field parameter
drho3 = (density_liquid - density_gas) / 3
# interface thickness
-W = 5
+W = sp.Symbol("interface_thickness")
# coefficients related to surface tension
beta = 12.0 * (surface_tension / W)
kappa = 1.5 * surface_tension * W
@@ -58,6 +63,7 @@ kappa = 1.5 * surface_tension * W
u = fields(f"vel_field({dimensions}): [{dimensions}D]", layout='fzyx')
# phase-field
C = fields(f"phase_field: [{dimensions}D]", layout='fzyx')
+C_tmp = fields(f"phase_field_tmp: [{dimensions}D]", layout='fzyx')
flag = fields(f"flag_field: uint8[{dimensions}D]", layout='fzyx')
# phase-field distribution functions
@@ -93,12 +99,17 @@ relaxation_rate_cutoff = sp.Piecewise((1 / (0.5 + relaxation_time_liquid), C.cen
# LBM METHODS #
###############
-method_phase = create_lb_method(stencil=stencil_phase, method='srt',
- relaxation_rate=relaxation_rate_allen_cahn, compressible=True)
+# method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
+# relaxation_rates=[1, 1.5, 1, 1.5, 1, 1.5])
+method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
+ relaxation_rates=[1, 1, 1, 1, 1, 1])
+
+method_phase.set_conserved_moments_relaxation_rate(relaxation_rate_allen_cahn)
method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1])
+
# create the kernels for the initialization of the g and h field
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W, fd_stencil=get_stencil("D3Q27"))
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
@@ -109,84 +120,84 @@ force_model_h = MultiphaseForceModel(force=force_h)
force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force,
fd_stencil=get_stencil("D3Q27"))
+force_model_g = MultiphaseForceModel(force=force_g, rho=density)
+
####################
# LBM UPDATE RULES #
####################
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
+phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
+ velocity_input=u,
+ output={'density': C_tmp},
+ force_model=force_model_h,
+ symbolic_fields={"symbolic_field": h,
+ "symbolic_temporary_field": h_tmp},
+ kernel_type='stream_pull_collide')
-method_phase.set_force_model(force_model_h)
-
-phase_field_LB_step = create_lb_update_rule(lb_method=method_phase,
- velocity_input=u,
- compressible=True,
- optimization={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
-
-phase_field_LB_step.set_main_assignments_from_dict({**phase_field_LB_step.main_assignments_dict, **{C.center: sum_h}})
-subexp = [ps.Assignment(a.lhs, float(a.rhs)) if a.rhs == 0 else a for a in phase_field_LB_step.subexpressions]
-phase_field_LB_step = AssignmentCollection(main_assignments=phase_field_LB_step.main_assignments,
- subexpressions=subexp)
phase_field_LB_step = sympy_cse(phase_field_LB_step)
-# ---------------------------------------------------------------------------------------------------------
hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
- force=force_g,
+ force_model=force_model_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
- kernel_type='collide_only')
+ kernel_type='collide_stream_push')
hydro_LB_step.set_sub_expressions_from_dict({**{relaxation_rate: relaxation_rate_cutoff},
**hydro_LB_step.subexpressions_dict})
-stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
- optimization={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='stream_pull_only')
+contact_angle = ContactAngle(contact_angle_in_degrees, W)
+
###################
# GENERATE SWEEPS #
###################
-cpu_vec = {'assume_inner_stride_one': True, 'nontemporal': True}
-
-vp = [('int32_t', 'cudaBlockSize0'),
- ('int32_t', 'cudaBlockSize1')]
-
info_header = f"""
+using namespace walberla;
#include "stencil/D3Q{q_phase}.h"\nusing Stencil_phase_T = walberla::stencil::D3Q{q_phase};
#include "stencil/D3Q{q_hydro}.h"\nusing Stencil_hydro_T = walberla::stencil::D3Q{q_hydro};
+using PdfField_phase_T = GhostLayerField<real_t, {q_phase}>;
+using PdfField_hydro_T = GhostLayerField<real_t, {q_hydro}>;
+using VelocityField_T = GhostLayerField<real_t, {dimensions}>;
+using PhaseField_T = GhostLayerField<real_t, 1>;
+#ifndef UTIL_H
+#define UTIL_H
+const char * stencil_phase_name = "{stencil_phase_name}";
+const char * stencil_hydro_name = "{stencil_hydro_name}";
+#endif
"""
with CodeGeneration() as ctx:
- if not ctx.optimize_for_localhost:
- cpu_vec = {'instruction_set': None}
-
generate_sweep(ctx, 'initialize_phase_field_distributions', h_updates, target='cpu')
generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates, target='cpu')
generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
- field_swaps=[(h, h_tmp)],
+ field_swaps=[(h, h_tmp), (C, C_tmp)],
inner_outer_split=True,
- cpu_vectorize_info=cpu_vec, target='cpu')
- generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target='cpu')
+ target='cpu')
+ generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target='cpu', streaming_pattern='pull')
generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
- inner_outer_split=True,
- cpu_vectorize_info=cpu_vec, target='cpu')
- generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target='cpu')
-
- generate_sweep(ctx, 'stream_hydro', stream_hydro,
field_swaps=[(g, g_tmp)],
inner_outer_split=True,
- cpu_vectorize_info=cpu_vec, target='cpu')
+ target='cpu')
+ generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target='cpu', streaming_pattern='push')
- ctx.write_file("GenDefines.h", info_header)
+ # communication
+
+ generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
+ streaming_pattern='pull', target='cpu')
- # TODO: generate communication (PackInfo)
+ generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
+ streaming_pattern='push', target='cpu')
+
+ generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='cpu')
+
+ pystencils_walberla.boundary.generate_boundary(ctx, 'ContactAngle', contact_angle,
+ C.name, stencil_hydro, index_shape=[], target='cpu')
+
+ ctx.write_file("GenDefines.h", info_header)
print("finished code generation successfully")
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
index d424d7c21c5255ec85ed28e493f54176d12ff4be..c040fc2f1df72793085b745787bd9dff6d00b9a5 100755
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
+++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py
@@ -14,8 +14,8 @@ class Scenario:
# simulation parameters
self.timesteps = 10000
- self.cells = (64, 32, 64)
- self.blocks = (1, 4, 1)
+ self.cells = (32, 32, 64)
+ self.blocks = (2, 4, 1)
self.periodic = (0, 0, 0)
self.size = (self.cells[0] * self.blocks[0],
self.cells[1] * self.blocks[1],
@@ -39,6 +39,8 @@ class Scenario:
density_ratio=1000,
viscosity_ratio=100)
+ self.interface_thickness = 5
+
# everything else
self.dbFile = "risingBubble3D.db"
@@ -72,6 +74,7 @@ class Scenario:
'gravitational_acceleration': self.parameters["gravitational_acceleration"],
'relaxation_time_liquid': self.parameters.get("relaxation_time_heavy"),
'relaxation_time_gas': self.parameters.get("relaxation_time_light"),
+ 'interface_thickness': self.interface_thickness
},
'Boundaries': {
'Border': [
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
index 51c6922032a44744f513890ac6642185930ba3af..a8ff39721f322b6bdd5deed047984176ed570c15 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt
@@ -10,14 +10,14 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenGPU
phase_field_LB_NoSlip.cu phase_field_LB_NoSlip.h
hydro_LB_step.cu hydro_LB_step.h
hydro_LB_NoSlip.cu hydro_LB_NoSlip.h
- stream_hydro.cu stream_hydro.h
PackInfo_phase_field_distributions.cu PackInfo_phase_field_distributions.h
PackInfo_phase_field.cu PackInfo_phase_field.h
PackInfo_velocity_based_distributions.cu PackInfo_velocity_based_distributions.h
+ ContactAngle.cu ContactAngle.h
GenDefines.h)
waLBerla_add_executable(NAME multiphaseGPU
- FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp CalculateNormals.cpp contact.cu multiphase_codegen.py
+ FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp util.cpp multiphase_codegen.py
DEPENDS blockforest core cuda field postprocessing lbm geometry timeloop gui PhaseFieldCodeGenGPU)
set_target_properties(multiphaseGPU PROPERTIES CXX_VISIBILITY_PRESET hidden)
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.cpp
deleted file mode 100644
index e20053e38130262fdf2bb8c556f1fd665888a945..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.cpp
+++ /dev/null
@@ -1,83 +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 CalculateNormals.cpp
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-#include "CalculateNormals.h"
-
-#include "core/Environment.h"
-#include "core/logging/Initialization.h"
-
-#include "field/FlagField.h"
-
-namespace walberla
-{
-using FlagField_T = FlagField< uint8_t >;
-using NormalsField_T = GhostLayerField< int8_t, 3 >;
-
-void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
- ConstBlockDataID flagFieldID, FlagUID domainFlagUID, FlagUID boundaryFlagUID)
-{
- for (auto& block : *blocks)
- {
- CellInterval globalCellBB = blocks->getBlockCellBB(block);
- CellInterval blockLocalCellBB;
- blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, block, globalCellBB);
-
- auto* normalsField = block.getData< NormalsField_T >(normalsFieldID);
- auto* flagField = block.getData< FlagField_T >(flagFieldID);
- auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
- auto domainFlag = flagField->getFlag(domainFlagUID);
-
- // clang-format off
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(normalsField,
-
- if( x < blockLocalCellBB.xMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x + 1, y, z) == domainFlag)
- normalsField->get(x, y, z, 0) = 1;
- }
-
- if( x > blockLocalCellBB.xMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x - 1, y, z) == domainFlag)
- normalsField->get(x, y, z, 0) = - 1;
- }
-
- if( y < blockLocalCellBB.yMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y + 1, z) == domainFlag)
- normalsField->get(x, y, z, 1) = 1;
- }
-
- if( y > blockLocalCellBB.yMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y - 1, z) == domainFlag)
- normalsField->get(x, y, z, 1) = - 1;
- }
-
- if( z < blockLocalCellBB.zMax() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z + 1) == domainFlag)
- normalsField->get(x, y, z, 2) = 1;
- }
-
- if( z > blockLocalCellBB.zMin() ){
- if(flagField->get(x, y, z) == boundaryFlag && flagField->get(x, y, z - 1) == domainFlag)
- normalsField->get(x, y, z, 2) = - 1;
- }
-
- )
- // clang-format on
- }
-}
-} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.h b/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.h
deleted file mode 100644
index 901db8544732cea2953e651fc10c78943326bf56..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/CalculateNormals.h
+++ /dev/null
@@ -1,32 +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 CalculateNormals.h
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-
-#include "blockforest/StructuredBlockForest.h"
-
-#include "domain_decomposition/BlockDataID.h"
-#include "domain_decomposition/IBlock.h"
-
-#include "field/FlagField.h"
-
-namespace walberla
-{
-void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
- ConstBlockDataID flagFieldID, FlagUID fluidFlagUID, FlagUID boundaryFlagUID);
-}
\ No newline at end of file
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
index 0fd64cb73cefcfefb25a565fee9c80b4c127b155..979c87ff330f669567df87c2e051a75271280b6b 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.cpp
@@ -53,24 +53,196 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B
}
}
+void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
+ const real_t D = 5, const real_t H = 2, const real_t DT = 20, const real_t Donut_x0 = 40)
+{
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ 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 Ri = D * sqrt(pow(H, 2) - pow(DT - sqrt(pow(globalCell[0] - Mx, 2) + pow(globalCell[2] - Mz, 2)), 2));
+
+ real_t shifter = atan2((globalCell[2] - Mz), (globalCell[0] - Mx));
+ if (shifter < 0) shifter = shifter + 2 * math::pi;
+ if ((globalCell[1] < Donut_x0 + Ri * sin(shifter / 2.0)) && (globalCell[1] > Donut_x0 - Ri)) {
+ phaseField->get(x, y, z) = 0.0;
+ } else { phaseField->get(x, y, z) = 1.0; })
+ }
+}
+
+void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+ real_t W = 5)
+{
+ Vector3< real_t > bubbleMidPoint;
+
+ auto X = blocks->getDomainCellBB().xMax();
+ auto Y = blocks->getDomainCellBB().yMax();
+ auto Z = blocks->getDomainCellBB().zMax();
+
+ // 20 percent from the top are filled with the gas phase
+ real_t gas_top = Y - Y / 5.0;
+
+ // Diameter of the bubble
+ real_t D = R * 2;
+
+ // distance in between the bubbles
+ int dist = 4;
+ auto nx = static_cast< unsigned int >(floor(X / (D + dist * W)));
+ auto nz = static_cast< unsigned int >(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));
+ std::vector< std::vector< real_t > > fluctuation_pos(nx, std::vector< real_t >(nz, 0.0));
+
+ real_t max_fluctuation_radius = R / 5;
+ real_t max_fluctuation_pos = (dist * W) / 3.0;
+ for (unsigned int i = 0; i < nx; ++i)
+ {
+ for (unsigned int j = 0; j < nz; ++j)
+ {
+ fluctuation_radius[i][j] = math::realRandom< real_t >(-max_fluctuation_radius, max_fluctuation_radius);
+ fluctuation_pos[i][j] = math::realRandom< real_t >(-max_fluctuation_pos, max_fluctuation_pos);
+ }
+ }
+
+ 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) {
+ for (unsigned int j = 0; j < nz; ++j)
+ {
+ bubbleMidPoint[0] = (i + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+ bubbleMidPoint[1] = R + W + 4;
+ bubbleMidPoint[2] = (j + 1) * (D + (dist * W)) - (D + (dist * W)) / 2.0 + fluctuation_pos[i][j];
+
+ real_t Ri = sqrt((globalCell[0] - bubbleMidPoint[0]) * (globalCell[0] - bubbleMidPoint[0]) +
+ (globalCell[1] - bubbleMidPoint[1]) * (globalCell[1] - bubbleMidPoint[1]) +
+ (globalCell[2] - bubbleMidPoint[2]) * (globalCell[2] - bubbleMidPoint[2]));
+ if (globalCell[0] >= i * (D + dist * W) && globalCell[0] <= (i + 1) * (D + dist * W) &&
+ globalCell[2] >= j * (D + dist * W) && globalCell[2] <= (j + 1) * (D + dist * W))
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (Ri - (R - fluctuation_radius[i][j])) / W);
+
+ if (globalCell[0] > nx * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ if (globalCell[2] > nz * (D + dist * W)) phaseField->get(x, y, z) = 1.0;
+ }
+ }
+
+ if (globalCell[1] > gas_top) {
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(2.0 * (gas_top + 10 - globalCell[1]) / W);
+ })
+ }
+}
+
void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID,
- const real_t W = 5)
+ const real_t W = 5, const bool pipe = true)
{
auto X = blocks->getDomainCellBB().xMax();
+ auto Z = blocks->getDomainCellBB().zMax();
auto halfY = (blocks->getDomainCellBB().yMax()) / 2.0;
- double perturbation = 0.05;
+ real_t perturbation = 0.05;
- for (auto& block : *blocks)
+ if (pipe)
{
- auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
- // clang-format off
- 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((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
- phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));
- )
- // clang-format on
+ 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((globalCell[0] - X / 2) * (globalCell[0] - X / 2) +
+ (globalCell[2] - Z / 2) * (globalCell[2] - Z / 2));
+ if (R > X) R = X; real_t tmp = perturbation * X * cos((2.0 * math::pi * R) / X);
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) + tmp) / (W / 2.0));)
+ }
+ }
+ else
+ {
+ 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 tmp = perturbation * X *
+ (cos((2.0 * math::pi * globalCell[0]) / X) + cos((2.0 * math::pi * globalCell[2]) / X));
+ phaseField->get(x, y, z) = 0.5 + 0.5 * tanh(((globalCell[1] - halfY) - tmp) / (W / 2.0));)
+ }
+ }
+}
+
+void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID flagFieldID,
+ field::FlagUID boundaryFlagUID, real_t const R_in, real_t const eccentricity,
+ real_t const start_transition, real_t const length_transition,
+ bool const eccentricity_or_pipe_ratio)
+{
+ if (eccentricity_or_pipe_ratio)
+ {
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+
+ real_t const shift = eccentricity * Mx / 2;
+
+ for (auto& block : *blocks)
+ {
+ 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));
+ real_t R1;
+ if (globalCell[1] <= start_transition) {
+ R1 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
+ real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
+ real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ R1 = sqrt((globalCell[0] - Mx - shift_tmp) * (globalCell[0] - Mx - shift_tmp) +
+ (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ } else {
+ R1 = sqrt((globalCell[0] - Mx - shift) * (globalCell[0] - Mx - shift) +
+ (globalCell[2] - Mz) * (globalCell[2] - Mz));
+ }
+
+ real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (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);)
+ }
+ }
+ else
+ {
+ auto Mx = blocks->getDomainCellBB().xMax() / 2.0;
+ auto Mz = blocks->getDomainCellBB().zMax() / 2.0;
+
+ auto R_outer = blocks->getDomainCellBB().xMax() / 2.0 + 1.0;
+
+ real_t const shift = eccentricity * R_in;
+ real_t R_tmp;
+
+ for (auto& block : *blocks)
+ {
+ 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 (globalCell[1] <= start_transition) {
+ R_tmp = R_in;
+ } else if (globalCell[1] > start_transition && globalCell[1] < start_transition + length_transition) {
+ real_t tmp = math::pi * (globalCell[1] - start_transition) / (length_transition);
+ real_t shift_tmp = shift * 0.5 * (1 - cos(tmp));
+ R_tmp = R_in + shift_tmp;
+ } else {
+ R_tmp = R_in + shift;
+ }
+
+ real_t R2 = sqrt((globalCell[0] - Mx) * (globalCell[0] - Mx) + (globalCell[2] - Mz) * (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);)
+ }
}
}
} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.h b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.h
index 04bac53f2e75d17ddd393da0e247b63df05d0704..585639ac9078beb048a8ce694ea97552e59b833a 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.h
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/InitializerFunctions.h
@@ -34,6 +34,17 @@ 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_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5);
+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_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
+ real_t W = 5);
+
+void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W = 5,
+ const bool pipe = true);
+
+void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID flagFieldID,
+ field::FlagUID boundaryFlagUID, real_t R_in, real_t eccentricity, real_t start_transition,
+ real_t length_transition, bool const eccentricity_or_pipe_ratio);
} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/contact.cu b/apps/showcases/PhaseFieldAllenCahn/GPU/contact.cu
deleted file mode 100644
index 48d26d1c5d7ec8e9353e48cd9a81731179ce4de4..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/contact.cu
+++ /dev/null
@@ -1,130 +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 contact.cu
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-
-#include "core/DataTypes.h"
-#include "core/Macros.h"
-
-#include <cmath>
-
-#include "contact.h"
-
-#define FUNC_PREFIX __global__
-
-namespace walberla
-{
-namespace lbm
-{
-#ifdef __GNUC__
-# pragma GCC diagnostic push
-#endif
-
-#ifdef __CUDACC__
-# pragma push
-#endif
-
-namespace internal_boundary_contact
-{
-static FUNC_PREFIX void contact_angle_treatment(uint8_t* WALBERLA_RESTRICT const _data_indexVector, double* WALBERLA_RESTRICT _data_phase,
- int64_t const _stride_phase_0, int64_t const _stride_phase_1,
- int64_t const _stride_phase_2, int64_t indexVectorSize, double alpha)
-{
- if (blockDim.x * blockIdx.x + threadIdx.x < indexVectorSize)
- {
- uint8_t* WALBERLA_RESTRICT _data_indexVector_10 = _data_indexVector;
- const int32_t x = *((int32_t*) (&_data_indexVector_10[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- uint8_t* WALBERLA_RESTRICT _data_indexVector_14 = _data_indexVector + 4;
- const int32_t y = *((int32_t*) (&_data_indexVector_14[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- uint8_t* WALBERLA_RESTRICT _data_indexVector_18 = _data_indexVector + 8;
- const int32_t z = *((int32_t*) (&_data_indexVector_18[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- uint8_t* WALBERLA_RESTRICT _data_indexVector_112 = _data_indexVector + 12;
- const int32_t nx = *((int32_t*) (&_data_indexVector_112[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- const int32_t x1 = x + nx;
- uint8_t* WALBERLA_RESTRICT _data_indexVector_116 = _data_indexVector + 16;
- const int32_t ny = *((int32_t*) (&_data_indexVector_116[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- const int32_t y1 = y + ny;
- uint8_t* WALBERLA_RESTRICT _data_indexVector_200 = _data_indexVector + 20;
- const int32_t nz = *((int32_t*) (&_data_indexVector_200[24 * blockDim.x * blockIdx.x + 24 * threadIdx.x]));
- const int32_t z1 = z + nz;
-
- const double h = 0.5 * sqrt((float) (nx * nx + ny * ny + nz * nz));
- const double a = cos(alpha);
- const double W = 5;
-
- double* WALBERLA_RESTRICT _phase_wall = _data_phase + _stride_phase_1 * y + _stride_phase_2 * z;
- double* WALBERLA_RESTRICT _phase_interior = _data_phase + _stride_phase_1 * y1 + _stride_phase_2 * z1;
- if (h < 0.001) { _phase_wall[_stride_phase_0 * x] = 1.0; }
- else if (a > 1e-8 || a < -1e-8)
- {
- const double var = -h * (4.0 / W) * a;
- _phase_wall[_stride_phase_0 * x] =
- (1 + var - sqrt((1 + var) * (1 + var) - 4 * var * _phase_interior[_stride_phase_0 * x1])) / (var + 1e-12) -
- _phase_interior[_stride_phase_0 * x1];
- }
- else
- {
- _phase_wall[_stride_phase_0 * x] = _phase_interior[_stride_phase_0 * x1];
- }
- }
-}
-} // namespace internal_boundary_contact
-
-#ifdef __GNUC__
-# pragma GCC diagnostic pop
-#endif
-
-#ifdef __CUDACC__
-# pragma pop
-#endif
-
-void contact::run(IBlock* block, IndexVectors::Type type, cudaStream_t stream)
-{
- auto* indexVectors = block->getData< IndexVectors >(indexVectorID);
- int64_t indexVectorSize = int64_c(indexVectors->indexVector(type).size());
- if (indexVectorSize == 0) return;
-
- auto pointer = indexVectors->pointerGpu(type);
-
- auto* _data_indexVector = reinterpret_cast< uint8_t* >(pointer);
-
- auto phaseField = block->getData< cuda::GPUField< double > >(phaseFieldID);
-
- auto& alpha = this->alpha_;
- WALBERLA_ASSERT_GREATER_EQUAL(0, -int_c(phaseField->nrOfGhostLayers()))
- double* WALBERLA_RESTRICT _data_phase = phaseField->dataAt(0, 0, 0, 0);
- const auto _stride_pdfs_0 = int64_t(phaseField->xStride());
- const auto _stride_pdfs_1 = int64_t(phaseField->yStride());
- const auto _stride_pdfs_2 = int64_t(phaseField->zStride());
- dim3 _block(int(((256 < indexVectorSize) ? 256 : indexVectorSize)), int(1), int(1));
- dim3 _grid(int(((indexVectorSize) % (((256 < indexVectorSize) ? 256 : indexVectorSize)) == 0 ?
- (int64_t)(indexVectorSize) / (int64_t)(((256 < indexVectorSize) ? 256 : indexVectorSize)) :
- ((int64_t)(indexVectorSize) / (int64_t)(((256 < indexVectorSize) ? 256 : indexVectorSize))) + 1)),
- int(1), int(1));
- internal_boundary_contact::contact_angle_treatment<<< _grid, _block, 0, stream >>>(
- _data_indexVector, _data_phase, _stride_pdfs_0, _stride_pdfs_1, _stride_pdfs_2, indexVectorSize, alpha);
-}
-
-void contact::operator()(IBlock* block, cudaStream_t stream) { run(block, IndexVectors::ALL, stream); }
-
-void contact::inner(IBlock* block, cudaStream_t stream) { run(block, IndexVectors::INNER, stream); }
-
-void contact::outer(IBlock* block, cudaStream_t stream) { run(block, IndexVectors::OUTER, stream); }
-
-} // namespace lbm
-} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/contact.h b/apps/showcases/PhaseFieldAllenCahn/GPU/contact.h
deleted file mode 100644
index f2cd84cd0c67577ab7a40d8ac720fda46f6e7fe3..0000000000000000000000000000000000000000
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/contact.h
+++ /dev/null
@@ -1,158 +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 contact.h
-//! \author Markus Holzer <markus.holzer@fau.de>
-//
-//======================================================================================================================
-
-#include "blockforest/StructuredBlockForest.h"
-
-#include "core/DataTypes.h"
-
-#include "cuda/GPUField.h"
-
-#include "domain_decomposition/BlockDataID.h"
-#include "domain_decomposition/IBlock.h"
-
-#include "field/FlagField.h"
-
-#include <set>
-#include <vector>
-
-namespace walberla
-{
-namespace lbm
-{
-class contact
-{
- public:
- struct IndexInfo
- {
- int32_t x1;
- int32_t y1;
- int32_t z1;
- int32_t x2;
- int32_t y2;
- int32_t z2;
- IndexInfo(int32_t x1_, int32_t y1_, int32_t z1_, int32_t x2_, int32_t y2_, int32_t z2_)
- : x1(x1_), y1(y1_), z1(z1_), x2(x2_), y2(y2_), z2(z2_)
- {}
- bool operator==(const IndexInfo& o) const
- {
- return x1 == o.x1 && y1 == o.y1 && z1 == o.z1 && x2 == o.x2 && y2 == o.y2 && z2 == o.z2;
- }
- };
-
- class IndexVectors
- {
- public:
- using CpuIndexVector = std::vector< IndexInfo >;
-
- enum Type { ALL = 0, INNER = 1, OUTER = 2, NUM_TYPES = 3 };
-
- IndexVectors() : cpuVectors_(NUM_TYPES) {}
- bool operator==(IndexVectors& other) { return other.cpuVectors_ == cpuVectors_; }
-
- ~IndexVectors()
- {
- for (auto& gpuVec : gpuVectors_)
- cudaFree(gpuVec);
- }
-
- CpuIndexVector& indexVector(Type t) { return cpuVectors_[t]; }
- IndexInfo* pointerCpu(Type t) { return &(cpuVectors_[t][0]); }
-
- IndexInfo* pointerGpu(Type t) { return gpuVectors_[t]; }
-
- void syncGPU()
- {
- gpuVectors_.resize(cpuVectors_.size());
- for (size_t i = 0; i < size_t(NUM_TYPES); ++i)
- {
- auto& gpuVec = gpuVectors_[i];
- auto& cpuVec = cpuVectors_[i];
- cudaFree(gpuVec);
- cudaMalloc(&gpuVec, sizeof(IndexInfo) * cpuVec.size());
- cudaMemcpy(gpuVec, &cpuVec[0], sizeof(IndexInfo) * cpuVec.size(), cudaMemcpyHostToDevice);
- }
- }
-
- private:
- std::vector< CpuIndexVector > cpuVectors_;
-
- using GpuIndexVector = IndexInfo*;
- std::vector< GpuIndexVector > gpuVectors_;
- };
-
- contact(const shared_ptr< StructuredBlockForest >& blocks, BlockDataID phaseFieldID_, double alpha)
- : phaseFieldID(phaseFieldID_), alpha_(alpha)
- {
- auto createIdxVector = [](IBlock* const, StructuredBlockStorage* const) { return new IndexVectors(); };
- indexVectorID = blocks->addStructuredBlockData< IndexVectors >(createIdxVector, "IndexField_hydro_NoSlip_gpu");
- };
-
- void operator()(IBlock* block, cudaStream_t stream = nullptr);
- void inner(IBlock* block, cudaStream_t stream = nullptr);
- void outer(IBlock* block, cudaStream_t stream = nullptr);
-
- template< typename NormalField_T >
- void fillFromNormalField(const shared_ptr< StructuredBlockForest >& blocks, ConstBlockDataID normalFieldID)
- {
- for (auto& block : *blocks)
- {
- auto* indexVectors = block.getData< IndexVectors >(indexVectorID);
- auto& indexVectorAll = indexVectors->indexVector(IndexVectors::ALL);
- auto& indexVectorInner = indexVectors->indexVector(IndexVectors::INNER);
- auto& indexVectorOuter = indexVectors->indexVector(IndexVectors::OUTER);
-
- auto* normalField = block.getData< NormalField_T >(normalFieldID);
-
- auto inner = normalField->xyzSize();
- inner.expand(cell_idx_t(-1));
-
- indexVectorAll.clear();
- indexVectorInner.clear();
- indexVectorOuter.clear();
- // clang-format off
- WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(normalField, Cell globalCell;
- blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z));
- if (normalField->get(x, y, z, 0) != 0 || normalField->get(x, y, z, 1) != 0 || normalField->get(x, y, z, 2) != 0)
- {
- auto element = IndexInfo(x, y, z, normalField->get(x, y, z, 0), normalField->get(x, y, z, 1), normalField->get(x, y, z, 2));
- indexVectorAll.push_back(element);
- if (inner.contains(x, y, z))
- indexVectorInner.push_back(element);
- else
- indexVectorOuter.push_back(element);
- }
- )
- // clang-format off
- indexVectors->syncGPU();
- }
- }
-
- private:
- void run(IBlock* block, IndexVectors::Type type, cudaStream_t stream = nullptr);
-
- BlockDataID indexVectorID;
-
- public:
- BlockDataID phaseFieldID;
- double alpha_;
-};
-
-} // namespace lbm
-} // namespace walberla
\ No newline at end of file
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/droplet_contact_angle.py b/apps/showcases/PhaseFieldAllenCahn/GPU/droplet_contact_angle.py
index 1b3e9d6a8e15642ea4855d6be13eb378bcb8cc1e..34bdb9fdceb1432638a89d5670c7c8ce537158d5 100755
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/droplet_contact_angle.py
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/droplet_contact_angle.py
@@ -18,25 +18,24 @@ class Scenario:
self.overlappingWidth = (8, 1, 1)
self.timeStepStrategy = 'normal'
- self.contactAngle = 22
-
# bubble parameters
self.dropletRadius = 24.0
self.dropletMidPoint = (64, 24, 64)
# everything else
- self.scenario = 1 # 1 rising bubble, 2 RTI
+ self.scenario = 1 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up
self.counter = 0
self.yPositions = []
@wlb.member_callback
- def config(self, **kwargs):
+ def config(self):
return {
'DomainSetup': {
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
+ 'tube': False
},
'Parameters': {
'timesteps': self.timesteps,
@@ -45,7 +44,6 @@ class Scenario:
'overlappingWidth': self.overlappingWidth,
'remainingTimeLoggerFrequency': 10.0,
'scenario': self.scenario,
- 'contactAngle': self.contactAngle
},
'PhysicalParameters': {
'density_liquid': 1.0,
@@ -55,6 +53,7 @@ class Scenario:
'gravitational_acceleration': 0.0,
'relaxation_time_liquid': 3 * 0.166,
'relaxation_time_gas': 3 * 0.0166,
+ 'interface_thickness': 5
},
'Boundaries': {
'Border': [
@@ -69,7 +68,7 @@ class Scenario:
'Bubble': {
'bubbleMidPoint': self.dropletMidPoint,
'bubbleRadius': self.dropletRadius,
- 'bubble': False
+ 'bubble': False # this means we are simulating a droplet rather than a bubble
},
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
index de54dfe91b39cce943fa15ea1338e0131842866c..a97ea043fd4facc1010f2838b99e5269c91a6023 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp
@@ -27,6 +27,7 @@
#include "cuda/AddGPUFieldToStorage.h"
#include "cuda/DeviceSelectMPI.h"
+#include "cuda/NVTX.h"
#include "cuda/ParallelStreams.h"
#include "cuda/communication/UniformGPUScheme.h"
@@ -36,6 +37,11 @@
#include "field/vtk/VTKWriter.h"
#include "geometry/InitBoundaryHandling.h"
+#include "geometry/mesh/TriangleMeshIO.h"
+
+#include "lbm/vtk/QCriterion.h"
+
+#include "postprocessing/FieldToSurfaceMesh.h"
#include "python_coupling/CreateConfig.h"
#include "python_coupling/PythonCallback.h"
@@ -43,10 +49,9 @@
#include "timeloop/SweepTimeloop.h"
-#include "CalculateNormals.h"
#include "InitializerFunctions.h"
#include "PythonExports.h"
-#include "contact.h"
+#include "util.h"
//////////////////////////////
// INCLUDE GENERATED FILES //
@@ -62,7 +67,7 @@
#include "initialize_velocity_based_distributions.h"
#include "phase_field_LB_NoSlip.h"
#include "phase_field_LB_step.h"
-#include "stream_hydro.h"
+#include "ContactAngle.h"
////////////
// USING //
@@ -70,14 +75,30 @@
using namespace walberla;
-using PdfField_phase_T = GhostLayerField< real_t, Stencil_phase_T::Size >;
-using PdfField_hydro_T = GhostLayerField< real_t, Stencil_hydro_T::Size >;
-using VelocityField_T = GhostLayerField< real_t, Stencil_hydro_T::Dimension >;
-using NormalsField_T = GhostLayerField< int8_t, Stencil_hydro_T::Dimension >;
-using PhaseField_T = GhostLayerField< real_t, 1 >;
-using FlagField_T = FlagField< uint8_t >;
+using NormalsField_T = GhostLayerField< int8_t, Stencil_hydro_T::Dimension >;
+using FlagField_T = FlagField< uint8_t >;
+
+typedef cuda::GPUField< real_t > GPUField;
+typedef cuda::GPUField< uint8_t > GPUField_int;
+
+class Filter
+{
+ public:
+ explicit Filter(Vector3< uint_t > numberOfCells) : numberOfCells_(numberOfCells) {}
+
+ void operator()(const IBlock& /*block*/) {}
+
+ bool operator()(const cell_idx_t x, const cell_idx_t y, const cell_idx_t z) const
+ {
+ return x >= -1 && x <= cell_idx_t(numberOfCells_[0]) && y >= -1 && y <= cell_idx_t(numberOfCells_[1]) &&
+ z >= -1 && z <= cell_idx_t(numberOfCells_[2]);
+ }
-typedef cuda::GPUField< double > GPUField;
+ private:
+ Vector3< uint_t > numberOfCells_;
+};
+
+using FluidFilter_T = Filter;
int main(int argc, char** argv)
{
@@ -99,6 +120,7 @@ int main(int argc, char** argv)
auto domainSetup = config->getOneBlock("DomainSetup");
Vector3< uint_t > cellsPerBlock = domainSetup.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const bool tube = domainSetup.getParameter< bool >("tube", true);
////////////////////////////////////////
// ADD GENERAL SIMULATION PARAMETERS //
@@ -109,9 +131,7 @@ int main(int argc, char** argv)
const uint_t timesteps = parameters.getParameter< uint_t >("timesteps", uint_c(50));
const real_t remainingTimeLoggerFrequency =
parameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0);
- const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
- const real_t alpha = parameters.getParameter< real_t >("contactAngle", real_c(90));
- const real_t alpha_rad = alpha * (math::pi / 180);
+ const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
Vector3< int > overlappingWidth =
parameters.getParameter< Vector3< int > >("overlappingWidth", Vector3< int >(1, 1, 1));
Vector3< int > gpuBlockSize =
@@ -124,7 +144,6 @@ int main(int argc, char** argv)
// CPU fields
BlockDataID vel_field = field::addToStorage< VelocityField_T >(blocks, "vel", real_t(0), field::fzyx);
BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_t(0), field::fzyx);
- BlockDataID normals = field::addToStorage< NormalsField_T >(blocks, "normals", int8_t(0), field::fzyx);
// GPU fields
BlockDataID lb_phase_field_gpu = cuda::addGPUFieldToStorage< cuda::GPUField< real_t > >(
blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
@@ -135,71 +154,94 @@ int main(int argc, char** argv)
BlockDataID phase_field_gpu =
cuda::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
// Flag field
- BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+ BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
+ BlockDataID flagFieldID_gpu = cuda::addGPUFieldToStorage< FlagField_T >(blocks, flagFieldID, "flag on GPU", true);
+
+ auto physical_parameters = config->getOneBlock("PhysicalParameters");
+ const real_t density_liquid = physical_parameters.getParameter< real_t >("density_liquid", real_c(1.0));
+ const real_t density_gas = physical_parameters.getParameter< real_t >("density_gas");
+ const real_t surface_tension = physical_parameters.getParameter< real_t >("surface_tension");
+ const real_t mobility = physical_parameters.getParameter< real_t >("mobility");
+ const real_t gravitational_acceleration =
+ physical_parameters.getParameter< real_t >("gravitational_acceleration");
+ const real_t relaxation_time_liquid = physical_parameters.getParameter< real_t >("relaxation_time_liquid");
+ const real_t relaxation_time_gas = physical_parameters.getParameter< real_t >("relaxation_time_gas");
+ const real_t interface_thickness = physical_parameters.getParameter< real_t >("interface_thickness");
+
+ std::array< real_t, 3 > center_of_mass = { 0.0, 0.0, 0.0 };
- WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field")
+ WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field")
if (scenario == 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", 20.0);
const bool bubble = bubbleParameters.getParameter< bool >("bubble", true);
- initPhaseField_sphere(blocks, phase_field, bubbleRadius, bubbleMidPoint, bubble);
+ initPhaseField_sphere(blocks, phase_field, bubbleRadius, bubbleMidPoint, bubble, interface_thickness);
+
}
else if (scenario == 2)
{
- initPhaseField_RTI(blocks, phase_field);
+ initPhaseField_RTI(blocks, phase_field, interface_thickness, tube);
+ }
+ else if (scenario == 3)
+ {
+ auto bubbleParameters = config->getOneBlock("Bubble");
+ const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", 20.0);
+ init_bubble_field(blocks, phase_field, bubbleRadius);
+ }
+ else if (scenario == 4)
+ {
+ auto TorusParameters = config->getOneBlock("Torus");
+ const real_t midpoint = TorusParameters.getParameter< real_t >("Donut_midpoint");
+ const real_t height = TorusParameters.getParameter< real_t >("Donut_h");
+ const real_t diameter = TorusParameters.getParameter< real_t >("Donut_D");
+ const real_t donutTime = TorusParameters.getParameter< real_t >("DonutTime");
+ init_Taylor_bubble(blocks, phase_field, diameter, height, donutTime, midpoint);
+ center_of_mass[0] = real_t(cellsPerBlock[0]);
+ center_of_mass[1] = real_t(midpoint);
+ center_of_mass[2] = real_t(cellsPerBlock[2]);
}
- WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done")
+ WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field done")
/////////////////
// ADD SWEEPS //
///////////////
- auto physical_parameters = config->getOneBlock("PhysicalParameters");
- const real_t density_liquid = physical_parameters.getParameter< real_t >("density_liquid", real_c(1.0));
- const real_t density_gas = physical_parameters.getParameter< real_t >("density_gas");
- const real_t surface_tension = physical_parameters.getParameter< real_t >("surface_tension");
- const real_t mobility = physical_parameters.getParameter< real_t >("mobility");
- const real_t gravitational_acceleration =
- physical_parameters.getParameter< real_t >("gravitational_acceleration");
- const real_t relaxation_time_liquid = physical_parameters.getParameter< real_t >("relaxation_time_liquid");
- const real_t relaxation_time_gas = physical_parameters.getParameter< real_t >("relaxation_time_gas");
-
- pystencils::initialize_phase_field_distributions init_h(lb_phase_field_gpu, phase_field_gpu, vel_field_gpu);
+ pystencils::initialize_phase_field_distributions init_h(lb_phase_field_gpu, phase_field_gpu, vel_field_gpu,
+ interface_thickness);
pystencils::initialize_velocity_based_distributions init_g(lb_velocity_field_gpu, vel_field_gpu);
pystencils::phase_field_LB_step phase_field_LB_step(
- lb_phase_field_gpu, phase_field_gpu, vel_field_gpu, mobility, gpuBlockSize[0], gpuBlockSize[1],
- Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
-
- pystencils::hydro_LB_step hydro_LB_step(
- lb_velocity_field_gpu, phase_field_gpu, vel_field_gpu, gravitational_acceleration, density_liquid, density_gas,
- surface_tension, relaxation_time_liquid, relaxation_time_gas, gpuBlockSize[0], gpuBlockSize[1],
+ flagFieldID_gpu, lb_phase_field_gpu, phase_field_gpu, vel_field_gpu, interface_thickness, mobility,
+ gpuBlockSize[0], gpuBlockSize[1], gpuBlockSize[2],
Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
- pystencils::stream_hydro stream_hydro(lb_velocity_field_gpu, gpuBlockSize[0], gpuBlockSize[1],
- Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
+ pystencils::hydro_LB_step hydro_LB_step(flagFieldID_gpu, lb_velocity_field_gpu, phase_field_gpu, vel_field_gpu,
+ gravitational_acceleration, interface_thickness, density_liquid,
+ density_gas, surface_tension, relaxation_time_liquid, relaxation_time_gas,
+ gpuBlockSize[0], gpuBlockSize[1], gpuBlockSize[2],
+ Cell(overlappingWidth[0], overlappingWidth[1], overlappingWidth[2]));
////////////////////////
// ADD COMMUNICATION //
//////////////////////
auto Comm_velocity_based_distributions =
- make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
+ make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 1);
auto generatedPackInfo_velocity_based_distributions =
- make_shared< pystencils::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
+ make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
Comm_velocity_based_distributions->addPackInfo(generatedPackInfo_velocity_based_distributions);
- auto Comm_phase_field = make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
+ auto Comm_phase_field = make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 1);
auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field_gpu);
Comm_phase_field->addPackInfo(generatedPackInfo_phase_field);
auto Comm_phase_field_distributions =
- make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 0);
+ make_shared< cuda::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, 1);
auto generatedPackInfo_phase_field_distributions =
- make_shared< pystencils::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
+ make_shared< lbm::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
Comm_phase_field_distributions->addPackInfo(generatedPackInfo_phase_field_distributions);
////////////////////////
@@ -213,17 +255,27 @@ int main(int argc, char** argv)
if (boundariesConfig)
{
geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
+ if (tube)
+ {
+ const real_t inner_radius = domainSetup.getParameter< real_t >("inner_radius", real_c(0));
+ const real_t eccentricity = domainSetup.getParameter< real_t >("ratio", real_c(0));
+ const real_t start_transition = domainSetup.getParameter< real_t >("start_transition", real_c(60));
+ const real_t length_transition = domainSetup.getParameter< real_t >("length_transition", real_c(10));
+ const bool eccentricity_or_pipe_ration = domainSetup.getParameter< bool >("eccentricity_or_pipe_ration", true);
+ initTubeWithCylinder(blocks, flagFieldID, wallFlagUID, inner_radius, eccentricity, start_transition,
+ length_transition, eccentricity_or_pipe_ration);
+ }
geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID);
}
+ cuda::fieldCpy< GPUField_int, FlagField_T >(blocks, flagFieldID_gpu, flagFieldID);
- calculate_normals(blocks, normals, flagFieldID, fluidFlagUID, wallFlagUID);
lbm::phase_field_LB_NoSlip phase_field_LB_NoSlip(blocks, lb_phase_field_gpu);
lbm::hydro_LB_NoSlip hydro_LB_NoSlip(blocks, lb_velocity_field_gpu);
- lbm::contact contact_angle(blocks, phase_field_gpu, alpha_rad);
+ pystencils::ContactAngle contact_angle(blocks, phase_field_gpu, interface_thickness);
phase_field_LB_NoSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldID, wallFlagUID, fluidFlagUID);
hydro_LB_NoSlip.fillFromFlagField< FlagField_T >(blocks, flagFieldID, wallFlagUID, fluidFlagUID);
- contact_angle.fillFromNormalField< NormalsField_T >(blocks, normals);
+ contact_angle.fillFromFlagField< FlagField_T >(blocks, flagFieldID, wallFlagUID, fluidFlagUID);
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the normals-field done")
@@ -251,7 +303,6 @@ int main(int argc, char** argv)
Comm_velocity_based_distributions->communicate(nullptr);
hydro_LB_NoSlip(&block);
- stream_hydro(&block);
}
};
auto simpleOverlapTimeStep = [&]() {
@@ -277,13 +328,6 @@ int main(int argc, char** argv)
for (auto& block : *blocks)
hydro_LB_NoSlip(&block);
-
- Comm_velocity_based_distributions->startCommunication(defaultStream);
- for (auto& block : *blocks)
- stream_hydro.inner(&block, defaultStream);
- Comm_velocity_based_distributions->wait(defaultStream);
- for (auto& block : *blocks)
- stream_hydro.outer(&block, defaultStream);
};
std::function< void() > timeStep;
if (timeStepStrategy == "overlap")
@@ -299,6 +343,15 @@ int main(int argc, char** argv)
timeLoop->add() << BeforeFunction(timeStep) << Sweep([](IBlock*) {}, "time step");
+ if (scenario == 4)
+ {
+ python_coupling::PythonCallback smear_interface("interface_diffusion");
+ if (smear_interface.isCallable())
+ {
+ smear_interface.data().exposeValue("blocks", blocks);
+ smear_interface();
+ }
+ }
cuda::fieldCpy< GPUField, PhaseField_T >(blocks, phase_field_gpu, phase_field);
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the PDFs")
@@ -309,23 +362,100 @@ int main(int argc, char** argv)
}
WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the PDFs done")
uint_t dbWriteFrequency = parameters.getParameter< uint_t >("dbWriteFrequency", 10000000);
+ int targetRank = 0;
timeLoop->addFuncAfterTimeStep(
[&]() {
if (timeLoop->getCurrentTimeStep() % dbWriteFrequency == 0)
{
cuda::fieldCpy< PhaseField_T, GPUField >(blocks, phase_field, phase_field_gpu);
- python_coupling::PythonCallback callback("at_end_of_time_step");
- if (callback.isCallable())
+ cuda::fieldCpy< VelocityField_T, GPUField >(blocks, vel_field, vel_field_gpu);
+ if (scenario == 4)
+ {
+ std::array< real_t, 4 > total_velocity = { 0.0, 0.0, 0.0, 0.0 };
+ real_t volume;
+ uint_t nrCells;
+ PhaseField_T gatheredPhaseField(0, 0, 0, 0);
+ VelocityField_T gatheredVelocityField(0, 0, 0, 0);
+
+ CellInterval boundingBox = blocks->getDomainCellBB();
+ if (cell_idx_t(center_of_mass[1] - cell_idx_t(cellsPerBlock[0]) * 1.5) >= 0)
+ boundingBox.min()[1] = cell_idx_t(center_of_mass[1] - cell_idx_t(cellsPerBlock[0]) * 1.5);
+ if (cell_idx_t(center_of_mass[1] + cell_idx_t(cellsPerBlock[0]) * 1.5) <= boundingBox.max()[1])
+ boundingBox.max()[1] = cell_idx_t(center_of_mass[1] + cell_idx_t(cellsPerBlock[0]) * 1.5);
+
+ field::gather< PhaseField_T >(gatheredPhaseField, blocks, phase_field, boundingBox, targetRank);
+ field::gather< VelocityField_T >(gatheredVelocityField, blocks, vel_field, boundingBox, targetRank);
+
+ WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank)
+ {
+ flood_fill(gatheredPhaseField, gatheredVelocityField, boundingBox, volume, nrCells, center_of_mass,
+ total_velocity);
+ }
+ WALBERLA_MPI_SECTION() { walberla::mpi::broadcastObject(center_of_mass, targetRank); }
+
+ python_coupling::PythonCallback callback("at_end_of_time_step");
+ if (callback.isCallable())
+ {
+ callback.data().exposeValue("blocks", blocks);
+ callback.data().exposeValue("timeStep", timeLoop->getCurrentTimeStep());
+ callback.data().exposeValue("target_rank", targetRank);
+ callback.data().exposeValue("bounding_box_min", boundingBox.min()[1]);
+ callback.data().exposeValue("bounding_box_max", boundingBox.max()[1]);
+ callback.data().exposeValue("total_velocity", total_velocity[0]);
+ callback.data().exposeValue("total_velocity_X", total_velocity[1]);
+ callback.data().exposeValue("total_velocity_Y", total_velocity[2]);
+ callback.data().exposeValue("total_velocity_Z", total_velocity[3]);
+ callback.data().exposeValue("center_of_mass_X", center_of_mass[0]);
+ callback.data().exposeValue("center_of_mass_Y", center_of_mass[1]);
+ callback.data().exposeValue("center_of_mass_Z", center_of_mass[2]);
+ callback.data().exposeValue("sum_inv_phi", volume);
+ callback.data().exposeValue("gas_cells_of_the_taylor_bubble", nrCells);
+ callback.data().exposeValue("stencil_phase", stencil_phase_name);
+ callback.data().exposeValue("stencil_hydro", stencil_hydro_name);
+ callback();
+ }
+ }
+ else
{
- callback.data().exposeValue("blocks", blocks);
- callback.data().exposeValue( "timeStep", timeLoop->getCurrentTimeStep());
- callback();
+ python_coupling::PythonCallback callback("at_end_of_time_step");
+ if (callback.isCallable())
+ {
+ callback.data().exposeValue("blocks", blocks);
+ callback.data().exposeValue("timeStep", timeLoop->getCurrentTimeStep());
+ callback.data().exposeValue("stencil_phase", stencil_phase_name);
+ callback.data().exposeValue("stencil_hydro", stencil_hydro_name);
+ callback();
+ }
}
}
},
"Python callback");
+ int meshWriteFrequency = parameters.getParameter< int >("meshWriteFrequency", 0);
+ int counter = 0;
+ if (meshWriteFrequency > 0)
+ {
+ timeLoop->addFuncAfterTimeStep(
+ [&]() {
+ if (timeLoop->getCurrentTimeStep() % uint_t(meshWriteFrequency) == 0)
+ {
+ auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field, 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 + "taylor_bubble_D_" + std::to_string(cellsPerBlock[0]) + "_" + out.str() + ".obj", *mesh);
+ counter++;
+ }
+ }
+ },
+ "Mesh writer");
+ }
+
// remaining time logger
timeLoop->addFuncAfterTimeStep(
timing::RemainingTimeLogger(timeLoop->getNrOfTimeSteps(), remainingTimeLoggerFrequency),
@@ -339,19 +469,22 @@ int main(int argc, char** argv)
"simulation_step", false, true, true, false, 0);
vtkOutput->addBeforeFunction([&]() {
cuda::fieldCpy< PhaseField_T, GPUField >(blocks, phase_field, phase_field_gpu);
- // cuda::fieldCpy<VelocityField_T, GPUField>( blocks, vel_field, vel_field_gpu );
+ cuda::fieldCpy<VelocityField_T, GPUField>( blocks, vel_field, vel_field_gpu );
});
- auto phaseWriter = make_shared< field::VTKWriter< PhaseField_T > >(phase_field, "PhaseField");
+ auto phaseWriter = make_shared< field::VTKWriter< PhaseField_T, float > >(phase_field, "PhaseField");
vtkOutput->addCellDataWriter(phaseWriter);
- // auto normlasWriter = make_shared<field::VTKWriter<NormalsField_T>>(normals, "Normals");
- // vtkOutput->addCellDataWriter(normlasWriter);
+ auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(flagFieldID, "flag");
+ vtkOutput->addCellDataWriter(flagWriter);
+
+ auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float > >(vel_field, "Velocity");
+ vtkOutput->addCellDataWriter(velWriter);
- // auto flagWriter = make_shared<field::VTKWriter<FlagField_T>>(flagFieldID, "flag");
- // vtkOutput->addCellDataWriter(flagWriter);
+ FluidFilter_T filter(cellsPerBlock);
- // auto velWriter = make_shared<field::VTKWriter<VelocityField_T>>(vel_field, "Velocity");
- // vtkOutput->addCellDataWriter(velWriter);
+ auto QCriterionWriter = make_shared< lbm::QCriterionVTKWriter< VelocityField_T, FluidFilter_T, float > >(
+ blocks, filter, vel_field, "Q-Criterion");
+ vtkOutput->addCellDataWriter(QCriterionWriter);
timeLoop->addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
}
@@ -371,6 +504,5 @@ int main(int argc, char** argv)
WALBERLA_LOG_RESULT_ON_ROOT("MLUPS per process: " << mlupsPerProcess)
WALBERLA_LOG_RESULT_ON_ROOT("Time per time step: " << time / real_c(timesteps) << " s")
}
-
return EXIT_SUCCESS;
}
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_RTI_3D.py b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_RTI_3D.py
index c9798504de79d23c5bbaaa686b209c917d8a63f5..024b66a9e15c4291f045d500d4493de2448f137b 100755
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_RTI_3D.py
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_RTI_3D.py
@@ -1,8 +1,11 @@
+import os
+
import waLBerla as wlb
-import waLBerla.tools.sqlitedb as wlbSqlite
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
from waLBerla.core_extension import makeSlice
import numpy as np
+import pandas as pd
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
@@ -10,12 +13,11 @@ class Scenario:
def __init__(self):
# output frequencies
self.vtkWriteFrequency = 1000
- self.dbWriteFrequency = 200
# simulation parameters
- self.timesteps = 27001
+ self.time = 2 # physical time in seconds
- self.cells = (64, 256, 64)
+ self.cells = (128, 512, 128)
self.blocks = (1, 1, 1)
self.periodic = (1, 0, 1)
self.size = (self.cells[0] * self.blocks[0],
@@ -24,25 +26,40 @@ class Scenario:
# physical parameters
self.density_heavy = 1.0
- self.reference_time = 6000
- self.parameters = calculate_parameters_rti(reference_length=128,
+ self.reference_time = 4000
+ self.dbWriteFrequency = self.reference_time // 20
+ self.timesteps = int(self.reference_time * self.time) + 1
+
+ self.capillary_number = 8.7
+ self.reynolds_number = 3000
+ self.atwood_number = 1
+ self.peclet_number = 744
+ self.density_ratio = 1000
+ self.viscosity_ratio = 100
+
+ self.parameters = calculate_parameters_rti(reference_length=self.cells[0],
reference_time=self.reference_time,
density_heavy=self.density_heavy,
- capillary_number=9.1,
- reynolds_number=128,
- atwood_number=1.0,
- peclet_number=140,
- density_ratio=3,
- viscosity_ratio=3)
+ capillary_number=self.capillary_number,
+ reynolds_number=self.reynolds_number,
+ atwood_number=self.atwood_number,
+ peclet_number=self.peclet_number,
+ density_ratio=self.density_ratio,
+ viscosity_ratio=self.viscosity_ratio)
+
+ self.interface_thickness = 5
+ self.tube = False
# everything else
- self.dbFile = "RTI.db"
+ self.dbFile = "RTI.csv"
- self.scenario = 2 # 1 rising bubble, 2 RTI
+ self.scenario = 2 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up
self.counter = 0
self.yPositions = []
+ self.config_dict = self.config()
+
@wlb.member_callback
def config(self):
return {
@@ -50,23 +67,24 @@ class Scenario:
'blocks': self.blocks,
'cellsPerBlock': self.cells,
'periodic': self.periodic,
+ 'tube': self.tube
},
'Parameters': {
'timesteps': self.timesteps,
'vtkWriteFrequency': self.vtkWriteFrequency,
'dbWriteFrequency': self.dbWriteFrequency,
- 'useGui': 0,
'remainingTimeLoggerFrequency': 10.0,
'scenario': self.scenario,
},
'PhysicalParameters': {
'density_liquid': self.density_heavy,
- 'density_gas': self.parameters["density_light"],
- 'surface_tension': self.parameters["surface_tension"],
- 'mobility': self.parameters.get("mobility", 0.1),
- 'gravitational_acceleration': self.parameters["gravitational_acceleration"],
+ 'density_gas': self.parameters.get("density_light"),
+ 'surface_tension': self.parameters.get("surface_tension"),
+ 'mobility': self.parameters.get("mobility"),
+ 'gravitational_acceleration': self.parameters.get("gravitational_acceleration"),
'relaxation_time_liquid': self.parameters.get("relaxation_time_heavy"),
'relaxation_time_gas': self.parameters.get("relaxation_time_light"),
+ 'interface_thickness': self.interface_thickness
},
'Boundaries': {
'Border': [
@@ -79,55 +97,79 @@ class Scenario:
@wlb.member_callback
def at_end_of_time_step(self, blocks, **kwargs):
t = kwargs['timeStep']
- ny = self.size[1]
- l0 = self.size[0]
if t % self.dbWriteFrequency == 0:
- location_of_spike = -100
- location_of_bubble = -100
- location_of_saddle = -100
- mass = -100
- spike_data = wlb.field.gather(blocks, 'phase', makeSlice[self.size[0] // 2, :, self.size[2] // 2])
- if spike_data:
- spike_field = np.asarray(spike_data).squeeze()
- location_of_spike = (np.argmax(spike_field > 0.5) - ny // 2) / l0
-
- bubble_data = wlb.field.gather(blocks, 'phase', makeSlice[0, :, 0])
- if bubble_data:
- bubble_field = np.asarray(bubble_data).squeeze()
- location_of_bubble = (np.argmax(bubble_field > 0.5) - ny // 2) / l0
-
- saddle_data = wlb.field.gather(blocks, 'phase', makeSlice[0, :, self.size[2] // 2])
- if saddle_data:
- saddle_field = np.asarray(saddle_data).squeeze()
- location_of_saddle = (np.argmax(saddle_field > 0.5) - ny // 2) / l0
-
phase = wlb.field.gather(blocks, 'phase', makeSlice[:, :, :])
if phase:
+ data = {'timestep': t}
+ data.update(self.config_dict['PhysicalParameters'])
+ data.update({'total_timesteps': self.timesteps})
+ data.update({'normalized_time': t / self.reference_time})
+ data.update({'tube_setup': self.tube})
+ data.update({'interface_thickness': self.interface_thickness})
+ data.update({'capillary_number': self.capillary_number})
+ data.update({'reynolds_number': self.reynolds_number})
+ data.update({'atwood_number': self.atwood_number})
+ data.update({'peclet_number': self.peclet_number})
+ data.update({'density_ratio': self.density_ratio})
+ data.update({'viscosity_ratio': self.viscosity_ratio})
+ data.update({'reference_time': self.reference_time})
+ data.update(kwargs)
+
phase_field = np.asarray(phase).squeeze()
+ stable = np.isfinite(np.sum(phase_field))
mass = np.sum(phase_field)
+ rel_max = np.max(phase_field) - 1
+ rel_min = np.min(phase_field)
+ data.update({'mass': mass})
+ data.update({'rel_max': rel_max})
+ data.update({'rel_min': rel_min})
+ data.update({'stable': stable})
+
+ if self.tube:
+ location_of_spike = self.get_interface_location(
+ phase_field[self.size[0] // 2, :, self.size[2] // 2])
+ a = np.where(phase_field < 0.5)
+ value = np.argmax(a[1])
+ location_of_bubble = self.get_interface_location(
+ phase_field[a[0][value], a[1][value] - 10:a[1][value] + 10, a[2][value]], a[1][value] - 10)
+
+ data.update({'location_of_spike': location_of_spike})
+ data.update({'location_of_bubble': location_of_bubble})
+ else:
+ location_of_spike = self.get_interface_location(
+ phase_field[self.size[0] // 2, :, self.size[2] // 2])
+ location_of_bubble = self.get_interface_location(phase_field[0, :, 0])
+ location_of_saddle = self.get_interface_location(phase_field[0, :, self.size[2] // 2])
+
+ data.update({'location_of_spike': location_of_spike})
+ data.update({'location_of_bubble': location_of_bubble})
+ data.update({'location_of_saddle': location_of_saddle})
+
+ sequenceValuesToScalars(data)
+
+ csv_file = f"RTI_{data['stencil_phase']}_{data['stencil_hydro']}_Re_{self.reynolds_number}_tube.csv"
+
+ df = pd.DataFrame.from_records([data])
+ if not os.path.isfile(csv_file):
+ df.to_csv(csv_file, index=False)
+ else:
+ df.to_csv(csv_file, index=False, mode='a', header=False)
+
+ def get_interface_location(self, one_dimensional_array, shift=None):
+ ny = self.size[1]
+ l0 = self.size[0]
- self.write_result_to_database(t, location_of_spike, location_of_bubble, location_of_saddle, mass)
-
- def write_result_to_database(self, t, spike, bubble, saddle, mass):
- """Writes the simulation result stored in the global variables shapeFactors and angles to
- an sqlite3 database, and resets the global variables."""
- result = {'waLBerlaVersion': wlb.build_info.version,
- 'xCells': self.size[0],
- 'yCells': self.size[1],
- 'zCells': self.size[2],
- 'spike': spike,
- 'bubble': bubble,
- 'saddle': saddle,
- 'mass': mass,
- 'timesteps': t,
- 'normalized_time': t / self.reference_time,
- 'processes': wlb.mpi.numProcesses(),
- }
- try:
- wlbSqlite.checkAndUpdateSchema(result, 'interface_location', self.dbFile)
- wlbSqlite.storeSingle(result, 'interface_location', self.dbFile)
- except Exception as e:
- wlb.log_warning("Failed to store run in database " + str(e) + "\n" + str(result))
+ index = np.argmax(one_dimensional_array > 0.5)
+
+ if index > 0:
+ zw1 = one_dimensional_array[index]
+ zw2 = one_dimensional_array[index - 1]
+ absolute_location = (index - 1) + (zw2 - 0.5) / (zw2 - zw1)
+ if shift:
+ absolute_location += shift
+ return (absolute_location - ny // 2) / l0
+ else:
+ return -100
scenarios = wlb.ScenarioManager()
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
index 48a19621db6f0e62be8c6bc21f6ecc4c65faa90a..0abf47375eaa5d43924d6582761307278ecb52ea 100644
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py
@@ -1,24 +1,33 @@
+from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from pystencils import fields, TypedSymbol
from pystencils.simp import sympy_cse
-from pystencils import AssignmentCollection
+from pystencils import Assignment
+from pystencils.astnodes import Block, Conditional
from lbmpy.boundaries import NoSlip
-from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
+from lbmpy.creationfunctions import create_lb_method
from lbmpy.stencils import get_stencil
-from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_from_kernel
-from lbmpy_walberla import generate_boundary
+import pystencils_walberla
+from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field
+from lbmpy_walberla import generate_boundary, generate_lb_pack_info
from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_phase_field_lb,\
- initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, get_collision_assignments_hydro
+ initializer_kernel_hydro_lb, interface_tracking_force, hydrodynamic_force, get_collision_assignments_hydro,\
+ get_collision_assignments_phase
from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
import numpy as np
import sympy as sp
-stencil_phase = get_stencil("D3Q19")
-stencil_hydro = get_stencil("D3Q27")
+stencil_phase_name = "D3Q27"
+stencil_hydro_name = "D3Q27"
+
+contact_angle_in_degrees = 22
+
+stencil_phase = get_stencil(stencil_phase_name)
+stencil_hydro = get_stencil(stencil_hydro_name)
q_phase = len(stencil_phase)
q_hydro = len(stencil_hydro)
@@ -43,7 +52,7 @@ relaxation_time_gas = sp.Symbol("tau_L")
# phase-field parameter
drho3 = (density_liquid - density_gas) / 3
# interface thickness
-W = 5
+W = sp.Symbol("interface_thickness")
# coefficients related to surface tension
beta = 12.0 * (surface_tension / W)
kappa = 1.5 * surface_tension * W
@@ -56,6 +65,7 @@ kappa = 1.5 * surface_tension * W
u = fields(f"vel_field({dimensions}): [{dimensions}D]", layout='fzyx')
# phase-field
C = fields(f"phase_field: [{dimensions}D]", layout='fzyx')
+C_tmp = fields(f"phase_field_tmp: [{dimensions}D]", layout='fzyx')
flag = fields(f"flag_field: uint8[{dimensions}D]", layout='fzyx')
# phase-field distribution functions
@@ -91,12 +101,17 @@ relaxation_rate_cutoff = sp.Piecewise((1 / (0.5 + relaxation_time_liquid), C.cen
# LBM METHODS #
###############
-method_phase = create_lb_method(stencil=stencil_phase, method='srt',
- relaxation_rate=relaxation_rate_allen_cahn, compressible=True)
+# method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
+# relaxation_rates=[1, 1.5, 1, 1.5, 1, 1.5])
+method_phase = create_lb_method(stencil=stencil_phase, method="mrt", compressible=True, weighted=True,
+ relaxation_rates=[1, 1, 1, 1, 1, 1])
+
+method_phase.set_conserved_moments_relaxation_rate(relaxation_rate_allen_cahn)
method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1])
+
# create the kernels for the initialization of the g and h field
h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W, fd_stencil=get_stencil("D3Q27"))
g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
@@ -107,60 +122,70 @@ force_model_h = MultiphaseForceModel(force=force_h)
force_g = hydrodynamic_force(g, C, method_hydro, relaxation_time, density_liquid, density_gas, kappa, beta, body_force,
fd_stencil=get_stencil("D3Q27"))
+force_model_g = MultiphaseForceModel(force=force_g, rho=density)
+
####################
# LBM UPDATE RULES #
####################
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
+phase_field_LB_step = get_collision_assignments_phase(lb_method=method_phase,
+ velocity_input=u,
+ output={'density': C_tmp},
+ force_model=force_model_h,
+ symbolic_fields={"symbolic_field": h,
+ "symbolic_temporary_field": h_tmp},
+ kernel_type='stream_pull_collide')
-method_phase.set_force_model(force_model_h)
-
-phase_field_LB_step = create_lb_update_rule(lb_method=method_phase,
- velocity_input=u,
- compressible=True,
- optimization={"symbolic_field": h,
- "symbolic_temporary_field": h_tmp},
- kernel_type='stream_pull_collide')
-
-phase_field_LB_step.set_main_assignments_from_dict({**phase_field_LB_step.main_assignments_dict, **{C.center: sum_h}})
-phase_field_LB_step = AssignmentCollection(main_assignments=phase_field_LB_step.main_assignments,
- subexpressions=phase_field_LB_step.subexpressions)
phase_field_LB_step = sympy_cse(phase_field_LB_step)
+phase_field_LB_step = [Conditional(sp.Eq(flag.center(), 2),
+ Block(phase_field_LB_step),
+ Block([Assignment(C_tmp.center, C.center)]))]
# ---------------------------------------------------------------------------------------------------------
hydro_LB_step = get_collision_assignments_hydro(lb_method=method_hydro,
density=density,
velocity_input=u,
- force=force_g,
+ force_model=force_model_g,
sub_iterations=2,
symbolic_fields={"symbolic_field": g,
"symbolic_temporary_field": g_tmp},
- kernel_type='collide_only')
+ kernel_type='collide_stream_push')
hydro_LB_step.set_sub_expressions_from_dict({**{relaxation_rate: relaxation_rate_cutoff},
**hydro_LB_step.subexpressions_dict})
-stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
- optimization={"symbolic_field": g,
- "symbolic_temporary_field": g_tmp},
- kernel_type='stream_pull_only')
+hydro_LB_step = [Conditional(sp.Eq(flag.center(), 2), Block(hydro_LB_step))]
+
+contact_angle = ContactAngle(contact_angle_in_degrees, W)
+
###################
# GENERATE SWEEPS #
###################
vp = [('int32_t', 'cudaBlockSize0'),
- ('int32_t', 'cudaBlockSize1')]
+ ('int32_t', 'cudaBlockSize1'),
+ ('int32_t', 'cudaBlockSize2')]
sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int32),
TypedSymbol("cudaBlockSize1", np.int32),
- 1)
+ TypedSymbol("cudaBlockSize2", np.int32))
+
sweep_params = {'block_size': sweep_block_size}
info_header = f"""
+using namespace walberla;
#include "stencil/D3Q{q_phase}.h"\nusing Stencil_phase_T = walberla::stencil::D3Q{q_phase};
#include "stencil/D3Q{q_hydro}.h"\nusing Stencil_hydro_T = walberla::stencil::D3Q{q_hydro};
+using PdfField_phase_T = GhostLayerField<real_t, {q_phase}>;
+using PdfField_hydro_T = GhostLayerField<real_t, {q_hydro}>;
+using VelocityField_T = GhostLayerField<real_t, {dimensions}>;
+using PhaseField_T = GhostLayerField<real_t, 1>;
+#ifndef UTIL_H
+#define UTIL_H
+const char * stencil_phase_name = "{stencil_phase_name}";
+const char * stencil_hydro_name = "{stencil_hydro_name}";
+#endif
"""
with CodeGeneration() as ctx:
@@ -168,35 +193,33 @@ with CodeGeneration() as ctx:
generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates, target='gpu')
generate_sweep(ctx, 'phase_field_LB_step', phase_field_LB_step,
- field_swaps=[(h, h_tmp)],
+ field_swaps=[(h, h_tmp), (C, C_tmp)],
inner_outer_split=True,
target='gpu',
gpu_indexing_params=sweep_params,
varying_parameters=vp)
- generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target='gpu')
+ generate_boundary(ctx, 'phase_field_LB_NoSlip', NoSlip(), method_phase, target='gpu', streaming_pattern='pull')
generate_sweep(ctx, 'hydro_LB_step', hydro_LB_step,
- inner_outer_split=True,
- target='gpu',
- gpu_indexing_params=sweep_params,
- varying_parameters=vp)
- generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target='gpu')
-
- generate_sweep(ctx, 'stream_hydro', stream_hydro,
field_swaps=[(g, g_tmp)],
inner_outer_split=True,
target='gpu',
gpu_indexing_params=sweep_params,
varying_parameters=vp)
+ generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target='gpu', streaming_pattern='push')
# communication
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field_distributions',
- phase_field_LB_step.main_assignments, target='gpu')
- generate_pack_info_from_kernel(ctx, 'PackInfo_phase_field',
- hydro_LB_step.all_assignments, target='gpu', kind='pull')
- generate_pack_info_from_kernel(ctx, 'PackInfo_velocity_based_distributions',
- stream_hydro.all_assignments, target='gpu', kind='pull')
+ generate_lb_pack_info(ctx, 'PackInfo_phase_field_distributions', stencil_phase, h,
+ streaming_pattern='pull', target='gpu')
+
+ generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
+ streaming_pattern='push', target='gpu')
+
+ generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target='gpu')
+
+ pystencils_walberla.boundary.generate_boundary(ctx, 'ContactAngle', contact_angle,
+ C.name, stencil_hydro, index_shape=[], target='gpu')
ctx.write_file("GenDefines.h", info_header)
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_rising_bubble.py b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_rising_bubble.py
index 9f5f08323daf8ee885fd35d510260eda43512daf..186590774aa43e8624eaf891f94616aad081360f 100755
--- a/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_rising_bubble.py
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_rising_bubble.py
@@ -39,6 +39,8 @@ class Scenario:
density_ratio=1000,
viscosity_ratio=100)
+ self.interface_thickness = 5
+
# everything else
self.dbFile = "risingBubble3D.db"
@@ -72,6 +74,7 @@ class Scenario:
'gravitational_acceleration': self.parameters["gravitational_acceleration"],
'relaxation_time_liquid': self.parameters.get("relaxation_time_heavy"),
'relaxation_time_gas': self.parameters.get("relaxation_time_light"),
+ 'interface_thickness': self.interface_thickness
},
'Boundaries': {
'Border': [
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/parameters_taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/GPU/parameters_taylor_bubble.py
new file mode 100644
index 0000000000000000000000000000000000000000..55db7aa693400b12f71d097470c57bed2609ee8e
--- /dev/null
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/parameters_taylor_bubble.py
@@ -0,0 +1,69 @@
+import math
+
+
+def calculate_parameters_taylor_bubble(reference_length=128,
+ reference_time=16000,
+ density_heavy=1.0,
+ d1=0.0254,
+ d2=0.0127):
+ r"""
+ Calculate the simulation parameters for a rising Taylor bubble. The calculation can be found in
+ 10.1016/S0009-2509(97)00210-8 by G. Das
+
+ Args:
+ reference_length: chosen reference length
+ reference_time: chosen reference time
+ density_heavy: density of the heavier fluid
+ d1: diameter of the outer tube
+ d2: diameter of the inner cylinder
+ """
+
+ water_rho = 998 # kg/m3
+ air_rho = 1.2047 # kg/m3
+ surface_tension = 0.07286 # kg/s2
+ water_mu = 1.002e-3 # kg/ms
+
+ water_nu = water_mu / water_rho # m2/s
+ air_mu = 1.8205e-5 # kg/ms
+ air_nu = air_mu / air_rho # m2/s
+ gravity = 9.81 # m/s2
+
+ dh = d1 - d2
+ dr = d1 / d2
+ de = d1 + d2
+ # ur = 0.1695 # (0.28913, 0.23882, 0.1695)
+
+ inverse_viscosity_number = math.sqrt((water_rho - air_rho) * water_rho * gravity * dh ** 3) / water_mu
+ bond_number = (water_rho - air_rho) * gravity * dh ** 2 / surface_tension
+ morton_number = gravity * water_mu ** 4 * (water_rho - air_rho) / (water_rho ** 2 * surface_tension ** 3)
+
+ d = reference_length / dr
+
+ density_light = 1.0 / (water_rho / air_rho)
+ dh = reference_length - d
+ g = dh / reference_time ** 2
+
+ mu_h = math.sqrt((density_heavy - density_light) * density_heavy * g * dh ** 3) / inverse_viscosity_number
+ mu_l = mu_h / (water_mu / air_mu)
+
+ dynamic_viscosity_heavy = mu_h / density_heavy
+ dynamic_viscosity_light = mu_l / density_light
+
+ relaxation_time_heavy = 3 * dynamic_viscosity_heavy
+ relaxation_time_light = 3 * dynamic_viscosity_light
+
+ sigma = (density_heavy - density_light) * g * dh ** 2 / bond_number
+
+ parameters = {
+ "inverse_viscosity_number": inverse_viscosity_number,
+ "bond_number": bond_number,
+ "morton_number": morton_number,
+ "density_light": density_light,
+ "dynamic_viscosity_heavy": dynamic_viscosity_heavy,
+ "dynamic_viscosity_light": dynamic_viscosity_light,
+ "relaxation_time_heavy": relaxation_time_heavy,
+ "relaxation_time_light": relaxation_time_light,
+ "gravitational_acceleration": -g,
+ "surface_tension": sigma
+ }
+ return parameters
diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/GPU/taylor_bubble.py
new file mode 100755
index 0000000000000000000000000000000000000000..0afb25234ac119bff33d6876e23bfec7cd79391b
--- /dev/null
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/taylor_bubble.py
@@ -0,0 +1,325 @@
+import math
+import os
+import pickle as pl
+import waLBerla as wlb
+from matplotlib import pyplot as plt
+import numpy as np
+import pandas as pd
+
+from waLBerla.core_extension import makeSlice
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+from scipy.ndimage.filters import gaussian_filter
+import scipy.spatial as spatial
+
+from parameters_taylor_bubble import calculate_parameters_taylor_bubble
+
+
+def intersection(points1, points2, eps):
+ tree = spatial.KDTree(points1)
+ distances, indices = tree.query(points2, k=1, distance_upper_bound=eps)
+ intersection_points = tree.data[indices[np.isfinite(distances)]]
+ return intersection_points
+
+
+def contour_points(contour, steps=1):
+ return np.row_stack([path.interpolated(steps).vertices
+ for linecol in contour.collections
+ for path in linecol.get_paths()])
+
+
+def test_line(center_x, center_z, angle, contour, tol):
+ contact = -1
+
+ line_size = 200
+
+ points1 = contour_points(contour)
+ points2 = np.zeros((line_size, 2))
+ points2[:, 0] = center_x + np.linspace(0, center_x, line_size) * np.cos(np.radians(angle))
+ points2[:, 1] = center_z + np.linspace(0, center_z, line_size) * np.sin(np.radians(angle))
+
+ intersection_points = intersection(points1, points2, tol)
+
+ if len(intersection_points) != 0:
+ contact = 1
+
+ return contact
+
+
+def get_circle(midpoint, radius):
+ theta = np.linspace(0, 2 * np.pi, 100)
+
+ a = midpoint[0] + radius * np.cos(theta)
+ b = midpoint[1] + radius * np.sin(theta)
+
+ return a, b
+
+
+class Scenario:
+ def __init__(self):
+ self.density_liquid = 1.0
+ self.reference_time = 16000
+ self.reference_length = 128
+ d1 = 0.0254 # (0.0508, 0.0381, 0.0254)
+ d2 = 0.0127 # (0.0254, 0.0127, 0.0127)
+ gpus = 1
+ self.interface_width = 5
+ self.mobility = 0.05
+
+ # output frequencies
+ self.vtkWriteFrequency = self.reference_time
+ self.dbWriteFrequency = self.reference_time // 25
+ self.meshWriteFrequency = self.reference_time
+ self.pngWriteFrequency = self.reference_time
+
+ # simulation parameters
+ self.diameter = self.reference_length
+ self.timesteps = self.reference_time * 15 + 1
+ self.cells = (self.diameter, (self.reference_length * 15) // gpus, self.diameter)
+ self.blocks = (1, gpus, 1)
+ self.periodic = (0, 0, 0)
+ self.size = (self.cells[0] * self.blocks[0],
+ self.cells[1] * self.blocks[1],
+ self.cells[2] * self.blocks[2])
+ self.inner_radius = self.diameter // 4
+
+ self.center_x = self.size[0] / 2
+ self.center_y = self.size[1] / 2
+ self.center_z = self.size[2] / 2
+
+ self.overlappingWidth = (8, 1, 1)
+ self.timeStepStrategy = 'normal'
+
+ self.scenario = 4 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 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.5
+
+ self.start_transition = (self.size[1] // 2) - 2 * self.diameter
+ self.length_transition = 4 * self.diameter
+
+ setup = "eccentricity" if self.eccentricity_or_pipe_ratio else "ratio"
+
+ self.csv_file = f"Taylor_bubble_D_{self.diameter}_DasC_{setup}_{self.ratio}_W_" \
+ f"{self.interface_width}_M_{self.mobility}.csv"
+
+ d = self.diameter / 2
+ dh = self.diameter - d
+
+ resolution = self.diameter / 128
+
+ self.Donut_D = 0.1 * self.diameter / resolution
+ self.Donut_h = dh / 6
+ self.DonutTime = 0.5 * (self.diameter + d) / 2
+
+ parameters = calculate_parameters_taylor_bubble(reference_length=self.reference_length,
+ reference_time=self.reference_time,
+ density_heavy=self.density_liquid,
+ d1=d1,
+ d2=d2)
+
+ self.density_gas = parameters["density_light"]
+ self.surface_tension = parameters["surface_tension"]
+
+ self.gravitational_acceleration = parameters["gravitational_acceleration"]
+
+ self.relaxation_time_liquid = parameters.get("relaxation_time_heavy")
+ self.relaxation_time_gas = parameters.get("relaxation_time_light")
+
+ self.config_dict = self.config()
+
+ @wlb.member_callback
+ def config(self):
+ return {
+ 'DomainSetup': {
+ 'blocks': self.blocks,
+ 'cellsPerBlock': self.cells,
+ '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_liquid,
+ 'density_gas': self.density_gas,
+ 'surface_tension': self.surface_tension,
+ 'mobility': self.mobility,
+ 'gravitational_acceleration': self.gravitational_acceleration,
+ 'relaxation_time_liquid': self.relaxation_time_liquid,
+ 'relaxation_time_gas': self.relaxation_time_gas,
+ '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'},
+ ],
+ },
+ 'Torus': {
+ 'Donut_midpoint': 5 * self.diameter,
+ 'Donut_h': self.Donut_h,
+ 'Donut_D': self.Donut_D,
+ 'DonutTime': self.DonutTime
+ }
+ }
+
+ @wlb.member_callback
+ def interface_diffusion(self, blocks):
+ for block in blocks:
+ phase_field_array = wlb.field.toArray(block['phase'])
+ phase_field_array[:, :, :] = gaussian_filter(phase_field_array[:, :, :], sigma=2)
+
+ @wlb.member_callback
+ def at_end_of_time_step(self, blocks, **kwargs):
+ t = kwargs["timeStep"]
+ target_rank = kwargs["target_rank"]
+ bounding_box_min = kwargs["bounding_box_min"]
+ bounding_box_max = kwargs["bounding_box_max"]
+ center_of_mass = [kwargs["center_of_mass_X"], kwargs["center_of_mass_Y"], kwargs["center_of_mass_Z"]]
+ if t % self.dbWriteFrequency == 0:
+ wlb_field = wlb.field.gather(blocks, 'phase', makeSlice[:, bounding_box_min:bounding_box_max, :],
+ target_rank)
+ if wlb_field:
+ data = {'timestep': t}
+ data.update(self.config_dict['Parameters'])
+ data.update(self.config_dict['DomainSetup'])
+ data.update(self.config_dict['Torus'])
+ data.update(kwargs)
+ del data["bounding_box_min"]
+ del data["bounding_box_max"]
+ data['total_velocity_Y / sum_inv_phi'] = kwargs['total_velocity_Y'] / kwargs['sum_inv_phi']
+
+ phase_field = np.asarray(wlb_field).squeeze()
+ assert np.isfinite(np.sum(phase_field)), "NaN detected in bounding Box"
+ location_of_gas = np.where(phase_field < 0.5)
+ bubble_tip = np.max(location_of_gas, axis=1)
+
+ fig_handle = plt.figure()
+
+ plt.axis('equal')
+ ax = plt.gca()
+ ax.set(ylim=(0, self.diameter + 1))
+ my_contour = plt.contour(np.rot90(phase_field[:, math.floor(center_of_mass[1]) - bounding_box_min, :]),
+ [0.5])
+
+ # For eccentricity test cases
+ center_x = self.center_x
+ center_z = self.center_z
+ new_radius = self.inner_radius
+
+ if self.eccentricity_or_pipe_ratio:
+ shift = self.ratio * self.center_x / 2
+
+ if center_of_mass[1] < self.start_transition:
+ center_x = self.center_x
+ elif self.start_transition < center_of_mass[1] < self.start_transition + self.length_transition:
+ tmp = math.pi * (center_of_mass[1] - self.start_transition) / self.length_transition
+ shift_tmp = shift * 0.5 * (1 - math.cos(tmp))
+ center_x = self.center_x + shift_tmp
+ else:
+ center_x = self.center_x + shift
+
+ else:
+ shift = self.ratio * self.center_x / 2
+
+ if center_of_mass[1] < self.start_transition:
+ new_radius = self.inner_radius
+ elif self.start_transition < center_of_mass[1] < self.start_transition + self.length_transition:
+ tmp = math.pi * (center_of_mass[1] - self.start_transition) / self.length_transition
+ shift_tmp = shift * 0.5 * (1 - math.cos(tmp))
+ new_radius = self.inner_radius + shift_tmp
+ else:
+ new_radius = self.inner_radius + shift
+
+ start_angle = 0
+ tol = 0.5
+
+ # Search for two lines where one intersects and one does not:
+ contact = test_line(center_x, center_z, start_angle, my_contour, tol)
+ angle = 0
+ angle1 = 0
+
+ if contact == -1:
+ num = np.linspace(0, 180, 500)
+ for i in range(0, len(num)):
+ test = test_line(center_x, center_z, num[i], my_contour, tol)
+ if test != -1:
+ angle = num[i]
+ break
+
+ num = np.linspace(0, -180, 500)
+ for i in range(0, len(num)):
+ test = test_line(center_x, center_z, num[i], my_contour, tol)
+ if test != -1:
+ angle1 = num[i]
+ break
+
+ theta = 360 - (angle - angle1)
+ else:
+ theta = 360
+
+ if t % self.pngWriteFrequency == 0:
+ plt.plot(center_x + np.linspace(0, center_x) * np.cos(np.radians(angle)),
+ center_z + np.linspace(0, center_z) * np.sin(np.radians(angle)), 'b-')
+ plt.plot(center_x + np.linspace(0, center_x) * np.cos(np.radians(angle1)),
+ center_z + np.linspace(0, center_z) * np.sin(np.radians(angle1)), 'r-')
+
+ radius = self.diameter // 2
+ circle1 = get_circle([radius, radius], radius)
+ plt.plot(circle1[0], circle1[1], 'k--', linewidth=2)
+
+ circle2 = get_circle([center_x, center_z], new_radius)
+ plt.fill(circle2[0], circle2[1], 'lightgrey')
+
+ plt.savefig(f"angle_measurement_ratio_{self.ratio}_{self.counter:06d}.png", dpi=600)
+
+ outfile = open(f"angle_measurement_ratio_{self.ratio}_{self.counter:06d}.pkl", 'wb')
+
+ pl.dump(fig_handle, outfile)
+ outfile.close()
+ self.counter += 1
+ plt.cla()
+
+ data['bubble_tip_y'] = bubble_tip[1]
+ 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['theta'] = theta
+ if self.eccentricity_or_pipe_ratio:
+ data['eccentricity'] = self.ratio
+ else:
+ data['pipe_ratio'] = self.ratio
+
+ 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/util.cpp b/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9912edb76e510947883101e72ca60af2b97959d6
--- /dev/null
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/util.cpp
@@ -0,0 +1,147 @@
+//======================================================================================================================
+//
+// 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 util.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "core/cell/Cell.h"
+#include "core/math/Constants.h"
+#include "core/math/Random.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/VTKWriter.h"
+
+#include <queue>
+namespace walberla
+{
+using flag_t = walberla::uint8_t;
+using PhaseField_T = GhostLayerField< real_t, 1 >;
+using VelocityField_T = GhostLayerField< real_t, 3 >;
+using FlagField_T = FlagField< flag_t >;
+
+void calc_total_velocity(const shared_ptr< StructuredBlockStorage >& blocks, std::array< real_t, 5 >& total_velocity,
+ BlockDataID phaseFieldID, BlockDataID velocityFieldID, ConstBlockDataID flagFieldID,
+ FlagUID fluidFlagUID)
+{
+ for (auto& block : *blocks)
+ {
+ auto phaseField = block.getData< PhaseField_T >(phaseFieldID);
+ auto velField = block.getData< VelocityField_T >(velocityFieldID);
+ auto flagField = block.getData< FlagField_T >(flagFieldID);
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(
+ phaseField,
+
+ auto fluidFlag = flagField->getFlag(fluidFlagUID); if (flagField->get(x, y, z) == fluidFlag) {
+ if (phaseField->get(x, y, z) < 0.5)
+ {
+ real_t invC = 1 - phaseField->get(x, y, z);
+ real_t U = velField->get(x, y, z, 0);
+ real_t V = velField->get(x, y, z, 1);
+ real_t W = velField->get(x, y, z, 2);
+ total_velocity[0] += sqrt((U * invC) * (U * invC) + (V * invC) * (V * invC) + (W * invC) * (W * invC));
+ total_velocity[1] += U * invC;
+ total_velocity[2] += V * invC;
+ total_velocity[3] += W * invC;
+ total_velocity[4] += invC;
+ }
+ })
+ }
+}
+
+void flood_fill(PhaseField_T& phaseField, VelocityField_T& velocityField, CellInterval boundingBox, real_t& volume,
+ uint_t& nrOfCells, std::array< real_t, 3 >& center_of_mass, std::array< real_t, 4 >& total_velocity)
+{
+ Cell startCell(boundingBox.xSize() / 2, boundingBox.ySize() / 2, boundingBox.zSize() / 2);
+ Field< bool, 1 > visit(phaseField.xSize(), phaseField.ySize(), phaseField.zSize(), false, field::fzyx);
+ using namespace stencil;
+
+ volume = 0;
+ nrOfCells = 0;
+ center_of_mass[0] = 0.0;
+ center_of_mass[1] = 0.0;
+ center_of_mass[2] = 0.0;
+
+ while (phaseField.get(startCell) > 0.5 && startCell.x() > 0)
+ --startCell.x();
+
+ if (phaseField.get(startCell) > 0.5) WALBERLA_ABORT("startCell for flood fill was not suitable")
+
+ std::queue< Cell > cellQueue;
+ cellQueue.push(startCell);
+ visit.get(startCell) = true;
+
+ real_t invC = 1 - phaseField.get(startCell);
+ real_t v_U = velocityField.get(startCell, 0);
+ real_t v_V = velocityField.get(startCell, 1);
+ real_t v_W = velocityField.get(startCell, 2);
+
+ nrOfCells++;
+ volume += invC;
+
+ total_velocity[0] += sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC));
+ total_velocity[1] += v_U * invC;
+ total_velocity[2] += v_V * invC;
+ total_velocity[3] += v_W * invC;
+
+ center_of_mass[0] += (startCell.x() + boundingBox.xMin());
+ center_of_mass[1] += (startCell.y() + boundingBox.yMin());
+ center_of_mass[2] += (startCell.z() + boundingBox.xMin());
+
+ const int DIRS[6] = { N, S, E, W, T, B };
+
+ CellInterval sizeInterval = phaseField.xyzSize();
+ while (!cellQueue.empty())
+ {
+ Cell& cell = cellQueue.front();
+ cellQueue.pop();
+
+ for (int i : DIRS)
+ {
+ Cell neighborCell(cell.x() + cx[i], cell.y() + cy[i], cell.z() + cz[i]);
+
+ if (!sizeInterval.contains(neighborCell)) { continue; }
+
+ if (phaseField.get(neighborCell) < 0.5 && !visit.get(neighborCell))
+ {
+ invC = 1 - phaseField.get(neighborCell);
+ v_U = velocityField.get(neighborCell, 0);
+ v_V = velocityField.get(neighborCell, 1);
+ v_W = velocityField.get(neighborCell, 2);
+
+ nrOfCells++;
+ volume += invC;
+
+ total_velocity[0] +=
+ sqrt((v_U * invC) * (v_U * invC) + (v_V * invC) * (v_V * invC) + (v_W * invC) * (v_W * invC));
+ total_velocity[1] += v_U * invC;
+ total_velocity[2] += v_V * invC;
+ total_velocity[3] += v_W * invC;
+
+ center_of_mass[0] += (neighborCell.x() + boundingBox.xMin());
+ center_of_mass[1] += (neighborCell.y() + boundingBox.yMin());
+ center_of_mass[2] += (neighborCell.z() + boundingBox.xMin());
+
+ visit.get(neighborCell) = true;
+ cellQueue.push(neighborCell);
+ }
+ }
+ }
+ center_of_mass[0] = center_of_mass[0] / real_t(nrOfCells);
+ center_of_mass[1] = center_of_mass[1] / real_t(nrOfCells);
+ center_of_mass[2] = center_of_mass[2] / real_t(nrOfCells);
+}
+} // namespace walberla
diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.h b/apps/showcases/PhaseFieldAllenCahn/GPU/util.h
similarity index 55%
rename from apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.h
rename to apps/showcases/PhaseFieldAllenCahn/GPU/util.h
index 901db8544732cea2953e651fc10c78943326bf56..cefd57d38736d9f8370a86c78e07215fa947e851 100644
--- a/apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.h
+++ b/apps/showcases/PhaseFieldAllenCahn/GPU/util.h
@@ -13,20 +13,31 @@
// 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 CalculateNormals.h
+//! \file util.h
//! \author Markus Holzer <markus.holzer@fau.de>
//
//======================================================================================================================
-#include "blockforest/StructuredBlockForest.h"
+#include "python_coupling/DictWrapper.h"
-#include "domain_decomposition/BlockDataID.h"
-#include "domain_decomposition/IBlock.h"
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Constants.h"
+#include "field/communication/PackInfo.h"
#include "field/FlagField.h"
+#include "field/vtk/VTKWriter.h"
+#include "GenDefines.h"
+#pragma once
+
+namespace walberla {
+
+ void calc_total_velocity(const shared_ptr <StructuredBlockStorage> &blocks, std::array<real_t, 5> &total_velocity,
+ BlockDataID phaseFieldID, BlockDataID velocityFieldID, ConstBlockDataID flagFieldID, FlagUID fluidFlagUID);
+
+ void flood_fill(PhaseField_T &phaseField, VelocityField_T &velocityField, CellInterval boundingBox,
+ real_t &volume, uint_t &nrOfCells,
+ std::array<real_t, 3> ¢er_of_mass, std::array<real_t, 4> &total_velocity);
+
+} // namespace walberla
-namespace walberla
-{
-void calculate_normals(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID normalsFieldID,
- ConstBlockDataID flagFieldID, FlagUID fluidFlagUID, FlagUID boundaryFlagUID);
-}
\ No newline at end of file
diff --git a/python/pystencils_walberla/boundary.py b/python/pystencils_walberla/boundary.py
index c81df201f43951031035c48771e3c891f3c36e33..22d3635dbb1143b7c52c0c4e1dc7339aeae16507 100644
--- a/python/pystencils_walberla/boundary.py
+++ b/python/pystencils_walberla/boundary.py
@@ -82,6 +82,7 @@ def generate_boundary(generation_context,
'target': target,
'namespace': namespace,
'inner_or_boundary': boundary_object.inner_or_boundary,
+ 'single_link': boundary_object.single_link,
'additional_data_handler': additional_data_handler
}
diff --git a/python/pystencils_walberla/templates/Boundary.tmpl.h b/python/pystencils_walberla/templates/Boundary.tmpl.h
index 57630039fd2f5b6b623b6b68fd081e25a3f289f0..33c355a530e9313f52d9c623a5e71da71d0a9c7c 100644
--- a/python/pystencils_walberla/templates/Boundary.tmpl.h
+++ b/python/pystencils_walberla/templates/Boundary.tmpl.h
@@ -197,6 +197,7 @@ public:
indexVectorInner.clear();
indexVectorOuter.clear();
+ {% if inner_or_boundary -%}
for( auto it = flagField->begin(); it != flagField->end(); ++it )
{
if( ! isFlagSet(it, domainFlag) )
@@ -204,11 +205,7 @@ public:
{%- for dirIdx, dirVec, offset in additional_data_handler.stencil_info %}
if ( isFlagSet( it.neighbor({{offset}} {%if dim == 3%}, 0 {%endif %}), boundaryFlag ) )
{
- {% if inner_or_boundary -%}
auto element = {{StructName}}(it.x(), it.y(), {%if dim == 3%} it.z(), {%endif %} {{dirIdx}} );
- {% else -%}
- auto element = {{StructName}}(it.x() + cell_idx_c({{dirVec[0]}}), it.y() + cell_idx_c({{dirVec[1]}}), {%if dim == 3%} it.z() + cell_idx_c({{dirVec[2]}}), {%endif %} {{additional_data_handler.inverse_directions[dirIdx]}} );
- {% endif -%}
{{additional_data_handler.data_initialisation(dirIdx)|indent(16)}}
indexVectorAll.push_back( element );
if( inner.contains( it.x(), it.y(), it.z() ) )
@@ -218,6 +215,61 @@ public:
}
{% endfor %}
}
+ {%else%}
+ auto flagWithGLayers = flagField->xyzSizeWithGhostLayer();
+ real_t dot = 0.0; real_t maxn = 0.0;
+ cell_idx_t calculated_idx = 0;
+ cell_idx_t dx = 0; cell_idx_t dy = 0; {%if dim == 3%} cell_idx_t dz = 0; {% endif %}
+ cell_idx_t sum_x = 0; cell_idx_t sum_y = 0; {%if dim == 3%} cell_idx_t sum_z = 0; {%endif %}
+ for( auto it = flagField->beginWithGhostLayerXYZ(); it != flagField->end(); ++it )
+ {
+ {% if single_link -%}
+ sum_x = 0; sum_y = 0; {%if dim == 3%} sum_z = 0; {%endif %}
+ {% endif %}
+ if( ! isFlagSet(it, boundaryFlag) )
+ continue;
+ {%- for dirIdx, dirVec, offset in additional_data_handler.stencil_info %}
+ if ( flagWithGLayers.contains(it.x() + cell_idx_c({{dirVec[0]}}), it.y() + cell_idx_c({{dirVec[1]}}), it.z() + cell_idx_c({{dirVec[2]}})) && isFlagSet( it.neighbor({{offset}} {%if dim == 3%}, 0 {%endif %}), domainFlag ) )
+ {
+ {% if single_link -%}
+ sum_x += cell_idx_c({{dirVec[0]}}); sum_y += cell_idx_c({{dirVec[1]}}); {%if dim == 3%} sum_z += cell_idx_c({{dirVec[2]}}); {%endif %}
+ {% else %}
+ auto element = {{StructName}}(it.x(), it.y(), {%if dim == 3%} it.z(), {%endif %} {{dirIdx}} );
+ {{additional_data_handler.data_initialisation(dirIdx)|indent(16)}}
+ indexVectorAll.push_back( element );
+ if( inner.contains( it.x(), it.y(), it.z() ) )
+ indexVectorInner.push_back( element );
+ else
+ indexVectorOuter.push_back( element );
+ {% endif %}
+ }
+ {% endfor %}
+
+ {% if single_link %}
+ dot = 0.0; maxn = 0.0; calculated_idx = 0;
+ if(sum_x != 0 or sum_y !=0 {%if dim == 3%} or sum_z !=0 {%endif %})
+ {
+ {%- for dirIdx, dirVec, offset in additional_data_handler.stencil_info %}
+ dx = {{dirVec[0]}}; dy = {{dirVec[1]}}; {%if dim == 3%} dz = {{dirVec[2]}}; {% endif %}
+ dot = dx*sum_x + dy*sum_y {%if dim == 3%} + dz*sum_z {% endif %};
+ if (dot > maxn)
+ {
+ maxn = dot;
+ calculated_idx = {{dirIdx}};
+ }
+ {% endfor %}
+ auto element = {{StructName}}(it.x(), it.y(), {%if dim == 3%} it.z(), {%endif %} calculated_idx );
+ indexVectorAll.push_back( element );
+ if( inner.contains( it.x(), it.y(), it.z() ) )
+ indexVectorInner.push_back( element );
+ else
+ indexVectorOuter.push_back( element );
+ }
+ {% endif -%}
+
+ }
+ {% endif %}
+
indexVectors->syncGPU();
}
diff --git a/tests/cuda/codegen/MicroBenchmarkGpuLbm.py b/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
index 298727b46c428384eeef7f755e8bfe4881d53d60..45bdc303c9f9983c7c4120748e4b1fba1f23ecf7 100644
--- a/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
+++ b/tests/cuda/codegen/MicroBenchmarkGpuLbm.py
@@ -1,5 +1,5 @@
import pystencils as ps
-from lbmpy.updatekernels import create_stream_pull_only_kernel
+from lbmpy.updatekernels import create_stream_only_kernel
from lbmpy.stencils import get_stencil
from pystencils_walberla import CodeGeneration, generate_sweep
@@ -8,15 +8,14 @@ with CodeGeneration() as ctx:
dtype = 'float64' if ctx.double_accuracy else 'float32'
# Copy sweep
- src, dst = ps.fields("src({f_size}), dst({f_size}) : {dtype}[3D]".format(dtype=dtype, f_size=f_size),
- layout='fzyx')
+ src, dst = ps.fields(f"src({f_size}), dst({f_size}) : {dtype}[3D]", layout='fzyx')
+
copy_only = [ps.Assignment(dst(i), src(i)) for i in range(f_size)]
generate_sweep(ctx, 'MicroBenchmarkCopyKernel', copy_only,
target='gpu', gpu_indexing_params={'block_size': (128, 1, 1)})
# Stream-only sweep
stencil = get_stencil("D3Q19")
- stream_only = create_stream_pull_only_kernel(stencil, src_field_name='src', dst_field_name='dst',
- generic_field_type=dtype, generic_layout='fzyx')
+ stream_only = create_stream_only_kernel(stencil, src_field=src, dst_field=dst)
generate_sweep(ctx, 'MicroBenchmarkStreamKernel', stream_only,
target='gpu', gpu_indexing_params={'block_size': (128, 1, 1)})