Introduce BoundingBox, skip simulation test

.gitlab-ci.yml

@@ -220,6 +220,7 @@ gcc_9_python_bindings_test:
     - export PYTHONPATH="Python"
     - export VF_WHEEL=$(find dist/*.whl)
     - pip3 install $VF_WHEEL
+    - pip3 install -r Python/requirements.txt
 
   script:
     - python3 -m unittest discover -s Python -v

Python/poiseuille/simulation.py

@@ -5,36 +5,31 @@ from pyfluids.kernel import LBMKernel, KernelType
 from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParameters
 from pyfluids.writer import Writer, OutputFormat
 
-grid_params = GridParameters()
-grid_params.node_distance = 1
-grid_params.number_of_nodes_per_direction = [1, 1, 10]
-grid_params.blocks_per_direction = [1, 1, 1]
-grid_params.periodic_boundary_in_x1 = True
-grid_params.periodic_boundary_in_x2 = True
+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.periodic_boundary_in_x1 = True
+default_grid_params.periodic_boundary_in_x2 = True
 
-physical_params = PhysicalParameters()
-physical_params.lattice_viscosity = 0.005
+default_physical_params = PhysicalParameters()
+default_physical_params.lattice_viscosity = 0.005
 
-runtime_params = RuntimeParameters()
-runtime_params.number_of_threads = 4
-runtime_params.number_of_timesteps = 1000
-runtime_params.timestep_log_interval = 100
+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
 
 
-def run_simulation(physical_params=physical_params, grid_params=grid_params, runtime_params=runtime_params):
+def run_simulation(physical_params=default_physical_params,
+                   grid_params=default_grid_params,
+                   runtime_params=default_runtime_params):
     simulation = Simulation()
 
     kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
     kernel.use_forcing = True
     kernel.forcing_in_x1 = 1e-6
 
-    node_distance = grid_params.node_distance
-    min_x1, min_x2, min_x3 = 0, 0, 0
-    max_x1, max_x2, max_x3 = tuple(map(
-        lambda n: n * node_distance,
-        grid_params.number_of_nodes_per_direction
-    ))
-
     writer = Writer()
     writer.output_path = "./output"
     writer.output_format = OutputFormat.BINARY
@@ -47,24 +42,26 @@ def run_simulation(physical_params=physical_params, grid_params=grid_params, run
 
     no_slip_bc = NoSlipBoundaryCondition()
 
-    block_width = 3 * node_distance
+    block_width = 3 * grid_params.node_distance
     simulation.add_object(
-        GbCuboid3D(min_x1 - block_width,
-                   min_x2 - block_width,
-                   min_x3 - block_width,
-                   max_x1 + block_width,
-                   max_x2 + block_width,
-                   min_x3),
+        GbCuboid3D(
+            grid_params.bounding_box.min_x1 - block_width,
+            grid_params.bounding_box.min_x2 - block_width,
+            grid_params.bounding_box.min_x3 - block_width,
+            grid_params.bounding_box.max_x1 + block_width,
+            grid_params.bounding_box.max_x2 + block_width,
+            grid_params.bounding_box.min_x3),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMin")
 
     simulation.add_object(
-        GbCuboid3D(min_x1 - block_width,
-                   min_x2 - block_width,
-                   max_x3,
-                   max_x1 + block_width,
-                   max_x2 + block_width,
-                   max_x3 + block_width),
+        GbCuboid3D(
+            grid_params.bounding_box.min_x1 - block_width,
+            grid_params.bounding_box.min_x2 - block_width,
+            grid_params.bounding_box.max_x3,
+            grid_params.bounding_box.max_x1 + block_width,
+            grid_params.bounding_box.max_x2 + block_width,
+            grid_params.bounding_box.max_x3 + block_width),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMax")
 

Python/poiseuille/test_poiseuille_l2.py

+import shutil
 import unittest
 
-import shutil
-import math
 import pyvista as pv
-from norms import root_mean_squared_error, mean_squared_error, mean_absolute_error
 from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+from norms import root_mean_squared_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
 
-from matplotlib.pyplot import plot, show, legend
-
 
 class TestPoiseuilleFlow(unittest.TestCase):
 
@@ -18,6 +16,7 @@ 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")
 
         physical_params = PhysicalParameters()
         physical_params.lattice_viscosity = 0.0005
@@ -27,57 +26,78 @@ class TestPoiseuilleFlow(unittest.TestCase):
         runtime_params.timestep_log_interval = 1000
 
         runtime_params.number_of_timesteps = 10000
-        grid_params = create_grid_params_with_nodes_in_column(nodes_in_column=5, height=10)
-        l2_norm_result_100 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, 11)
+        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)
 
         runtime_params.number_of_timesteps *= 2
         physical_params.lattice_viscosity *= 2
-        grid_params = create_grid_params_with_nodes_in_column(nodes_in_column=10, height=10)
-        l2_norm_result_200 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, 11)
+        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)
 
         runtime_params.number_of_timesteps *= 2
         physical_params.lattice_viscosity *= 2
-        grid_params = create_grid_params_with_nodes_in_column(nodes_in_column=20, height=10)
-        l2_norm_result_400 = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, 11)
+        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)
 
-        # nodes = [5, 10, 20]
-        # l2_norms = [l2_norm_result_100, l2_norm_result_200, l2_norm_result_400]
-        # plot(nodes, l2_norms)
-        # show()
-        #
-        # self.assertTrue(l2_norm_result_200 <= l2_norm_result_100)
-        # self.assertTrue(l2_norm_result_400 <= l2_norm_result_200)
 
+        self.assertTrue(l2_norm_result_200 <= l2_norm_result_100)
+        self.assertTrue(l2_norm_result_400 <= l2_norm_result_200)
 
-def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, total_height):
-    output_folder = "./output"
+
+def get_delta_x(number_of_nodes, height):
+    return height / number_of_nodes
+
+
+def run_simulation_with_settings(grid_params, physical_params, runtime_params, output_folder):
     remove_existing_output_directory(output_folder)
     run_simulation(physical_params, grid_params, runtime_params)
+
+
+def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, channel_height):
+    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)
+
+    return root_mean_squared_error(analytical_results, numerical_results)
+
+
+def get_heights(output_folder, runtime_params):
     mesh_of_last_timestep = get_mesh_for_last_timestep(output_folder, runtime_params)
     column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
-    velocities_in_x_direction = mesh_of_last_timestep.get_array("Vx")
-    numerical_results = get_values_from_indices(velocities_in_x_direction, column_indices)
     heights = get_heights_from_indices(mesh_of_last_timestep, column_indices)
-    analytical_results = get_analytical_results(physical_params, total_height, heights)
+    return heights
 
-    print("")
-    # print(f"Numerical results: {numerical_results}")
-    # print(f"Analytical results: {analytical_results}")
-    print(f"Heights: {(min(heights), max(heights))}")
-    plot(heights, numerical_results)
-    plot(heights, analytical_results)
-    legend(["numerical", "analytical"])
-    show()
 
-    return root_mean_squared_error(analytical_results, numerical_results)
+def get_numerical_results(runtime_params, output_folder, heights):
+    mesh_of_last_timestep = get_mesh_for_last_timestep(output_folder, runtime_params)
+    velocities_in_x_direction = mesh_of_last_timestep.get_array("Vx")
+    column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
+    numerical_results = get_values_from_indices(velocities_in_x_direction, column_indices)
 
+    return numerical_results
 
-def get_analytical_results(physical_params, total_height, height_values):
-    settings = get_analytical_poiseuille_settings(total_height, physical_params)
+
+def calculate_analytical_results(physical_params, height_values, channel_height):
+    settings = get_analytical_poiseuille_settings(channel_height, physical_params)
     analytical_results = poiseuille_at_heights(settings, height_values)
     return analytical_results
 
 
+def get_analytical_results(physical_params, heights, channel_height):
+    analytical_results = calculate_analytical_results(physical_params, heights, channel_height)
+    return analytical_results
+
+
 def get_mesh_for_last_timestep(output_folder, runtime_params):
     file_name_of_last_timestep = get_output_file_name(output_folder, runtime_params)
     mesh_of_last_timestep = pv.read(file_name_of_last_timestep)
@@ -94,14 +114,13 @@ def get_analytical_poiseuille_settings(height, physical_params):
     settings.viscosity = physical_params.lattice_viscosity
     settings.density = 1
     settings.force = 1e-6
-
-    print(f"Poiseuille height {settings.height}")
     return settings
 
 
 def get_output_file_name(output_folder, runtime_params):
     timesteps = runtime_params.number_of_timesteps
     file_name = f"{output_folder}/mq/mq{timesteps}/mq0_{timesteps}.bin.vtu"
+
     return file_name
 
 
@@ -109,12 +128,16 @@ def get_heights_from_indices(mesh, indices):
     return [mesh.points[index][2] for index in indices]
 
 
-def create_grid_params_with_nodes_in_column(nodes_in_column, height):
+def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x):
     grid_params = GridParameters()
-    grid_params.node_distance = height / (nodes_in_column - 1)
-    grid_params.number_of_nodes_per_direction = [nodes_in_column, nodes_in_column, nodes_in_column]
-    grid_params.blocks_per_direction = [1, 1, 1]
+    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.periodic_boundary_in_x1 = True
     grid_params.periodic_boundary_in_x2 = True
+    grid_params.periodic_boundary_in_x3 = False
+
+    print(f"GridParameters.node_distance = {grid_params.node_distance}")
+    print(f"GridParameters.number_of_nodes_per_direction = {grid_params.number_of_nodes_per_direction}")
 
     return grid_params

Python/requirements.txt

+appdirs==1.4.4
+attrs==20.3.0
+cycler==0.10.0
+imageio==2.9.0
+iniconfig==1.1.1
+kiwisolver==1.3.1
+meshio==4.3.8
+numpy==1.19.5
+packaging==20.8
+Pillow==8.1.0
+pluggy==0.13.1
+py==1.10.0
+pyparsing==2.4.7
+pytest==6.2.1
+python-dateutil==2.8.1
+pyvista==0.27.4
+scooby==0.5.6
+six==1.15.0
+toml==0.10.2
+vtk==9.0.1

src/cpu/pythonbindings/src/submodules/simulationparameters.cpp

 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <complex>
 #include <simulationconfig/SimulationParameters.h>
 
 namespace py = pybind11;
@@ -22,7 +23,29 @@ void makeParametersModule(py::module_ &parentModule)
             .def_readwrite("blocks_per_direction", &GridParameters::blocksPerDirection)
             .def_readwrite("periodic_boundary_in_x1", &GridParameters::periodicBoundaryInX1)
             .def_readwrite("periodic_boundary_in_x2", &GridParameters::periodicBoundaryInX2)
-            .def_readwrite("periodic_boundary_in_x3", &GridParameters::periodicBoundaryInX3);
+            .def_readwrite("periodic_boundary_in_x3", &GridParameters::periodicBoundaryInX3)
+            .def_property_readonly("bounding_box", &GridParameters::boundingBox);
+
+    py::class_<BoundingBox, std::shared_ptr<BoundingBox>>(parametersModule, "BoundingBox")
+            .def_readonly("min_x1", &BoundingBox::minX1)
+            .def_readonly("min_x2", &BoundingBox::minX2)
+            .def_readonly("min_x3", &BoundingBox::minX3)
+            .def_readonly("max_x1", &BoundingBox::maxX1)
+            .def_readonly("max_x2", &BoundingBox::maxX2)
+            .def_readonly("max_x3", &BoundingBox::maxX3)
+            .def("__repr__", [](BoundingBox &self)
+            {
+                std::ostringstream stream;
+                stream << "<BoundingBox" << std::endl
+                       << "min x1: " << self.minX1 << std::endl
+                       << "min x2: " << self.minX2 << std::endl
+                       << "min x3: " << self.minX3 << std::endl
+                       << "max x1: " << self.maxX1 << std::endl
+                       << "max x2: " << self.maxX2 << std::endl
+                       << "max x3: " << self.maxX3 << std::endl << ">";
+
+                return stream.str();
+            });
 
     py::class_<RuntimeParameters, std::shared_ptr<RuntimeParameters>>(parametersModule, "RuntimeParameters")
             .def(py::init())

src/cpu/simulationconfig/include/simulationconfig/Simulation.h

@@ -63,7 +63,7 @@ public:
     void run();
 
 private:
-    std::shared_ptr<GbObject3D> makeSimulationBoundingBox(const int &nodesInX1, const int &nodesInX2, const int &nodesInX3) const;
+    std::shared_ptr<GbObject3D> makeSimulationBoundingBox() const;
 
     void writeBlocksToFile() const;
 

src/cpu/simulationconfig/include/simulationconfig/SimulationParameters.h

@@ -2,14 +2,36 @@
 #define VIRTUALFLUIDSPYTHONBINDINGS_SIMULATIONPARAMETERS_H
 
 #include <array>
+#include <memory>
 #include <geometry3d/GbPoint3D.h>
 
-struct PhysicalParameters {
+struct PhysicalParameters
+{
     double latticeViscosity{};
     double bulkViscosityFactor{1};
 };
 
-struct GridParameters {
+struct BoundingBox
+{
+    const double minX1;
+    const double minX2;
+    const double minX3;
+    const double maxX1;
+    const double maxX2;
+    const double maxX3;
+
+    BoundingBox(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3) :
+            minX1(minX1),
+            minX2(minX2),
+            minX3(minX3),
+            maxX1(maxX1),
+            maxX2(maxX2),
+            maxX3(maxX3)
+    {}
+};
+
+struct GridParameters
+{
     std::array<int, 3> numberOfNodesPerDirection{1, 1, 1};
     std::array<int, 3> blocksPerDirection{1, 1, 1};
     int referenceDirectionIndex{};
@@ -17,9 +39,20 @@ struct GridParameters {
     bool periodicBoundaryInX1{};
     bool periodicBoundaryInX2{};
     bool periodicBoundaryInX3{};
+
+    std::shared_ptr<BoundingBox> boundingBox()
+    {
+        return std::make_shared<BoundingBox>(
+                0, 0, 0,
+                numberOfNodesPerDirection[0] * nodeDistance,
+                numberOfNodesPerDirection[1] * nodeDistance,
+                numberOfNodesPerDirection[2] * nodeDistance
+        );
+    }
 };
 
-struct RuntimeParameters {
+struct RuntimeParameters
+{
     int numberOfTimeSteps{};
     int timeStepLogInterval{};
     int numberOfThreads{};

src/cpu/simulationconfig/src/Simulation.cpp

@@ -108,7 +108,7 @@ void Simulation::run()
     logSimulationData(nodesInX1, nodesInX2, nodesInX3);
 
     setBlockSize(nodesInX1, nodesInX2, nodesInX3);
-    auto gridCube = makeSimulationBoundingBox(nodesInX1, nodesInX2, nodesInX3);
+    auto gridCube = makeSimulationBoundingBox();
 
     generateBlockGrid(gridCube);
 
@@ -254,18 +254,15 @@ void Simulation::writeBlocksToFile() const
 }
 
 std::shared_ptr<GbObject3D>
-Simulation::makeSimulationBoundingBox(const int &nodesInX1, const int &nodesInX2,
-                                      const int &nodesInX3) const
+Simulation::makeSimulationBoundingBox() const
 {
-    double minX1 = 0, minX2 = 0, minX3 = 0;
-    const double maxX1 = minX1 + gridParameters->nodeDistance * nodesInX1;
-    const double maxX2 = minX2 + gridParameters->nodeDistance * nodesInX2;
-    const double maxX3 = minX3 + gridParameters->nodeDistance * nodesInX3;
+    auto box = gridParameters->boundingBox();
+
     UBLOG(logINFO, "Bounding box dimensions = [("
-            << minX1 << ", " << minX2 << ", " << minX3 << "); ("
-            << maxX1 << ", " << maxX2 << ", " << maxX3 << ")]")
+            << box->minX1 << ", " << box->minX2 << ", " << box->minX3 << "); ("
+            << box->maxX1 << ", " << box->maxX2 << ", " << box->maxX3 << ")]")
 
-    auto gridCube = std::make_shared<GbCuboid3D>(minX1, minX2, minX3, maxX1, maxX2, maxX3);
+    auto gridCube = std::make_shared<GbCuboid3D>(box->minX1, box->minX2, box->minX3, box->maxX1, box->maxX2, box->maxX3);
     GbSystem3D::writeGeoObject(gridCube.get(), writerConfig.outputPath + "/geo/gridCube", writerConfig.getWriter());
     return gridCube;
 }