diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 8b1886363..72359ffb4 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -18,14 +18,17 @@ include(addGoogleTest) include(addGoogleBench) find_package(Threads) +find_package(CLI11) set(SAMURAI_BENCHMARKS benchmark_celllist_construction.cpp benchmark_search.cpp - benchmark_set.cpp + # benchmark_set.cpp # buggy for now. main.cpp ) + + # foreach(filename IN LISTS SAMURAI_BENCHMARKS) # string(REPLACE ".cpp" "" targetname ${filename}) # add_executable(${targetname} ${COMMON_BASE} ${filename} ${SAMURAI_HEADERS}) @@ -35,8 +38,16 @@ set(SAMURAI_BENCHMARKS # endforeach() add_executable(bench_samurai ${SAMURAI_BENCHMARKS}) + +add_executable(bench_advection_2d benchmark_advection_2d.cpp) +add_executable(bench_linear_convection benchmark_linear_convection.cpp) + + + # target_include_directories(bench_samurai PRIVATE ${SAMURAI_INCLUDE_DIR}) -target_link_libraries(bench_samurai samurai benchmark::benchmark) +target_link_libraries(bench_samurai samurai CLI11::CLI11 benchmark::benchmark) +target_link_libraries(bench_advection_2d samurai CLI11::CLI11 benchmark::benchmark) +target_link_libraries(bench_linear_convection samurai CLI11::CLI11 benchmark::benchmark) # target_include_directories(bench_samurai_lib PRIVATE ${SAMURAI_INCLUDE_DIR}) # # if(DOWNLOAD_GTEST OR GTEST_SRC_DIR) diff --git a/benchmark/benchmark_advection_2d.cpp b/benchmark/benchmark_advection_2d.cpp new file mode 100644 index 000000000..e4302a4d5 --- /dev/null +++ b/benchmark/benchmark_advection_2d.cpp @@ -0,0 +1,286 @@ +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +template +void init(Field& u) +{ + auto& mesh = u.mesh(); + u.resize(); + + samurai::for_each_cell(mesh, + [&](auto& cell) + { + auto center = cell.center(); + const double radius = .2; + const double x_center = 0.3; + // const double y_center = 0.3; + if (std::abs(center[0] - x_center) <= radius * radius) + { + u[cell] = 1; + } + else + { + u[cell] = 0; + } + }); +} + +template +void flux_correction(double dt, const std::array& a, const Field& u, Field& unp1) +{ + using mesh_t = typename Field::mesh_t; + using mesh_id_t = typename mesh_t::mesh_id_t; + using interval_t = typename mesh_t::interval_t; + constexpr std::size_t dim = Field::dim; + + auto mesh = u.mesh(); + + for (std::size_t level = mesh.min_level(); level < mesh.max_level(); ++level) + { + xt::xtensor_fixed> stencil; + + stencil = { + {-1, 0} + }; + + auto subset_right = samurai::intersection(samurai::translate(mesh[mesh_id_t::cells][level + 1], stencil), + mesh[mesh_id_t::cells][level]) + .on(level); + + subset_right( + [&](const auto& i, const auto& index) + { + auto j = index[0]; + const double dx = mesh.cell_length(level); + + unp1(level, i, j) = unp1(level, i, j) + + dt / dx + * (samurai::upwind_op(level, i, j).right_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i + 1, 2 * j).right_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i + 1, 2 * j + 1).right_flux(a, u)); + }); + + stencil = { + {1, 0} + }; + + auto subset_left = samurai::intersection(samurai::translate(mesh[mesh_id_t::cells][level + 1], stencil), + mesh[mesh_id_t::cells][level]) + .on(level); + + subset_left( + [&](const auto& i, const auto& index) + { + auto j = index[0]; + const double dx = mesh.cell_length(level); + + unp1(level, i, j) = unp1(level, i, j) + - dt / dx + * (samurai::upwind_op(level, i, j).left_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i, 2 * j).left_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i, 2 * j + 1).left_flux(a, u)); + }); + + stencil = { + {0, -1} + }; + + auto subset_up = samurai::intersection(samurai::translate(mesh[mesh_id_t::cells][level + 1], stencil), mesh[mesh_id_t::cells][level]) + .on(level); + + subset_up( + [&](const auto& i, const auto& index) + { + auto j = index[0]; + const double dx = mesh.cell_length(level); + + unp1(level, i, j) = unp1(level, i, j) + + dt / dx + * (samurai::upwind_op(level, i, j).up_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i, 2 * j + 1).up_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i + 1, 2 * j + 1).up_flux(a, u)); + }); + + stencil = { + {0, 1} + }; + + auto subset_down = samurai::intersection(samurai::translate(mesh[mesh_id_t::cells][level + 1], stencil), + mesh[mesh_id_t::cells][level]) + .on(level); + + subset_down( + [&](const auto& i, const auto& index) + { + auto j = index[0]; + const double dx = mesh.cell_length(level); + + unp1(level, i, j) = unp1(level, i, j) + - dt / dx + * (samurai::upwind_op(level, i, j).down_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i, 2 * j).down_flux(a, u) + - .5 * samurai::upwind_op(level + 1, 2 * i + 1, 2 * j).down_flux(a, u)); + }); + } +} + +template +void ADVECTION_2D(benchmark::State& state) +{ + for (auto _ : state) + { + state.PauseTiming(); + + double Tf = 0.1; + + constexpr std::size_t dim = 2; + using Config = samurai::MRConfig; + + // Simulation parameters + xt::xtensor_fixed> min_corner = {0., 0.}; + xt::xtensor_fixed> max_corner = {1., 1.}; + std::array a{ + {1, 0} + }; + double cfl = 0.5; + double t = 0.; + std::string restart_file; + + // Multiresolution parameters + double mr_epsilon = 2.e-4; // Threshold used by multiresolution + double mr_regularity = 1.; // Regularity guess for multiresolution + // bool correction = false; + + const samurai::Box box(min_corner, max_corner); + samurai::MRMesh mesh; + auto u = samurai::make_scalar_field("u", mesh); + + mesh = {box, min_level, max_level}; + init(u); + + samurai::make_bc>(u, 0.); + + double dt = cfl * mesh.cell_length(max_level); + + auto unp1 = samurai::make_scalar_field("unp1", mesh); + + auto MRadaptation = samurai::make_MRAdapt(u); + MRadaptation(mr_epsilon, mr_regularity); + +#ifdef SAMURAI_WITH_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + state.ResumeTiming(); + int ITER_STEPS = 10; + for (int i = 0; i < ITER_STEPS; i++) + { + MRadaptation(mr_epsilon, mr_regularity); + t += dt; + if (t > Tf) + { + dt += Tf - t; + t = Tf; + } + samurai::update_ghost_mr(u); + unp1.resize(); + unp1 = u - dt * samurai::upwind(a, u); + // if (correction) + // { + // flux_correction(dt, a, u, unp1); + // } + std::swap(u.array(), unp1.array()); + } +#ifdef SAMURAI_WITH_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + } +} + +// BENCHMARK(ADVECTION_2D); + +// MRA with min_level = 5 +BENCHMARK_TEMPLATE(ADVECTION_2D, 5, 8)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 5, 10)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 5, 12)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 5, 14)->Unit(benchmark::kMillisecond)->Iterations(1); + +// MRA with max_level - min-level = 2 +BENCHMARK_TEMPLATE(ADVECTION_2D, 6, 8)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 8, 10)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 10, 12)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 12, 14)->Unit(benchmark::kMillisecond)->Iterations(1); + +// Uniform +BENCHMARK_TEMPLATE(ADVECTION_2D, 6, 6)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 8, 8)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 10, 10)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 12, 12)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(ADVECTION_2D, 14, 14)->Unit(benchmark::kMillisecond)->Iterations(1); + +/** SOURCE : https://gist.github.com/mdavezac/eb16de7e8fc08e522ff0d420516094f5 + **/ + +// This reporter does nothing. +// We can use it to disable output from all but the root process +class NullReporter : public ::benchmark::BenchmarkReporter +{ + public: + + NullReporter() + { + } + + virtual bool ReportContext(const Context&) + { + return true; + } + + virtual void ReportRuns(const std::vector&) + { + } + + virtual void Finalize() + { + } +}; + +int main(int argc, char** argv) +{ +#ifdef SAMURAI_WITH_MPI + MPI_Init(&argc, &argv); + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + ::benchmark::Initialize(&argc, argv); + if (rank == 0) + { +#endif + ::benchmark::RunSpecifiedBenchmarks(); +#ifdef SAMURAI_WITH_MPI + } + else + { + NullReporter null; + ::benchmark::RunSpecifiedBenchmarks(&null); + } +#endif + +#ifdef SAMURAI_WITH_MPI + MPI_Finalize(); +#endif + return 0; +} diff --git a/benchmark/benchmark_linear_convection.cpp b/benchmark/benchmark_linear_convection.cpp new file mode 100644 index 000000000..1432dc9cb --- /dev/null +++ b/benchmark/benchmark_linear_convection.cpp @@ -0,0 +1,260 @@ +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +template +void LINEAR_CONVECTION(benchmark::State& state) +{ + for (auto _ : state) + { + state.PauseTiming(); + + // static constexpr std::size_t dim = 2; + using Config = samurai::MRConfig; + using Box = samurai::Box; + using point_t = typename Box::point_t; + + //--------------------// + // Program parameters // + //--------------------// + + // Simulation parameters + double left_box = -1; + double right_box = 1; + + // Time integration + double Tf = 3; + double dt = 0; + double cfl = 0.95; + double t = 0.; + std::string restart_file; + + // Multiresolution parameters + double mr_epsilon = 1e-4; // Threshold used by multiresolution + double mr_regularity = 1.; // Regularity guess for multiresolution + + //--------------------// + // Problem definition // + //--------------------// + + point_t box_corner1, box_corner2; + box_corner1.fill(left_box); + box_corner2.fill(right_box); + Box box(box_corner1, box_corner2); + std::array periodic; + periodic.fill(false); + samurai::MRMesh mesh; + auto u = samurai::make_scalar_field("u", mesh); + + samurai::make_bc>(u, 0.); + + mesh = {box, min_level, max_level, periodic}; + // Initial solution + u = samurai::make_scalar_field( + "u", + mesh, + [](const auto& coords) + { + if constexpr (dim == 1) + { + const auto& x = coords(0); + return (x >= -0.8 && x <= -0.3) ? 1. : 0.; + } + else if constexpr (dim == 2) + { + const auto& x = coords(0); + const auto& y = coords(1); + return (x >= -0.8 && x <= -0.3 && y >= 0.3 && y <= 0.8) ? 1. : 0.; + } + else + { + const auto& x = coords(0); + const auto& y = coords(1); + const auto& z = coords(2); + return (x >= -0.8 && x <= -0.3 && y >= 0.3 && y <= 0.8 && z >= -0.8 && z <= -0.3) ? 1. : 0.; + } + }); + + auto unp1 = samurai::make_scalar_field("unp1", mesh); + // Intermediary fields for the RK3 scheme + auto u1 = samurai::make_scalar_field("u1", mesh); + auto u2 = samurai::make_scalar_field("u2", mesh); + + unp1.fill(0); + u1.fill(0); + u2.fill(0); + + // Convection operator + samurai::VelocityVector velocity; + velocity.fill(1); + if constexpr (dim == 2) + { + velocity(1) = -1; + } + if constexpr (dim == 3) + { + velocity(1) = -1; + velocity(2) = 1; + } + + auto conv = samurai::make_convection_weno5(velocity); + + //--------------------// + // Time iteration // + //--------------------// + + if (dt == 0) + { + double dx = mesh.cell_length(max_level); + auto a = xt::abs(velocity); + double sum_velocities = xt::sum(xt::abs(velocity))(); + dt = cfl * dx / sum_velocities; + } + + auto MRadaptation = samurai::make_MRAdapt(u); + MRadaptation(mr_epsilon, mr_regularity); + +#ifdef SAMURAI_WITH_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + state.ResumeTiming(); + int ITER_STEPS = 10; + for (int i = 0; i < ITER_STEPS; i++) + { + // Move to next timestep + t += dt; + if (t > Tf) + { + dt += Tf - t; + t = Tf; + } + + // Mesh adaptation + MRadaptation(mr_epsilon, mr_regularity); + samurai::update_ghost_mr(u); + unp1.resize(); + u1.resize(); + u2.resize(); + u1.fill(0); + u2.fill(0); + + // unp1 = u - dt * conv(u); + + // TVD-RK3 (SSPRK3) + u1 = u - dt * conv(u); + samurai::update_ghost_mr(u1); + u2 = 3. / 4 * u + 1. / 4 * (u1 - dt * conv(u1)); + samurai::update_ghost_mr(u2); + unp1 = 1. / 3 * u + 2. / 3 * (u2 - dt * conv(u2)); + + // u <-- unp1 + std::swap(u.array(), unp1.array()); + } +#ifdef SAMURAI_WITH_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif + } +} + +//// 2D + +// MRA with min_level = 5 +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 5, 8, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 5, 10, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 5, 12, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +// BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 5, 14, 2)->Unit(benchmark::kMillisecond)->Iterations(1); + +// MRA with max_level - min-level = 2 +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 6, 8, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 8, 10, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 10, 12, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +// BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 12, 14, 2)->Unit(benchmark::kMillisecond)->Iterations(1); + +// Uniform +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 6, 6, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 8, 8, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 10, 10, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 12, 12, 2)->Unit(benchmark::kMillisecond)->Iterations(1); +// BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 14, 14, 2)->Unit(benchmark::kMillisecond)->Iterations(1); + +//// 3D + +// MRA with min_level = +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 3, 5, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 3, 6, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 3, 7, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 3, 8, 3)->Unit(benchmark::kMillisecond)->Iterations(1); + +// Uniform +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 4, 4, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 5, 5, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 6, 6, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 7, 7, 3)->Unit(benchmark::kMillisecond)->Iterations(1); +BENCHMARK_TEMPLATE(LINEAR_CONVECTION, 8, 8, 3)->Unit(benchmark::kMillisecond)->Iterations(1); + +/** SOURCE : https://gist.github.com/mdavezac/eb16de7e8fc08e522ff0d420516094f5 + **/ + +// This reporter does nothing. +// We can use it to disable output from all but the root process +class NullReporter : public ::benchmark::BenchmarkReporter +{ + public: + + NullReporter() + { + } + + virtual bool ReportContext(const Context&) + { + return true; + } + + virtual void ReportRuns(const std::vector&) + { + } + + virtual void Finalize() + { + } +}; + +int main(int argc, char** argv) +{ +#ifdef SAMURAI_WITH_MPI + MPI_Init(&argc, &argv); + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + ::benchmark::Initialize(&argc, argv); + if (rank == 0) + { +#endif + ::benchmark::RunSpecifiedBenchmarks(); +#ifdef SAMURAI_WITH_MPI + } + else + { + NullReporter null; + ::benchmark::RunSpecifiedBenchmarks(&null); + } +#endif + +#ifdef SAMURAI_WITH_MPI + MPI_Finalize(); +#endif + return 0; +} diff --git a/benchmark/benchmark_search.cpp b/benchmark/benchmark_search.cpp index 9bb3cbf67..ce5a7b881 100644 --- a/benchmark/benchmark_search.cpp +++ b/benchmark/benchmark_search.cpp @@ -73,7 +73,8 @@ class MyFixture : public ::benchmark::Fixture for (std::size_t s = 0; s < state.range(0); ++s) { auto level = std::experimental::randint(min_level, max_level); - std::array coord; + // std::array coord; + xt::xtensor_fixed> coord; for (auto& c : coord) { c = std::experimental::randint(-bound << level, (bound << level) - 1); diff --git a/demos/FiniteVolume/linear_convection.cpp b/demos/FiniteVolume/linear_convection.cpp index aca48e50a..78150110f 100644 --- a/demos/FiniteVolume/linear_convection.cpp +++ b/demos/FiniteVolume/linear_convection.cpp @@ -110,22 +110,30 @@ int main(int argc, char* argv[]) { mesh = {box, min_level, max_level, periodic}; // Initial solution - u = samurai::make_scalar_field("u", - mesh, - [](const auto& coords) - { - if constexpr (dim == 1) - { - const auto& x = coords(0); - return (x >= -0.8 && x <= -0.3) ? 1. : 0.; - } - else - { - const auto& x = coords(0); - const auto& y = coords(1); - return (x >= -0.8 && x <= -0.3 && y >= 0.3 && y <= 0.8) ? 1. : 0.; - } - }); + u = samurai::make_scalar_field( + "u", + mesh, + [](const auto& coords) + { + if constexpr (dim == 1) + { + const auto& x = coords(0); + return (x >= -0.8 && x <= -0.3) ? 1. : 0.; + } + else if constexpr (dim == 2) + { + const auto& x = coords(0); + const auto& y = coords(1); + return (x >= -0.8 && x <= -0.3 && y >= 0.3 && y <= 0.8) ? 1. : 0.; + } + else if constexpr (dim == 3) + { + const auto& x = coords(0); + const auto& y = coords(1); + const auto& z = coords(2); + return (x >= -0.8 && x <= -0.3 && y >= 0.3 && y <= 0.8 && z >= -0.8 && z <= -0.3) ? 1. : 0.; + } + }); } else { @@ -148,6 +156,11 @@ int main(int argc, char* argv[]) { velocity(1) = -1; } + if constexpr (dim == 3) + { + velocity(1) = -1; + velocity(2) = 1; + } auto conv = samurai::make_convection_weno5(velocity); //--------------------//