Skip to content

GitLab

Commit 03b9f95f authored by Christoph Schwarzmeier's avatar Christoph Schwarzmeier
Browse files

Merge branch 'FixAllenCahn' into 'master' 

Update Allen Cahn model

See merge request walberla/walberla!473
parents 8f7272ba 076ebb5c
Hide whitespace changes
Inline Side-by-side
apps/benchmarks/PhaseFieldAllenCahn/benchmark_multiphase.cpp
View file @ 03b9f95f
......@@ -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
......
apps/benchmarks/PhaseFieldAllenCahn/multiphase_codegen.py
View file @ 03b9f95f
......@@ -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)
......
apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt
View file @ 03b9f95f
......@@ -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)
apps/showcases/PhaseFieldAllenCahn/CPU/CalculateNormals.cpp deleted 100644 → 0
View file @ 8f7272ba
//======================================================================================================================
//
// 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
apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp
View file @ 03b9f95f
......@@ -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
apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h
View file @ 03b9f95f
......@@ -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
apps/showcases/PhaseFieldAllenCahn/CPU/contact.cpp deleted 100644 → 0
View file @ 8f7272ba
//======================================================================================================================
//
// 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)