Merge branch 'AdaptPhaseField' into 'master'

Add drop impact and Taylor bubble test case for PFLBM

See merge request !556

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)

apps/showcases/PhaseFieldAllenCahn/CPU/InitializerFunctions.h

@@ -31,12 +31,24 @@
 
 namespace walberla
 {
-void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
-                           Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5);
+void initPhaseField_sphere(const shared_ptr< StructuredBlockStorage >& blocks,
+                           BlockDataID phaseFieldID, BlockDataID velocityFieldID,
+                           real_t R,
+                           Vector3< real_t > bubbleMidPoint, bool bubble = true, real_t W = 5,
+                           const Vector3< real_t >& initialVelocity = Vector3< real_t >(real_c(0)));
+
+void initPhaseField_pool(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t W,
+                         real_t poolDepth);
+
+void init_hydrostatic_pressure(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID densityFieldID,
+                               real_t GravitationalAcceleration, real_t poolDepth);
 
 void init_Taylor_bubble(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t D = 5,
                         real_t H = 2, real_t DT = 20, real_t Donut_x0 = 40);
 
+void init_Taylor_bubble_cylindric(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& phaseFieldID,
+                                  real_t R, real_t H, real_t L, real_t W);
+
 void init_bubble_field(const shared_ptr< StructuredBlockStorage >& blocks, BlockDataID phaseFieldID, real_t R,
                        real_t W = 5);
 

apps/showcases/PhaseFieldAllenCahn/CPU/multiphase.cpp

@@ -30,9 +30,12 @@
 #include "field/vtk/VTKWriter.h"
 
 #include "geometry/InitBoundaryHandling.h"
+#include "geometry/mesh/TriangleMeshIO.h"
 
 #include "lbm/PerformanceEvaluation.h"
 
+#include "postprocessing/FieldToSurfaceMesh.h"
+
 #include "python_coupling/CreateConfig.h"
 #include "python_coupling/PythonCallback.h"
 
@@ -68,8 +71,9 @@ int main(int argc, char** argv)
       // ADD DOMAIN PARAMETERS //
       //////////////////////////
 
-      auto domainSetup                = config->getOneBlock("DomainSetup");
-      const bool tube                 = domainSetup.getParameter< bool >("tube", true);
+      auto domainSetup             = config->getOneBlock("DomainSetup");
+      Vector3< uint_t > domainSize = domainSetup.getParameter< Vector3< uint_t > >("domainSize");
+      const bool tube              = domainSetup.getParameter< bool >("tube", false);
 
       ////////////////////////////////////////
       // ADD GENERAL SIMULATION PARAMETERS //
@@ -92,7 +96,8 @@ int main(int argc, char** argv)
          field::addToStorage< PdfField_hydro_T >(blocks, "LB velocity field", real_c(0.0), field::fzyx);
       BlockDataID velocity_field_ID =
          field::addToStorage< VelocityField_T >(blocks, "velocity", real_c(0.0), field::fzyx);
-      BlockDataID phase_field_ID = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
+      BlockDataID density_field_ID = field::addToStorage< PhaseField_T >(blocks, "density", real_c(1.0), field::fzyx);
+      BlockDataID phase_field_ID   = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
 
       BlockDataID flagFieldID = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
 
@@ -108,15 +113,45 @@ int main(int argc, char** argv)
       const real_t interface_thickness    = physical_parameters.getParameter< real_t >("interface_thickness");
 
       WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field")
