Skip to content
Snippets Groups Projects
Commit 65876888 authored by Sven Marcus's avatar Sven Marcus
Browse files

Fix l2 norm test

parent a8f3ca3b
No related branches found
No related tags found
1 merge request!6Proper Python simulation tests
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
+ 6
− 0
View file @ 65876888
...@@ -172,6 +172,11 @@ gcc_9_python: ...@@ -172,6 +172,11 @@ gcc_9_python:
needs: ["gcc_9"] needs: ["gcc_9"]
cache:
key: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG"
paths:
- $BUILD_FOLDER
artifacts: artifacts:
paths: paths:
- build/ - build/
...@@ -219,6 +224,7 @@ gcc_9_python_bindings_test: ...@@ -219,6 +224,7 @@ gcc_9_python_bindings_test:
before_script: before_script:
- export PYTHONPATH="Python" - export PYTHONPATH="Python"
- export VF_WHEEL=$(find dist/*.whl) - export VF_WHEEL=$(find dist/*.whl)
- apt-get update && apt-get install -y libglx0
- pip3 install $VF_WHEEL - pip3 install $VF_WHEEL
- pip3 install -r Python/requirements.txt - pip3 install -r Python/requirements.txt
......
Python/norms.py
+ 7
− 1
View file @ 65876888
...@@ -15,7 +15,6 @@ def root_mean_squared_error(real_values, numerical_values): ...@@ -15,7 +15,6 @@ def root_mean_squared_error(real_values, numerical_values):
raise ValueError("Real and numerical value lists must be same length") 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_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) return math.sqrt(sum_of_squared_distances / num_values)
...@@ -41,3 +40,10 @@ def mean_squared_error(real_values, numerical_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) sum_of_squared_distances = get_sum_of_squared_distances(real_values, numerical_values)
return sum_of_squared_distances / num_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
Python/poiseuille/analytical.py
+ 18
− 14
View file @ 65876888
from dataclasses import dataclass from dataclasses import dataclass
@dataclass
class PoiseuilleSettings: class PoiseuilleSettings:
density = 1
viscosity = 0.005
height = 10
length = 1
pressure_in = 0 def __init__(self):
pressure_out = 0 self.density = 1
self.viscosity = 0.005
force = 0 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): def poiseuille_at_z(settings: PoiseuilleSettings, z: float):
...@@ -27,9 +26,14 @@ def poiseuille_at_heights(settings: PoiseuilleSettings, heights): ...@@ -27,9 +26,14 @@ def poiseuille_at_heights(settings: PoiseuilleSettings, heights):
if __name__ == '__main__': if __name__ == '__main__':
h1 = [1, 2, 3, 4, 5, 6, 7, 8, 9] # 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] # 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 = PoiseuilleSettings()
settings.force = 1e-6 settings.force = 1e-8
print(max(poiseuille_at_heights(settings, h1))) settings.height = 32
print(max(poiseuille_at_heights(settings, h2)))
# 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
Python/poiseuille/simulation.py
+ 6
− 6
View file @ 65876888
...@@ -6,9 +6,9 @@ from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParam ...@@ -6,9 +6,9 @@ from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParam
from pyfluids.writer import Writer, OutputFormat from pyfluids.writer import Writer, OutputFormat
default_grid_params = GridParameters() default_grid_params = GridParameters()
default_grid_params.node_distance = 1 default_grid_params.node_distance = 10 / 32
default_grid_params.number_of_nodes_per_direction = [1, 1, 10] default_grid_params.number_of_nodes_per_direction = [8, 8, 32]
default_grid_params.blocks_per_direction = [1, 1, 1] default_grid_params.blocks_per_direction = [1, 1, 4]
default_grid_params.periodic_boundary_in_x1 = True default_grid_params.periodic_boundary_in_x1 = True
default_grid_params.periodic_boundary_in_x2 = True default_grid_params.periodic_boundary_in_x2 = True
...@@ -17,8 +17,8 @@ default_physical_params.lattice_viscosity = 0.005 ...@@ -17,8 +17,8 @@ default_physical_params.lattice_viscosity = 0.005
default_runtime_params = RuntimeParameters() default_runtime_params = RuntimeParameters()
default_runtime_params.number_of_threads = 4 default_runtime_params.number_of_threads = 4
default_runtime_params.number_of_timesteps = 1000 default_runtime_params.number_of_timesteps = 100000
default_runtime_params.timestep_log_interval = 100 default_runtime_params.timestep_log_interval = 10000
def run_simulation(physical_params=default_physical_params, def run_simulation(physical_params=default_physical_params,
...@@ -28,7 +28,7 @@ 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 = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
kernel.use_forcing = True kernel.use_forcing = True
kernel.forcing_in_x1 = 1e-6 kernel.forcing_in_x1 = 1e-8
writer = Writer() writer = Writer()
writer.output_path = "./output" writer.output_path = "./output"
......
Python/poiseuille/test_poiseuille_l2.py
+ 33
− 15
View file @ 65876888
import shutil import shutil
import unittest import unittest
import math
import pyvista as pv import pyvista as pv
import matplotlib.pyplot as plt
from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters 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.analytical import poiseuille_at_heights, PoiseuilleSettings
from poiseuille.simulation import run_simulation from poiseuille.simulation import run_simulation
from vtk_utilities import vertical_column_from_mesh, get_values_from_indices from vtk_utilities import vertical_column_from_mesh, get_values_from_indices
...@@ -16,36 +19,41 @@ class TestPoiseuilleFlow(unittest.TestCase): ...@@ -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 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 = PhysicalParameters()
physical_params.lattice_viscosity = 0.0005 physical_params.lattice_viscosity = 0.005
runtime_params = RuntimeParameters() runtime_params = RuntimeParameters()
runtime_params.number_of_threads = 4 runtime_params.number_of_threads = 4
runtime_params.timestep_log_interval = 1000 runtime_params.timestep_log_interval = 1000
runtime_params.number_of_timesteps = 10000 runtime_params.number_of_timesteps = 10000
channel_height = 10
nodes_in_column = 8 nodes_in_column = 8
grid_params = create_grid_params_with_nodes_in_column(nodes_in_column, grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
delta_x=channel_height / 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 runtime_params.number_of_timesteps *= 2
physical_params.lattice_viscosity *= 2 physical_params.lattice_viscosity *= 2
nodes_in_column *= 2 nodes_in_column *= 2
grid_params = create_grid_params_with_nodes_in_column(nodes_in_column, grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
delta_x=channel_height / 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 physical_params.lattice_viscosity *= 2
nodes_in_column *= 2 nodes_in_column *= 2
grid_params = create_grid_params_with_nodes_in_column(nodes_in_column, grid_params = create_grid_params_with_nodes_in_column(nodes_in_column,
delta_x=channel_height / 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_200 <= l2_norm_result_100)
self.assertTrue(l2_norm_result_400 <= l2_norm_result_200) 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 ...@@ -60,15 +68,20 @@ def run_simulation_with_settings(grid_params, physical_params, runtime_params, o
run_simulation(physical_params, grid_params, runtime_params) 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" 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, output_folder)
heights = get_heights(output_folder, runtime_params) heights = get_heights(output_folder, runtime_params)
numerical_results = get_numerical_results(runtime_params, output_folder, heights) 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): def get_heights(output_folder, runtime_params):
...@@ -89,6 +102,8 @@ def get_numerical_results(runtime_params, output_folder, heights): ...@@ -89,6 +102,8 @@ def get_numerical_results(runtime_params, output_folder, heights):
def calculate_analytical_results(physical_params, height_values, channel_height): def calculate_analytical_results(physical_params, height_values, channel_height):
settings = get_analytical_poiseuille_settings(channel_height, physical_params) 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) analytical_results = poiseuille_at_heights(settings, height_values)
return analytical_results return analytical_results
...@@ -113,7 +128,10 @@ def get_analytical_poiseuille_settings(height, physical_params): ...@@ -113,7 +128,10 @@ def get_analytical_poiseuille_settings(height, physical_params):
settings.height = height settings.height = height
settings.viscosity = physical_params.lattice_viscosity settings.viscosity = physical_params.lattice_viscosity
settings.density = 1 settings.density = 1
settings.force = 1e-6 settings.force = 1e-8
# print(settings)
return settings return settings
...@@ -131,8 +149,8 @@ def get_heights_from_indices(mesh, indices): ...@@ -131,8 +149,8 @@ def get_heights_from_indices(mesh, indices):
def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x): def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x):
grid_params = GridParameters() grid_params = GridParameters()
grid_params.node_distance = delta_x grid_params.node_distance = delta_x
grid_params.number_of_nodes_per_direction = [8, 8, nodes_in_column] grid_params.number_of_nodes_per_direction = [2, 2, nodes_in_column]
grid_params.blocks_per_direction = [1, 1, 4] grid_params.blocks_per_direction = [1, 1, 6]
grid_params.periodic_boundary_in_x1 = True grid_params.periodic_boundary_in_x1 = True
grid_params.periodic_boundary_in_x2 = True grid_params.periodic_boundary_in_x2 = True
grid_params.periodic_boundary_in_x3 = False grid_params.periodic_boundary_in_x3 = False
......
src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+ 56
− 36
View file @ 65876888
...@@ -47,26 +47,31 @@ ...@@ -47,26 +47,31 @@
#include "basics/writer/WbWriterVtkXmlASCII.h" #include "basics/writer/WbWriterVtkXmlASCII.h"
WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor() = default; WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor() = default;
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
const std::string &path, const std::string &path,
WbWriter *const writer, WbWriter *const writer,
SPtr<LBMUnitConverter> conv, SPtr<LBMUnitConverter> conv,
SPtr<Communicator> comm) 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(); minInitLevel = this->grid->getCoarsestInitializedLevel();
maxInitLevel = this->grid->getFinestInitializedLevel(); maxInitLevel = this->grid->getFinestInitializedLevel();
blockVector.resize(maxInitLevel + 1); 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]); grid->getBlocks(level, gridRank, true, blockVector[level]);
} }
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::init() {} void WriteMacroscopicQuantitiesCoProcessor::init()
{}
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::process(double step) void WriteMacroscopicQuantitiesCoProcessor::process(double step)
{ {
...@@ -75,14 +80,18 @@ void WriteMacroscopicQuantitiesCoProcessor::process(double step) ...@@ -75,14 +80,18 @@ void WriteMacroscopicQuantitiesCoProcessor::process(double step)
UBLOG(logDEBUG3, "WriteMacroscopicQuantitiesCoProcessor::update:" << step); UBLOG(logDEBUG3, "WriteMacroscopicQuantitiesCoProcessor::update:" << step);
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::collectData(double step) void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
{ {
int istep = static_cast<int>(step); int istep = static_cast<int>(step);
for (int level = minInitLevel; level <= maxInitLevel; level++) { for (int level = minInitLevel; level <= maxInitLevel; level++)
for (SPtr<Block3D> block : blockVector[level]) { {
if (block) { for (SPtr<Block3D> block : blockVector[level])
{
if (block)
{
addDataMQ(block); addDataMQ(block);
} }
} }
...@@ -93,26 +102,29 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step) ...@@ -93,26 +102,29 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
subfolder = "mq" + UbSystem::toString(istep); subfolder = "mq" + UbSystem::toString(istep);
pfilePath = path + "/mq/" + subfolder; pfilePath = path + "/mq/" + subfolder;
cfilePath = path + "/mq/mq_collection"; 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); std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
size_t found = partName.find_last_of("/"); size_t found = partName.find_last_of("/");
std::string piece = partName.substr(found + 1); std::string piece = partName.substr(found + 1);
piece = subfolder + "/" + piece; piece = subfolder + "/" + piece;
std::vector<std::string> cellDataNames; std::vector<std::string> cellDataNames;
std::vector<std::string> pieces = comm->gather(piece); std::vector<std::string> pieces = comm->gather(piece);
if (comm->getProcessID() == comm->getRoot()) { if (comm->getProcessID() == comm->getRoot())
{
std::string pname = std::string pname =
WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames); WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
found = pname.find_last_of("/"); found = pname.find_last_of("/");
piece = pname.substr(found + 1); piece = pname.substr(found + 1);
std::vector<std::string> filenames; std::vector<std::string> filenames;
filenames.push_back(piece); filenames.push_back(piece);
if (step == CoProcessor::scheduler->getMinBegin()) { if (step == CoProcessor::scheduler->getMinBegin())
{
WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false); WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
} else { } else
{
WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false); WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
} }
UBLOG(logINFO, "WriteMacroscopicQuantitiesCoProcessor step: " << istep); UBLOG(logINFO, "WriteMacroscopicQuantitiesCoProcessor step: " << istep);
...@@ -120,6 +132,7 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step) ...@@ -120,6 +132,7 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
clearData(); clearData();
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::clearData() void WriteMacroscopicQuantitiesCoProcessor::clearData()
{ {
...@@ -128,10 +141,11 @@ void WriteMacroscopicQuantitiesCoProcessor::clearData() ...@@ -128,10 +141,11 @@ void WriteMacroscopicQuantitiesCoProcessor::clearData()
datanames.clear(); datanames.clear();
data.clear(); data.clear();
} }
////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
{ {
double level = (double)block->getLevel(); double level = (double) block->getLevel();
// double blockID = (double)block->getGlobalID(); // double blockID = (double)block->getGlobalID();
// Diese Daten werden geschrieben: // Diese Daten werden geschrieben:
...@@ -146,18 +160,20 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) ...@@ -146,18 +160,20 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
data.resize(datanames.size()); data.resize(datanames.size());
SPtr<ILBMKernel> kernel = block->getKernel(); SPtr<ILBMKernel> kernel = block->getKernel();
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray(); SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions(); SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
LBMReal f[D3Q27System::ENDF + 1]; LBMReal f[D3Q27System::ENDF + 1];
LBMReal vx1, vx2, vx3, rho; LBMReal vx1, vx2, vx3, rho;
// knotennummerierung faengt immer bei 0 an! // 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; calcMacros = &D3Q27System::calcCompMacroscopicValues;
} else { } else
{
calcMacros = &D3Q27System::calcIncompMacroscopicValues; calcMacros = &D3Q27System::calcIncompMacroscopicValues;
} }
...@@ -165,9 +181,9 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) ...@@ -165,9 +181,9 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
int minX2 = 0; int minX2 = 0;
int minX3 = 0; int minX3 = 0;
int maxX1 = (int)(distributions->getNX1()); int maxX1 = (int) (distributions->getNX1());
int maxX2 = (int)(distributions->getNX2()); int maxX2 = (int) (distributions->getNX2());
int maxX3 = (int)(distributions->getNX3()); int maxX3 = (int) (distributions->getNX3());
// int minX1 = 1; // int minX1 = 1;
// int minX2 = 1; // int minX2 = 1;
...@@ -178,21 +194,25 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) ...@@ -178,21 +194,25 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
// int maxX3 = (int)(distributions->getNX3()); // int maxX3 = (int)(distributions->getNX3());
// nummern vergeben und node vector erstellen + daten sammeln // 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; maxX1 -= 2;
maxX2 -= 2; maxX2 -= 2;
maxX3 -= 2; maxX3 -= 2;
// D3Q27BoundaryConditionPtr bcPtr; // D3Q27BoundaryConditionPtr bcPtr;
int nr = (int)nodes.size(); int nr = (int) nodes.size();
for (int ix3 = minX3; ix3 <= maxX3; ix3++) { for (int ix3 = minX3; ix3 <= maxX3; ix3++)
for (int ix2 = minX2; ix2 <= maxX2; ix2++) { {
for (int ix1 = minX1; ix1 <= maxX1; ix1++) { for (int ix2 = minX2; ix2 <= maxX2; ix2++)
if (!bcArray->isUndefined(ix1, ix2, ix3) && !bcArray->isSolid(ix1, ix2, ix3)) { {
int index = 0; 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++; 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]), nodes.emplace_back(float(worldCoordinates[0]), float(worldCoordinates[1]),
float(worldCoordinates[2])); float(worldCoordinates[2]));
...@@ -202,7 +222,7 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) ...@@ -202,7 +222,7 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
if (UbMath::isNaN(rho) || UbMath::isInfinity(rho)) if (UbMath::isNaN(rho) || UbMath::isInfinity(rho))
UB_THROW(UbException( 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) + "," + block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
UbSystem::toString(ix2) + "," + UbSystem::toString(ix3))); UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
// rho=999.0; // rho=999.0;
...@@ -213,19 +233,19 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block) ...@@ -213,19 +233,19 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
// press=999.0; // press=999.0;
if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1)) if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
UB_THROW(UbException( 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) + "," + block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
UbSystem::toString(ix2) + "," + UbSystem::toString(ix3))); UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
// vx1=999.0; // vx1=999.0;
if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2)) if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
UB_THROW(UbException( 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) + "," + block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
UbSystem::toString(ix2) + "," + UbSystem::toString(ix3))); UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
// vx2=999.0; // vx2=999.0;
if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3)) if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
UB_THROW(UbException( 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) + "," + block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
UbSystem::toString(ix2) + "," + UbSystem::toString(ix3))); UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
......
src/cpu/simulationconfig/src/Simulation.cpp
+ 4
− 3
View file @ 65876888
...@@ -117,7 +117,7 @@ void Simulation::run() ...@@ -117,7 +117,7 @@ void Simulation::run()
auto metisVisitor = std::make_shared<MetisPartitioningGridVisitor>(communicator, auto metisVisitor = std::make_shared<MetisPartitioningGridVisitor>(communicator,
MetisPartitioningGridVisitor::LevelBased, MetisPartitioningGridVisitor::LevelBased,
D3Q27System::B); D3Q27System::B, MetisPartitioner::RECURSIVE);
InteractorsHelper intHelper(grid, metisVisitor); InteractorsHelper intHelper(grid, metisVisitor);
for (auto const &interactor : interactors) for (auto const &interactor : interactors)
...@@ -125,8 +125,7 @@ void Simulation::run() ...@@ -125,8 +125,7 @@ void Simulation::run()
intHelper.selectBlocks(); intHelper.selectBlocks();
// important: run this after metis & intHelper.selectBlocks()
writeBlocksToFile();
int numberOfProcesses = communicator->getNumberOfProcesses(); int numberOfProcesses = communicator->getNumberOfProcesses();
SetKernelBlockVisitor kernelVisitor(lbmKernel, physicalParameters->latticeViscosity, SetKernelBlockVisitor kernelVisitor(lbmKernel, physicalParameters->latticeViscosity,
numberOfProcesses); numberOfProcesses);
...@@ -147,6 +146,8 @@ void Simulation::run() ...@@ -147,6 +146,8 @@ void Simulation::run()
grid->accept(bcVisitor); grid->accept(bcVisitor);
writeBoundaryConditions(); writeBoundaryConditions();
// important: run this after metis & intHelper.selectBlocks()
writeBlocksToFile();
auto visualizationScheduler = std::make_shared<UbScheduler>(simulationParameters->timeStepLogInterval); auto visualizationScheduler = std::make_shared<UbScheduler>(simulationParameters->timeStepLogInterval);
auto mqCoProcessor = makeMacroscopicQuantitiesCoProcessor(converter, auto mqCoProcessor = makeMacroscopicQuantitiesCoProcessor(converter,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment