From ce5b4cd7988a8edbe2cf7805b8ec75d8530bfdad Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 24 Feb 2024 16:20:15 +0100 Subject: [PATCH 01/13] Added subprogram to calculate the density contrast PDF for turbulence simulations --- main/src/analytical_solutions/CMakeLists.txt | 1 + .../turbulence/CMakeLists.txt | 12 ++ .../turbulence/density_pdf.cpp | 104 ++++++++++++++++++ .../turbulence/density_pdf.hpp | 48 ++++++++ 4 files changed, 165 insertions(+) create mode 100644 main/src/analytical_solutions/turbulence/CMakeLists.txt create mode 100644 main/src/analytical_solutions/turbulence/density_pdf.cpp create mode 100644 main/src/analytical_solutions/turbulence/density_pdf.hpp diff --git a/main/src/analytical_solutions/CMakeLists.txt b/main/src/analytical_solutions/CMakeLists.txt index e088a7d97..b75efad96 100644 --- a/main/src/analytical_solutions/CMakeLists.txt +++ b/main/src/analytical_solutions/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(sedov_solution) +add_subdirectory(turbulence) install(FILES compare_solutions.py TYPE BIN) install(FILES compare_noh.py TYPE BIN) diff --git a/main/src/analytical_solutions/turbulence/CMakeLists.txt b/main/src/analytical_solutions/turbulence/CMakeLists.txt new file mode 100644 index 000000000..de0251b4e --- /dev/null +++ b/main/src/analytical_solutions/turbulence/CMakeLists.txt @@ -0,0 +1,12 @@ +set(SOURCES density_pdf.hpp density_pdf.cpp) +add_executable(density_pdf ${SOURCES}) +target_include_directories(density_pdf PRIVATE ${CSTONE_DIR} ${PROJECT_SOURCE_DIR}/main/src ${CMAKE_BINARY_DIR}/main/src) +target_link_libraries(density_pdf io stdc++fs OpenMP::OpenMP_CXX ${MPI_CXX_LIBRARIES}) + +install(TARGETS density_pdf RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) + +include(setup_gitinfo) +configure_file( + ${PROJECT_SOURCE_DIR}/main/src/version.h.in + ${CMAKE_BINARY_DIR}/main/src/version.h +) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp new file mode 100644 index 000000000..dffed8e92 --- /dev/null +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -0,0 +1,104 @@ +/* + * MIT License + * + * Copyright (c) 2024 CSCS, ETH Zurich + * 2024 University of Basel + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ + +#include +#include + +#include "density_pdf.hpp" +#include "io/arg_parser.hpp" +#include "io/factory.hpp" +#include "util/utils.hpp" +#include "cstone/primitives/mpi_wrappers.hpp" + +using namespace sphexa; +using T = float; + +int main(int argc, char** argv) +{ + auto [rank, numRanks] = initMpi(); + const ArgParser parser(argc, (const char**)argv); + + const std::string inputFile = parser.get("--file"); + const size_t nBins = parser.get("-n", 50); + const T referenceDensity = parser.get("--rho0", 1.0); + std::string outputFile = parser.get("-o", std::string("density_pdf.txt")); + const std::string sph_type = parser.get("--sph", std::string("std")); + const T minValue = parser.get("--min", -6.0); + const T maxValue = parser.get("--max", 8.0); + + if (!std::filesystem::exists(inputFile)) + { + printf("Please provide a existing file: no file found at %s\n", inputFile.c_str()); + return exitSuccess(); + } + + auto h5reader = fileReaderFactory(false, MPI_COMM_WORLD); + h5reader->setStep(inputFile, -1, FileMode::collective); + + const size_t localNumParticles = h5reader->localNumParticles(); + const size_t globalNumParticles = h5reader->globalNumParticles(); + std::vector rho(localNumParticles); + std::vector bins; + + if (sph_type == "std") { h5reader->readField("rho", rho.data()); } + else + { + std::vector m(localNumParticles); + std::vector xm(localNumParticles); + h5reader->readField("kx", rho.data()); + h5reader->readField("xm", xm.data()); + h5reader->readField("m", m.data()); +#pragma omp for schedule(static) + for (size_t i = 0; i < localNumParticles; ++i) + { + rho[i] = rho[i] * m[i] / xm[i]; + } + } + + h5reader->closeStep(); + + bins = computeProbabilityDistribution(rho, referenceDensity, nBins, minValue, maxValue); + std::vector reduced_bins(nBins, 0.0); + MPI_Reduce(bins.data(), reduced_bins.data(), nBins, MpiType{}, MPI_SUM, 0, MPI_COMM_WORLD); + + if (rank == 0) + { + std::for_each(reduced_bins.begin(), reduced_bins.end(), + [globalNumParticles](T& i) { i /= globalNumParticles; }); + std::ofstream outFile(std::filesystem::path(outputFile), std::ofstream::out); + + T binSize = (maxValue - minValue) / nBins; + T firstMiddle = minValue + 0.5 * binSize; + + for (size_t i = 0; i < nBins; i++) + { + T binCenter = binSize * i + firstMiddle; + outFile << binCenter << ' ' << reduced_bins[i] << std::endl; + } + printf("Calculated PDF for %lu particles in %lu bins.\n", globalNumParticles, nBins); + } + + exitSuccess(); +} diff --git a/main/src/analytical_solutions/turbulence/density_pdf.hpp b/main/src/analytical_solutions/turbulence/density_pdf.hpp new file mode 100644 index 000000000..b85f7f1c4 --- /dev/null +++ b/main/src/analytical_solutions/turbulence/density_pdf.hpp @@ -0,0 +1,48 @@ +/* + * MIT License + * + * Copyright (c) 2024 CSCS, ETH Zurich + * 2024 University of Basel + * + * Permission is hereby granted, free of charge, to any person obtaining a copy + * of this software and associated documentation files (the "Software"), to deal + * in the Software without restriction, including without limitation the rights + * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included in all + * copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + * SOFTWARE. + */ +#include +#include +#include + +template +std::vector computeProbabilityDistribution(std::vector& data, const T referenceValue, size_t binCount, T binStart, + T binEnd) +{ + std::vector bins(binCount); + +#pragma omp for schedule(static) + for (size_t i = 0; i < data.size(); ++i) + { + data[i] = std::log(data[i] / referenceValue); + } + + T binSize = (binEnd - binStart) / binCount; + for (size_t bin = 0; bin < binCount; bin++) + { + bins[bin] = std::count_if(data.begin(), data.end(), [bin, binStart, binSize](T i) + { return i > binSize * bin + binStart && i <= binSize * (bin + 1) + binStart; }); + } + return bins; +} From d16a22bb35bb86d811f7d5faf22e13d93e9aef7d Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Tue, 27 Feb 2024 22:07:35 +0100 Subject: [PATCH 02/13] use average density to calculate density-contrast PDF --- .../turbulence/density_pdf.cpp | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index dffed8e92..f48e64aca 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -25,6 +25,7 @@ #include #include +#include #include "density_pdf.hpp" #include "io/arg_parser.hpp" @@ -40,13 +41,12 @@ int main(int argc, char** argv) auto [rank, numRanks] = initMpi(); const ArgParser parser(argc, (const char**)argv); - const std::string inputFile = parser.get("--file"); - const size_t nBins = parser.get("-n", 50); - const T referenceDensity = parser.get("--rho0", 1.0); - std::string outputFile = parser.get("-o", std::string("density_pdf.txt")); - const std::string sph_type = parser.get("--sph", std::string("std")); - const T minValue = parser.get("--min", -6.0); - const T maxValue = parser.get("--max", 8.0); + const std::string inputFile = parser.get("--file"); + const size_t nBins = parser.get("-n", 50); + std::string outputFile = parser.get("-o", std::string("density_pdf.txt")); + const std::string sph_type = parser.get("--sph", std::string("std")); + const T minValue = parser.get("--min", -8.0); + const T maxValue = parser.get("--max", 6.0); if (!std::filesystem::exists(inputFile)) { @@ -76,9 +76,13 @@ int main(int argc, char** argv) rho[i] = rho[i] * m[i] / xm[i]; } } - h5reader->closeStep(); + T localTotalDensity = std::reduce(rho.begin(), rho.end()); + T referenceDensity; + MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); + referenceDensity /= globalNumParticles; + bins = computeProbabilityDistribution(rho, referenceDensity, nBins, minValue, maxValue); std::vector reduced_bins(nBins, 0.0); MPI_Reduce(bins.data(), reduced_bins.data(), nBins, MpiType{}, MPI_SUM, 0, MPI_COMM_WORLD); @@ -92,6 +96,9 @@ int main(int argc, char** argv) T binSize = (maxValue - minValue) / nBins; T firstMiddle = minValue + 0.5 * binSize; + // header line containing metadata + outFile << nBins << ' ' << binSize << ' ' << referenceDensity << std::endl; + for (size_t i = 0; i < nBins; i++) { T binCenter = binSize * i + firstMiddle; From 4531d9d721d83cb58f02854e2c4e82e96a45d206 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 20 Apr 2024 18:28:49 +0200 Subject: [PATCH 03/13] Density PDF: adressed review, formatting --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 5 +++-- main/src/analytical_solutions/turbulence/density_pdf.hpp | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index f48e64aca..aea33fe72 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -41,6 +41,7 @@ int main(int argc, char** argv) auto [rank, numRanks] = initMpi(); const ArgParser parser(argc, (const char**)argv); + // The default pdf range was taken in comparison to Federrath et al. 2021, DOI 10.1038/s41550-020-01282-z const std::string inputFile = parser.get("--file"); const size_t nBins = parser.get("-n", 50); std::string outputFile = parser.get("-o", std::string("density_pdf.txt")); @@ -89,11 +90,11 @@ int main(int argc, char** argv) if (rank == 0) { + T binSize = (maxValue - minValue) / nBins; std::for_each(reduced_bins.begin(), reduced_bins.end(), - [globalNumParticles](T& i) { i /= globalNumParticles; }); + [globalNumParticles, binSize](T& i) { i /= globalNumParticles * binSize; }); std::ofstream outFile(std::filesystem::path(outputFile), std::ofstream::out); - T binSize = (maxValue - minValue) / nBins; T firstMiddle = minValue + 0.5 * binSize; // header line containing metadata diff --git a/main/src/analytical_solutions/turbulence/density_pdf.hpp b/main/src/analytical_solutions/turbulence/density_pdf.hpp index b85f7f1c4..afd5ffec3 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.hpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.hpp @@ -41,7 +41,8 @@ std::vector computeProbabilityDistribution(std::vector& data, const T refe T binSize = (binEnd - binStart) / binCount; for (size_t bin = 0; bin < binCount; bin++) { - bins[bin] = std::count_if(data.begin(), data.end(), [bin, binStart, binSize](T i) + bins[bin] = std::count_if(data.begin(), data.end(), + [bin, binStart, binSize](T i) { return i > binSize * bin + binStart && i <= binSize * (bin + 1) + binStart; }); } return bins; From 7e5d326a8d1e56bbfd4b03e3e0dab389f49e1899 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Fri, 26 Apr 2024 19:04:38 +0200 Subject: [PATCH 04/13] density pdf: add step as a runtime argument, check count of loaded particles --- .../turbulence/density_pdf.cpp | 26 ++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index aea33fe72..0b141744d 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -44,6 +44,7 @@ int main(int argc, char** argv) // The default pdf range was taken in comparison to Federrath et al. 2021, DOI 10.1038/s41550-020-01282-z const std::string inputFile = parser.get("--file"); const size_t nBins = parser.get("-n", 50); + const int step = parser.get("-s", -1); std::string outputFile = parser.get("-o", std::string("density_pdf.txt")); const std::string sph_type = parser.get("--sph", std::string("std")); const T minValue = parser.get("--min", -8.0); @@ -55,13 +56,20 @@ int main(int argc, char** argv) return exitSuccess(); } + std::vector rho; + std::vector bins; + auto h5reader = fileReaderFactory(false, MPI_COMM_WORLD); - h5reader->setStep(inputFile, -1, FileMode::collective); + h5reader->setStep(inputFile, step, FileMode::collective); - const size_t localNumParticles = h5reader->localNumParticles(); - const size_t globalNumParticles = h5reader->globalNumParticles(); - std::vector rho(localNumParticles); - std::vector bins; + const size_t localNumParticles = h5reader->localNumParticles(); + const size_t globalNumParticles = h5reader->globalNumParticles(); + + rho.resize(localNumParticles); + if (rho.size() != localNumParticles) + { + throw new std::runtime_error("rho length doesn't match local count: " + rho.size() + localNumParticles); + } if (sph_type == "std") { h5reader->readField("rho", rho.data()); } else @@ -79,11 +87,17 @@ int main(int argc, char** argv) } h5reader->closeStep(); - T localTotalDensity = std::reduce(rho.begin(), rho.end()); + T localTotalDensity = std::reduce(rho.begin(), rho.end(), 0); T referenceDensity; MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; + if (rank == 0) + { + printf("starting PDF calculation with %lu global, %lu local particles and reference density %f\n", + globalNumParticles, localNumParticles, referenceDensity); + } + bins = computeProbabilityDistribution(rho, referenceDensity, nBins, minValue, maxValue); std::vector reduced_bins(nBins, 0.0); MPI_Reduce(bins.data(), reduced_bins.data(), nBins, MpiType{}, MPI_SUM, 0, MPI_COMM_WORLD); From 427d3b562fdcfa8ef8bb460dea7a40417973d9c0 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 27 Apr 2024 11:23:56 +0200 Subject: [PATCH 05/13] density pdf: print more info at runtime --- .../analytical_solutions/turbulence/density_pdf.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index 0b141744d..1e84b14be 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -64,6 +64,10 @@ int main(int argc, char** argv) const size_t localNumParticles = h5reader->localNumParticles(); const size_t globalNumParticles = h5reader->globalNumParticles(); + if (rank == 0) + { + printf("Density-PDF: local particles: %lu \t global particles: %lu\n", localNumParticles, globalNumParticles); + } rho.resize(localNumParticles); if (rho.size() != localNumParticles) @@ -92,11 +96,7 @@ int main(int argc, char** argv) MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; - if (rank == 0) - { - printf("starting PDF calculation with %lu global, %lu local particles and reference density %f\n", - globalNumParticles, localNumParticles, referenceDensity); - } + if (rank == 0) { printf("starting PDF calculation with reference density %f\n", referenceDensity); } bins = computeProbabilityDistribution(rho, referenceDensity, nBins, minValue, maxValue); std::vector reduced_bins(nBins, 0.0); From b075696f5ee72d36ae933f4a431420728356e7d6 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 27 Apr 2024 11:36:39 +0200 Subject: [PATCH 06/13] density pdf: parallelize local density reduction --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index 1e84b14be..aca6007d0 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -91,7 +91,13 @@ int main(int argc, char** argv) } h5reader->closeStep(); - T localTotalDensity = std::reduce(rho.begin(), rho.end(), 0); + T localTotalDensity = 0.0; +#pragma omp parallel for reduction(+ : localTotalDensity) + for (size_t i = 0; i < localNumParticles; ++i) + { + localTotalDensity += rho[i]; + } + T referenceDensity; MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; From 2d5b7575a9b883a664e11132b0f2337fb430b681 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 27 Apr 2024 11:46:50 +0200 Subject: [PATCH 07/13] density pdf: fix appending int to string --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index aca6007d0..54661a080 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -72,7 +72,8 @@ int main(int argc, char** argv) rho.resize(localNumParticles); if (rho.size() != localNumParticles) { - throw new std::runtime_error("rho length doesn't match local count: " + rho.size() + localNumParticles); + throw new std::runtime_error("rho length doesn't match local count: " + std::to_string(rho.size()) + "\t" + + std::to_string(localNumParticles)); } if (sph_type == "std") { h5reader->readField("rho", rho.data()); } From 9128c22f9c38fab851fba8755da6fc8779288eb2 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 27 Apr 2024 20:00:56 +0200 Subject: [PATCH 08/13] density pdf: additional debug info --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index 54661a080..dbf27b7db 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -72,8 +72,8 @@ int main(int argc, char** argv) rho.resize(localNumParticles); if (rho.size() != localNumParticles) { - throw new std::runtime_error("rho length doesn't match local count: " + std::to_string(rho.size()) + "\t" + - std::to_string(localNumParticles)); + throw std::runtime_error("rho length doesn't match local count: " + std::to_string(rho.size()) + "\t" + + std::to_string(localNumParticles)); } if (sph_type == "std") { h5reader->readField("rho", rho.data()); } @@ -99,7 +99,8 @@ int main(int argc, char** argv) localTotalDensity += rho[i]; } - T referenceDensity; + printf("rank %i, local average density: %f ", rank, localTotalDensity / localNumParticles); + T referenceDensity = 0.0; MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; From 9e44f3cbd7475fb10276928975c9d85f5fecbfb5 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sat, 27 Apr 2024 20:36:35 +0200 Subject: [PATCH 09/13] density pdf: test with mpi barrier --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index dbf27b7db..7a9d6c38d 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -99,7 +99,8 @@ int main(int argc, char** argv) localTotalDensity += rho[i]; } - printf("rank %i, local average density: %f ", rank, localTotalDensity / localNumParticles); + MPI_Barrier(MPI_COMM_WORLD); + printf("rank %i, local average density: %f \n", rank, localTotalDensity / localNumParticles); T referenceDensity = 0.0; MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; From 156f7601f2099c8196e586051deed3143544414a Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Sun, 28 Apr 2024 15:54:56 +0200 Subject: [PATCH 10/13] density pdf: check size after load --- .../analytical_solutions/turbulence/density_pdf.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index 7a9d6c38d..df50edc58 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -70,11 +70,6 @@ int main(int argc, char** argv) } rho.resize(localNumParticles); - if (rho.size() != localNumParticles) - { - throw std::runtime_error("rho length doesn't match local count: " + std::to_string(rho.size()) + "\t" + - std::to_string(localNumParticles)); - } if (sph_type == "std") { h5reader->readField("rho", rho.data()); } else @@ -90,6 +85,13 @@ int main(int argc, char** argv) rho[i] = rho[i] * m[i] / xm[i]; } } + rho.shrink_to_fit(); + if (rho.size() != localNumParticles) + { + throw std::runtime_error("rho length doesn't match local count: " + std::to_string(rho.size()) + "\t" + + std::to_string(localNumParticles)); + } + h5reader->closeStep(); T localTotalDensity = 0.0; From 272353413e8d64c660a5f5d88de10d6e94201ef0 Mon Sep 17 00:00:00 2001 From: Sebastian Keller Date: Mon, 29 Apr 2024 18:24:29 +0200 Subject: [PATCH 11/13] added missing MPI_CXX_INCLUDE_DIR in CMake --- main/src/analytical_solutions/turbulence/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/src/analytical_solutions/turbulence/CMakeLists.txt b/main/src/analytical_solutions/turbulence/CMakeLists.txt index de0251b4e..ae748b535 100644 --- a/main/src/analytical_solutions/turbulence/CMakeLists.txt +++ b/main/src/analytical_solutions/turbulence/CMakeLists.txt @@ -1,6 +1,6 @@ set(SOURCES density_pdf.hpp density_pdf.cpp) add_executable(density_pdf ${SOURCES}) -target_include_directories(density_pdf PRIVATE ${CSTONE_DIR} ${PROJECT_SOURCE_DIR}/main/src ${CMAKE_BINARY_DIR}/main/src) +target_include_directories(density_pdf PRIVATE ${CSTONE_DIR} ${PROJECT_SOURCE_DIR}/main/src ${CMAKE_BINARY_DIR}/main/src ${MPI_CXX_INCLUDE_PATH}) target_link_libraries(density_pdf io stdc++fs OpenMP::OpenMP_CXX ${MPI_CXX_LIBRARIES}) install(TARGETS density_pdf RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) From e658e8dd1e8db0947a2cb7fce9e7a09c2bdd5a63 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Thu, 2 May 2024 15:53:01 +0200 Subject: [PATCH 12/13] density pdf: convert analysis type to double --- .../turbulence/density_pdf.cpp | 15 +++++++-------- .../turbulence/density_pdf.hpp | 4 ++-- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index df50edc58..fe7e41fb1 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -57,7 +57,7 @@ int main(int argc, char** argv) } std::vector rho; - std::vector bins; + std::vector bins; auto h5reader = fileReaderFactory(false, MPI_COMM_WORLD); h5reader->setStep(inputFile, step, FileMode::collective); @@ -94,30 +94,29 @@ int main(int argc, char** argv) h5reader->closeStep(); - T localTotalDensity = 0.0; + double localTotalDensity = 0.0; #pragma omp parallel for reduction(+ : localTotalDensity) for (size_t i = 0; i < localNumParticles; ++i) { localTotalDensity += rho[i]; } - MPI_Barrier(MPI_COMM_WORLD); printf("rank %i, local average density: %f \n", rank, localTotalDensity / localNumParticles); - T referenceDensity = 0.0; - MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); + double referenceDensity = 0.0; + MPI_Allreduce(&localTotalDensity, &referenceDensity, 1, MpiType{}, MPI_SUM, MPI_COMM_WORLD); referenceDensity /= globalNumParticles; if (rank == 0) { printf("starting PDF calculation with reference density %f\n", referenceDensity); } bins = computeProbabilityDistribution(rho, referenceDensity, nBins, minValue, maxValue); - std::vector reduced_bins(nBins, 0.0); - MPI_Reduce(bins.data(), reduced_bins.data(), nBins, MpiType{}, MPI_SUM, 0, MPI_COMM_WORLD); + std::vector reduced_bins(nBins, 0.0); + MPI_Reduce(bins.data(), reduced_bins.data(), nBins, MpiType{}, MPI_SUM, 0, MPI_COMM_WORLD); if (rank == 0) { T binSize = (maxValue - minValue) / nBins; std::for_each(reduced_bins.begin(), reduced_bins.end(), - [globalNumParticles, binSize](T& i) { i /= globalNumParticles * binSize; }); + [globalNumParticles, binSize](double& i) { i /= globalNumParticles * binSize; }); std::ofstream outFile(std::filesystem::path(outputFile), std::ofstream::out); T firstMiddle = minValue + 0.5 * binSize; diff --git a/main/src/analytical_solutions/turbulence/density_pdf.hpp b/main/src/analytical_solutions/turbulence/density_pdf.hpp index afd5ffec3..a3ee102b8 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.hpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.hpp @@ -27,10 +27,10 @@ #include template -std::vector computeProbabilityDistribution(std::vector& data, const T referenceValue, size_t binCount, T binStart, +std::vector computeProbabilityDistribution(std::vector& data, const double referenceValue, size_t binCount, T binStart, T binEnd) { - std::vector bins(binCount); + std::vector bins(binCount); #pragma omp for schedule(static) for (size_t i = 0; i < data.size(); ++i) From 80250a67ab6413c86b89872d40f67cfabac35520 Mon Sep 17 00:00:00 2001 From: luke_sch02 Date: Thu, 2 May 2024 16:58:20 +0200 Subject: [PATCH 13/13] density pdf: adjust formatting --- main/src/analytical_solutions/turbulence/density_pdf.cpp | 2 +- main/src/analytical_solutions/turbulence/density_pdf.hpp | 7 +++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/main/src/analytical_solutions/turbulence/density_pdf.cpp b/main/src/analytical_solutions/turbulence/density_pdf.cpp index fe7e41fb1..2b5ad09eb 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.cpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.cpp @@ -56,7 +56,7 @@ int main(int argc, char** argv) return exitSuccess(); } - std::vector rho; + std::vector rho; std::vector bins; auto h5reader = fileReaderFactory(false, MPI_COMM_WORLD); diff --git a/main/src/analytical_solutions/turbulence/density_pdf.hpp b/main/src/analytical_solutions/turbulence/density_pdf.hpp index a3ee102b8..6149393e7 100644 --- a/main/src/analytical_solutions/turbulence/density_pdf.hpp +++ b/main/src/analytical_solutions/turbulence/density_pdf.hpp @@ -27,8 +27,8 @@ #include template -std::vector computeProbabilityDistribution(std::vector& data, const double referenceValue, size_t binCount, T binStart, - T binEnd) +std::vector computeProbabilityDistribution(std::vector& data, const double referenceValue, size_t binCount, + T binStart, T binEnd) { std::vector bins(binCount); @@ -41,8 +41,7 @@ std::vector computeProbabilityDistribution(std::vector& data, const d T binSize = (binEnd - binStart) / binCount; for (size_t bin = 0; bin < binCount; bin++) { - bins[bin] = std::count_if(data.begin(), data.end(), - [bin, binStart, binSize](T i) + bins[bin] = std::count_if(data.begin(), data.end(), [bin, binStart, binSize](T i) { return i > binSize * bin + binStart && i <= binSize * (bin + 1) + binStart; }); } return bins;