Merge branch 'python-simulation-testing' into 'develop'

Adds a convergence analysis test case that runs on Phoenix

See merge request irmb/VirtualFluids_dev!48

.gitlab-ci.yml

@@ -3,6 +3,7 @@ image: irmb/virtualfluids-python-deps:latest
 stages:
   - build
   - build_python
+  - container_upload
   - test
   - benchmark
   - analyze
@@ -141,26 +142,6 @@ msvc_16:
     paths:
       - $CI_PROJECT_DIR/$env:BUILD_FOLDER/
 
-###############################################################################
-build_singularity_image:
-  stage: build
-
-  tags:
-    - priviliged
-    - linux
-
-  rules:
-    - if: $CI_COMMIT_TAG
-
-  artifacts:
-    expire_in: 1 hrs
-    paths:
-      - Containers/VirtualFluidsOpenMPI.sif
-
-  script:
-    - singularity build --fakeroot Containers/VirtualFluidsOpenMPI.sif Containers/VirtualFluidsOpenMPI.def
-    - ls -sh Containers/VirtualFluidsOpenMPI.sif
-
 
 ###############################################################################
 ##                             Build Python                                  ##
@@ -189,6 +170,26 @@ gcc_9_python:
   script:
     - python3 setup.py bdist_wheel
 
+###############################################################################
+##                            Container Upload                               ##
+###############################################################################
+build_singularity_image:
+  stage: container_upload
+
+  needs:
+    - gcc_9_python
+
+  tags:
+    - linux
+    - privileged
+
+  rules:
+    - if: $CI_COMMIT_TAG
+
+  script:
+    - singularity build Containers/VirtualFluidsPython.sif Containers/VirtualFluidsPython.def
+    - singularity push --docker-username "${CI_REGISTRY_USER}" --docker-password "${CI_REGISTRY_PASSWORD}" Containers/VirtualFluidsPython.sif oras://"$CI_REGISTRY_IMAGE"/"$CI_PROJECT_NAME":"$CI_COMMIT_TAG"
+
 ###############################################################################
 ##                                Tests                                      ##
 ###############################################################################
@@ -236,6 +237,47 @@ gcc_9_python_bindings_test:
   script:
     - python3 -m unittest discover -s Python -v
 
+###############################################################################
+gcc_9_python_slurm_test:
+  stage: test
+
+  needs: ["gcc_9_python"]
+
+  rules:
+    - if: $PHOENIX_PRIVATE_KEY
+
+  tags:
+    - linux
+    - privileged
+
+  variables:
+    SSH_KEY: "$PHOENIX_PRIVATE_KEY"
+    HOST: "$PHOENIX_HOSTNAME"
+    USER: "$PHOENIX_USER"
+
+  before_script:
+    - 'command -v ssh-agent >/dev/null || ( apt-get update -y && apt-get install openssh-client -y )'
+    - apt-get install -y rsync
+    - mkdir -p ~/.ssh
+    - chmod 700 ~/.ssh
+    - eval $(ssh-agent -s)
+    - echo "$SSH_KEY" | tr -d '\r' | ssh-add -
+    - echo $SSH_KEY >> ansible/private_key
+    - ssh-keyscan -t rsa $HOST >> ~/.ssh/known_hosts
+    - ssh $USER@$HOST "rm -rf output-*"
+    - ssh $USER@$HOST "rm -f *.out"
+    - pip3 install git+git://github.com/SvenMarcus/ssh-slurm-runner
+
+  script:
+    - singularity build PoiseuilleTestContainer.sif Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def
+    - scp PoiseuilleTestContainer.sif $USER@$HOST:PoiseuilleTestContainer.sif
+    - scp Python/SlurmTests/poiseuille/slurm.job $USER@$HOST:slurm.job
+    - python3 -m ssh_slurm_runner slurm.job --host $HOST --user $USER --keyfile ansible/private_key
+    - ssh $USER@$HOST "rm -rf output-*"
+    - ssh $USER@$HOST "rm -f *.out"
+    - ssh $USER@$HOST "rm PoiseuilleTestContainer.sif"
+    - ssh $USER@$HOST "rm slurm.job"
+
 ###############################################################################
 ##                            Benchmark                                      ##
 ###############################################################################
@@ -549,10 +591,10 @@ vf_wheel_to_jupyterhub:
   needs: ["gcc_9_python", "gcc_9_unit_tests", "gcc_9_python_bindings_test"]
 
   variables:
-    HOST: "finrod.irmb.bau.tu-bs.de"
-    SSH_KEY: "$SSH_PRIVATE_KEY_JUPYTER_HOST_AT_FINROD"
-    REMOTE_USER: "jupyter_host"
-    jupyter_host: "jupyter_host"
+    HOST: "gitlab-runner01.irmb.bau.tu-bs.de"
+    SSH_KEY: "$SSH_PRIVATE_KEY"
+    REMOTE_USER: "runner"
+    jupyter_host: "runner"
 
   script:
     - ansible-playbook -i ansible/hosts.cfg -u $REMOTE_USER ansible/playbook_jupyter_update.yml

