From fa05055b844c0260c5e3b4af144085331f3733e3 Mon Sep 17 00:00:00 2001
From: Sven Marcus <s.marcus@outlook.de>
Date: Fri, 11 Dec 2020 19:45:22 +0100
Subject: [PATCH] Experimenting with bounding box
---
Python/norms.py | 39 ++++++-
Python/poiseuille/analytical.py | 9 ++
Python/poiseuille/simulation.py | 2 +-
Python/poiseuille/test_poiseuille_l2.py | 106 ++++++++++++++++----
src/cpu/simulationconfig/src/Simulation.cpp | 10 +-
5 files changed, 136 insertions(+), 30 deletions(-)
diff --git a/Python/norms.py b/Python/norms.py
index 40018ac01..cd4666ed8 100644
--- a/Python/norms.py
+++ b/Python/norms.py
@@ -1,14 +1,43 @@
import math
-def l2_norm(real_values, numerical_values):
+def get_sum_of_squared_distances(real_values, numerical_values):
+ combined_values = zip(real_values, numerical_values)
+ sum_of_squared_distances = sum((numerical_value - real_value) ** 2
+ for real_value, numerical_value
+ in combined_values)
+ return sum_of_squared_distances
+
+
+def root_mean_squared_error(real_values, numerical_values):
+ num_values = len(real_values)
+ if num_values != len(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)
+
+
+def mean_absolute_error(real_values, numerical_values):
num_values = len(real_values)
if num_values != len(numerical_values):
raise ValueError("Real and numerical value lists must be same length")
combined_values = zip(real_values, numerical_values)
- sum_of_squared_distances = sum((real_value - numerical_value) ** 2
- for real_value, numerical_value
- in combined_values)
+ sum_of_absolute_distances = sum(abs(numerical_value - real_value)
+ for real_value, numerical_value
+ in combined_values)
+
+ return sum_of_absolute_distances / num_values
+
+
+def mean_squared_error(real_values, numerical_values):
+ num_values = len(real_values)
+ if num_values != len(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)
- return math.sqrt(1 / num_values * sum_of_squared_distances)
+ return sum_of_squared_distances / num_values
diff --git a/Python/poiseuille/analytical.py b/Python/poiseuille/analytical.py
index 4ee428a4a..fb5ea9ad0 100644
--- a/Python/poiseuille/analytical.py
+++ b/Python/poiseuille/analytical.py
@@ -24,3 +24,12 @@ def poiseuille_at_z(settings: PoiseuilleSettings, z: float):
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-6
+ print(max(poiseuille_at_heights(settings, h1)))
+ print(max(poiseuille_at_heights(settings, h2)))
diff --git a/Python/poiseuille/simulation.py b/Python/poiseuille/simulation.py
index 40c40698a..13966fc9e 100644
--- a/Python/poiseuille/simulation.py
+++ b/Python/poiseuille/simulation.py
@@ -7,7 +7,7 @@ from pyfluids.writer import Writer, OutputFormat
grid_params = GridParameters()
grid_params.node_distance = 1
-grid_params.number_of_nodes_per_direction = [2, 2, 10]
+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
diff --git a/Python/poiseuille/test_poiseuille_l2.py b/Python/poiseuille/test_poiseuille_l2.py
index e9391f0c3..852d458f4 100644
--- a/Python/poiseuille/test_poiseuille_l2.py
+++ b/Python/poiseuille/test_poiseuille_l2.py
@@ -1,11 +1,16 @@
import unittest
+import shutil
+import math
import pyvista as pv
-from norms import l2_norm
+from norms import root_mean_squared_error, mean_squared_error, mean_absolute_error
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
from poiseuille.analytical import poiseuille_at_heights, PoiseuilleSettings
-from poiseuille.simulation import run_simulation, grid_params, physical_params, runtime_params
+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):
@@ -14,39 +19,100 @@ class TestPoiseuilleFlow(unittest.TestCase):
WHEN comparing the simulation results to the analytical solution THEN the L2-Norm should be less than 1e-4
"""
- run_simulation()
- file_name_of_last_timestep = get_output_file_name(runtime_params)
- mesh_of_last_timestep = pv.read(file_name_of_last_timestep)
- column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
- numerical_results_from_single_column = get_values_from_indices(mesh_of_last_timestep.get_array("Vx"), column_indices)
+ physical_params = PhysicalParameters()
+ physical_params.lattice_viscosity = 0.0005
+
+ runtime_params = RuntimeParameters()
+ runtime_params.number_of_threads = 4
+ 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)
+
+ runtime_params.number_of_timesteps = 20000
+ 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)
+
+ # runtime_params.number_of_timesteps = 40000
+ # 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 = [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)
+
+
+def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params, total_height):
+ output_folder = "./output"
+ remove_existing_output_directory(output_folder)
+ run_simulation(physical_params, grid_params, 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)
- heights = get_heights_from_indices(mesh_of_last_timestep, column_indices)
- settings = get_analytical_poiseuille_settings(grid_params, physical_params)
- analytical_results = poiseuille_at_heights(settings, 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()
- l2_norm_result = l2_norm(analytical_results, numerical_results_from_single_column)
+ return root_mean_squared_error(analytical_results, numerical_results)
- max_acceptable_error = 1e-4
- self.assertLessEqual(l2_norm_result, max_acceptable_error)
- print(f"The L2-Norm is: {l2_norm_result}")
+def get_analytical_results(physical_params, total_height, height_values):
+ settings = get_analytical_poiseuille_settings(total_height, physical_params)
+ analytical_results = poiseuille_at_heights(settings, height_values)
+ return analytical_results
-def get_analytical_poiseuille_settings(grid_params, physical_params):
+
+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 remove_existing_output_directory(output_dir):
+ shutil.rmtree(output_dir, ignore_errors=True)
+
+
+def get_analytical_poiseuille_settings(height, physical_params):
settings = PoiseuilleSettings()
- settings.length = grid_params.number_of_nodes_per_direction[0]
- settings.height = grid_params.number_of_nodes_per_direction[2]
+ settings.height = height
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(rumtime_params):
- timesteps = rumtime_params.number_of_timesteps
- file_name = f"output/mq/mq{timesteps}/mq0_{timesteps}.bin.vtu"
+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_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):
+ 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.periodic_boundary_in_x1 = True
+ grid_params.periodic_boundary_in_x2 = True
+
+ return grid_params
diff --git a/src/cpu/simulationconfig/src/Simulation.cpp b/src/cpu/simulationconfig/src/Simulation.cpp
index 558f6d3c2..a3574fedc 100644
--- a/src/cpu/simulationconfig/src/Simulation.cpp
+++ b/src/cpu/simulationconfig/src/Simulation.cpp
@@ -258,10 +258,12 @@ std::shared_ptr<GbObject3D>
Simulation::makeSimulationBoundingBox(const int &nodesInX1, const int &nodesInX2,
const int &nodesInX3) 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;
+
+ double halfDx = -gridParameters->nodeDistance / 2.0;
+ double minX1 = halfDx, minX2 = halfDx, minX3 = halfDx;
+ const double maxX1 = minX1 + gridParameters->nodeDistance * (nodesInX1 - 1) + halfDx;
+ const double maxX2 = minX2 + gridParameters->nodeDistance * (nodesInX2 - 1) + halfDx;
+ const double maxX3 = minX3 + gridParameters->nodeDistance * (nodesInX3 - 1) + halfDx;
UBLOG(logINFO, "Bounding box dimensions = [("
<< minX1 << ", " << minX2 << ", " << minX3 << "); ("
<< maxX1 << ", " << maxX2 << ", " << maxX3 << ")]")
--
GitLab