From 5fcf3d9c4413d7e4ede97dc847ed5a113d469fb0 Mon Sep 17 00:00:00 2001 From: Markus Holzer Date: Fri, 15 Jul 2022 16:33:48 +0200 Subject: [PATCH 01/16] Added dropplet --- .../CPU/InitializerFunctions.cpp | 52 ++++- .../CPU/InitializerFunctions.h | 13 +- .../PhaseFieldAllenCahn/CPU/multiphase.cpp | 19 +- .../CPU/multiphase_codegen.py | 3 +- .../CPU/multiphase_drop.py | 86 ++++++++ .../PhaseFieldAllenCahn/CPU/taylor_bubble.py | 189 ++++++++++++++++++ 6 files changed, 354 insertions(+), 8 deletions(-) create mode 100755 apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py create mode 100755 apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index d7cb146a0..951babf4c 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -31,15 +31,19 @@ namespace walberla { using PhaseField_T = GhostLayerField< real_t, 1 >; using FlagField_T = FlagField< uint8_t >; +using VelocityField_T = GhostLayerField< real_t, 3 >; -void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, +void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, + BlockDataID phaseFieldID, BlockDataID velocityFieldID, const real_t R = real_c(10.0), const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0), - const bool bubble = true, const real_t W = real_c(5.0)) + const bool bubble = true, const real_t W = real_c(5.0), + const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0))) { for (auto& block : *blocks) { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); + auto velocityField = block.getData< VelocityField_T >(velocityFieldID); // clang-format off WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); @@ -47,12 +51,54 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); if (bubble) phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); - else phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); + else + { + phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); + + if (Ri < R + W) + { + velocityField->get(x, y, z, 0) = initialVelocity[0]; + velocityField->get(x, y, z, 1) = initialVelocity[1]; + velocityField->get(x, y, z, 2) = initialVelocity[2]; + } + } ) // clang-format on } } +void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W, + real_t poolDepth) +{ + for (auto& block : *blocks) + { + auto phaseField = block.getData< PhaseField_T >(phaseFieldID); + + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + + phaseField->get(x, y, z) += 0.5 - 0.5 * tanh(2.0 * (real_c(globalCell[1]) - poolDepth) / W); + if (phaseField->get(x, y, z) > 1) { phaseField->get(x, y, z) = 1; } + }) + } +} + +void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID densityFieldID, + real_t GravitationalAcceleration, real_t poolDepth) +{ + for (auto& block : *blocks) + { + auto densityField = block.getData< PhaseField_T >(densityFieldID); + // clang-format off + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(densityField, Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + real_t pressure = 3.0 * GravitationalAcceleration * (globalCell[1] - poolDepth); + densityField->get(x, y, z) += pressure;) + // clang-format on + } +} + void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0), const real_t Donut_x0 = real_c(40.0)) { diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h index 585639ac9..6d403d124 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h @@ -31,8 +31,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_sphere(const shared_ptr< StructuredBlockStorage >& blocks, + BlockDataID phaseFieldID, BlockDataID velocityFieldID, + real_t R, + Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5, + const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0))); + +void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W, + real_t poolDepth); + +void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID densityFieldID, + real_t GravitationalAcceleration, real_t poolDepth); void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t D = 5, real_t H = 2, real_t DT = 20, real_t Donut_x0 = 40); diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index 8e0f2cb62..969732d61 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -92,6 +92,7 @@ int main(int argc, char** argv) field::addToStorage< PdfField_hydro_T >(blocks, "LB velocity field", real_c(0.0), field::fzyx); BlockDataID velocity_field_ID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx); + BlockDataID density_field_ID = field::addToStorage< PhaseField_T >(blocks, "density", real_c(1.0), field::fzyx); BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx); BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field"); @@ -114,9 +115,23 @@ int main(int argc, char** argv) const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint"); const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0)); const bool bubble = bubbleParameters.getParameter< bool >("bubble", true); - initPhaseField_sphere(blocks, phase_field_ID, bubbleRadius, bubbleMidPoint, bubble); + initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble); } else if (scenario == 2) { initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); } + else if (scenario == 3) + { + auto dropParameters = config->getOneBlock("Drop"); + const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("dropMidPoint"); + const real_t dropRadius = dropParameters.getParameter< real_t >("dropRadius"); + const real_t poolDepth = dropParameters.getParameter< real_t >("poolDepth"); + const Vector3< real_t > impact_velocity = + dropParameters.getParameter< Vector3< real_t > >("impact_velocity"); + init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth); + + initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false, + interface_thickness, impact_velocity); + initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth); + } WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done") ///////////////// @@ -125,7 +140,7 @@ int main(int argc, char** argv) pystencils::initialize_phase_field_distributions init_h(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID, interface_thickness); - pystencils::initialize_velocity_based_distributions init_g(hydrodynamic_PDFs_ID, velocity_field_ID); + pystencils::initialize_velocity_based_distributions init_g(density_field_ID,hydrodynamic_PDFs_ID, velocity_field_ID); pystencils::phase_field_LB_step phase_field_LB_step(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID, mobility, interface_thickness); diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py index 1b5113454..1e2031057 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py @@ -37,6 +37,7 @@ with CodeGeneration() as ctx: # velocity field u = fields(f"vel_field({stencil_hydro.D}): {field_type}[{stencil_hydro.D}D]", layout='fzyx') + density_field = fields(f"density_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx') # phase-field C = fields(f"phase_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx') C_tmp = fields(f"phase_field_tmp: {field_type}[{stencil_hydro.D}D]", layout='fzyx') @@ -87,7 +88,7 @@ with CodeGeneration() as ctx: # create the kernels for the initialization of the g and h field h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters) - g_updates = initializer_kernel_hydro_lb(method_hydro, 1.0, u, g) + g_updates = initializer_kernel_hydro_lb(method_hydro, density_field, u, g) force_h = interface_tracking_force(C, stencil_phase, parameters) hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py new file mode 100755 index 000000000..071ecd940 --- /dev/null +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -0,0 +1,86 @@ +import waLBerla as wlb + + +class Scenario: + def __init__(self): + # output frequencies + self.dropRadius = 20 + + self.surface_tension = 0.0001 + self.gravitational_acceleration = -2.48e-10 + self.impact_velocity = (0.02574, -0.0474, 0) + + self.density_ratio = 1000 + self.kin_visc_ratio = 100 + self.density_heavy = 1 + self.density_light = 0.0012 # self.density_heavy / self.density_ratio + + self.relaxation_time_heavy = 0.003 # 1 / 1.9889 + self.relaxation_time_light = 0.03 # 1 / self.kin_visc_ratio * 1 / self.density_ratio * (self.relaxation_time_heavy - 0.5) + 0.5 + + self.vtkWriteFrequency = 10000 + self.meshWriteFrequency = 100 + + # simulation parameters + self.timesteps = 15000 + self.cells = (self.dropRadius * 2, self.dropRadius * 2, self.dropRadius * 4) + self.blocks = (10, 10, 5) + self.periodic = (1, 0, 1) + self.size = (self.cells[0] * self.blocks[0], + self.cells[1] * self.blocks[1], + self.cells[2] * self.blocks[2]) + self.poolDepth = self.size[1] * 0.5 + + # physical parameters + self.dropMidPoint = (self.size[0] / 2, self.poolDepth + self.dropRadius, self.size[2] / 2) + + self.interface_thickness = 5 + self.mobility = 0.1 + + self.scenario = 3 # 1 rising bubble, 2 RTI + + self.counter = 0 + self.yPositions = [] + + @wlb.member_callback + def config(self): + return { + 'DomainSetup': { + 'blocks': self.blocks, + 'cellsPerBlock': self.cells, + 'periodic': self.periodic, + }, + 'Parameters': { + 'timesteps': self.timesteps, + 'vtkWriteFrequency': self.vtkWriteFrequency, + 'meshWriteFrequency': self.meshWriteFrequency, + 'remainingTimeLoggerFrequency': 10.0, + 'scenario': self.scenario, + }, + 'PhysicalParameters': { + 'density_liquid': self.density_heavy, + 'density_gas': self.density_light, + 'surface_tension': self.surface_tension, + 'mobility': self.mobility, + 'gravitational_acceleration': self.gravitational_acceleration, + 'relaxation_time_liquid': self.relaxation_time_heavy, + 'relaxation_time_gas': self.relaxation_time_light, + 'interface_thickness': self.interface_thickness, + }, + 'Boundaries': { + 'Border': [ + {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'}, + ] + }, + 'Drop': { + 'dropMidPoint': self.dropMidPoint, + 'dropRadius': self.dropRadius, + 'poolDepth': self.poolDepth, + 'impact_velocity': self.impact_velocity, + }, + } + + +scenarios = wlb.ScenarioManager() +scenarios.add(Scenario()) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py new file mode 100755 index 000000000..4f1183581 --- /dev/null +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py @@ -0,0 +1,189 @@ +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 + + +class Scenario: + def __init__(self): + self.diameter = 64 + + self.bond_number = 100 + self.morton_number = 0.015 + + self.relaxation_time_liquid = 0.56797 + self.kinematic_viscosity_liquid = 1 / 3 * (self.relaxation_time_liquid - 0.5) + + self.density_liquid = 1.0 + self.density_gas = self.density_liquid / 744 + + self.dynamic_viscosity_liquid = self.density_liquid * self.kinematic_viscosity_liquid + self.dynamic_viscosity_gas = self.dynamic_viscosity_liquid / 4236 + self.kinematic_viscosity_gas = self.dynamic_viscosity_gas / self.density_gas + self.relaxation_time_gas = 3 * self.kinematic_viscosity_gas + 0.5 + + self.gravitational_acceleration = - (self.bond_number**3 * self.dynamic_viscosity_liquid**4 / self.density_liquid**4 / + self.morton_number / self.diameter**6)**0.5 + + self.surface_tension = self.density_liquid * abs(self.gravitational_acceleration) * self.diameter**2 / self.bond_number + + self.interface_width = 5 + self.mobility = 0.1 + + self.timesteps = 20 * (self.diameter / abs(self.gravitational_acceleration))**0.5 + + # output frequencies + self.vtkWriteFrequency = self.timesteps // 20 + self.dbWriteFrequency = self.timesteps // 100 + self.meshWriteFrequency = self.timesteps // 20 + + # simulation parameters + self.size = (self.diameter, self.diameter * 10, self.diameter) + self.blocks = (1, 1, 1) + self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) + self.periodic = (0, 0, 0) + self.inner_radius = self.diameter + + self.center_x = self.size[0] / 2 + self.center_y = self.size[1] / 2 + self.center_z = self.size[2] / 2 + + self.scenario = 4 # 1 rising bubble 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 + + 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}_C_{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 + + self.config_dict = self.config() + + @wlb.member_callback + def config(self): + return { + 'DomainSetup': { + 'blocks': self.blocks, + 'domainSize': self.size, + 'cellsPerBlock': self.cells, + 'diameter': self.diameter, + 'periodic': self.periodic, + 'inner_radius': self.inner_radius, + 'ratio': self.ratio, + 'start_transition': self.start_transition, + 'length_transition': self.length_transition, + 'eccentricity_or_pipe_ration': self.eccentricity_or_pipe_ratio, + 'tube': True + }, + 'Parameters': { + 'timesteps': self.timesteps, + 'vtkWriteFrequency': self.vtkWriteFrequency, + 'dbWriteFrequency': self.dbWriteFrequency, + 'meshWriteFrequency': self.meshWriteFrequency, + 'remainingTimeLoggerFrequency': 60.0, + 'scenario': self.scenario, + }, + 'PhysicalParameters': { + 'density_liquid': self.density_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 - 0.5, + 'relaxation_time_gas': self.relaxation_time_gas - 0.5, + 'interface_thickness': self.interface_width + }, + 'Boundaries': { + 'Border': [ + {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'W', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'E', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'}, + {'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'}, + ], + }, + 'TaylorBubble': { + 'BubbleRadius': 0.75 * 0.5 * self.diameter, + 'InitialHeight': 1 * self.diameter, + 'Length': 3 * self.diameter, + } + } + + # @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"] + if t % self.dbWriteFrequency == 0: + wlb_field = wlb.field.gather(blocks, 'phase', makeSlice[:, :, self.size[2] // 2]) + if wlb_field: + data = {'timestep': t} + data.update(self.config_dict['Parameters']) + data.update(self.config_dict['DomainSetup']) + data.update(self.config_dict['PhysicalParameters']) + data.update(self.config_dict['TaylorBubble']) + data.update(kwargs) + + phase_field = np.asarray(wlb_field).squeeze() + location_of_gas = np.where(phase_field < 0.5) + center_of_mass = np.mean(location_of_gas, axis=1) + + assert np.isfinite(np.sum(phase_field)), "NaN detected in the phasefield" + + self.yPositions.append(center_of_mass[1]) + if len(self.yPositions) > 1: + speed = self.yPositions[-1] - self.yPositions[-2] + else: + speed = 0 + + data['center_of_mass_x'] = center_of_mass[0] + data['center_of_mass_y'] = center_of_mass[1] + # data['center_of_mass_z'] = center_of_mass[2] + + data['xCells'] = self.size[0] + data['yCells'] = self.size[1] + data['zCells'] = self.size[2] + + data['rising_velocity'] = speed + + data['StencilHydroLB'] = kwargs["stencil_hydro"] + data['StencilPhaseLB'] = kwargs["stencil_phase"] + sequenceValuesToScalars(data) + + df = pd.DataFrame.from_records([data]) + if not os.path.isfile(self.csv_file): + df.to_csv(self.csv_file, index=False) + else: + df.to_csv(self.csv_file, index=False, mode='a', header=False) + + +scenarios = wlb.ScenarioManager() +scenarios.add(Scenario()) -- GitLab From 3af05b5173fa7afd96b42270b7773a0431451247 Mon Sep 17 00:00:00 2001 From: Markus Holzer Date: Sat, 16 Jul 2022 13:00:46 +0200 Subject: [PATCH 02/16] Added Taylor bubble --- .../PhaseFieldAllenCahn/CPU/CMakeLists.txt | 2 +- .../CPU/InitializerFunctions.cpp | 45 ++++++++++ .../CPU/InitializerFunctions.h | 3 + .../PhaseFieldAllenCahn/CPU/multiphase.cpp | 90 ++++++++++++++++--- .../PhaseFieldAllenCahn/CPU/taylor_bubble.py | 8 +- .../PhaseFieldAllenCahn/GPU/CMakeLists.txt | 2 +- 6 files changed, 131 insertions(+), 19 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt index abce2c997..383c9d31a 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/CMakeLists.txt @@ -18,4 +18,4 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenCPU waLBerla_add_executable(NAME multiphaseCPU FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp multiphase_codegen.py - DEPENDS blockforest core field postprocessing python_coupling lbm geometry timeloop gui PhaseFieldCodeGenCPU) + DEPENDS blockforest core field postprocessing python_coupling lbm geometry timeloop PhaseFieldCodeGenCPU) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index 951babf4c..707a04d2d 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -121,6 +121,51 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc } } +void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& phaseFieldID, + real_t R, real_t H, real_t L, real_t W) +{ + auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0); + auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0); + + for (auto& block : *blocks) + { + PhaseField_T* phaseField = block.getData< PhaseField_T >(phaseFieldID); + + // initialize top and bottom walls of cylinder + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + + const real_t Ri = real_c(sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2))); + + const real_t cylinderMidY = real_c(H) + real_c(0.5) * real_c(L); + + const real_t distance = std::abs(real_c(globalCell[1]) - cylinderMidY); + const real_t lHalf = real_c(0.5) * real_c(L); + + // cylinder side walls + if (real_c(globalCell[1]) >= H - W && real_c(globalCell[1]) < H + L + W) + { + phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - R) / W); + + // cylinder top and bottom walls (must be added to not overwrite the side wall initialization) + if (Ri < R + W) + { + phaseField->get(x, y, z) += real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (distance - lHalf) / W); + + // limit maximum values of phase field + if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); } + if (phaseField->get(x, y, z) < real_c(0)) { phaseField->get(x, y, z) = real_c(0); } + } + } + else + { + phaseField->get(x, y, z) = real_c(1); + } + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ + } +} + void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R, real_t W = real_c(5.0)) { diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h index 6d403d124..ae99d3ab0 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h @@ -46,6 +46,9 @@ void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& block void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t D = 5, real_t H = 2, real_t DT = 20, real_t Donut_x0 = 40); +void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& phaseFieldID, + real_t R, real_t H, real_t L, real_t W); + void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R, real_t W = 5); diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index 969732d61..48f7e9011 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -30,9 +30,12 @@ #include "field/vtk/VTKWriter.h" #include "geometry/InitBoundaryHandling.h" +#include "geometry/mesh/TriangleMeshIO.h" #include "lbm/PerformanceEvaluation.h" +#include "postprocessing/FieldToSurfaceMesh.h" + #include "python_coupling/CreateConfig.h" #include "python_coupling/PythonCallback.h" @@ -69,6 +72,7 @@ int main(int argc, char** argv) ////////////////////////// auto domainSetup = config->getOneBlock("DomainSetup"); + Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize"); const bool tube = domainSetup.getParameter< bool >("tube", true); //////////////////////////////////////// @@ -109,17 +113,21 @@ int main(int argc, char** argv) const real_t interface_thickness = physical_parameters.getParameter< real_t >("interface_thickness"); WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field") - if (scenario == 1) - { - auto bubbleParameters = config->getOneBlock("Bubble"); - const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint"); - const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0)); - const bool bubble = bubbleParameters.getParameter< bool >("bubble", true); - initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble); - } - else if (scenario == 2) { initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); } - else if (scenario == 3) + switch (scenario) { + case 1: { + auto bubbleParameters = config->getOneBlock("Bubble"); + const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint"); + const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0)); + const bool bubble = bubbleParameters.getParameter< bool >("bubble", true); + initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble); + break; + } + case 2: { + initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); + break; + } + case 3:{ auto dropParameters = config->getOneBlock("Drop"); const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("dropMidPoint"); const real_t dropRadius = dropParameters.getParameter< real_t >("dropRadius"); @@ -131,6 +139,19 @@ int main(int argc, char** argv) initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false, interface_thickness, impact_velocity); initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth); + break; + } + case 4:{ + auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble"); + const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("BubbleRadius"); + const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("InitialHeight"); + const real_t Length = TaylorBubbleParameters.getParameter< real_t >("Length"); + init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length, + real_c(interface_thickness)); + break; + } + default: + WALBERLA_ABORT("Scenario is not defined.") } WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done") @@ -195,6 +216,29 @@ int main(int argc, char** argv) if (boundariesConfig) { geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig); + if (tube) + { + // initialize cylindrical domain walls + const Vector3< real_t > domainCylinderBottomEnd = + Vector3< real_t >(real_c(domainSize[0]) * real_c(0.5), real_c(0), real_c(domainSize[2]) * real_c(0.5)); + const Vector3< real_t > domainCylinderTopEnd = Vector3< real_t >( + real_c(domainSize[0]) * real_c(0.5), real_c(domainSize[1]), real_c(domainSize[2]) * real_c(0.5)); + const real_t radius = real_c(domainSize[0]) * real_c(0.5); + const geometry::Cylinder cylinderTube(domainCylinderBottomEnd, domainCylinderTopEnd, radius); + + for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt) + { + FlagField_T* flagField = blockIt->template getData< FlagField_T >(flagFieldID); + auto wallFlag = flagField->getOrRegisterFlag(wallFlagUID); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, { + Cell globalCell = Cell(x, y, z); + blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt); + const Vector3< real_t > globalPoint = blocks->getCellCenter(globalCell); + + if (!geometry::contains(cylinderTube, globalPoint)) { addFlag(flagField->get(x, y, z), wallFlag); } + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ + } + } geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID); } @@ -263,6 +307,32 @@ int main(int argc, char** argv) }, "Python callback"); + int meshWriteFrequency = parameters.getParameter< int >("meshWriteFrequency", 0); + int counter = 0; + int targetRank = 0; + if (meshWriteFrequency > 0) + { + timeloop.addFuncAfterTimeStep( + [&]() { + if (timeloop.getCurrentTimeStep() % uint_t(meshWriteFrequency) == 0) + { + auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field_ID, 0.5, 0, true, + targetRank, MPI_COMM_WORLD); + WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank) + { + std::string path = ""; + std::ostringstream out; + out << std::internal << std::setfill('0') << std::setw(6) << counter; + geometry::writeMesh( + path + "taylor_bubble_D_" + std::to_string(domainSize[0]) + "_" + out.str() + ".obj", *mesh); + counter++; + } + } + }, + "Mesh writer"); + } + + // remaining time logger timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency), diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py index 4f1183581..32ce51f8a 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py @@ -30,7 +30,7 @@ class Scenario: self.relaxation_time_gas = 3 * self.kinematic_viscosity_gas + 0.5 self.gravitational_acceleration = - (self.bond_number**3 * self.dynamic_viscosity_liquid**4 / self.density_liquid**4 / - self.morton_number / self.diameter**6)**0.5 + self.morton_number / self.diameter**6)**0.5 self.surface_tension = self.density_liquid * abs(self.gravitational_acceleration) * self.diameter**2 / self.bond_number @@ -133,12 +133,6 @@ class Scenario: } } - # @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"] diff --git a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt index ac4feda33..9116f0b19 100644 --- a/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt +++ b/apps/showcases/PhaseFieldAllenCahn/GPU/CMakeLists.txt @@ -18,4 +18,4 @@ waLBerla_generate_target_from_python(NAME PhaseFieldCodeGenGPU waLBerla_add_executable(NAME multiphaseGPU FILES multiphase.cpp PythonExports.cpp InitializerFunctions.cpp util.cpp multiphase_codegen.py - DEPENDS blockforest core cuda field postprocessing python_coupling lbm geometry timeloop gui PhaseFieldCodeGenGPU) + DEPENDS blockforest core cuda field postprocessing python_coupling lbm geometry timeloop PhaseFieldCodeGenGPU) -- GitLab From 4a532034cfcb8fd4c2a4eed589381ab6572c5004 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Wed, 20 Jul 2022 16:08:33 +0200 Subject: [PATCH 03/16] Fix floating point conversion warnings --- .../CPU/InitializerFunctions.cpp | 111 ++++++++++-------- 1 file changed, 62 insertions(+), 49 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index 707a04d2d..32a667456 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -29,20 +29,19 @@ namespace walberla { -using PhaseField_T = GhostLayerField< real_t, 1 >; -using FlagField_T = FlagField< uint8_t >; +using PhaseField_T = GhostLayerField< real_t, 1 >; +using FlagField_T = FlagField< uint8_t >; using VelocityField_T = GhostLayerField< real_t, 3 >; -void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, - BlockDataID phaseFieldID, BlockDataID velocityFieldID, - const real_t R = real_c(10.0), +void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, + BlockDataID velocityFieldID, const real_t R = real_c(10.0), const Vector3< real_t > bubbleMidPoint = Vector3< real_t >(0.0, 0.0, 0.0), const bool bubble = true, const real_t W = real_c(5.0), const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0))) { for (auto& block : *blocks) { - auto phaseField = block.getData< PhaseField_T >(phaseFieldID); + auto phaseField = block.getData< PhaseField_T >(phaseFieldID); auto velocityField = block.getData< VelocityField_T >(velocityFieldID); // clang-format off WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell; @@ -78,8 +77,9 @@ void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, Blo Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - phaseField->get(x, y, z) += 0.5 - 0.5 * tanh(2.0 * (real_c(globalCell[1]) - poolDepth) / W); - if (phaseField->get(x, y, z) > 1) { phaseField->get(x, y, z) = 1; } + phaseField->get(x, y, z) += + real_c(0.5) - real_c(0.5) * real_c(tanh(real_c(2.0) * (real_c(globalCell[1]) - poolDepth) / W)); + if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); } }) } } @@ -93,14 +93,15 @@ void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& block // clang-format off WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(densityField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t pressure = 3.0 * GravitationalAcceleration * (globalCell[1] - poolDepth); + real_t pressure = real_c(3.0) * GravitationalAcceleration * (real_c(globalCell[1]) - poolDepth); densityField->get(x, y, z) += pressure;) // clang-format on } } void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, - const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0), const real_t Donut_x0 = real_c(40.0)) + const real_t D = real_c(5.0), const real_t H = real_c(2.0), const real_t DT = real_c(20.0), + const real_t Donut_x0 = real_c(40.0)) { auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0); auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0); @@ -111,7 +112,10 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t Ri = D * real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2))); + real_t Ri = + D * + real_c(sqrt(pow(H, 2) - + pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2))); real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx)); if (shifter < 0) shifter = shifter + 2 * math::pi; @@ -158,10 +162,7 @@ void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& bl if (phaseField->get(x, y, z) < real_c(0)) { phaseField->get(x, y, z) = real_c(0); } } } - else - { - phaseField->get(x, y, z) = real_c(1); - } + else { phaseField->get(x, y, z) = real_c(1); } }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -183,8 +184,8 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block // distance in between the bubbles real_t dist = real_c(4.0); - auto nx = uint_c(floor(X / (D + dist * W))); - auto nz = uint_c(floor(Z / (D + dist * W))); + auto nx = uint_c(floor(X / (D + dist * W))); + auto nz = uint_c(floor(Z / (D + dist * W))); // fluctuation of the bubble radii std::vector< std::vector< real_t > > fluctuation_radius(nx, std::vector< real_t >(nz, 0.0)); @@ -209,16 +210,22 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block for (unsigned int i = 0; i < nx; ++i) { for (unsigned int j = 0; j < nz; ++j) { - bubbleMidPoint[0] = real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j]; + bubbleMidPoint[0] = + real_c(i + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j]; bubbleMidPoint[1] = R + W + 4; - bubbleMidPoint[2] = real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j]; - - real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + - (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + - (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); - if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) && - real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W)) - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W); + bubbleMidPoint[2] = + real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j]; + + real_t Ri = + sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + + (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + + (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); + if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && + real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) && + real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && + real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W)) + phaseField->get(x, y, z) = + real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W); if (real_c(globalCell[0]) > real_c(nx) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0); if (real_c(globalCell[2]) > real_c(nz) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0); @@ -226,7 +233,8 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block } if (real_c(globalCell[1]) > gas_top) { - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W); + phaseField->get(x, y, z) = + real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W); }) } } @@ -234,9 +242,9 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, const real_t W = real_c(5.0), const bool pipe = true) { - auto X = real_c(blocks->getDomainCellBB().xMax()); - auto Z = real_c(blocks->getDomainCellBB().zMax()); - auto halfY = real_c(blocks->getDomainCellBB().yMax()) / real_c(2.0); + auto X = real_c(blocks->getDomainCellBB().xMax()); + auto Z = real_c(blocks->getDomainCellBB().zMax()); + auto halfY = real_c(blocks->getDomainCellBB().yMax()) / real_c(2.0); auto perturbation = real_c(0.05); if (pipe) @@ -247,9 +255,10 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) + - (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)); + (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)); if (R > X) R = X; real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X); - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)));) + phaseField->get(x, y, z) = + real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)));) } } else @@ -260,8 +269,10 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t tmp = perturbation * X * - (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X)); - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)));) + (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + + cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X)); + phaseField->get(x, y, z) = + real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)));) } } } @@ -284,22 +295,24 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl { auto flagField = block.template getData< FlagField_T >(flagFieldID); auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, Cell globalCell; - blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t R1; - if (real_c(globalCell[1]) <= start_transition) { - R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); - } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) { + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( + flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + real_t R1; if (real_c(globalCell[1]) <= start_transition) { + R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + } else if (real_c(globalCell[1]) > start_transition && + real_c(globalCell[1]) < start_transition + length_transition) { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); - R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) + + R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); } else { R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); } - real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag); if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);) } @@ -308,7 +321,7 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl { auto Mx = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0); auto Mz = real_c(blocks->getDomainCellBB().zMax()) / real_c(2.0); - + auto R_outer = real_c(blocks->getDomainCellBB().xMax()) / real_c(2.0) + real_c(1.0); auto const shift = real_c(eccentricity * R_in); @@ -322,15 +335,15 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); if (real_c(globalCell[1]) <= start_transition) { R_tmp = R_in; - } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) { + } else if (real_c(globalCell[1]) > start_transition && + real_c(globalCell[1]) < start_transition + length_transition) { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); - R_tmp = R_in + shift_tmp; - } else { - R_tmp = R_in + shift; - } + R_tmp = R_in + shift_tmp; + } else { R_tmp = R_in + shift; } - real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag); if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);) } -- GitLab From 185d7d7d5f0fd50ebcc2c93585e947a8efbab303 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Wed, 20 Jul 2022 16:13:22 +0200 Subject: [PATCH 04/16] Improve formatting of iterator macros --- .../CPU/InitializerFunctions.cpp | 122 +++++++++++------- 1 file changed, 73 insertions(+), 49 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index 32a667456..bdad06c50 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -43,16 +43,17 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); auto velocityField = block.getData< VelocityField_T >(velocityFieldID); - // clang-format off - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, Cell globalCell; + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); - if (bubble) phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); + if (bubble) + phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); else { - phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); + phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); if (Ri < R + W) { @@ -61,8 +62,7 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B velocityField->get(x, y, z, 2) = initialVelocity[2]; } } - ) - // clang-format on + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -80,7 +80,7 @@ void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, Blo phaseField->get(x, y, z) += real_c(0.5) - real_c(0.5) * real_c(tanh(real_c(2.0) * (real_c(globalCell[1]) - poolDepth) / W)); if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); } - }) + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -90,12 +90,12 @@ void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& block for (auto& block : *blocks) { auto densityField = block.getData< PhaseField_T >(densityFieldID); - // clang-format off - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(densityField, Cell globalCell; - blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t pressure = real_c(3.0) * GravitationalAcceleration * (real_c(globalCell[1]) - poolDepth); - densityField->get(x, y, z) += pressure;) - // clang-format on + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(densityField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + real_t pressure = real_c(3.0) * GravitationalAcceleration * (real_c(globalCell[1]) - poolDepth); + densityField->get(x, y, z) += pressure; + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -109,19 +109,23 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc for (auto& block : *blocks) { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t Ri = D * real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2))); - real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx)); + real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx)); if (shifter < 0) shifter = shifter + 2 * math::pi; - if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) { + if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) + { phaseField->get(x, y, z) = 0.0; - } else { phaseField->get(x, y, z) = real_c(1.0); }) + } + else { phaseField->get(x, y, z) = real_c(1.0); } + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -205,9 +209,11 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block for (auto& block : *blocks) { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - for (unsigned int i = 0; i < nx; ++i) { + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + for (unsigned int i = 0; i < nx; ++i) + { for (unsigned int j = 0; j < nz; ++j) { bubbleMidPoint[0] = @@ -232,10 +238,12 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block } } - if (real_c(globalCell[1]) > gas_top) { + if (real_c(globalCell[1]) > gas_top) + { phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W); - }) + } + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -252,13 +260,16 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc for (auto& block : *blocks) { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) + - (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)); - if (R > X) R = X; real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) + + (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)); + if (R > X) R = X; + real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X); phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)));) + real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0))); + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } else @@ -266,13 +277,15 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc for (auto& block : *blocks) { auto phaseField = block.getData< PhaseField_T >(phaseFieldID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - phaseField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t tmp = perturbation * X * (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X)); phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)));) + real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0))); + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } } @@ -295,26 +308,34 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl { auto flagField = block.template getData< FlagField_T >(flagFieldID); auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t R1; if (real_c(globalCell[1]) <= start_transition) { + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + real_t R1; + if (real_c(globalCell[1]) <= start_transition) + { R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); - } else if (real_c(globalCell[1]) > start_transition && - real_c(globalCell[1]) < start_transition + length_transition) { + } + else if (real_c(globalCell[1]) > start_transition && + real_c(globalCell[1]) < start_transition + length_transition) + { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); - } else { + } + else + { R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); } - real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag); - if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);) + if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag); + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } else @@ -331,21 +352,24 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl { auto flagField = block.template getData< FlagField_T >(flagFieldID); auto boundaryFlag = flagField->getOrRegisterFlag(boundaryFlagUID); - WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ( - flagField, Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - if (real_c(globalCell[1]) <= start_transition) { - R_tmp = R_in; - } else if (real_c(globalCell[1]) > start_transition && - real_c(globalCell[1]) < start_transition + length_transition) { + WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, { + Cell globalCell; + blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); + if (real_c(globalCell[1]) <= start_transition) { R_tmp = R_in; } + else if (real_c(globalCell[1]) > start_transition && + real_c(globalCell[1]) < start_transition + length_transition) + { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); R_tmp = R_in + shift_tmp; - } else { R_tmp = R_in + shift; } + } + else { R_tmp = R_in + shift; } real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag); - if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag);) + if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag); + }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } } -- GitLab From 1b06aa554989f25949813415697b1bcb0990d7a0 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Wed, 20 Jul 2022 16:14:33 +0200 Subject: [PATCH 05/16] Remove unused includes --- .../showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index bdad06c50..b0b83abb9 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -23,9 +23,6 @@ #include "field/FlagField.h" #include "field/communication/PackInfo.h" -#include "field/vtk/VTKWriter.h" - -#include "python_coupling/DictWrapper.h" namespace walberla { -- GitLab From 5daff681e4fdebadc288eb1b748114062d05339f Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Wed, 20 Jul 2022 16:46:18 +0200 Subject: [PATCH 06/16] Update parameters in Taylor bubble showcase --- .../PhaseFieldAllenCahn/CPU/multiphase_codegen.py | 2 +- .../{taylor_bubble.py => multiphase_taylor_bubble.py} | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) rename apps/showcases/PhaseFieldAllenCahn/CPU/{taylor_bubble.py => multiphase_taylor_bubble.py} (96%) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py index 1e2031057..0b117fc23 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py @@ -28,7 +28,7 @@ with CodeGeneration() as ctx: stencil_hydro = LBStencil(Stencil.D3Q27) assert (stencil_phase.D == stencil_hydro.D) - # In the codegneration skript we only need the symbolic parameters + # In the code-generation skript, we only need the symbolic parameters parameters = AllenCahnParameters(0, 0, 0, 0, 0) ######################## diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py similarity index 96% rename from apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py rename to apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py index 32ce51f8a..0f1450dbe 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/taylor_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py @@ -18,7 +18,7 @@ class Scenario: self.bond_number = 100 self.morton_number = 0.015 - self.relaxation_time_liquid = 0.56797 + self.relaxation_time_liquid = 1 / 1.76 self.kinematic_viscosity_liquid = 1 / 3 * (self.relaxation_time_liquid - 0.5) self.density_liquid = 1.0 @@ -34,14 +34,14 @@ class Scenario: self.surface_tension = self.density_liquid * abs(self.gravitational_acceleration) * self.diameter**2 / self.bond_number - self.interface_width = 5 - self.mobility = 0.1 + self.interface_width = 3 + self.mobility = 0.08 - self.timesteps = 20 * (self.diameter / abs(self.gravitational_acceleration))**0.5 + self.timesteps = 20 / (abs(self.gravitational_acceleration) / self.diameter)**0.5 # output frequencies self.vtkWriteFrequency = self.timesteps // 20 - self.dbWriteFrequency = self.timesteps // 100 + self.dbWriteFrequency = self.timesteps // 20 self.meshWriteFrequency = self.timesteps // 20 # simulation parameters -- GitLab From 305a813c18ac52bea2965014068e50d52646d059 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Wed, 20 Jul 2022 16:47:33 +0200 Subject: [PATCH 07/16] Fix more floating point conversion warnings --- .../CPU/InitializerFunctions.cpp | 65 ++++++++++--------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp index b0b83abb9..da5d24f34 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.cpp @@ -43,11 +43,11 @@ void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, B WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t Ri = sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + - (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + - (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); - if (bubble) - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); + real_t Ri = + real_c(sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + + (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + + (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]))); + if (bubble) { phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); } else { phaseField->get(x, y, z) = real_c(0.5) - real_c(0.5) * real_c(tanh(2.0 * (Ri - R) / W)); @@ -115,9 +115,9 @@ void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, Bloc real_c(sqrt(pow(H, 2) - pow(DT - sqrt(pow(real_c(globalCell[0]) - Mx, 2) + pow(real_c(globalCell[2]) - Mz, 2)), 2))); - real_t shifter = atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx)); + real_t shifter = real_c(atan2((real_c(globalCell[2]) - Mz), (real_c(globalCell[0]) - Mx))); if (shifter < 0) shifter = shifter + 2 * math::pi; - if ((real_c(globalCell[1]) < Donut_x0 + Ri * sin(shifter / 2.0)) && (real_c(globalCell[1]) > Donut_x0 - Ri)) + if ((real_c(globalCell[1]) < Donut_x0 + Ri * real_c(sin(shifter / 2.0))) && (real_c(globalCell[1]) > Donut_x0 - Ri)) { phaseField->get(x, y, z) = 0.0; } @@ -151,12 +151,13 @@ void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& bl // cylinder side walls if (real_c(globalCell[1]) >= H - W && real_c(globalCell[1]) < H + L + W) { - phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - R) / W); + phaseField->get(x, y, z) = real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (Ri - R) / W)); // cylinder top and bottom walls (must be added to not overwrite the side wall initialization) if (Ri < R + W) { - phaseField->get(x, y, z) += real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (distance - lHalf) / W); + phaseField->get(x, y, z) += + real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (distance - lHalf) / W)); // limit maximum values of phase field if (phaseField->get(x, y, z) > real_c(1)) { phaseField->get(x, y, z) = real_c(1); } @@ -219,16 +220,16 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block bubbleMidPoint[2] = real_c(j + 1) * (D + (dist * W)) - (D + (dist * W)) / real_c(2.0) + fluctuation_pos[i][j]; - real_t Ri = + real_t Ri = real_c( sqrt((real_c(globalCell[0]) - bubbleMidPoint[0]) * (real_c(globalCell[0]) - bubbleMidPoint[0]) + (real_c(globalCell[1]) - bubbleMidPoint[1]) * (real_c(globalCell[1]) - bubbleMidPoint[1]) + - (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2])); + (real_c(globalCell[2]) - bubbleMidPoint[2]) * (real_c(globalCell[2]) - bubbleMidPoint[2]))); if (real_c(globalCell[0]) >= real_c(i) * (D + dist * W) && real_c(globalCell[0]) <= real_c(i + 1) * (D + dist * W) && real_c(globalCell[2]) >= real_c(j) * (D + dist * W) && real_c(globalCell[2]) <= real_c(j + 1) * (D + dist * W)) phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W); + real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (Ri - (R - fluctuation_radius[i][j])) / W)); if (real_c(globalCell[0]) > real_c(nx) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0); if (real_c(globalCell[2]) > real_c(nz) * (D + dist * W)) phaseField->get(x, y, z) = real_c(1.0); @@ -238,7 +239,7 @@ void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, Block if (real_c(globalCell[1]) > gas_top) { phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W); + real_c(0.5) + real_c(0.5) * real_c(tanh(real_c(2.0) * (gas_top + 10 - real_c(globalCell[1])) / W)); } }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } @@ -260,12 +261,12 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(phaseField, { Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); - real_t R = sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) + - (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2)); + real_t R = real_c(sqrt((real_c(globalCell[0]) - X / 2) * (real_c(globalCell[0]) - X / 2) + + (real_c(globalCell[2]) - Z / 2) * (real_c(globalCell[2]) - Z / 2))); if (R > X) R = X; - real_t tmp = perturbation * X * cos((real_c(2.0) * math::pi * R) / X); + real_t tmp = perturbation * X * real_c(cos((real_c(2.0) * math::pi * R) / X)); phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0))); + real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) + tmp) / (W / real_c(2.0)))); }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -278,10 +279,10 @@ void initPhaseField_RTI(const shared_ptr< StructuredBlockStorage >& blocks, Bloc Cell globalCell; blocks->transformBlockLocalToGlobalCell(globalCell, block, Cell(x, y, z)); real_t tmp = perturbation * X * - (cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X) + - cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X)); + (real_c(cos((real_c(2.0) * math::pi * real_c(globalCell[0])) / X)) + + real_c(cos((real_c(2.0) * math::pi * real_c(globalCell[2])) / X))); phaseField->get(x, y, z) = - real_c(0.5) + real_c(0.5) * tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0))); + real_c(0.5) + real_c(0.5) * real_c(tanh(((real_c(globalCell[1]) - halfY) - tmp) / (W / real_c(2.0)))); }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ } } @@ -311,25 +312,25 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl real_t R1; if (real_c(globalCell[1]) <= start_transition) { - R1 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + R1 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz))); } else if (real_c(globalCell[1]) > start_transition && real_c(globalCell[1]) < start_transition + length_transition) { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); - real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); - R1 = sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t shift_tmp = shift * real_c(0.5) * (1 - real_c(cos(tmp))); + R1 = real_c(sqrt((real_c(globalCell[0]) - Mx - shift_tmp) * (real_c(globalCell[0]) - Mx - shift_tmp) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz))); } else { - R1 = sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + R1 = real_c(sqrt((real_c(globalCell[0]) - Mx - shift) * (real_c(globalCell[0]) - Mx - shift) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz))); } - real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t R2 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz))); if (R1 < R_in) addFlag(flagField->get(x, y, z), boundaryFlag); if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag); }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ @@ -357,13 +358,13 @@ void initTubeWithCylinder(const shared_ptr< StructuredBlockStorage >& blocks, Bl real_c(globalCell[1]) < start_transition + length_transition) { real_t tmp = math::pi * (real_c(globalCell[1]) - start_transition) / (length_transition); - real_t shift_tmp = shift * real_c(0.5) * (1 - cos(tmp)); + real_t shift_tmp = shift * real_c(0.5) * (1 - real_c(cos(tmp))); R_tmp = R_in + shift_tmp; } else { R_tmp = R_in + shift; } - real_t R2 = sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + - (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz)); + real_t R2 = real_c(sqrt((real_c(globalCell[0]) - Mx) * (real_c(globalCell[0]) - Mx) + + (real_c(globalCell[2]) - Mz) * (real_c(globalCell[2]) - Mz))); if (R2 < R_tmp) addFlag(flagField->get(x, y, z), boundaryFlag); if (R2 > R_outer) addFlag(flagField->get(x, y, z), boundaryFlag); }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ -- GitLab From 9069a056b22e0504693e422a53c607b4dabc7500 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 10:00:09 +0200 Subject: [PATCH 08/16] Update comment for scenario setting in setup files --- apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py | 2 +- apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py | 2 +- .../PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py | 2 +- .../PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py index d7e2f1959..58ac2fbac 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_RTI_3D.py @@ -61,7 +61,7 @@ class Scenario: # everything else self.dbFile = "RTI.csv" - self.scenario = 2 # 1 rising bubble or droplet, 2 RTI + self.scenario = 2 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up self.config_dict = self.config() @wlb.member_callback diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index 071ecd940..a1a3c251f 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -37,7 +37,7 @@ class Scenario: self.interface_thickness = 5 self.mobility = 0.1 - self.scenario = 3 # 1 rising bubble, 2 RTI + self.scenario = 3 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up self.counter = 0 self.yPositions = [] diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py index 5caf4b343..213dd9b1b 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_rising_bubble.py @@ -50,7 +50,7 @@ class Scenario: # everything else self.dbFile = "risingBubble3D.db" - self.scenario = 1 # 1 rising bubble, 2 RTI + self.scenario = 1 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up self.counter = 0 self.yPositions = [] diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py index 0f1450dbe..473f331b7 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py @@ -55,7 +55,7 @@ class Scenario: self.center_y = self.size[1] / 2 self.center_z = self.size[2] / 2 - self.scenario = 4 # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up + self.scenario = 4 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up self.counter = 0 self.yPositions = [] -- GitLab From 4b77a97edfb0133fd12b7f1775fd03d365d4b497 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 11:15:08 +0200 Subject: [PATCH 09/16] Update parametrization of drop showcase --- .../PhaseFieldAllenCahn/CPU/multiphase.cpp | 6 +- .../CPU/multiphase_drop.py | 73 ++++++++++++------- 2 files changed, 48 insertions(+), 31 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index 48f7e9011..dae39b88c 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -129,9 +129,9 @@ int main(int argc, char** argv) } case 3:{ auto dropParameters = config->getOneBlock("Drop"); - const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("dropMidPoint"); - const real_t dropRadius = dropParameters.getParameter< real_t >("dropRadius"); - const real_t poolDepth = dropParameters.getParameter< real_t >("poolDepth"); + const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("drop_mid_point"); + const real_t dropRadius = dropParameters.getParameter< real_t >("drop_radius"); + const real_t poolDepth = dropParameters.getParameter< real_t >("pool_depth"); const Vector3< real_t > impact_velocity = dropParameters.getParameter< Vector3< real_t > >("impact_velocity"); init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth); diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index a1a3c251f..fc88d02c5 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -1,41 +1,57 @@ import waLBerla as wlb +import math class Scenario: def __init__(self): - # output frequencies - self.dropRadius = 20 + # physical parameters + self.drop_diameter = 20 + self.poolDepth = self.drop_diameter * 0.5 - self.surface_tension = 0.0001 - self.gravitational_acceleration = -2.48e-10 - self.impact_velocity = (0.02574, -0.0474, 0) + self.bond_number = 3.18 + self.weber_number = 2010 + self.ohnesorge_number = 0.0384 - self.density_ratio = 1000 - self.kin_visc_ratio = 100 + self.relaxation_rate_heavy = 1.988 self.density_heavy = 1 - self.density_light = 0.0012 # self.density_heavy / self.density_ratio + self.density_ratio = 1000 + self.dynamic_viscosity_ratio = 100 - self.relaxation_time_heavy = 0.003 # 1 / 1.9889 - self.relaxation_time_light = 0.03 # 1 / self.kin_visc_ratio * 1 / self.density_ratio * (self.relaxation_time_heavy - 0.5) + 0.5 + self.impact_angle_degree = 0 # drop impact angle in degree - self.vtkWriteFrequency = 10000 - self.meshWriteFrequency = 100 + self.interface_thickness = 4 + self.mobility = 0.03 - # simulation parameters - self.timesteps = 15000 - self.cells = (self.dropRadius * 2, self.dropRadius * 2, self.dropRadius * 4) - self.blocks = (10, 10, 5) - self.periodic = (1, 0, 1) - self.size = (self.cells[0] * self.blocks[0], - self.cells[1] * self.blocks[1], - self.cells[2] * self.blocks[2]) - self.poolDepth = self.size[1] * 0.5 + self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy + self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5) + self.dynamic_viscosity_heavy = self.kinematic_viscosity_heavy * self.density_heavy - # physical parameters - self.dropMidPoint = (self.size[0] / 2, self.poolDepth + self.dropRadius, self.size[2] / 2) + self.density_light = self.density_heavy / self.density_ratio + self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio + self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.dynamic_viscosity_ratio + self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5 - self.interface_thickness = 5 - self.mobility = 0.1 + self.surface_tension = (self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter / self.density_heavy + self.gravitational_acceleration = - self.bond_number * self.surface_tension / self.drop_diameter**2 / self.density_heavy + self.impact_velocity_magnitude = (self.weber_number * self.surface_tension / self.drop_diameter / self.density_heavy)**0.5 + + self.impact_velocity = (self.impact_velocity_magnitude * math.sin(self.impact_angle_degree * math.pi / 180), + -self.impact_velocity_magnitude * math.cos(self.impact_angle_degree * math.pi / 180), + 0) + + self.reference_time = self.drop_diameter / self.impact_velocity_magnitude + self.timesteps = 15 * self.reference_time + self.vtkWriteFrequency = self.reference_time + self.meshWriteFrequency = self.reference_time + + # domain parameters + self.size = (self.drop_diameter * 10, + self.drop_diameter * 5, + self.drop_diameter * 10) + self.blocks = (1, 1, 1) + self.periodic = (1, 0, 1) + self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) + self.drop_mid_point = (self.size[0] / 2, self.poolDepth + self.drop_diameter / 2, self.size[2] / 2) self.scenario = 3 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up @@ -47,6 +63,7 @@ class Scenario: return { 'DomainSetup': { 'blocks': self.blocks, + 'domainSize': self.size, 'cellsPerBlock': self.cells, 'periodic': self.periodic, }, @@ -74,9 +91,9 @@ class Scenario: ] }, 'Drop': { - 'dropMidPoint': self.dropMidPoint, - 'dropRadius': self.dropRadius, - 'poolDepth': self.poolDepth, + 'drop_mid_point': self.drop_mid_point, + 'drop_radius': self.drop_diameter / 2, + 'pool_depth': self.poolDepth, 'impact_velocity': self.impact_velocity, }, } -- GitLab From 06774b0aa0aee2939c29e8b0c5a76362f08b5d64 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 12:09:15 +0200 Subject: [PATCH 10/16] Update parametrization of Taylor bubble showcase --- .../PhaseFieldAllenCahn/CPU/multiphase.cpp | 6 +- .../CPU/multiphase_taylor_bubble.py | 79 ++++++++++--------- 2 files changed, 45 insertions(+), 40 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index dae39b88c..8cae979f6 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -143,9 +143,9 @@ int main(int argc, char** argv) } case 4:{ auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble"); - const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("BubbleRadius"); - const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("InitialHeight"); - const real_t Length = TaylorBubbleParameters.getParameter< real_t >("Length"); + const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("bubble_radius"); + const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("initial_height"); + const real_t Length = TaylorBubbleParameters.getParameter< real_t >("length"); init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length, real_c(interface_thickness)); break; diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py index 473f331b7..42c7093b3 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py @@ -13,43 +13,48 @@ from scipy.ndimage.filters import gaussian_filter class Scenario: def __init__(self): - self.diameter = 64 + # physical parameters + self.tube_diameter = 64 self.bond_number = 100 self.morton_number = 0.015 - self.relaxation_time_liquid = 1 / 1.76 - self.kinematic_viscosity_liquid = 1 / 3 * (self.relaxation_time_liquid - 0.5) + self.relaxation_rate_heavy = 1.76 + self.density_heavy = 1 + self.density_ratio = 744 + self.dynamic_viscosity_ratio = 4236 - self.density_liquid = 1.0 - self.density_gas = self.density_liquid / 744 + self.interface_width = 3 + self.mobility = 0.08 - self.dynamic_viscosity_liquid = self.density_liquid * self.kinematic_viscosity_liquid - self.dynamic_viscosity_gas = self.dynamic_viscosity_liquid / 4236 - self.kinematic_viscosity_gas = self.dynamic_viscosity_gas / self.density_gas - self.relaxation_time_gas = 3 * self.kinematic_viscosity_gas + 0.5 + self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy + self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5) + self.dynamic_viscosity_heavy = self.density_heavy * self.kinematic_viscosity_heavy - self.gravitational_acceleration = - (self.bond_number**3 * self.dynamic_viscosity_liquid**4 / self.density_liquid**4 / - self.morton_number / self.diameter**6)**0.5 + self.density_light = self.density_heavy / self.density_ratio + self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio + self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light + self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5 - self.surface_tension = self.density_liquid * abs(self.gravitational_acceleration) * self.diameter**2 / self.bond_number + self.surface_tension = (self.bond_number * self.dynamic_viscosity_heavy**4 / self.morton_number / + self.density_heavy**2 / self.tube_diameter**2)**0.5 - self.interface_width = 3 - self.mobility = 0.08 + self.gravitational_acceleration = - (self.morton_number * self.density_heavy * self.surface_tension**3 / + self.dynamic_viscosity_heavy**4) - self.timesteps = 20 / (abs(self.gravitational_acceleration) / self.diameter)**0.5 + self.reference_time = (self.tube_diameter / abs(self.gravitational_acceleration))**0.5 + self.timesteps = 20 * self.reference_time - # output frequencies - self.vtkWriteFrequency = self.timesteps // 20 - self.dbWriteFrequency = self.timesteps // 20 - self.meshWriteFrequency = self.timesteps // 20 + self.vtkWriteFrequency = self.reference_time + self.dbWriteFrequency = self.reference_time + self.meshWriteFrequency = self.reference_time # simulation parameters - self.size = (self.diameter, self.diameter * 10, self.diameter) + self.size = (self.tube_diameter, self.tube_diameter * 10, self.tube_diameter) self.blocks = (1, 1, 1) self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) self.periodic = (0, 0, 0) - self.inner_radius = self.diameter + self.inner_radius = self.tube_diameter self.center_x = self.size[0] / 2 self.center_y = self.size[1] / 2 @@ -63,22 +68,22 @@ class Scenario: self.eccentricity_or_pipe_ratio = False # if True eccentricity is conducted otherwise pipe ratio self.ratio = 0 - self.start_transition = (self.size[1] // 2) - 2 * self.diameter - self.length_transition = 4 * self.diameter + self.start_transition = (self.size[1] // 2) - 2 * self.tube_diameter + self.length_transition = 4 * self.tube_diameter setup = "eccentricity" if self.eccentricity_or_pipe_ratio else "ratio" - self.csv_file = f"Taylor_bubble_D_{self.diameter}_C_{setup}_{self.ratio}_W_" \ + self.csv_file = f"Taylor_bubble_D_{self.tube_diameter}_C_{setup}_{self.ratio}_W_" \ f"{self.interface_width}_M_{self.mobility}.csv" - d = self.diameter / 2 - dh = self.diameter - d + d = self.tube_diameter / 2 + dh = self.tube_diameter - d - resolution = self.diameter / 128 + resolution = self.tube_diameter / 128 - self.Donut_D = 0.1 * self.diameter / resolution + self.Donut_D = 0.1 * self.tube_diameter / resolution self.Donut_h = dh / 6 - self.DonutTime = 0.5 * (self.diameter + d) / 2 + self.DonutTime = 0.5 * (self.tube_diameter + d) / 2 self.config_dict = self.config() @@ -89,7 +94,7 @@ class Scenario: 'blocks': self.blocks, 'domainSize': self.size, 'cellsPerBlock': self.cells, - 'diameter': self.diameter, + 'diameter': self.tube_diameter, 'periodic': self.periodic, 'inner_radius': self.inner_radius, 'ratio': self.ratio, @@ -107,13 +112,13 @@ class Scenario: 'scenario': self.scenario, }, 'PhysicalParameters': { - 'density_liquid': self.density_liquid, - 'density_gas': self.density_gas, + 'density_liquid': self.density_heavy, + 'density_gas': self.density_light, 'surface_tension': self.surface_tension, 'mobility': self.mobility, 'gravitational_acceleration': self.gravitational_acceleration, - 'relaxation_time_liquid': self.relaxation_time_liquid - 0.5, - 'relaxation_time_gas': self.relaxation_time_gas - 0.5, + 'relaxation_time_liquid': self.relaxation_time_heavy - 0.5, + 'relaxation_time_gas': self.relaxation_time_light - 0.5, 'interface_thickness': self.interface_width }, 'Boundaries': { @@ -127,9 +132,9 @@ class Scenario: ], }, 'TaylorBubble': { - 'BubbleRadius': 0.75 * 0.5 * self.diameter, - 'InitialHeight': 1 * self.diameter, - 'Length': 3 * self.diameter, + 'bubble_radius': 0.75 * 0.5 * self.tube_diameter, + 'initial_height': 1 * self.tube_diameter, + 'length': 3 * self.tube_diameter, } } -- GitLab From 134a1227f45a54d0053b9b9fec34877c20690131 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 16:43:40 +0200 Subject: [PATCH 11/16] Fix error in drop parameterization --- .../PhaseFieldAllenCahn/CPU/multiphase.cpp | 90 +++++++++---------- .../CPU/multiphase_drop.py | 14 ++- 2 files changed, 56 insertions(+), 48 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index 8cae979f6..36c6aaed2 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -71,9 +71,9 @@ int main(int argc, char** argv) // ADD DOMAIN PARAMETERS // ////////////////////////// - auto domainSetup = config->getOneBlock("DomainSetup"); - Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize"); - const bool tube = domainSetup.getParameter< bool >("tube", true); + auto domainSetup = config->getOneBlock("DomainSetup"); + Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize"); + const bool tube = domainSetup.getParameter< bool >("tube", true); //////////////////////////////////////// // ADD GENERAL SIMULATION PARAMETERS // @@ -97,7 +97,7 @@ int main(int argc, char** argv) BlockDataID velocity_field_ID = field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx); BlockDataID density_field_ID = field::addToStorage< PhaseField_T >(blocks, "density", real_c(1.0), field::fzyx); - BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx); + BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx); BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field"); @@ -115,43 +115,42 @@ int main(int argc, char** argv) WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field") switch (scenario) { - case 1: { - auto bubbleParameters = config->getOneBlock("Bubble"); - const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint"); - const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0)); - const bool bubble = bubbleParameters.getParameter< bool >("bubble", true); - initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble); - break; - } - case 2: { - initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); - break; - } - case 3:{ - auto dropParameters = config->getOneBlock("Drop"); - const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("drop_mid_point"); - const real_t dropRadius = dropParameters.getParameter< real_t >("drop_radius"); - const real_t poolDepth = dropParameters.getParameter< real_t >("pool_depth"); - const Vector3< real_t > impact_velocity = - dropParameters.getParameter< Vector3< real_t > >("impact_velocity"); - init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth); - - initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false, - interface_thickness, impact_velocity); - initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth); - break; - } - case 4:{ - auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble"); - const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("bubble_radius"); - const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("initial_height"); - const real_t Length = TaylorBubbleParameters.getParameter< real_t >("length"); - init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length, + case 1: { + auto bubbleParameters = config->getOneBlock("Bubble"); + const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint"); + const real_t bubbleRadius = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0)); + const bool bubble = bubbleParameters.getParameter< bool >("bubble", true); + initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble); + break; + } + case 2: { + initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); + break; + } + case 3: { + auto dropParameters = config->getOneBlock("Drop"); + const Vector3< real_t > dropMidPoint = dropParameters.getParameter< Vector3< real_t > >("drop_mid_point"); + const real_t dropRadius = dropParameters.getParameter< real_t >("drop_radius"); + const real_t poolDepth = dropParameters.getParameter< real_t >("pool_depth"); + const Vector3< real_t > impact_velocity = dropParameters.getParameter< Vector3< real_t > >("impact_velocity"); + init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth); + + initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false, + interface_thickness, impact_velocity); + initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth); + break; + } + case 4: { + auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble"); + const real_t BubbleRadius = TaylorBubbleParameters.getParameter< real_t >("bubble_radius"); + const real_t InitialHeight = TaylorBubbleParameters.getParameter< real_t >("initial_height"); + const real_t Length = TaylorBubbleParameters.getParameter< real_t >("length"); + init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length, real_c(interface_thickness)); - break; - } - default: - WALBERLA_ABORT("Scenario is not defined.") + break; + } + default: + WALBERLA_ABORT("Scenario is not defined.") } WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done") @@ -161,7 +160,8 @@ int main(int argc, char** argv) pystencils::initialize_phase_field_distributions init_h(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID, interface_thickness); - pystencils::initialize_velocity_based_distributions init_g(density_field_ID,hydrodynamic_PDFs_ID, velocity_field_ID); + pystencils::initialize_velocity_based_distributions init_g(density_field_ID, hydrodynamic_PDFs_ID, + velocity_field_ID); pystencils::phase_field_LB_step phase_field_LB_step(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID, mobility, interface_thickness); @@ -316,15 +316,14 @@ int main(int argc, char** argv) [&]() { if (timeloop.getCurrentTimeStep() % uint_t(meshWriteFrequency) == 0) { - auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field_ID, 0.5, 0, true, - targetRank, MPI_COMM_WORLD); + auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field_ID, 0.5, 0, + true, targetRank, MPI_COMM_WORLD); WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank) { std::string path = ""; std::ostringstream out; out << std::internal << std::setfill('0') << std::setw(6) << counter; - geometry::writeMesh( - path + "taylor_bubble_D_" + std::to_string(domainSize[0]) + "_" + out.str() + ".obj", *mesh); + geometry::writeMesh(path + std::to_string(domainSize[0]) + "_" + out.str() + ".obj", *mesh); counter++; } } @@ -332,7 +331,6 @@ int main(int argc, char** argv) "Mesh writer"); } - // remaining time logger timeloop.addFuncAfterTimeStep( timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency), diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index fc88d02c5..cb84869d9 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -28,7 +28,7 @@ class Scenario: self.density_light = self.density_heavy / self.density_ratio self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio - self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.dynamic_viscosity_ratio + self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5 self.surface_tension = (self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter / self.density_heavy @@ -48,7 +48,7 @@ class Scenario: self.size = (self.drop_diameter * 10, self.drop_diameter * 5, self.drop_diameter * 10) - self.blocks = (1, 1, 1) + self.blocks = (2, 2, 1) self.periodic = (1, 0, 1) self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) self.drop_mid_point = (self.size[0] / 2, self.poolDepth + self.drop_diameter / 2, self.size[2] / 2) @@ -58,6 +58,16 @@ class Scenario: self.counter = 0 self.yPositions = [] + print('meshWriteFrequency', self.meshWriteFrequency) + print('density_liquid', self.density_heavy) + print('density_gas', self.density_light) + print('surface_tension', self.surface_tension) + print('mobility', self.mobility) + print('gravitational_acceleration', self.gravitational_acceleration) + print('relaxation_time_liquid', self.relaxation_time_heavy) + print('relaxation_time_gas', self.relaxation_time_light) + print('interface_thickness', self.interface_thickness) + @wlb.member_callback def config(self): return { -- GitLab From 6f4ed1ee08f5a0db647bb9e5d0ad58c9a6f58acf Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 16:56:42 +0200 Subject: [PATCH 12/16] Remove debuggin output in drop parameter file --- .../PhaseFieldAllenCahn/CPU/multiphase_drop.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index cb84869d9..38bfc42b7 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -58,16 +58,6 @@ class Scenario: self.counter = 0 self.yPositions = [] - print('meshWriteFrequency', self.meshWriteFrequency) - print('density_liquid', self.density_heavy) - print('density_gas', self.density_light) - print('surface_tension', self.surface_tension) - print('mobility', self.mobility) - print('gravitational_acceleration', self.gravitational_acceleration) - print('relaxation_time_liquid', self.relaxation_time_heavy) - print('relaxation_time_gas', self.relaxation_time_light) - print('interface_thickness', self.interface_thickness) - @wlb.member_callback def config(self): return { -- GitLab From 6a061ea34b845cdad6e4f08db46b748e8af038d7 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 17:39:44 +0200 Subject: [PATCH 13/16] Change default value for tube to false in multiphase.cpp --- apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp | 4 ++-- apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp index 36c6aaed2..2ed8706f3 100644 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp @@ -73,7 +73,7 @@ int main(int argc, char** argv) auto domainSetup = config->getOneBlock("DomainSetup"); Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize"); - const bool tube = domainSetup.getParameter< bool >("tube", true); + const bool tube = domainSetup.getParameter< bool >("tube", false); //////////////////////////////////////// // ADD GENERAL SIMULATION PARAMETERS // @@ -323,7 +323,7 @@ int main(int argc, char** argv) std::string path = ""; std::ostringstream out; out << std::internal << std::setfill('0') << std::setw(6) << counter; - geometry::writeMesh(path + std::to_string(domainSize[0]) + "_" + out.str() + ".obj", *mesh); + geometry::writeMesh(path + out.str() + ".obj", *mesh); counter++; } } diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index 38bfc42b7..2695a3f48 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -5,7 +5,7 @@ import math class Scenario: def __init__(self): # physical parameters - self.drop_diameter = 20 + self.drop_diameter = 10 self.poolDepth = self.drop_diameter * 0.5 self.bond_number = 3.18 @@ -48,7 +48,7 @@ class Scenario: self.size = (self.drop_diameter * 10, self.drop_diameter * 5, self.drop_diameter * 10) - self.blocks = (2, 2, 1) + self.blocks = (1, 1, 1) self.periodic = (1, 0, 1) self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) self.drop_mid_point = (self.size[0] / 2, self.poolDepth + self.drop_diameter / 2, self.size[2] / 2) -- GitLab From d73bb658240b2f593e1cee9f6a21b52e670dd597 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Thu, 21 Jul 2022 17:40:20 +0200 Subject: [PATCH 14/16] Change default drop diameter to 20 --- apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index 2695a3f48..ebe9d5193 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -5,7 +5,7 @@ import math class Scenario: def __init__(self): # physical parameters - self.drop_diameter = 10 + self.drop_diameter = 20 self.poolDepth = self.drop_diameter * 0.5 self.bond_number = 3.18 -- GitLab From 32888d718c8eb2d06d645e0ed3f665d5fea437f9 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Fri, 22 Jul 2022 10:39:08 +0200 Subject: [PATCH 15/16] Use relaxation time for drop correctly --- .../CPU/multiphase_drop.py | 23 +++++++++++-------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py index ebe9d5193..2845605b1 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py @@ -6,7 +6,7 @@ class Scenario: def __init__(self): # physical parameters self.drop_diameter = 20 - self.poolDepth = self.drop_diameter * 0.5 + self.pool_depth = self.drop_diameter * 0.5 self.bond_number = 3.18 self.weber_number = 2010 @@ -17,7 +17,7 @@ class Scenario: self.density_ratio = 1000 self.dynamic_viscosity_ratio = 100 - self.impact_angle_degree = 0 # drop impact angle in degree + self.impact_angle_degree = 0 # drop impact angle in degree self.interface_thickness = 4 self.mobility = 0.03 @@ -31,9 +31,14 @@ class Scenario: self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5 - self.surface_tension = (self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter / self.density_heavy - self.gravitational_acceleration = - self.bond_number * self.surface_tension / self.drop_diameter**2 / self.density_heavy - self.impact_velocity_magnitude = (self.weber_number * self.surface_tension / self.drop_diameter / self.density_heavy)**0.5 + self.surface_tension = ((self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter + / self.density_heavy) + + self.gravitational_acceleration = (- self.bond_number * self.surface_tension / self.drop_diameter**2 / + self.density_heavy) + + self.impact_velocity_magnitude = (self.weber_number * self.surface_tension / self.drop_diameter / + self.density_heavy)**0.5 self.impact_velocity = (self.impact_velocity_magnitude * math.sin(self.impact_angle_degree * math.pi / 180), -self.impact_velocity_magnitude * math.cos(self.impact_angle_degree * math.pi / 180), @@ -51,7 +56,7 @@ class Scenario: self.blocks = (1, 1, 1) self.periodic = (1, 0, 1) self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2]) - self.drop_mid_point = (self.size[0] / 2, self.poolDepth + self.drop_diameter / 2, self.size[2] / 2) + self.drop_mid_point = (self.size[0] / 2, self.pool_depth + self.drop_diameter / 2, self.size[2] / 2) self.scenario = 3 # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up @@ -80,8 +85,8 @@ class Scenario: 'surface_tension': self.surface_tension, 'mobility': self.mobility, 'gravitational_acceleration': self.gravitational_acceleration, - 'relaxation_time_liquid': self.relaxation_time_heavy, - 'relaxation_time_gas': self.relaxation_time_light, + 'relaxation_time_liquid': self.relaxation_time_heavy - 0.5, + 'relaxation_time_gas': self.relaxation_time_light - 0.5, 'interface_thickness': self.interface_thickness, }, 'Boundaries': { @@ -93,7 +98,7 @@ class Scenario: 'Drop': { 'drop_mid_point': self.drop_mid_point, 'drop_radius': self.drop_diameter / 2, - 'pool_depth': self.poolDepth, + 'pool_depth': self.pool_depth, 'impact_velocity': self.impact_velocity, }, } -- GitLab From 5f640b5df4c7f4bd3da9efa290999b196ee48903 Mon Sep 17 00:00:00 2001 From: Christoph Schwarzmeier Date: Fri, 22 Jul 2022 10:40:04 +0200 Subject: [PATCH 16/16] Remove unused imports in Taylor bubble setup --- .../PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py index 42c7093b3..f07adfe57 100755 --- a/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py +++ b/apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py @@ -1,14 +1,10 @@ -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 class Scenario: -- GitLab