-      if (scenario == 1)
+      switch (scenario)
       {
+      case 1: {
          auto bubbleParameters                  = config->getOneBlock("Bubble");
          const Vector3< real_t > bubbleMidPoint = bubbleParameters.getParameter< Vector3< real_t > >("bubbleMidPoint");
          const real_t bubbleRadius              = bubbleParameters.getParameter< real_t >("bubbleRadius", real_c(20.0));
          const bool bubble                      = bubbleParameters.getParameter< bool >("bubble", true);
-         initPhaseField_sphere(blocks, phase_field_ID, bubbleRadius, bubbleMidPoint, bubble);
+         initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, bubbleRadius, bubbleMidPoint, bubble);
+         break;
+      }
+      case 2: {
+         initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube);
+         break;
+      }
+      case 3: {
+         auto dropParameters                     = config->getOneBlock("Drop");
+         const Vector3< real_t > dropMidPoint    = dropParameters.getParameter< Vector3< real_t > >("drop_mid_point");
+         const real_t dropRadius                 = dropParameters.getParameter< real_t >("drop_radius");
+         const real_t poolDepth                  = dropParameters.getParameter< real_t >("pool_depth");
+         const Vector3< real_t > impact_velocity = dropParameters.getParameter< Vector3< real_t > >("impact_velocity");
+         init_hydrostatic_pressure(blocks, density_field_ID, gravitational_acceleration, poolDepth);
+
+         initPhaseField_sphere(blocks, phase_field_ID, velocity_field_ID, dropRadius, dropMidPoint, false,
+                               interface_thickness, impact_velocity);
+         initPhaseField_pool(blocks, phase_field_ID, interface_thickness, poolDepth);
+         break;
+      }
+      case 4: {
+         auto TaylorBubbleParameters = config->getOneBlock("TaylorBubble");
+         const real_t BubbleRadius   = TaylorBubbleParameters.getParameter< real_t >("bubble_radius");
+         const real_t InitialHeight  = TaylorBubbleParameters.getParameter< real_t >("initial_height");
+         const real_t Length         = TaylorBubbleParameters.getParameter< real_t >("length");
+         init_Taylor_bubble_cylindric(blocks, phase_field_ID, BubbleRadius, InitialHeight, Length,
+                                      real_c(interface_thickness));
+         break;
+      }
+      default:
+         WALBERLA_ABORT("Scenario is not defined.")
       }
-      else if (scenario == 2) { initPhaseField_RTI(blocks, phase_field_ID, interface_thickness, tube); }
       WALBERLA_LOG_INFO_ON_ROOT("Initialisation of the phase-field done")
 
       /////////////////
@@ -125,7 +160,8 @@ int main(int argc, char** argv)
 
       pystencils::initialize_phase_field_distributions init_h(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID,
                                                               interface_thickness);
-      pystencils::initialize_velocity_based_distributions init_g(hydrodynamic_PDFs_ID, velocity_field_ID);
+      pystencils::initialize_velocity_based_distributions init_g(density_field_ID, hydrodynamic_PDFs_ID,
+                                                                 velocity_field_ID);
 
       pystencils::phase_field_LB_step phase_field_LB_step(allen_cahn_PDFs_ID, phase_field_ID, velocity_field_ID,
                                                           mobility, interface_thickness);
@@ -180,6 +216,29 @@ int main(int argc, char** argv)
       if (boundariesConfig)
       {
          geometry::initBoundaryHandling< FlagField_T >(*blocks, flagFieldID, boundariesConfig);
+         if (tube)
+         {
+            // initialize cylindrical domain walls
+            const Vector3< real_t > domainCylinderBottomEnd =
+               Vector3< real_t >(real_c(domainSize[0]) * real_c(0.5), real_c(0), real_c(domainSize[2]) * real_c(0.5));
+            const Vector3< real_t > domainCylinderTopEnd = Vector3< real_t >(
+               real_c(domainSize[0]) * real_c(0.5), real_c(domainSize[1]), real_c(domainSize[2]) * real_c(0.5));
+            const real_t radius = real_c(domainSize[0]) * real_c(0.5);
+            const geometry::Cylinder cylinderTube(domainCylinderBottomEnd, domainCylinderTopEnd, radius);
+
+            for (auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt)
+            {
+               FlagField_T* flagField = blockIt->template getData< FlagField_T >(flagFieldID);
+               auto wallFlag          = flagField->getOrRegisterFlag(wallFlagUID);
+               WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, {
+                  Cell globalCell = Cell(x, y, z);
+                  blocks->transformBlockLocalToGlobalCell(globalCell, *blockIt);
+                  const Vector3< real_t > globalPoint = blocks->getCellCenter(globalCell);
+
+                  if (!geometry::contains(cylinderTube, globalPoint)) { addFlag(flagField->get(x, y, z), wallFlag); }
+               }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+            }
+         }
          geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, flagFieldID, fluidFlagUID);
       }
 
@@ -248,6 +307,30 @@ int main(int argc, char** argv)
          },
          "Python callback");
 
