From 0f7c485bcf974399eb544af4907cfc4c3f4ee9ee Mon Sep 17 00:00:00 2001
From: Sven Marcus <s.marcus@outlook.de>
Date: Mon, 25 Jan 2021 20:41:33 +0100
Subject: [PATCH] Fix l2 norm test

---
 .gitlab-ci.yml                                |   6 +
 Python/norms.py                               |   8 +-
 Python/poiseuille/analytical.py               |  32 ++---
 Python/poiseuille/simulation.py               |  12 +-
 Python/poiseuille/test_poiseuille_l2.py       |  48 +++++---
 .../WriteMacroscopicQuantitiesCoProcessor.cpp | 109 +++++++++++-------
 src/cpu/simulationconfig/src/Simulation.cpp   |   7 +-
 7 files changed, 142 insertions(+), 80 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 9dc5df9ea..d3a44fcdf 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -173,6 +173,11 @@ gcc_9_python:
 
   needs: ["gcc_9"]
 
+  cache:
+    key: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG"
+    paths:
+      - $BUILD_FOLDER
+
   artifacts:
     paths:
       - build/
@@ -221,6 +226,7 @@ gcc_9_python_bindings_test:
   before_script:
     - export PYTHONPATH="Python"
     - export VF_WHEEL=$(find dist/*.whl)
+    - apt-get update && apt-get install -y libglx0
     - pip3 install $VF_WHEEL
     - pip3 install -r Python/requirements.txt
 
diff --git a/Python/norms.py b/Python/norms.py
index cd4666ed8..78ae34459 100644
--- a/Python/norms.py
+++ b/Python/norms.py
@@ -15,7 +15,6 @@ def root_mean_squared_error(real_values, numerical_values):
         raise ValueError("Real and numerical value lists must be same length")
 
     sum_of_squared_distances = get_sum_of_squared_distances(real_values, numerical_values)
-    sum_of_squared_real_values = sum(real_value ** 2 for real_value in real_values)
 
     return math.sqrt(sum_of_squared_distances / num_values)
 
@@ -41,3 +40,10 @@ def mean_squared_error(real_values, numerical_values):
     sum_of_squared_distances = get_sum_of_squared_distances(real_values, numerical_values)
 
     return sum_of_squared_distances / num_values
+
+
+def l2_norm_error(real_values, numerical_values):
+    sum_of_squared_distances = get_sum_of_squared_distances(real_values, numerical_values)
+    sum_of_squared_real_values = sum(real_value ** 2 for real_value in real_values)
+
+    return math.sqrt(sum_of_squared_distances / sum_of_squared_real_values)
\ No newline at end of file
diff --git a/Python/poiseuille/analytical.py b/Python/poiseuille/analytical.py
index fb5ea9ad0..bca1a3ff9 100644
--- a/Python/poiseuille/analytical.py
+++ b/Python/poiseuille/analytical.py
@@ -1,17 +1,16 @@
 from dataclasses import dataclass
 
 
-@dataclass
 class PoiseuilleSettings:
-    density = 1
-    viscosity = 0.005
-    height = 10
-    length = 1
 
-    pressure_in = 0
-    pressure_out = 0
-
-    force = 0
+    def __init__(self):
+        self.density = 1
+        self.viscosity = 0.005
+        self.height = 10
+        self.length = 1
+        self.pressure_in = 0
+        self.pressure_out = 0
+        self.force = 0
 
 
 def poiseuille_at_z(settings: PoiseuilleSettings, z: float):
@@ -27,9 +26,14 @@ def poiseuille_at_heights(settings: PoiseuilleSettings, heights):
 
 
 if __name__ == '__main__':
-    h1 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
-    h2 = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5]
+    # h1 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+    # h2 = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5]
     settings = PoiseuilleSettings()
-    settings.force = 1e-6
-    print(max(poiseuille_at_heights(settings, h1)))
-    print(max(poiseuille_at_heights(settings, h2)))
+    settings.force = 1e-8
+    settings.height = 32
+
+    # print(max(poiseuille_at_heights(settings, h1)))
+    # print(max(poiseuille_at_heights(settings, h2)))
+
+    v = poiseuille_at_z(settings, 16)
+    print(v)
\ No newline at end of file
diff --git a/Python/poiseuille/simulation.py b/Python/poiseuille/simulation.py
index f2abe8a6e..f7f6d468f 100644
--- a/Python/poiseuille/simulation.py
+++ b/Python/poiseuille/simulation.py
@@ -6,9 +6,9 @@ from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParam
 from pyfluids.writer import Writer, OutputFormat
 
 default_grid_params = GridParameters()
-default_grid_params.node_distance = 1
-default_grid_params.number_of_nodes_per_direction = [1, 1, 10]
-default_grid_params.blocks_per_direction = [1, 1, 1]
+default_grid_params.node_distance = 10 / 32
+default_grid_params.number_of_nodes_per_direction = [8, 8, 32]
+default_grid_params.blocks_per_direction = [1, 1, 4]
 default_grid_params.periodic_boundary_in_x1 = True
 default_grid_params.periodic_boundary_in_x2 = True
 
@@ -17,8 +17,8 @@ default_physical_params.lattice_viscosity = 0.005
 
 default_runtime_params = RuntimeParameters()
 default_runtime_params.number_of_threads = 4
-default_runtime_params.number_of_timesteps = 1000
-default_runtime_params.timestep_log_interval = 100
+default_runtime_params.number_of_timesteps = 100000
+default_runtime_params.timestep_log_interval = 10000
 
 
 def run_simulation(physical_params=default_physical_params,
@@ -28,7 +28,7 @@ def run_simulation(physical_params=default_physical_params,
 
     kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
     kernel.use_forcing = True
-    kernel.forcing_in_x1 = 1e-6
+    kernel.forcing_in_x1 = 1e-8
 
     writer = Writer()
     writer.output_path = "./output"
diff --git a/Python/poiseuille/test_poiseuille_l2.py b/Python/poiseuille/test_poiseuille_l2.py
index 0661f0c06..bbb975fd7 100644
--- a/Python/poiseuille/test_poiseuille_l2.py
+++ b/Python/poiseuille/test_poiseuille_l2.py
@@ -1,10 +1,13 @@
 import shutil
 import unittest
 
+import math
 import pyvista as pv
+import matplotlib.pyplot as plt
+
 from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
 
-from norms import root_mean_squared_error
+from norms import root_mean_squared_error, l2_norm_error
 from poiseuille.analytical import poiseuille_at_heights, PoiseuilleSettings
 from poiseuille.simulation import run_simulation
 from vtk_utilities import vertical_column_from_mesh, get_values_from_indices
@@ -16,36 +19,41 @@ class TestPoiseuilleFlow(unittest.TestCase):
         """
         WHEN comparing the simulation results to the analytical solution THEN the L2-Norm should be less than 1e-4
         """
-        self.skipTest("Skipping test! This test is not implemented correctly")
+        # self.skipTest("Skipping test! This test is not implemented correctly")
+
+        channel_height = 10
 
         physical_params = PhysicalParameters()
-        physical_params.lattice_viscosity = 0.0005
+        physical_params.lattice_viscosity = 0.005
 
         runtime_params = RuntimeParameters()
         runtime_params.number_of_threads = 4
         runtime_params.timestep_log_interval = 1000
-
         runtime_params.number_of_timesteps = 10000
-        channel_height = 10
+
         nodes_in_column = 8
         grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
                                                               delta_x=channel_height / nodes_in_column)
-        l2_norm_result_100 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, channel_height)
+        l2_norm_result_100 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params)
 
         runtime_params.number_of_timesteps *= 2
         physical_params.lattice_viscosity *= 2
         nodes_in_column *= 2
         grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
                                                               delta_x=channel_height / nodes_in_column)
-        l2_norm_result_200 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, channel_height)
+        l2_norm_result_200 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params)
 
-        runtime_params.number_of_timesteps *= 2
+        runtime_params.number_of_timesteps *= 12
         physical_params.lattice_viscosity *= 2
         nodes_in_column *= 2
         grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
                                                               delta_x=channel_height / nodes_in_column)
-        l2_norm_result_400 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, channel_height)
+        l2_norm_result_400 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params)
 
+        nodes = [8, 16, 32]
+        l2_results = [l2_norm_result_100, l2_norm_result_200, l2_norm_result_400]
+        plt.plot(nodes, l2_results)
+        plt.show()
 
         self.assertTrue(l2_norm_result_200 <= l2_norm_result_100)
         self.assertTrue(l2_norm_result_400 <= l2_norm_result_200)
@@ -60,15 +68,20 @@ def run_simulation_with_settings(grid_params, physical_params, runtime_params, o
     run_simulation(physical_params, grid_params, runtime_params)
 
 
-def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, channel_height):
+def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params):
     output_folder = "./output"
     run_simulation_with_settings(grid_params, physical_params, runtime_params, output_folder)
     heights = get_heights(output_folder, runtime_params)
 
     numerical_results = get_numerical_results(runtime_params, output_folder, heights)
-    analytical_results = get_analytical_results(physical_params, heights, channel_height)
+    analytical_results = get_analytical_results(physical_params, heights, grid_params.number_of_nodes_per_direction[2])
+
+    plt.plot(heights, numerical_results)
+    plt.plot(heights, analytical_results)
+    plt.legend(["numerical", "analytical"])
+    plt.show()
 
-    return root_mean_squared_error(analytical_results, numerical_results)
+    return l2_norm_error(analytical_results, numerical_results)
 
 
 def get_heights(output_folder, runtime_params):
@@ -89,6 +102,8 @@ def get_numerical_results(runtime_params, output_folder, heights):
 
 def calculate_analytical_results(physical_params, height_values, channel_height):
     settings = get_analytical_poiseuille_settings(channel_height, physical_params)
+    max_height = max(height_values)
+    height_values = [value / max_height * channel_height for value in height_values]
     analytical_results = poiseuille_at_heights(settings, height_values)
     return analytical_results
 
@@ -113,7 +128,10 @@ def get_analytical_poiseuille_settings(height, physical_params):
     settings.height = height
     settings.viscosity = physical_params.lattice_viscosity
     settings.density = 1
-    settings.force = 1e-6
+    settings.force = 1e-8
+
+    # print(settings)
+
     return settings
 
 
@@ -131,8 +149,8 @@ def get_heights_from_indices(mesh, indices):
 def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x):
     grid_params = GridParameters()
     grid_params.node_distance = delta_x
-    grid_params.number_of_nodes_per_direction = [8, 8, nodes_in_column]
-    grid_params.blocks_per_direction = [1, 1, 4]
+    grid_params.number_of_nodes_per_direction = [2, 2, nodes_in_column]
+    grid_params.blocks_per_direction = [1, 1, 6]
     grid_params.periodic_boundary_in_x1 = True
     grid_params.periodic_boundary_in_x2 = True
     grid_params.periodic_boundary_in_x3 = False
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
index c34dccf9f..50dc98d4f 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
@@ -47,26 +47,31 @@
 #include "basics/writer/WbWriterVtkXmlASCII.h"
 
 WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor() = default;
+
 //////////////////////////////////////////////////////////////////////////
 WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
                                                                              const std::string &path,
                                                                              WbWriter *const writer,
                                                                              SPtr<LBMUnitConverter> conv,
                                                                              SPtr<Communicator> comm)
-    : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
+        : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
 {
-    gridRank     = comm->getProcessID();
+    gridRank = comm->getProcessID();
     minInitLevel = this->grid->getCoarsestInitializedLevel();
     maxInitLevel = this->grid->getFinestInitializedLevel();
 
     blockVector.resize(maxInitLevel + 1);
 
-    for (int level = minInitLevel; level <= maxInitLevel; level++) {
+    for (int level = minInitLevel; level <= maxInitLevel; level++)
+    {
         grid->getBlocks(level, gridRank, true, blockVector[level]);
     }
 }
+
 //////////////////////////////////////////////////////////////////////////
-void WriteMacroscopicQuantitiesCoProcessor::init() {}
+void WriteMacroscopicQuantitiesCoProcessor::init()
+{}
+
 //////////////////////////////////////////////////////////////////////////
 void WriteMacroscopicQuantitiesCoProcessor::process(double step)
 {
@@ -75,14 +80,18 @@ void WriteMacroscopicQuantitiesCoProcessor::process(double step)
 
     UBLOG(logDEBUG3, "WriteMacroscopicQuantitiesCoProcessor::update:" << step);
 }
+
 //////////////////////////////////////////////////////////////////////////
 void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
 {
     int istep = static_cast<int>(step);
 
-    for (int level = minInitLevel; level <= maxInitLevel; level++) {
-        for (SPtr<Block3D> block : blockVector[level]) {
-            if (block) {
+    for (int level = minInitLevel; level <= maxInitLevel; level++)
+    {
+        for (SPtr<Block3D> block : blockVector[level])
+        {
+            if (block)
+            {
                 addDataMQ(block);
             }
         }
@@ -93,26 +102,29 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
     subfolder = "mq" + UbSystem::toString(istep);
     pfilePath = path + "/mq/" + subfolder;
     cfilePath = path + "/mq/mq_collection";
-    partPath  = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
+    partPath = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
 
     std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
-    size_t found         = partName.find_last_of("/");
-    std::string piece    = partName.substr(found + 1);
-    piece                = subfolder + "/" + piece;
+    size_t found = partName.find_last_of("/");
+    std::string piece = partName.substr(found + 1);
+    piece = subfolder + "/" + piece;
 
     std::vector<std::string> cellDataNames;
     std::vector<std::string> pieces = comm->gather(piece);
-    if (comm->getProcessID() == comm->getRoot()) {
+    if (comm->getProcessID() == comm->getRoot())
+    {
         std::string pname =
-            WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
+                WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
         found = pname.find_last_of("/");
         piece = pname.substr(found + 1);
 
         std::vector<std::string> filenames;
         filenames.push_back(piece);
-        if (step == CoProcessor::scheduler->getMinBegin()) {
+        if (step == CoProcessor::scheduler->getMinBegin())
+        {
             WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
-        } else {
+        } else
+        {
             WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
         }
         UBLOG(logINFO, "WriteMacroscopicQuantitiesCoProcessor step: " << istep);
@@ -120,6 +132,7 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
 
     clearData();
 }
+
 //////////////////////////////////////////////////////////////////////////
 void WriteMacroscopicQuantitiesCoProcessor::clearData()
 {
@@ -128,10 +141,11 @@ void WriteMacroscopicQuantitiesCoProcessor::clearData()
     datanames.clear();
     data.clear();
 }
+
 //////////////////////////////////////////////////////////////////////////
 void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 {
-    double level = (double)block->getLevel();
+    double level = (double) block->getLevel();
     //   double blockID = (double)block->getGlobalID();
 
     // Diese Daten werden geschrieben:
@@ -146,18 +160,20 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 
     data.resize(datanames.size());
 
-    SPtr<ILBMKernel> kernel                 = block->getKernel();
-    SPtr<BCArray3D> bcArray                 = kernel->getBCProcessor()->getBCArray();
+    SPtr<ILBMKernel> kernel = block->getKernel();
+    SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
     SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
     LBMReal f[D3Q27System::ENDF + 1];
     LBMReal vx1, vx2, vx3, rho;
 
     // knotennummerierung faengt immer bei 0 an!
-    unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+    int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
 
-    if (block->getKernel()->getCompressible()) {
+    if (block->getKernel()->getCompressible())
+    {
         calcMacros = &D3Q27System::calcCompMacroscopicValues;
-    } else {
+    } else
+    {
         calcMacros = &D3Q27System::calcIncompMacroscopicValues;
     }
 
@@ -165,9 +181,9 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
     int minX2 = 0;
     int minX3 = 0;
 
-    int maxX1 = (int)(distributions->getNX1());
-    int maxX2 = (int)(distributions->getNX2());
-    int maxX3 = (int)(distributions->getNX3());
+    int maxX1 = (int) (distributions->getNX1());
+    int maxX2 = (int) (distributions->getNX2());
+    int maxX3 = (int) (distributions->getNX3());
 
     // int minX1 = 1;
     // int minX2 = 1;
@@ -178,21 +194,25 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
     // int maxX3 = (int)(distributions->getNX3());
 
     // nummern vergeben und node vector erstellen + daten sammeln
-    CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
+    CbArray3D<int> nodeNumbers((int) maxX1, (int) maxX2, (int) maxX3, -1);
     maxX1 -= 2;
     maxX2 -= 2;
     maxX3 -= 2;
 
     // D3Q27BoundaryConditionPtr bcPtr;
-    int nr = (int)nodes.size();
+    int nr = (int) nodes.size();
 
-    for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
-        for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
-            for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
-                if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) {
-                    int index                  = 0;
+    for (int ix3 = minX3; ix3 <= maxX3; ix3++)
+    {
+        for (int ix2 = minX2; ix2 <= maxX2; ix2++)
+        {
+            for (int ix1 = minX1; ix1 <= maxX1; ix1++)
+            {
+                if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3))
+                {
+                    int index = 0;
                     nodeNumbers(ix1, ix2, ix3) = nr++;
-                    Vector3D worldCoordinates  = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+                    Vector3D worldCoordinates = grid->getNodeCoordinates(block, ix1, ix2, ix3);
                     nodes.emplace_back(float(worldCoordinates[0]), float(worldCoordinates[1]),
                                        float(worldCoordinates[2]));
 
@@ -202,7 +222,7 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
 
                     if (UbMath::isNaN(rho) || UbMath::isInfinity(rho))
                         UB_THROW(UbException(
-                            UB_EXARGS, "rho is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                UB_EXARGS, "rho is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
                                            UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
                     // rho=999.0;
@@ -213,19 +233,19 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
                     // press=999.0;
                     if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
                         UB_THROW(UbException(
-                            UB_EXARGS, "vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                UB_EXARGS, "vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
                                            UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
                     // vx1=999.0;
                     if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
                         UB_THROW(UbException(
-                            UB_EXARGS, "vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                UB_EXARGS, "vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
                                            UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
                     // vx2=999.0;
                     if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
                         UB_THROW(UbException(
-                            UB_EXARGS, "vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
+                                UB_EXARGS, "vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
                                            block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
                                            UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
 
@@ -249,15 +269,22 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
     maxX2 -= 1;
     maxX3 -= 1;
     // cell vector erstellen
-    for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
-        for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
-            for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
+    for (int ix3 = minX3; ix3 <= maxX3; ix3++)
+    {
+        for (int ix2 = minX2; ix2 <= maxX2; ix2++)
+        {
+            for (int ix1 = minX1; ix1 <= maxX1; ix1++)
+            {
                 if ((SWB = nodeNumbers(ix1, ix2, ix3)) >= 0 && (SEB = nodeNumbers(ix1 + 1, ix2, ix3)) >= 0 &&
                     (NEB = nodeNumbers(ix1 + 1, ix2 + 1, ix3)) >= 0 && (NWB = nodeNumbers(ix1, ix2 + 1, ix3)) >= 0 &&
                     (SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
                     (NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
-                    (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
-                    cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
+                    (NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0)
+                {
+                    cells.push_back(
+                            makeUbTuple((unsigned int) SWB, (unsigned int) SEB, (unsigned int) NEB, (unsigned int) NWB,
+                                        (unsigned int) SWT, (unsigned int) SET, (unsigned int) NET,
+                                        (unsigned int) NWT));
                 }
             }
         }
diff --git a/src/cpu/simulationconfig/src/Simulation.cpp b/src/cpu/simulationconfig/src/Simulation.cpp
index 978e7e2d5..c9a5884e8 100644
--- a/src/cpu/simulationconfig/src/Simulation.cpp
+++ b/src/cpu/simulationconfig/src/Simulation.cpp
@@ -117,7 +117,7 @@ void Simulation::run()
 
     auto metisVisitor = std::make_shared<MetisPartitioningGridVisitor>(communicator,
                                                                        MetisPartitioningGridVisitor::LevelBased,
-                                                                       D3Q27System::B);
+                                                                       D3Q27System::B, MetisPartitioner::RECURSIVE);
 
     InteractorsHelper intHelper(grid, metisVisitor);
     for (auto const &interactor : interactors)
@@ -125,8 +125,7 @@ void Simulation::run()
 
     intHelper.selectBlocks();
 
-    // important: run this after metis & intHelper.selectBlocks()
-    writeBlocksToFile();
+
     int numberOfProcesses = communicator->getNumberOfProcesses();
     SetKernelBlockVisitor kernelVisitor(lbmKernel, physicalParameters->latticeViscosity,
                                         numberOfProcesses);
@@ -147,6 +146,8 @@ void Simulation::run()
     grid->accept(bcVisitor);
 
     writeBoundaryConditions();
+    // important: run this after metis & intHelper.selectBlocks()
+    writeBlocksToFile();
 
     auto visualizationScheduler = std::make_shared<UbScheduler>(simulationParameters->timeStepLogInterval);
     auto mqCoProcessor = makeMacroscopicQuantitiesCoProcessor(converter,
-- 
GitLab