Containers/VirtualFluidsPython.def

+BootStrap: docker
+From: ubuntu:20.04
+
+%files
+    Python Python
+    dist dist
+
+
+%post
+    export DEBIAN_FRONTEND=noninteractive
+    apt-get update &&          \
+    apt-get install -y         \
+    build-essential            \
+    cmake=3.16.3-1ubuntu1      \
+    python3                    \
+    python3-dev                \
+    python3-pip                \
+    mpich                      \
+    libomp-dev
+
+    pip3 install setuptools wheel $(find dist/*.whl)
+
+%environment
+    export PYTHONPATH=/Python
+
+%runscript
+    python3 /Python/liddrivencavity/simulation.py
+
+%appenv poiseuille
+    export PYTHONPATH=Python
+
+%apprun poisueille
+    python3 /Python/poiseuille/poiseuille_hpc.py

Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def

+BootStrap: docker
+From: ubuntu:20.04
+
+%files
+    3rdParty 3rdParty
+    apps apps
+    CMake CMake
+    Python Python
+    src src
+    CMakeLists.txt CMakeLists.txt
+    cpu.cmake cpu.cmake
+    gpu.cmake gpu.cmake
+    setup.py setup.py
+    pyproject.toml pyproject.toml
+
+
+%post
+    export DEBIAN_FRONTEND=noninteractive
+    apt-get update &&          \
+    apt-get install -y         \
+    build-essential            \
+    cmake=3.16.3-1ubuntu1      \
+    python3                    \
+    python3-dev                \
+    python3-pip                \
+    mpich                      \
+    libomp-dev                 \
+    libgl1
+
+    pip3 install setuptools wheel numpy scipy pyvista
+
+    export PYTHONPATH=Python
+    python3 /setup.py install
+
+%environment
+    export PYTHONPATH=/Python
+
+%apprun testcase
+    python3 /Python/SlurmTests/poiseuille/simulation_runner.py
+
+%apprun evaluate
+    python3 /Python/SlurmTests/poiseuille/evaluator.py
\ No newline at end of file

Python/SlurmTests/poiseuille/evaluator.py

+import numpy as np
+import scipy.stats as stats
+import errors
+from SlurmTests.poiseuille.result_collector import collect_results
+from SlurmTests.poiseuille.settings import Scaling
+
+analytical_results, numerical_results = collect_results()
+normalized_l2_errors = [errors.normalized_l2_error(analytical, numerical)
+                        for analytical, numerical in zip(analytical_results, numerical_results)]
+
+nodes_in_x3_per_run = []
+for simulation_run in range(0, 3):
+    grid_params, _, _, _ = Scaling.configuration_for_scale_level(simulation_run)
+    nodes_in_x3_per_run.append(grid_params.number_of_nodes_per_direction[2])
+
+nodes_as_log = [np.log10(node) for node in nodes_in_x3_per_run]
+l2_norms_as_log = [np.log10(l2) for l2 in normalized_l2_errors]
+res = stats.linregress(nodes_as_log, l2_norms_as_log)
+
+assert res.slope <= -2, f"Expected slope of l2 error to be <= -2, but was {res.slope}"

Python/SlurmTests/poiseuille/result_collector.py

+from typing import Collection, List
+
+import pyvista as pv
+from poiseuille.analytical import PoiseuilleSettings, poiseuille_at_heights
+from vtk_utilities import vertical_column_from_mesh, get_values_from_indices
+from SlurmTests.poiseuille.settings import Scaling
+
+
+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
+
+
+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)
+    return mesh_of_last_timestep
+
+
+def get_heights_from_indices(mesh, indices):
+    return [mesh.points[index][2] for index in indices]
+
+
+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)
+    heights = get_heights_from_indices(mesh_of_last_timestep, column_indices)
+    return heights
+
+
+def get_numerical_results(runtime_params, output_folder):
+    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(grid_params, physical_params, kernel, height_values):
+    channel_height = grid_params.number_of_nodes_per_direction[2]
+    settings = get_analytical_poiseuille_settings(channel_height, physical_params, kernel)
+    max_grid_height = channel_height * grid_params.node_distance
+    adjusted_height_values = [value / max_grid_height * channel_height for value in height_values]
+    analytical_results = poiseuille_at_heights(settings, adjusted_height_values)
+    return analytical_results
+
+
+def get_analytical_poiseuille_settings(height, physical_params, kernel):
+    settings = PoiseuilleSettings()
+    settings.height = height
+    settings.viscosity = physical_params.lattice_viscosity
+    settings.density = 1
+    settings.force = kernel.forcing_in_x1
+
+    return settings
+
+
+def collect_results() -> (List[List[float]], List[List[float]]):
+    analytical_results = []
+    numerical_results = []
+
+    for simulation_run in range(0, 3):
+        output_folder = f"output-{simulation_run}"
+        grid_params, physical_params, runtime_params, kernel = Scaling.configuration_for_scale_level(simulation_run)
+        heights = get_heights(output_folder, runtime_params)
+        analytical_results.append(
+            get_analytical_results(grid_params, physical_params, kernel, heights))
+        numerical_results.append(get_numerical_results(runtime_params, output_folder))
+
+    return analytical_results, numerical_results

Python/SlurmTests/poiseuille/settings.py

+import os
+from acousticscaling import OneDirectionalAcousticScaling
+from pyfluids.kernel import LBMKernel, KernelType
+from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParameters
+
+
+grid_params = GridParameters()
+grid_params.node_distance = 1
+grid_params.number_of_nodes_per_direction = [1, 1, 16]
+grid_params.blocks_per_direction = [1, 1, 4]
+grid_params.periodic_boundary_in_x1 = True
+grid_params.periodic_boundary_in_x2 = True
+
+physical_params = PhysicalParameters()
+physical_params.lattice_viscosity = 1e-4
+
+runtime_params = RuntimeParameters()
+runtime_params.number_of_threads = int(os.environ["PYFLUIDS_NUM_THREADS"])
+runtime_params.number_of_timesteps = 4_000_000
+runtime_params.timestep_log_interval = 1_000_000
+
+kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+kernel.use_forcing = True
+kernel.forcing_in_x1 = 5e-10
+
+Scaling = OneDirectionalAcousticScaling(grid_params, physical_params, runtime_params, kernel)

Python/SlurmTests/poiseuille/simulation_runner.py

+import os
+
+from SlurmTests.poiseuille.settings import Scaling
+from poiseuille.simulation import run_simulation
+from pyfluids.writer import Writer, OutputFormat
+
+
+scale_level = int(os.environ["PYFLUIDS_SCALE_LEVEL"])
+grid_params, physical_params, runtime_params, kernel = Scaling.configuration_for_scale_level(scale_level)
+
+writer = Writer()
+writer.output_format = OutputFormat.BINARY
+writer.output_path = "./output-" + str(scale_level)
+
+run_simulation(grid_params=grid_params,
+               physical_params=physical_params,
+               runtime_params=runtime_params,
+               kernel=kernel,
+               writer=writer)

Python/SlurmTests/poiseuille/slurm.job

+#!/bin/bash
+#SBATCH -J PyFluidsTest
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=1
+#SBATCH --cpus-per-task=20
+
+#SBATCH --mem-per-cpu=3000
+#SBATCH --time=07:00:00
+#SBATCH --partition=standard
+
+source $HOME/.bashrc
+
+echo "PyFluids Poiseuille Test Case"
+echo "Number of tasks: ${SLURM_NTASKS}"
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=0
+export SINGULARITYENV_PYFLUIDS_NUM_THREADS=4
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=1
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=2
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+srun singularity run --app evaluate PoiseuilleTestContainer.sif

Python/acousticscaling.py

+from pyfluids.kernel import LBMKernel
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+
+class OneDirectionalAcousticScaling:
+
+    def __init__(self, grid_parameters: GridParameters,
+                 physical_parameters: PhysicalParameters,
+                 runtime_parameters: RuntimeParameters,
+                 kernel: LBMKernel):
+        self._grid_params = grid_parameters
+        self._physical_params = physical_parameters
+        self._runtime_params = runtime_parameters
+        self._kernel = kernel
+
+    def configuration_for_scale_level(self, level: int = 1) -> (GridParameters,
+                                                                PhysicalParameters,
+                                                                RuntimeParameters,
+                                                                LBMKernel):
+        if level < 0:
+            raise ValueError("level must be >= 0")
+
+        grid_params = self.clone_grid_params_for_level(level)
+        physical_params = self.clone_physical_parameters(level)
+        runtime_params = self.clone_runtime_params_for_level(level)
+        kernel = self.clone_kernel_for_level(level)
+
+        return grid_params, physical_params, runtime_params, kernel
+
+    def clone_grid_params_for_level(self, level) -> GridParameters:
+        grid_params = GridParameters()
+        grid_params.reference_direction_index = self._grid_params.reference_direction_index
+        grid_params.periodic_boundary_in_x1 = self._grid_params.periodic_boundary_in_x1
+        grid_params.periodic_boundary_in_x2 = self._grid_params.periodic_boundary_in_x2
+        grid_params.periodic_boundary_in_x3 = self._grid_params.periodic_boundary_in_x3
+
+        grid_params.number_of_nodes_per_direction = list(self._grid_params.number_of_nodes_per_direction)
+        grid_params.blocks_per_direction = list(self._grid_params.blocks_per_direction)
+        grid_params.node_distance = self._grid_params.node_distance
+
+        if level > 0:
+            grid_params.node_distance /= (level * 2)
+            grid_params.number_of_nodes_per_direction = [grid_params.number_of_nodes_per_direction[0],
+                                                         grid_params.number_of_nodes_per_direction[1],
+                                                         grid_params.number_of_nodes_per_direction[2] * (level * 2)]
+
+            grid_params.blocks_per_direction = [grid_params.blocks_per_direction[0],
+                                                grid_params.blocks_per_direction[1],
+                                                grid_params.blocks_per_direction[2] * (level * 2)]
+
+        return grid_params
+
+    def clone_physical_parameters(self, level):
+        physical_params = PhysicalParameters()
+        physical_params.lattice_viscosity = self._physical_params.lattice_viscosity
+
+        if level > 0:
+            physical_params.lattice_viscosity *= (level * 2)
+
+        return physical_params
+
+    def clone_runtime_params_for_level(self, level):
+        runtime_params = RuntimeParameters()
+        runtime_params.number_of_timesteps = self._runtime_params.number_of_timesteps
+        runtime_params.number_of_threads = self._runtime_params.number_of_threads
+        runtime_params.timestep_log_interval = self._runtime_params.timestep_log_interval
+
+        if level > 0:
+            runtime_params.number_of_timesteps *= (level * 2)
+
+        return runtime_params
+
+    def clone_kernel_for_level(self, level):
+        kernel = LBMKernel(self._kernel.type)
+        kernel.use_forcing = self._kernel.use_forcing
+        kernel.forcing_in_x1 = self._kernel.forcing_in_x1
+        kernel.forcing_in_x2 = self._kernel.forcing_in_x2
+        kernel.forcing_in_x3 = self._kernel.forcing_in_x3
+
+        if level > 0:
+            kernel.forcing_in_x1 /= (level * 2)
+            kernel.forcing_in_x2 /= (level * 2)
+            kernel.forcing_in_x3 /= (level * 2)
+
+        return kernel

Python/norms.py → Python/errors.py

@@ -42,8 +42,8 @@ def mean_squared_error(real_values, numerical_values):
     return sum_of_squared_distances / num_values
 
 
-def l2_norm_error(real_values, numerical_values):
+def normalized_l2_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
+    return math.sqrt(sum_of_squared_distances / sum_of_squared_real_values)

Python/poiseuille/analytical.py

 from dataclasses import dataclass
 
 
+@dataclass
 class PoiseuilleSettings:
-
-    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
+    density = 1
+    viscosity = 0.005
+    height = 10
+    length = 1
+    pressure_in = 0
+    pressure_out = 0
+    force = 0
 
 
 def poiseuille_at_z(settings: PoiseuilleSettings, z: float):
     pressure_grad = ((settings.pressure_out - settings.pressure_in) / settings.length)
 
-    return (1 / settings.viscosity
+    return ((1 / settings.viscosity)
             * (- pressure_grad + settings.density * settings.force)
-            * z / 2 * (settings.height - z))
+            * (z / 2)
+            * (settings.height - z))
 
 
 def poiseuille_at_heights(settings: PoiseuilleSettings, heights):
     return [poiseuille_at_z(settings, z) for z in 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]
-    settings = PoiseuilleSettings()
-    settings.force = 1e-8
-    settings.height = 32
+def reynolds_number(settings: PoiseuilleSettings):
+    max_v = poiseuille_at_z(settings, settings.height / 2)
+    return max_v * settings.height / settings.viscosity
 
-    # 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
+if __name__ == '__main__':
+    sim_settings = PoiseuilleSettings()
+
+    sim_settings.force = 2e-7
+    sim_settings.viscosity = 1e-3
+    sim_settings.height = 16
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))
+
+    sim_settings.viscosity *= 2
+    sim_settings.height *= 2
+    sim_settings.force /= 2
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))
+
+    sim_settings.viscosity *= 2
+    sim_settings.height *= 2
+    sim_settings.force /= 2
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))

Python/poiseuille/simulation.py

@@ -17,23 +17,30 @@ default_physical_params.lattice_viscosity = 0.005
 
 default_runtime_params = RuntimeParameters()
 default_runtime_params.number_of_threads = 4
-default_runtime_params.number_of_timesteps = 100000
-default_runtime_params.timestep_log_interval = 10000
+default_runtime_params.number_of_timesteps = 10000
+default_runtime_params.timestep_log_interval = 1000
+
+default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel.use_forcing = True
+default_kernel.forcing_in_x1 = 1e-8
+
+default_writer = Writer()
+default_writer.output_path = "./output"
+default_writer.output_format = OutputFormat.BINARY
+
+
+default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel.use_forcing = True
+default_kernel.forcing_in_x1 = 1e-8
 
 
 def run_simulation(physical_params=default_physical_params,
                    grid_params=default_grid_params,
-                   runtime_params=default_runtime_params):
+                   runtime_params=default_runtime_params,
+                   kernel=default_kernel,
+                   writer=default_writer):
     simulation = Simulation()
 
-    kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
-    kernel.use_forcing = True
-    kernel.forcing_in_x1 = 1e-8
-
-    writer = Writer()
-    writer.output_path = "./output"
-    writer.output_format = OutputFormat.BINARY
-
     simulation.set_kernel_config(kernel)
     simulation.set_physical_parameters(physical_params)
     simulation.set_grid_parameters(grid_params)
@@ -42,26 +49,26 @@ def run_simulation(physical_params=default_physical_params,
 
     no_slip_bc = NoSlipBoundaryCondition()
 
-    block_width = 3 * grid_params.node_distance
+    block_thickness = 3 * grid_params.node_distance
     simulation.add_object(
         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_x1 - block_thickness,
+            grid_params.bounding_box.min_x2 - block_thickness,
+            grid_params.bounding_box.min_x3 - block_thickness,
+            grid_params.bounding_box.max_x1 + block_thickness,
+            grid_params.bounding_box.max_x2 + block_thickness,
             grid_params.bounding_box.min_x3),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMin")
 
     simulation.add_object(
         GbCuboid3D(
-            grid_params.bounding_box.min_x1 - block_width,
-            grid_params.bounding_box.min_x2 - block_width,
+            grid_params.bounding_box.min_x1 - block_thickness,
+            grid_params.bounding_box.min_x2 - block_thickness,
             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),
+            grid_params.bounding_box.max_x1 + block_thickness,
+            grid_params.bounding_box.max_x2 + block_thickness,
+            grid_params.bounding_box.max_x3 + block_thickness),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMax")
 

Python/poiseuille/test_poiseuille_l2.py

+import os
 import shutil
 import unittest
 
 import matplotlib.pyplot as plt
+import numpy as np
 import pyvista as pv
+from pyfluids.kernel import LBMKernel, KernelType
 from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+from scipy import stats
 
-from norms import l2_norm_error
+from errors import normalized_l2_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
 
 
 class TestPoiseuilleFlow(unittest.TestCase):
+    node_distances = [1, 0.5, 0.25]
+    number_of_nodes = [16, 32, 64]
+    number_of_timesteps = [2_500_000, 5_000_000, 10_000_000]
+    forcings = [1e-9, 5e-10, 2.5e-10]
+    viscosities = [1e-3, 2e-3, 4e-3]
+
+    def zipped_settings(self):
+        return zip(self.node_distances,
+                   self.number_of_nodes,
+                   self.number_of_timesteps,
+                   self.forcings,
+                   self.viscosities)
 
     def test_poiseuille_flow(self):
-        self.skipTest("Skipping test! This test has not been implemented correctly yet")
+        self.skipTest("This test is not implemented correctly yet")
         plt.ion()
 
-        channel_height = 10
-        number_of_nodes = [8, 16, 32]
-        number_of_timesteps = [10_000, 20_000, 40_000]
-        viscosities = [5e-3, 1e-2, 2e-2]
-        l2_norm_results = []
-
         physical_params = PhysicalParameters()
 
         runtime_params = RuntimeParameters()
-        runtime_params.number_of_threads = 4
-        runtime_params.timestep_log_interval = 1000
-
-        for test_number, nodes_in_column in enumerate(number_of_nodes):
-            runtime_params.number_of_timesteps = number_of_timesteps[test_number]
-            physical_params.lattice_viscosity = viscosities[test_number]
-            delta_x = channel_height / nodes_in_column
-            grid_params = create_grid_params_with_nodes_in_column(nodes_in_column, delta_x)
-            l2_norm_result = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params)
-            l2_norm_results.append(l2_norm_result)
-
-        plt.plot(number_of_nodes, l2_norm_results)
+        runtime_params.number_of_threads = os.cpu_count()
+        runtime_params.timestep_log_interval = 10000
+
+        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel.use_forcing = True
+
+        normalized_l2_errors = []
+        for delta_x, nodes, timesteps, forcing, viscosity in self.zipped_settings():
+            physical_params.lattice_viscosity = viscosity
+            runtime_params.number_of_timesteps = timesteps
+            kernel.forcing_in_x1 = forcing
+
+            grid_params = create_grid_params_with_nodes_in_column(nodes, delta_x)
+            l2_error = get_l2_error_for_simulation(grid_params, physical_params, runtime_params, kernel)
+            normalized_l2_errors.append(l2_error)
+
+        nodes_as_log = [np.log10(node) for node in self.number_of_nodes]
+        l2_norms_as_log = [np.log10(l2) for l2 in normalized_l2_errors]
+        res = stats.linregress(nodes_as_log, l2_norms_as_log)
+
+        plt.xscale("log")
+        plt.yscale("log")
+        plt.plot(self.number_of_nodes, [np.power(10, res.intercept + res.slope * node) for node in nodes_as_log], 'r-')
+        plt.plot(self.number_of_nodes, normalized_l2_errors, "x:")
         plt.show()
 
-        self.assertTrue(l2_norm_results[1] <= l2_norm_results[0])
-        self.assertTrue(l2_norm_results[2] <= l2_norm_results[1])
-
-
-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)
+        print(normalized_l2_errors)
+        self.assertAlmostEqual(res.slope, -2, places=2)
 
 
-def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params):
+def get_l2_error_for_simulation(grid_params, physical_params, runtime_params, kernel):
     output_folder = "./output"
-    run_simulation_with_settings(grid_params, physical_params, runtime_params, output_folder)
+    run_simulation_with_settings(grid_params, physical_params, runtime_params, kernel, 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, grid_params.number_of_nodes_per_direction[2])
+    numerical_results = get_numerical_results(runtime_params, output_folder)
+    analytical_results = get_analytical_results(grid_params, physical_params, kernel, heights)
 
     plt.plot(heights, numerical_results)
     plt.plot(heights, analytical_results)
     plt.legend(["numerical", "analytical"])
     plt.show()
 
-    return l2_norm_error(analytical_results, numerical_results)
+    return normalized_l2_error(analytical_results, numerical_results)
+
+
+def run_simulation_with_settings(grid_params, physical_params, runtime_params, kernel, output_folder):
+    shutil.rmtree(output_folder, ignore_errors=True)
+    run_simulation(physical_params, grid_params, runtime_params, kernel)
 
 
 def get_heights(output_folder, runtime_params):
@@ -72,7 +94,7 @@ def get_heights(output_folder, runtime_params):
     return heights
 
 
-def get_numerical_results(runtime_params, output_folder, heights):
+def get_numerical_results(runtime_params, output_folder):
     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)
@@ -81,16 +103,12 @@ def get_numerical_results(runtime_params, output_folder, heights):
     return numerical_results
 
 
-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
-
-
-def get_analytical_results(physical_params, heights, channel_height):
-    analytical_results = calculate_analytical_results(physical_params, heights, channel_height)
+def get_analytical_results(grid_params, physical_params, kernel, height_values):
+    channel_height = grid_params.number_of_nodes_per_direction[2]
+    settings = get_analytical_poiseuille_settings(channel_height, physical_params, kernel)
+    max_grid_height = channel_height * grid_params.node_distance
+    adjusted_height_values = [value / max_grid_height * channel_height for value in height_values]
+    analytical_results = poiseuille_at_heights(settings, adjusted_height_values)
     return analytical_results
 
 
@@ -100,18 +118,12 @@ def get_mesh_for_last_timestep(output_folder, runtime_params):
     return mesh_of_last_timestep
 
 
-def remove_existing_output_directory(output_dir):
-    shutil.rmtree(output_dir, ignore_errors=True)
-
-
-def get_analytical_poiseuille_settings(height, physical_params):
+def get_analytical_poiseuille_settings(height, physical_params, kernel):
     settings = PoiseuilleSettings()
     settings.height = height
     settings.viscosity = physical_params.lattice_viscosity
     settings.density = 1
-    settings.force = 1e-8
-
-    # print(settings)
+    settings.force = kernel.forcing_in_x1
 
     return settings
 
@@ -130,13 +142,10 @@ 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 = [2, 2, nodes_in_column]
-    grid_params.blocks_per_direction = [1, 1, 6]
+    grid_params.number_of_nodes_per_direction = [1, 1, nodes_in_column]
+    grid_params.blocks_per_direction = [1, 1, 8]
     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

@@ -14,8 +14,10 @@ py==1.10.0
 pyparsing==2.4.7
 pytest==6.2.1
 python-dateutil==2.8.1
-pyvista==0.27.4
+pyvista==0.28.1
+scipy==1.6.1
 scooby==0.5.6
 six==1.15.0
 toml==0.10.2
+transforms3d==0.3.1
 vtk==9.0.1

Python/tests/test_acousticscaling.py

+import unittest
+from typing import List
+
+from pyfluids.kernel import LBMKernel, KernelType
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+from acousticscaling import OneDirectionalAcousticScaling
+
+
+class OneDirectionalAcousticScalingTest(unittest.TestCase):
+
+    def setUp(self) -> None:
+        self.grid_params = self.make_grid_params()
+        self.physical_params = self.make_physical_params()
+        self.runtime_params = self.make_runtime_params()
+        self.kernel = self.make_kernel()
+
+        self.sut = OneDirectionalAcousticScaling(self.grid_params, self.physical_params, self.runtime_params,
+                                                 self.kernel)
+
+    def test_given_sim_parameters__when_scaling_level_zero__should_return_equal_sim_parameters(self):
+        factor = 1
+        actual_params = self.sut.configuration_for_scale_level(0)
+        actual_grid_params = actual_params[0]
+        actual_physical_params = actual_params[factor]
+        actual_runtime_params = actual_params[2]
+        actual_kernel = actual_params[3]
+
+        self.assert_parameters_scaled_by_factor(actual_grid_params, actual_kernel,
+                                                actual_physical_params, actual_runtime_params, factor)
+
+    def test_given_sim_parameters__when_scaling_level_one__should_return_sim_parameters_scaled_by_two(self):
+        actual_params = self.sut.configuration_for_scale_level(1)
+        actual_grid_params = actual_params[0]
+        actual_physical_params = actual_params[1]
+        actual_runtime_params = actual_params[2]
+        actual_kernel = actual_params[3]
+
+        self.assert_parameters_scaled_by_factor(actual_grid_params, actual_kernel,
+                                                actual_physical_params, actual_runtime_params, 2)
+
+    def assert_parameters_scaled_by_factor(self, actual_grid_params, actual_kernel,
+                                           actual_physical_params, actual_runtime_params, factor):
+        self.assert_grid_params_scaled_by_factor(actual_grid_params, factor=factor)
+        self.assert_physical_params_scaled_by_factor(actual_physical_params, factor=factor)
+        self.assert_runtime_params_scaled_by_factor(actual_runtime_params, factor=factor)
+        self.assert_kernel_forcing_scaled_by_factor(actual_kernel, factor=factor)
+
+    def assert_grid_params_scaled_by_factor(self, actual_grid_params: GridParameters, factor: int):
+        expected_nodes_per_direction = self.scaled_list(self.grid_params.number_of_nodes_per_direction, factor)
+        expected_blocks_per_direction = self.scaled_list(self.grid_params.blocks_per_direction, factor)
+        expected_node_distance = self.grid_params.node_distance / factor
+        self.assertEqual(expected_node_distance, actual_grid_params.node_distance)
+        self.assertEqual(expected_nodes_per_direction, actual_grid_params.number_of_nodes_per_direction)
+        self.assertEqual(expected_blocks_per_direction, actual_grid_params.blocks_per_direction)
+        self.assertEqual(self.grid_params.reference_direction_index, actual_grid_params.reference_direction_index)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x1, actual_grid_params.periodic_boundary_in_x1)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x2, actual_grid_params.periodic_boundary_in_x2)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x3, actual_grid_params.periodic_boundary_in_x3)
+
+    def assert_physical_params_scaled_by_factor(self, actual_params: PhysicalParameters, factor: int):
+        self.assertEqual(self.physical_params.lattice_viscosity * factor, actual_params.lattice_viscosity)
+        self.assertEqual(self.physical_params.bulk_viscosity_factor, actual_params.bulk_viscosity_factor)
+
+    def assert_runtime_params_scaled_by_factor(self, actual_params: RuntimeParameters, factor: int):
+        self.assertEqual(self.runtime_params.number_of_timesteps * factor, actual_params.number_of_timesteps)
+        self.assertEqual(self.runtime_params.number_of_threads, actual_params.number_of_threads)
+        self.assertEqual(self.runtime_params.timestep_log_interval, actual_params.timestep_log_interval)
+
+    def assert_kernel_forcing_scaled_by_factor(self, actual_kernel: LBMKernel, factor: int):
+        self.assertEqual(self.kernel.type, actual_kernel.type)
+        self.assertEqual(self.kernel.use_forcing, actual_kernel.use_forcing)
+        self.assertAlmostEqual(self.kernel.forcing_in_x1 / factor, actual_kernel.forcing_in_x1)
+        self.assertAlmostEqual(self.kernel.forcing_in_x2, actual_kernel.forcing_in_x2)
+        self.assertAlmostEqual(self.kernel.forcing_in_x3, actual_kernel.forcing_in_x3)
+
+    @staticmethod
+    def scaled_list(list_to_scale: List[int], factor: int) -> List[int]:
+        return [list_to_scale[0], list_to_scale[1], list_to_scale[2] * factor]
+
+    @staticmethod
+    def make_kernel():
+        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel.use_forcing = True
+        kernel.forcing_in_x1 = 5e-10
+        return kernel
+
+    @staticmethod
+    def make_runtime_params():
+        runtime_params = RuntimeParameters()
+        runtime_params.number_of_threads = 4
+        runtime_params.number_of_timesteps = 4_000_000
+        runtime_params.timestep_log_interval = 1_000_000
+        return runtime_params
+
+    @staticmethod
+    def make_physical_params():
+        physical_params = PhysicalParameters()
+        physical_params.lattice_viscosity = 1e-4
+        return physical_params
+
+    @staticmethod
+    def make_grid_params():
+        grid_params = GridParameters()
+        grid_params.node_distance = 1
+        grid_params.number_of_nodes_per_direction = [1, 1, 16]
+        grid_params.blocks_per_direction = [1, 1, 16]
+        grid_params.periodic_boundary_in_x1 = True
+        grid_params.periodic_boundary_in_x2 = True
+
+        return grid_params
+
+
+if __name__ == '__main__':
+    unittest.main()

Python/tests/test_boundaryconditions.py

+import unittest
+from pyfluids.boundaryconditions import *
+
+
+class BoundaryConditionsTest(unittest.TestCase):
+
+    def test__can_create_no_slip_bc(self):
+        """
+        Should be able to create NoSlipBoundaryCondition
+        """
+        sut = NoSlipBoundaryCondition()
+
+    def test__can_create_velocity_bc(self):
+        """
+        Should be able to create VelocityBoundaryCondition
+        """
+        sut = VelocityBoundaryCondition()
+
+    def test__can_create_velocity_bc_with_directions_function_and_time(self):
+        """
+        Should be able to create VelocityBoundaryCondition with directions, function and start/end time
+        """
+        from pymuparser import Parser
+
+        parser = Parser()
+        parser.expression = "1"
+        sut = VelocityBoundaryCondition(True, True, True, parser, 0, 1)
+
+    def test__can_create_velocity_bc_with_directions__function_per_direction__and__time(self):
+        """
+        Should be able to create VelocityBoundaryCondition with directions, function per direction and start/end time
+        """
+        from pymuparser import Parser
+
+        f1 = Parser()
+        f1.expression = "1"
+
+        f2 = Parser()
+        f2.expression = "1"
+
+        f3 = Parser()
+        f3.expression = "1"
+        sut = VelocityBoundaryCondition(True, True, True, f1, f2, f3, 0, 1)
+
+    def test__can_create_velocity_bc_with_speeds_and_times_per_direction(self):
+        """
+        Should be able to create VelocityBoundaryCondition with speeds and start/end times per direction
+        """
+        vx1, vx2, vx3 = 1, 2, 3
+        start1, end1 = 0, 1
+        start2, end2 = 1, 2
+        start3, end3 = 2, 3
+
+        sut = VelocityBoundaryCondition(vx1, start1, end1, vx2, start2, end2, vx3, start3, end3)
+
+    def test__can_create_non_reflecting_outflow(self):
+        """
+        Should be able to create NonReflectingOutflow
+        """
+
+        sut = NonReflectingOutflow()

Python/test_kernel.py → Python/tests/test_kernel.py

@@ -51,3 +51,10 @@ class TestLBMKernel(unittest.TestCase):
         self.assertEqual(self.sut.forcing_in_x2, 8)
         self.assertEqual(self.sut.forcing_in_x3, 5)
 
+    def test_lbm_kernel__when_getting_type__should_equal_kernel_type_enum_value(self) -> None:
+        """
+        WHEN getting the kernel type IT should equal the corresponding KernelType enum value
+        """
+
+        actual = self.sut.type
+        self.assertEqual(KernelType.BGK, actual)

Python/vtk_utilities.py

 import math
+from typing import List
+
 import pyvista as pv
 
 
 def vertical_column_from_mesh(mesh):
     last_seen = math.inf
     relevant_indices = []
+    first_x = 0
+    first_y = 0
     for index, point in enumerate(mesh.points):
-        if point[2] == last_seen:
+        if index == 0:
+            first_x = point[0]
+            first_y = point[1]
+
+        if (point[0] != first_x or point[1] != first_y) and point[2] == last_seen:
             continue
+
         relevant_indices.append(index)
         last_seen = point[2]
 
     return relevant_indices
 
 
-def get_values_from_indices(array, indices):
+def get_values_from_indices(array, indices) -> List[float]:
     return [array[index] for index in indices]
 
 

setup.py

@@ -11,6 +11,7 @@ from distutils.version import LooseVersion
 vf_cmake_args = [
     "-DBUILD_VF_PYTHON_BINDINGS=ON",
     "-DBUILD_VF_CPU:BOOL=ON",
+    "-DBUILD_VF_GPU:BOOL=OFF",
     "-DUSE_METIS=ON",
     "-DUSE_MPI=ON",
     "-DBUILD_SHARED_LIBS=OFF",