+      int meshWriteFrequency = parameters.getParameter< int >("meshWriteFrequency", 0);
+      int counter            = 0;
+      int targetRank         = 0;
+      if (meshWriteFrequency > 0)
+      {
+         timeloop.addFuncAfterTimeStep(
+            [&]() {
+               if (timeloop.getCurrentTimeStep() % uint_t(meshWriteFrequency) == 0)
+               {
+                  auto mesh = postprocessing::realFieldToSurfaceMesh< PhaseField_T >(blocks, phase_field_ID, 0.5, 0,
+                                                                                     true, targetRank, MPI_COMM_WORLD);
+                  WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank)
+                  {
+                     std::string path = "";
+                     std::ostringstream out;
+                     out << std::internal << std::setfill('0') << std::setw(6) << counter;
+                     geometry::writeMesh(path + out.str() + ".obj", *mesh);
+                     counter++;
+                  }
+               }
+            },
+            "Mesh writer");
+      }
+
       // remaining time logger
       timeloop.addFuncAfterTimeStep(
          timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),

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

apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_codegen.py

@@ -28,7 +28,7 @@ with CodeGeneration() as ctx:
     stencil_hydro = LBStencil(Stencil.D3Q27)
     assert (stencil_phase.D == stencil_hydro.D)
 
-    # In the codegneration skript we only need the symbolic parameters
+    # In the code-generation skript, we only need the symbolic parameters
     parameters = AllenCahnParameters(0, 0, 0, 0, 0)
 
     ########################
@@ -37,6 +37,7 @@ with CodeGeneration() as ctx:
 
     # velocity field
     u = fields(f"vel_field({stencil_hydro.D}): {field_type}[{stencil_hydro.D}D]", layout='fzyx')
+    density_field = fields(f"density_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
     # phase-field
     C = fields(f"phase_field: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
     C_tmp = fields(f"phase_field_tmp: {field_type}[{stencil_hydro.D}D]", layout='fzyx')
@@ -87,7 +88,7 @@ with CodeGeneration() as ctx:
 
     # create the kernels for the initialization of the g and h field
     h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
-    g_updates = initializer_kernel_hydro_lb(method_hydro, 1.0, u, g)
+    g_updates = initializer_kernel_hydro_lb(method_hydro, density_field, u, g)
 
     force_h = interface_tracking_force(C, stencil_phase, parameters)
     hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)

apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_drop.py

+import waLBerla as wlb
+import math
+
+
+class Scenario:
+    def __init__(self):
+        # physical parameters
+        self.drop_diameter = 20
+        self.pool_depth = self.drop_diameter * 0.5
+
+        self.bond_number = 3.18
+        self.weber_number = 2010
+        self.ohnesorge_number = 0.0384
+
+        self.relaxation_rate_heavy = 1.988
+        self.density_heavy = 1
+        self.density_ratio = 1000
+        self.dynamic_viscosity_ratio = 100
+
+        self.impact_angle_degree = 0    # drop impact angle in degree
+
+        self.interface_thickness = 4
+        self.mobility = 0.03
+
+        self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy
+        self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5)
+        self.dynamic_viscosity_heavy = self.kinematic_viscosity_heavy * self.density_heavy
+
+        self.density_light = self.density_heavy / self.density_ratio
+        self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio
+        self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light
+        self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5
+
+        self.surface_tension = ((self.dynamic_viscosity_heavy / self.ohnesorge_number)**2 / self.drop_diameter
+                                / self.density_heavy)
+
+        self.gravitational_acceleration = (- self.bond_number * self.surface_tension / self.drop_diameter**2 /
+                                           self.density_heavy)
+
+        self.impact_velocity_magnitude = (self.weber_number * self.surface_tension / self.drop_diameter /
+                                          self.density_heavy)**0.5
+
+        self.impact_velocity = (self.impact_velocity_magnitude * math.sin(self.impact_angle_degree * math.pi / 180),
+                                -self.impact_velocity_magnitude * math.cos(self.impact_angle_degree * math.pi / 180),
+                                0)
+
+        self.reference_time = self.drop_diameter / self.impact_velocity_magnitude
+        self.timesteps = 15 * self.reference_time
+        self.vtkWriteFrequency = self.reference_time
+        self.meshWriteFrequency = self.reference_time
+
+        # domain parameters
+        self.size = (self.drop_diameter * 10,
+                     self.drop_diameter * 5,
+                     self.drop_diameter * 10)
+        self.blocks = (1, 1, 1)
+        self.periodic = (1, 0, 1)
+        self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2])
+        self.drop_mid_point = (self.size[0] / 2, self.pool_depth + self.drop_diameter / 2, self.size[2] / 2)
+
+        self.scenario = 3   # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
+
+        self.counter = 0
+        self.yPositions = []
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'domainSize': self.size,
+                'cellsPerBlock': self.cells,
+                'periodic': self.periodic,
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'meshWriteFrequency': self.meshWriteFrequency,
+                'remainingTimeLoggerFrequency': 10.0,
+                'scenario': self.scenario,
+            },
+            'PhysicalParameters': {
+                'density_liquid': self.density_heavy,
+                'density_gas': self.density_light,
+                'surface_tension': self.surface_tension,
+                'mobility': self.mobility,
+                'gravitational_acceleration': self.gravitational_acceleration,
+                'relaxation_time_liquid': self.relaxation_time_heavy - 0.5,
+                'relaxation_time_gas': self.relaxation_time_light - 0.5,
+                'interface_thickness': self.interface_thickness,
+            },
+            'Boundaries': {
+                'Border': [
+                    {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
+                ]
+            },
+            'Drop': {
+                'drop_mid_point': self.drop_mid_point,
+                'drop_radius': self.drop_diameter / 2,
+                'pool_depth': self.pool_depth,
+                'impact_velocity': self.impact_velocity,
+            },
+        }
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())

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 = []

apps/showcases/PhaseFieldAllenCahn/CPU/multiphase_taylor_bubble.py

+import os
+import waLBerla as wlb
+import numpy as np
+import pandas as pd
+
+from waLBerla.core_extension import makeSlice
+from waLBerla.tools.sqlitedb import sequenceValuesToScalars
+
+
+class Scenario:
+    def __init__(self):
+        # physical parameters
+        self.tube_diameter = 64
+
+        self.bond_number = 100
+        self.morton_number = 0.015
+
+        self.relaxation_rate_heavy = 1.76
+        self.density_heavy = 1
+        self.density_ratio = 744
+        self.dynamic_viscosity_ratio = 4236
+
+        self.interface_width = 3
+        self.mobility = 0.08
+
+        self.relaxation_time_heavy = 1 / self.relaxation_rate_heavy
+        self.kinematic_viscosity_heavy = 1 / 3 * (self.relaxation_time_heavy - 0.5)
+        self.dynamic_viscosity_heavy = self.density_heavy * self.kinematic_viscosity_heavy
+
+        self.density_light = self.density_heavy / self.density_ratio
+        self.dynamic_viscosity_light = self.dynamic_viscosity_heavy / self.dynamic_viscosity_ratio
+        self.kinematic_viscosity_light = self.dynamic_viscosity_light / self.density_light
+        self.relaxation_time_light = 3 * self.kinematic_viscosity_light + 0.5
+
+        self.surface_tension = (self.bond_number * self.dynamic_viscosity_heavy**4 / self.morton_number /
+                                self.density_heavy**2 / self.tube_diameter**2)**0.5
+
+        self.gravitational_acceleration = - (self.morton_number * self.density_heavy * self.surface_tension**3 /
+                                             self.dynamic_viscosity_heavy**4)
+
+        self.reference_time = (self.tube_diameter / abs(self.gravitational_acceleration))**0.5
+        self.timesteps = 20 * self.reference_time
+
+        self.vtkWriteFrequency = self.reference_time
+        self.dbWriteFrequency = self.reference_time
+        self.meshWriteFrequency = self.reference_time
+
+        # simulation parameters
+        self.size = (self.tube_diameter, self.tube_diameter * 10, self.tube_diameter)
+        self.blocks = (1, 1, 1)
+        self.cells = (self.size[0] // self.blocks[0], self.size[1] // self.blocks[1], self.size[2] // self.blocks[2])
+        self.periodic = (0, 0, 0)
+        self.inner_radius = self.tube_diameter
+
+        self.center_x = self.size[0] / 2
+        self.center_y = self.size[1] / 2
+        self.center_z = self.size[2] / 2
+
+        self.scenario = 4   # 1 rising bubble, 2 RTI, 3 drop, 4 taylor bubble set up
+
+        self.counter = 0
+        self.yPositions = []
+
+        self.eccentricity_or_pipe_ratio = False  # if True eccentricity is conducted otherwise pipe ratio
+        self.ratio = 0
+
+        self.start_transition = (self.size[1] // 2) - 2 * self.tube_diameter
+        self.length_transition = 4 * self.tube_diameter
+
+        setup = "eccentricity" if self.eccentricity_or_pipe_ratio else "ratio"
+
+        self.csv_file = f"Taylor_bubble_D_{self.tube_diameter}_C_{setup}_{self.ratio}_W_" \
+                        f"{self.interface_width}_M_{self.mobility}.csv"
+
+        d = self.tube_diameter / 2
+        dh = self.tube_diameter - d
+
+        resolution = self.tube_diameter / 128
+
+        self.Donut_D = 0.1 * self.tube_diameter / resolution
+        self.Donut_h = dh / 6
+        self.DonutTime = 0.5 * (self.tube_diameter + d) / 2
+
+        self.config_dict = self.config()
+
+    @wlb.member_callback
+    def config(self):
+        return {
+            'DomainSetup': {
+                'blocks': self.blocks,
+                'domainSize': self.size,
+                'cellsPerBlock': self.cells,
+                'diameter': self.tube_diameter,
+                'periodic': self.periodic,
+                'inner_radius': self.inner_radius,
+                'ratio': self.ratio,
+                'start_transition': self.start_transition,
+                'length_transition': self.length_transition,
+                'eccentricity_or_pipe_ration': self.eccentricity_or_pipe_ratio,
+                'tube': True
+            },
+            'Parameters': {
+                'timesteps': self.timesteps,
+                'vtkWriteFrequency': self.vtkWriteFrequency,
+                'dbWriteFrequency': self.dbWriteFrequency,
+                'meshWriteFrequency': self.meshWriteFrequency,
+                'remainingTimeLoggerFrequency': 60.0,
+                'scenario': self.scenario,
+            },
+            'PhysicalParameters': {
+                'density_liquid': self.density_heavy,
+                'density_gas': self.density_light,
+                'surface_tension': self.surface_tension,
+                'mobility': self.mobility,
+                'gravitational_acceleration': self.gravitational_acceleration,
+                'relaxation_time_liquid': self.relaxation_time_heavy - 0.5,
+                'relaxation_time_gas': self.relaxation_time_light - 0.5,
+                'interface_thickness': self.interface_width
+            },
+            'Boundaries': {
+                'Border': [
+                    {'direction': 'N', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'S', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'W', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'E', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'T', 'walldistance': -1, 'flag': 'NoSlip'},
+                    {'direction': 'B', 'walldistance': -1, 'flag': 'NoSlip'},
+                ],
+            },
+            'TaylorBubble': {
+                'bubble_radius': 0.75 * 0.5 * self.tube_diameter,
+                'initial_height': 1 * self.tube_diameter,
+                'length': 3 * self.tube_diameter,
+            }
+        }
+
+    @wlb.member_callback
+    def at_end_of_time_step(self, blocks, **kwargs):
+        t = kwargs["timeStep"]
+        if t % self.dbWriteFrequency == 0:
+            wlb_field = wlb.field.gather(blocks, 'phase', makeSlice[:, :, self.size[2] // 2])
+            if wlb_field:
+                data = {'timestep': t}
+                data.update(self.config_dict['Parameters'])
+                data.update(self.config_dict['DomainSetup'])
+                data.update(self.config_dict['PhysicalParameters'])
+                data.update(self.config_dict['TaylorBubble'])
+                data.update(kwargs)
+
+                phase_field = np.asarray(wlb_field).squeeze()
+                location_of_gas = np.where(phase_field < 0.5)
+                center_of_mass = np.mean(location_of_gas, axis=1)
+
+                assert np.isfinite(np.sum(phase_field)), "NaN detected in the phasefield"
+
+                self.yPositions.append(center_of_mass[1])
+                if len(self.yPositions) > 1:
+                    speed = self.yPositions[-1] - self.yPositions[-2]
+                else:
+                    speed = 0
+
+                data['center_of_mass_x'] = center_of_mass[0]
+                data['center_of_mass_y'] = center_of_mass[1]
+                # data['center_of_mass_z'] = center_of_mass[2]
+
+                data['xCells'] = self.size[0]
+                data['yCells'] = self.size[1]
+                data['zCells'] = self.size[2]
+
+                data['rising_velocity'] = speed
+
+                data['StencilHydroLB'] = kwargs["stencil_hydro"]
+                data['StencilPhaseLB'] = kwargs["stencil_phase"]
+                sequenceValuesToScalars(data)
+
+                df = pd.DataFrame.from_records([data])
+                if not os.path.isfile(self.csv_file):
+                    df.to_csv(self.csv_file, index=False)
+                else:
+                    df.to_csv(self.csv_file, index=False, mode='a', header=False)
+
+
+scenarios = wlb.ScenarioManager()
+scenarios.add(Scenario())

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)