From 53bd0590b020919ee95b3ca12b9a096e3f58ed19 Mon Sep 17 00:00:00 2001 From: Loic Gouarin Date: Sun, 13 Sep 2026 19:08:49 +0200 Subject: [PATCH] Measure the performance numbers of the article, and build its table The article this repository reproduces is first of all a performance paper, and its tables are made of three numbers this solver did not measure: the sparsity index of the adapted mesh, the cell updates per second, and the time to solution. Without them the repository reproduces pictures and compares nothing. euler/metrics.hpp measures them where they happen rather than leaving them to be reconstructed from a log, and the conventions are the article's: sparsity index cells of the adapted mesh over the cells of the uniform mesh at max-level, reported at the initial and at the final time. A Riemann problem fills its mesh up as the waves spread, and one number taken at one end would flatter or damn it. cell updates one cell advanced by one time step, counted once per cell per step whatever the integrator does inside. SSP-RK2 evaluates the flux twice and Strang sweeps 2*dim - 1 times; counting those would make the more expensive integrator look faster. time to solution the time loop. The mesh adaptation is part of it and is reported apart, which is the number the comparison of a multiresolution against a gradient criterion turns on. Writing files is not part of it: --nfiles is a choice of whoever runs the solver, and a table of times that moved with it would compare nothing. Everything is accumulated on the local subdomain and reduced once, at report time, which is a collective call: cells and updates summed over the ranks, times taken at their maximum. --metrics-file writes the same numbers as JSON, and python/performance.py makes the table out of one run per resolution. On the reference case of the article, configuration 3 of Lax & Liu to t_f = 0.8 at second order, on one core: l_min l_max resolution Mcu/s time (s) AMR cells ti/tf sparsity ti/tf 3 6 64^2 2.2 0.39 25.1% 1720 / 3784 42.0% / 92.4% 3 7 128^2 2.4 2.49 29.6% 3916 / 13522 23.9% / 82.5% 3 8 256^2 2.5 15.41 28.1% 8416 / 45394 12.8% / 69.3% 3 9 512^2 2.5 102.92 28.4% 17524 / 151480 6.7% / 57.8% 9 9 512^2 4.7 147.64 0.0% 262144 / 262144 100.0% / 100.0% The last row is the uniform mesh at the same resolution, the reference the article puts at the bottom of its own table: the adapted run does 2.7 times fewer cell updates and is 1.4 times faster, running at a little more than half the throughput of the uniform one. About a third of what is lost is the adaptation itself; the rest is what a cell of an adapted mesh costs. The sparsity column is where a comparison with the article goes wrong most easily. What the multiresolution keeps is a neighbourhood of the discontinuities, which are curves in a plane, so their share of the mesh falls with the resolution -- 92%, 83%, 69%, 58% over the four rows above. A sparsity index quoted without the max-level it was measured at compares nothing. The threshold is the other half: --mr-eps takes several values and measures each, and two decades of it take the final sparsity at level 9 from 58% to 16% and the time to solution from 115 s to 59 s, with the initial mesh unmoved, the details of a piecewise constant state being of order one at a discontinuity and far above every threshold in that range. Those rows say nothing about accuracy, which is what python/error_analysis.py measures. The metrics are held against what the run computed rather than against themselves: the sparsity index is 100% on a uniform mesh whose cell count is arithmetic on a unit domain, the final cell count is the cell count of the HDF5 file the run wrote, and the time to solution holds the adaptation and not the output. One slow test holds the fall of the sparsity index with the resolution, which is the reading every comparison with the article depends on. --- euler/metrics.hpp | 302 +++++++++++++++++++++++++++++++++++++++ main_1d.cpp | 29 +++- main_2d.cpp | 29 +++- main_3d.cpp | 29 +++- python/performance.py | 215 ++++++++++++++++++++++++++++ readme.md | 106 ++++++++++++++ tests/test_metrics.py | 76 ++++++++++ tests/test_validation.py | 43 +++++- tests/util.py | 15 ++ 9 files changed, 834 insertions(+), 10 deletions(-) create mode 100644 euler/metrics.hpp create mode 100644 python/performance.py create mode 100644 tests/test_metrics.py diff --git a/euler/metrics.hpp b/euler/metrics.hpp new file mode 100644 index 0000000..28b3343 --- /dev/null +++ b/euler/metrics.hpp @@ -0,0 +1,302 @@ +// Copyright 2025 the samurai team +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include + +#ifdef SAMURAI_WITH_MPI +#include +namespace mpi = boost::mpi; +#endif + +// ============================================================================= +// Performance metrics +// ----------------------------------------------------------------------------- +// The article this repository reproduces is first of all a performance paper, +// and it reports three numbers per run. They are cheap to measure and useless +// to guess, so the solver measures them itself rather than leaving them to be +// reconstructed from a log: +// +// sparsity index cells of the adapted mesh over the cells a uniform mesh +// at max-level would hold, as a percentage. 100% is a mesh +// that never coarsened. Reported at the initial and at the +// final time, as the article does: the mesh of a Riemann +// problem fills up as the waves spread, and one number +// taken at one end would flatter or damn it. +// +// Mcu/s millions of cell updates per second, a cell update being +// one cell advanced by one time step. Counted as one per +// cell per time step whatever the integrator does +// internally: SSP-RK2 evaluates the flux twice and Strang +// sweeps 2*dim - 1 times, and counting those would make a +// more expensive integrator look faster. It is the same +// convention as the article, whose scheme is the +// single-pass MUSCL-Hancock. +// +// time to solution the time loop alone. Initialization and output are +// measured too but kept out of it: nfiles is a choice of +// the person running, not a property of the solver. +// +// Plus the fraction of the time to solution spent adapting the mesh, which is +// what the comparison of a multiresolution against a gradient criterion is +// really about. +// +// All of it is accumulated on the local subdomain and reduced once, at print +// time: cells and updates are summed over the ranks, times are taken at their +// maximum, which is the rank that everyone else waits for. +// ============================================================================= + +class Metrics +{ + public: + + // The denominator of the sparsity index is fixed once, from the domain at + // max-level: it is the mesh the run would have used without adaptation. The + // levels and the dimension travel with it so that a metrics file says which + // run it describes without the command line that produced it. + template + explicit Metrics(const Mesh& mesh) + : m_uniform_cells(mesh.domain().nb_cells()) + , m_dim(Mesh::dim) + , m_min_level(mesh.min_level()) + , m_max_level(mesh.max_level()) + { + } + + // Starts the time to solution, and records the mesh the run starts from. + template + void start(const Mesh& mesh) + { + m_initial_cells = local_cells(mesh); + m_start = clock::now(); + } + + // One time step over the cells the mesh carries, whatever the integrator + // does inside: one cell advanced by one time step is one cell update. + template + void step(const Mesh& mesh) + { + m_cell_updates += local_cells(mesh); + ++m_steps; + } + + // The mesh adaptation, timed apart. It is part of the time to solution. + template + void adapt(Fn&& adaptation) + { + const auto begin = clock::now(); + adaptation(); + m_adapt_time += seconds_since(begin); + } + + // Writing a file, timed apart and subtracted from the time to solution. + template + void output(Fn&& write) + { + const auto begin = clock::now(); + write(); + m_output_time += seconds_since(begin); + } + + // Stops the time to solution and records the mesh the run ends on. + template + void stop(const Mesh& mesh) + { + m_run_time = seconds_since(m_start) - m_output_time; + m_final_cells = local_cells(mesh); + } + + // Everything above, reduced over the ranks. Call once, after stop(). + struct Summary + { + std::size_t dim = 0; + std::size_t min_level = 0; + std::size_t max_level = 0; + std::uint64_t uniform_cells = 0; + std::uint64_t initial_cells = 0; + std::uint64_t final_cells = 0; + std::uint64_t cell_updates = 0; + std::uint64_t steps = 0; + double run_time = 0.; + double adapt_time = 0.; + double output_time = 0.; + + double initial_sparsity() const + { + return sparsity(initial_cells); + } + + double final_sparsity() const + { + return sparsity(final_cells); + } + + // Millions of cell updates per second, averaged over the whole run. + double mcu_per_second() const + { + return run_time > 0. ? static_cast(cell_updates) / run_time * 1e-6 : 0.; + } + + double adapt_fraction() const + { + return run_time > 0. ? 100. * adapt_time / run_time : 0.; + } + + private: + + double sparsity(std::uint64_t cells) const + { + return uniform_cells > 0 ? 100. * static_cast(cells) / static_cast(uniform_cells) : 0.; + } + }; + + // Prints the block below and, when given a name, writes the same numbers as + // JSON. Collective: the reduction is inside, so every rank calls it and only + // the root one writes. + void report(const std::string& filename = "") const + { + const auto s = summary(); + if (!is_root()) + { + return; + } + print(s); + write(s, filename); + } + + // The reduced numbers, for a caller that wants them rather than a report. + // Collective as well. + Summary summary() const + { + Summary s{.dim = m_dim, + .min_level = m_min_level, + .max_level = m_max_level, + .uniform_cells = m_uniform_cells, + .initial_cells = m_initial_cells, + .final_cells = m_final_cells, + .cell_updates = m_cell_updates, + .steps = m_steps, + .run_time = m_run_time, + .adapt_time = m_adapt_time, + .output_time = m_output_time}; +#ifdef SAMURAI_WITH_MPI + mpi::communicator world; + // The uniform mesh is the whole domain and is already global; the rest + // is what this rank carried. + s.initial_cells = mpi::all_reduce(world, s.initial_cells, std::plus()); + s.final_cells = mpi::all_reduce(world, s.final_cells, std::plus()); + s.cell_updates = mpi::all_reduce(world, s.cell_updates, std::plus()); + s.run_time = mpi::all_reduce(world, s.run_time, mpi::maximum()); + s.adapt_time = mpi::all_reduce(world, s.adapt_time, mpi::maximum()); + s.output_time = mpi::all_reduce(world, s.output_time, mpi::maximum()); +#endif + return s; + } + + private: + + // The block printed at the end of a run, in the conventions above. + static void print(const Summary& s) + { + std::cout << "\nperformance" << std::endl; + std::cout << fmt::format(" cells {} -> {} of {} uniform", s.initial_cells, s.final_cells, s.uniform_cells) << std::endl; + std::cout << fmt::format(" sparsity index {:.2f}% -> {:.2f}%", s.initial_sparsity(), s.final_sparsity()) << std::endl; + std::cout << fmt::format(" cell updates {} over {} time steps", s.cell_updates, s.steps) << std::endl; + std::cout << fmt::format(" time to solution {:.2f} s, {:.1f}% of it adapting the mesh", s.run_time, s.adapt_fraction()) << std::endl; + std::cout << fmt::format(" throughput {:.2f} Mcu/s", s.mcu_per_second()) << std::endl; + } + + // The same numbers as JSON, for whoever builds a table out of several runs. + static void write(const Summary& s, const std::string& filename) + { + if (filename.empty()) + { + return; + } + + std::ofstream out(filename); + if (!out) + { + throw std::runtime_error("cannot write the metrics file " + filename); + } + + out << fmt::format("{{\n" + " \"dim\": {},\n" + " \"min_level\": {},\n" + " \"max_level\": {},\n" + " \"uniform_cells\": {},\n" + " \"initial_cells\": {},\n" + " \"final_cells\": {},\n" + " \"initial_sparsity\": {:.6e},\n" + " \"final_sparsity\": {:.6e},\n" + " \"cell_updates\": {},\n" + " \"steps\": {},\n" + " \"run_time\": {:.6e},\n" + " \"adapt_time\": {:.6e},\n" + " \"output_time\": {:.6e},\n" + " \"adapt_fraction\": {:.6e},\n" + " \"mcu_per_second\": {:.6e}\n" + "}}\n", + s.dim, + s.min_level, + s.max_level, + s.uniform_cells, + s.initial_cells, + s.final_cells, + s.initial_sparsity(), + s.final_sparsity(), + s.cell_updates, + s.steps, + s.run_time, + s.adapt_time, + s.output_time, + s.adapt_fraction(), + s.mcu_per_second()); + } + + using clock = std::chrono::steady_clock; + + // What this rank holds, the leaves only: ghosts are not cells of the mesh. + template + static std::uint64_t local_cells(const Mesh& mesh) + { + return mesh.nb_cells(Mesh::mesh_id_t::cells); + } + + static double seconds_since(const clock::time_point& begin) + { + return std::chrono::duration(clock::now() - begin).count(); + } + + static bool is_root() + { +#ifdef SAMURAI_WITH_MPI + return mpi::communicator().rank() == 0; +#else + return true; +#endif + } + + std::uint64_t m_uniform_cells = 0; + std::size_t m_dim = 0; + std::size_t m_min_level = 0; + std::size_t m_max_level = 0; + std::uint64_t m_initial_cells = 0; + std::uint64_t m_final_cells = 0; + std::uint64_t m_cell_updates = 0; + std::uint64_t m_steps = 0; + + clock::time_point m_start{}; + double m_run_time = 0.; + double m_adapt_time = 0.; + double m_output_time = 0.; +}; diff --git a/main_1d.cpp b/main_1d.cpp index 7357398..5eef940 100644 --- a/main_1d.cpp +++ b/main_1d.cpp @@ -14,6 +14,7 @@ #include "euler/config.hpp" #include "euler/eos.hpp" #include "euler/init/cases.hpp" +#include "euler/metrics.hpp" #include "euler/save.hpp" #include "euler/reconstruction.hpp" #include "euler/schemes.hpp" @@ -44,6 +45,7 @@ int main(int argc, char* argv[]) fs::path path = "results"; std::string filename; std::size_t nfiles = 1; + std::string metrics_file; auto available = test_case::TestCaseRegistry::instance().available_test_cases(); @@ -73,6 +75,7 @@ int main(int argc, char* argv[]) app.add_option("--path", path, "Output path")->capture_default_str()->group("Output"); app.add_option("--filename", filename, "File name prefix (defaults to _)")->group("Output"); app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output"); + app.add_option("--metrics-file", metrics_file, "Write the performance metrics of the run, as JSON, to this file")->group("Output"); // The cases that take a parameter of their own declare it here, before the // parse. All of them do, not only the selected one: which case runs is @@ -202,9 +205,18 @@ int main(int argc, char* argv[]) auto MRadaptation = samurai::make_MRAdapt(u); auto mra_config = samurai::mra_config().relative_detail(true); + // The three numbers the article reports per run, measured here rather than + // reconstructed afterwards from a log: see euler/metrics.hpp. + Metrics metrics(mesh); + metrics.start(mesh); + while (t != Tf) { - MRadaptation(mra_config); + metrics.adapt( + [&] + { + MRadaptation(mra_config); + }); double dt = cfl * dx / get_max_lambda(u, eos); t += dt; @@ -223,6 +235,8 @@ int main(int argc, char* argv[]) std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << std::endl; + metrics.step(mesh); + if (order == 1) { advance(u, unp1, unp2, first_order, first_order_sweeps, dt_for_flux, dt, integrator); @@ -235,11 +249,20 @@ int main(int argc, char* argv[]) if (t >= static_cast(nsave + 1) * dt_save || t == Tf) { const std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : ""; - save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); - samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + // Writing is not part of the time to solution: how many files a run + // produces is a choice of the person running it. + metrics.output( + [&] + { + save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); + samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + }); } } + metrics.stop(mesh); + metrics.report(metrics_file); + samurai::finalize(); return 0; } diff --git a/main_2d.cpp b/main_2d.cpp index 1b8e3bd..be8f9dc 100644 --- a/main_2d.cpp +++ b/main_2d.cpp @@ -17,6 +17,7 @@ #include "euler/config.hpp" #include "euler/eos.hpp" #include "euler/init/cases.hpp" +#include "euler/metrics.hpp" #include "euler/prediction.hpp" #include "euler/save.hpp" #include "euler/reconstruction.hpp" @@ -117,6 +118,7 @@ int main(int argc, char* argv[]) fs::path path = "results"; std::string filename; std::size_t nfiles = 1; + std::string metrics_file; auto available = test_case::TestCaseRegistry::instance().available_test_cases(); @@ -148,6 +150,7 @@ int main(int argc, char* argv[]) app.add_option("--path", path, "Output path")->capture_default_str()->group("Output"); app.add_option("--filename", filename, "File name prefix (defaults to _)")->group("Output"); app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output"); + app.add_option("--metrics-file", metrics_file, "Write the performance metrics of the run, as JSON, to this file")->group("Output"); // The cases that take a parameter of their own declare it here, before the // parse. All of them do, not only the selected one: which case runs is @@ -290,13 +293,22 @@ int main(int argc, char* argv[]) return make_second_order_scheme(scheme, eos, options); }); + // The three numbers the article reports per run, measured here rather than + // reconstructed afterwards from a log: see euler/metrics.hpp. + Metrics metrics(mesh); + samurai::times::timers.start("TimeLoop"); + metrics.start(mesh); bool done = false; while (!done) { double dt = cfl * dx / get_max_lambda(u, eos); - MRadaptation(mra_config); + metrics.adapt( + [&] + { + MRadaptation(mra_config); + }); if (check_positivity) { @@ -316,6 +328,8 @@ int main(int argc, char* argv[]) } std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << "\r"; + metrics.step(mesh); + if (order == 1) { advance(u, unp1, unp2, first_order, first_order_sweeps, dt_for_flux, dt, integrator); @@ -330,12 +344,21 @@ int main(int argc, char* argv[]) if (t >= static_cast(nsave + 1) * dt_save || t == Tf) { const std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : ""; - save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); - samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + // Writing is not part of the time to solution: how many files a run + // produces is a choice of the person running it. + metrics.output( + [&] + { + save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); + samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + }); } } + metrics.stop(mesh); samurai::times::timers.stop("TimeLoop"); + metrics.report(metrics_file); + samurai::finalize(); return 0; } diff --git a/main_3d.cpp b/main_3d.cpp index b33890f..28a4c49 100644 --- a/main_3d.cpp +++ b/main_3d.cpp @@ -19,6 +19,7 @@ #include "euler/eos.hpp" #include "euler/init/cases.hpp" #include "euler/limiter.hpp" +#include "euler/metrics.hpp" #include "euler/prediction.hpp" #include "euler/save.hpp" #include "euler/reconstruction.hpp" @@ -117,6 +118,7 @@ int main(int argc, char* argv[]) fs::path path = "results"; std::string filename; std::size_t nfiles = 1; + std::string metrics_file; auto available = test_case::TestCaseRegistry::instance().available_test_cases(); @@ -148,6 +150,7 @@ int main(int argc, char* argv[]) app.add_option("--path", path, "Output path")->capture_default_str()->group("Output"); app.add_option("--filename", filename, "File name prefix (defaults to _)")->group("Output"); app.add_option("--nfiles", nfiles, "Number of output files")->capture_default_str()->group("Output"); + app.add_option("--metrics-file", metrics_file, "Write the performance metrics of the run, as JSON, to this file")->group("Output"); // The cases that take a parameter of their own declare it here, before the // parse. All of them do, not only the selected one: which case runs is @@ -291,14 +294,23 @@ int main(int argc, char* argv[]) return make_second_order_scheme(scheme, eos, options); }); + // The three numbers the article reports per run, measured here rather than + // reconstructed afterwards from a log: see euler/metrics.hpp. + Metrics metrics(mesh); + samurai::times::timers.start("TimeLoop"); + metrics.start(mesh); std::size_t limited_cells = 0; bool done = false; while (!done) { double dt = cfl * dx / get_max_lambda(u, eos); - MRadaptation(mra_config); + metrics.adapt( + [&] + { + MRadaptation(mra_config); + }); if (check_positivity) { @@ -318,6 +330,8 @@ int main(int argc, char* argv[]) } std::cout << fmt::format("iteration {}: t = {}, dt = {}", nt++, t, dt) << "\r"; + metrics.step(mesh); + if (order == 1) { advance(u, unp1, unp2, first_order, first_order_sweeps, dt_for_flux, dt, integrator); @@ -338,15 +352,24 @@ int main(int argc, char* argv[]) if (t >= static_cast(nsave + 1) * dt_save || t == Tf) { const std::string suffix = (nfiles != 1) ? fmt::format("_ite_{}", nsave++) : ""; - save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); - samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + // Writing is not part of the time to solution: how many files a run + // produces is a choice of the person running it. + metrics.output( + [&] + { + save(path.string(), fmt::format("{}{}", filename, suffix), u, eos); + samurai::dump(path, fmt::format("{}_restart{}", filename, suffix), mesh, u); + }); } } + metrics.stop(mesh); samurai::times::timers.stop("TimeLoop"); std::cout << std::endl << fmt::format("positivity floor applied to {} cell updates over {} iterations", limited_cells, nt) << std::endl; + metrics.report(metrics_file); + samurai::finalize(); return 0; } diff --git a/python/performance.py b/python/performance.py new file mode 100644 index 0000000..94d22dd --- /dev/null +++ b/python/performance.py @@ -0,0 +1,215 @@ +# Copyright 2025 the samurai team +# SPDX-License-Identifier: BSD-3-Clause +""" +Build the performance table of a test case: sparsity index, throughput and time +to solution, one row per resolution. + + python performance.py --levels 6 7 8 9 + +Runs the solver once per max-level, reads back the metrics each run writes +(`--metrics-file`, see euler/metrics.hpp) and prints them in the column layout +of the article this repository reproduces, so that the two tables can be read +side by side. + +Reading the table +----------------- + l_min/l_max the levels the multiresolution was allowed to span + resolution the uniform mesh at l_max, i.e. the equivalent resolution + Mcu/s millions of cell updates per second, one update being one cell + advanced by one time step, averaged over the whole run + time time to solution in seconds: the time loop, output excluded + AMR the share of that time spent adapting the mesh + cells leaves at the initial and at the final time + sparsity those two counts over the uniform mesh, in percent + +What the rows are, and what they are not +---------------------------------------- +The article varies the number of cells per octree leaf at a FIXED equivalent +resolution; samurai carries one cell per leaf, so that axis does not exist here +and the table sweeps the resolution instead. The row to compare with a row of +the article is the one with the same equivalent resolution, and nothing else is +comparable: the sparsity index of a two-dimensional Riemann problem falls +roughly as 2^-l_max, the discontinuities being curves in a plane, so a table +read across rows says as much about the resolutions chosen as about the two +codes. + +The threshold of the multiresolution is the knob that trades cells for +accuracy, the way the block size is in the article, and --mr-eps sweeps it: pass +several values and each one is measured at each level. Sparsity alone is not a +figure of merit, though — a large enough threshold makes any mesh sparse and the +solution wrong — so a row is only worth reading next to the error it carries, +which is what python/error_analysis.py measures on the cases that have an exact +solution. + +--uniform adds the row the article puts last: the same run on a uniform mesh at +l_max, which is the reference both for the throughput (no adaptation to pay +for, no level interface to cross) and for the time to solution (the speed-up +adaptation buys). It costs as much as the whole sweep above it, hence opt-in. + +Examples +-------- +# the case the article uses for its performance table, over four resolutions +python performance.py --levels 6 7 8 9 + +# with the uniform reference run at the finest level +python performance.py --levels 6 7 8 9 --uniform + +# what the multiresolution threshold buys, at one resolution +python performance.py --levels 8 --mr-eps 1e-4 1e-3 1e-2 + +# the same in three dimensions, where only configuration 3 exists +python performance.py --dim 3 --levels 4 5 6 --min-level 2 + +# first order, to see what the reconstruction costs and what it saves +python performance.py --levels 6 7 8 --order 1 +""" + +import argparse +import json +import os +import subprocess +import sys +import tempfile + +# The reference case of the article: Lax & Liu configuration 3, to t_f = 0.8. +DEFAULT_CASE = "lax_liu" +DEFAULT_TF = 0.8 + + +def run(exe, workdir, level, min_level, eps, args): + """Run one simulation and return the metrics it wrote, as a dict.""" + workdir = os.path.abspath(workdir) + os.makedirs(workdir, exist_ok=True) + metrics_file = os.path.join(workdir, "metrics.json") + + cmd = [ + exe, + "--test-case", args.test_case, + "--scheme", args.scheme, + "--order", str(args.order), + "--Tf", str(args.Tf), + "--cfl", str(args.cfl), + "--min-level", str(min_level), + "--max-level", str(level), + "--nfiles", "1", + "--metrics-file", metrics_file, + ] + + if args.test_case == "lax_liu": + cmd += ["--riemann-config", str(args.riemann_config)] + if eps is not None: + cmd += ["--mr-eps", str(eps)] + cmd += args.extra + + print(f"# {' '.join(cmd)}", flush=True) + proc = subprocess.run(cmd, cwd=workdir, capture_output=True, text=True) + if proc.returncode != 0: + sys.stderr.write(proc.stdout[-2000:] + proc.stderr[-2000:]) + raise RuntimeError(f"{os.path.basename(exe)} failed at level {level}") + + with open(metrics_file) as handle: + metrics = json.load(handle) + # The levels and the dimension come back from the run itself; the threshold + # is ours, the solver having no reason to know it was swept. + metrics["mr_eps"] = eps + return metrics + + +def resolution(metrics): + """The equivalent uniform resolution, as the article writes it.""" + return f"{2 ** metrics['max_level']}^{metrics['dim']}" + + +def cells(n): + """Cell counts the way the article prints them: millions past a million.""" + return f"{n / 1e6:.2f} M" if n >= 1e6 else f"{n}" + + +def print_table(rows, header): + print() + print(header) + columns = ( + f"{'l_min':>6} {'l_max':>6} {'resolution':>12} {'mr-eps':>9} {'Mcu/s':>8}" + f" {'time (s)':>9} {'AMR':>6} {'cells ti/tf':>21} {'sparsity ti/tf':>17}" + ) + print(columns) + print("-" * len(columns)) + for m in rows: + counts = f"{cells(m['initial_cells'])} / {cells(m['final_cells'])}" + sparsity = f"{m['initial_sparsity']:.1f}% / {m['final_sparsity']:.1f}%" + eps = f"{m['mr_eps']:.0e}" if m["mr_eps"] is not None else "default" + print( + f"{m['min_level']:>6} {m['max_level']:>6} {resolution(m):>12} {eps:>9}" + f" {m['mcu_per_second']:>8.1f} {m['run_time']:>9.2f} {m['adapt_fraction']:>5.1f}%" + f" {counts:>21} {sparsity:>17}" + ) + + +def speedup(adapted, uniform): + """What the adaptation bought at equal resolution, printed under the table.""" + print() + print( + f"# uniform reference at l_max = {uniform['max_level']}: " + f"{uniform['run_time'] / adapted['run_time']:.2f}x the time to solution of the adapted run, " + f"{uniform['mcu_per_second'] / adapted['mcu_per_second']:.2f}x its throughput, " + f"{uniform['cell_updates'] / adapted['cell_updates']:.2f}x its cell updates" + ) + + +def main(argv=None): + parser = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) + parser.add_argument("--dim", type=int, default=2, choices=(1, 2, 3), help="which solver to drive") + parser.add_argument("--levels", type=int, nargs="+", default=[6, 7, 8], help="max-levels to measure") + parser.add_argument("--min-level", type=int, default=3, help="min-level, the same for every row") + parser.add_argument("--test-case", default=DEFAULT_CASE) + parser.add_argument("--riemann-config", type=int, default=3, help="Lax & Liu configuration, for the lax_liu case") + parser.add_argument("--Tf", type=float, default=DEFAULT_TF) + parser.add_argument("--cfl", type=float, default=0.4) + parser.add_argument("--scheme", default="hllc") + parser.add_argument("--order", type=int, default=2, choices=(1, 2), help="the article runs second order") + parser.add_argument("--mr-eps", type=float, nargs="+", default=[None], + help="multiresolution thresholds to measure; every one of them is run at every level") + parser.add_argument("--uniform", action="store_true", help="add the uniform reference run at the finest level") + parser.add_argument("--exe", default=None, help="path to the solver (default: build/euler_d)") + parser.add_argument("--workdir", default=None, help="where to run (default: a temporary directory)") + parser.add_argument("--json", default=None, help="also write every row to this file") + parser.add_argument("--extra", nargs=argparse.REMAINDER, default=[], help="passed on to the solver; must come last") + args = parser.parse_args(argv) + + exe = args.exe or os.path.join("build", f"euler_{args.dim}d") + exe = os.path.abspath(exe) + if not os.path.exists(exe): + parser.error(f"{exe} not found; build the project or pass --exe") + + workdir = args.workdir or tempfile.mkdtemp(prefix="performance_") + rows = [ + run(exe, os.path.join(workdir, f"level{level}_eps{eps}"), level, args.min_level, eps, args) + for level in args.levels + for eps in args.mr_eps + ] + + case = args.test_case + (f" #{args.riemann_config}" if args.test_case == "lax_liu" else "") + header = f"# {case}, order {args.order}, {args.scheme}, Tf = {args.Tf}, euler_{args.dim}d" + print_table(rows, header) + + if args.uniform: + # The threshold plays no part here: with one level there is nothing to + # coarsen. One reference run per sweep, at the finest resolution of it. + finest = max(args.levels) + reference = run(exe, os.path.join(workdir, "uniform"), finest, finest, None, args) + print_table([reference], "# uniform mesh, same resolution") + speedup([row for row in rows if row["max_level"] == finest][0], reference) + rows.append(reference) + + if args.json: + with open(args.json, "w") as handle: + json.dump(rows, handle, indent=2) + print(f"\n# rows written to {args.json}") + + return rows + + +if __name__ == "__main__": + main() diff --git a/readme.md b/readme.md index e3a2a56..5b00cc5 100644 --- a/readme.md +++ b/readme.md @@ -159,6 +159,7 @@ and costs a fraction of an order that says nothing about the scheme. | `--path` | Output directory path | `results` | | `--filename` | Output file name prefix | `_` | | `--nfiles` | Number of output files to generate | `1` | +| `--metrics-file` | Write the performance metrics of the run, as JSON, to this file | (none) | ### Example Usage @@ -194,6 +195,111 @@ centre: ./euler_2d --test-case lax_liu --riemann-config 5 --riemann-interface 0.5 --Tf 0.3 ``` +## Performance + +The article this repository reproduces is first of all a performance paper, and +its tables are made of three numbers. Every run reports them: + +``` +performance + cells 17524 -> 151480 of 262144 uniform + sparsity index 6.68% -> 57.79% + cell updates 262193561 over 2672 time steps + time to solution 102.92 s, 28.4% of it adapting the mesh + throughput 2.55 Mcu/s +``` + +| Metric | What it is | +| :----- | :--------- | +| sparsity index | cells of the adapted mesh over the cells of the uniform mesh at `--max-level`, in percent. 100% is a mesh that never coarsened. Given at the initial and at the final time, as the article gives it: a Riemann problem fills its mesh up as the waves spread, and one number taken at one end would flatter or damn it. | +| cell updates, Mcu/s | one cell advanced by one time step is one cell update; the throughput is millions of those per second over the whole run. Counted once per cell per step whatever the integrator does inside, so that the `2*dim - 1` sweeps of Strang do not read as more work done. | +| time to solution | the time loop. The mesh adaptation is part of it and is reported apart; writing files is not, `--nfiles` being a choice of whoever runs the solver. | + +`--metrics-file ` writes the same numbers as JSON, which is what +`python/performance.py` builds a table out of: + +```bash +python python/performance.py --levels 6 7 8 9 --uniform +``` + +One run per resolution, on the reference case of the article — configuration 3 +of Lax & Liu, to `t_f = 0.8`, second order — gives, on one core: + +``` + l_min l_max resolution mr-eps Mcu/s time (s) AMR cells ti/tf sparsity ti/tf + 3 6 64^2 default 2.2 0.39 25.1% 1720 / 3784 42.0% / 92.4% + 3 7 128^2 default 2.4 2.49 29.6% 3916 / 13522 23.9% / 82.5% + 3 8 256^2 default 2.5 15.41 28.1% 8416 / 45394 12.8% / 69.3% + 3 9 512^2 default 2.5 102.92 28.4% 17524 / 151480 6.7% / 57.8% + 9 9 512^2 default 4.7 147.64 0.0% 262144 / 262144 100.0% / 100.0% +``` + +The cell counts are reproducible to the cell; the times are wall clock on one +core and move by ten percent or so between runs, which is worth remembering +before reading anything into a small difference. + +The last row is the uniform mesh at the same resolution, which is the reference +the article puts at the bottom of its own table. Reading the two bottom rows +together is the whole point of the exercise: **the adapted run does 2.7 times +fewer cell updates and is 1.4 times faster**, because it runs at a little more +than half the throughput of the uniform one. Roughly a third of what is lost is +the adaptation itself, at 28% of the time to solution; the rest is what an +adapted mesh costs per cell — level interfaces, prediction, intervals that are +shorter than a uniform row. + +Comparing that with the table of the article takes some care, and the sparsity +column is where it goes wrong most easily: + +- **Equivalent resolution is the only fair pairing.** The article varies the + number of cells per octree leaf at a fixed equivalent resolution of 4096²; + samurai carries one cell per leaf, so that axis does not exist here and the + table above sweeps the resolution instead. A sparsity index quoted without the + max-level it was measured at compares nothing: what the adaptation keeps is a + neighbourhood of the discontinuities, which are curves in a plane, so their + share of the mesh falls as the resolution rises — 92%, 83%, 69%, 58% over the + four rows above. + +- **The refinement criterion is not the same one**, and this is the real + difference between the two codes rather than a defect of either. The article + refines on a Löhner criterion, a normalised second difference thresholded at + `r_refine = 0.4`; samurai refines on the details of the multiresolution + thresholded at `--mr-eps`. The multiresolution comes with an error estimate + that the gradient criterion has not, and it keeps more cells for it. The + threshold is the knob that trades the two against each other, and `--mr-eps` + sweeps it: + +```bash +python python/performance.py --levels 9 --mr-eps 1e-4 1e-3 1e-2 +``` + +``` + l_min l_max resolution mr-eps Mcu/s time (s) AMR cells ti/tf sparsity ti/tf + 3 9 512^2 1e-04 2.3 115.04 28.5% 17524 / 151480 6.7% / 57.8% + 3 9 512^2 1e-03 1.8 69.19 38.9% 17524 / 72616 6.7% / 27.7% + 3 9 512^2 1e-02 1.5 58.86 43.6% 17524 / 42139 6.7% / 16.1% +``` + + Two decades of threshold take the final sparsity from 58% to 16%, which is the + order of magnitude the article reports, and the time to solution from 115 s to + 59 s. The initial mesh does not move at all: the details of a piecewise + constant state are of order one at the discontinuities and far above every + threshold in this range, so the three runs start from the same cells and part + company as the solution develops structure. **These rows say nothing about + accuracy**, and a large enough threshold makes any mesh sparse and any + solution wrong; the error of an adapted run against a uniform one is what + `python/error_analysis.py` measures, on the cases that have an exact solution. + +- **The AMR share is not measured over the same cadence.** The article runs its + AMR cycle once every 10 time steps and this solver adapts at every one, which + is most of the distance between the 28% above and the few percent it reports + on CPU at its nominal block size. + +- **Throughput is architecture, not method.** The numbers of the article are + measured on 72 ARM cores or on a Hopper GPU, against one core here, and its + solver works on blocks of 16² cells where this one works on intervals. The + column worth comparing is the sparsity index; the Mcu/s column is worth + comparing against *itself*, between two runs of this solver. + ## Tests The suite drives the built binaries as subprocesses, so it checks what a user diff --git a/tests/test_metrics.py b/tests/test_metrics.py new file mode 100644 index 0000000..407a7bd --- /dev/null +++ b/tests/test_metrics.py @@ -0,0 +1,76 @@ +# Copyright 2025 the samurai team +# SPDX-License-Identifier: BSD-3-Clause +"""The performance numbers a run reports, checked against what it computed. + +A metric nobody checks drifts from what it claims to measure, and a wrong one is +worse than none: the whole point of the sparsity index and of the throughput is +to be compared with the numbers of another code. So they are held here against +quantities measured independently — the cells of the file the run wrote, and the +cell count of the uniform mesh, which on a unit domain is arithmetic. + +The timings themselves are not asserted on. What can be asserted is that they +add up: the mesh adaptation is part of the time to solution and the output is +not. +""" + +import pytest + +from util import read, run_case_with_metrics + +# Resolutions chosen so that each run takes a second or two. +DIMENSIONS = [("euler_1d", 1, 8), ("euler_2d", 2, 6), ("euler_3d", 3, 4)] + + +@pytest.mark.parametrize("binary,dim,level", DIMENSIONS) +def test_a_uniform_mesh_is_the_reference_of_the_sparsity_index(binary, dim, level, tmp_path): + """On a uniform mesh the sparsity index is 100% by definition. + + It is the denominator of every other run, so it is worth pinning: the domain + of sod_x is the unit cube, and the mesh at max-level holds 2^(dim * level) + cells whatever the adaptation does or does not do. + """ + metrics, _ = run_case_with_metrics(binary, tmp_path, "sod_x", min_level=level, max_level=level, Tf=0.05) + + assert metrics["uniform_cells"] == 2 ** (dim * level) + assert metrics["initial_cells"] == metrics["uniform_cells"] + assert metrics["final_cells"] == metrics["uniform_cells"] + assert metrics["initial_sparsity"] == pytest.approx(100.0) + assert metrics["final_sparsity"] == pytest.approx(100.0) + + # Nothing refines or coarsens, so every step updates the same cells. + assert metrics["steps"] > 0 + assert metrics["cell_updates"] == metrics["steps"] * metrics["uniform_cells"] + + +def test_the_sparsity_index_counts_the_cells_of_the_output(tmp_path): + """The mesh the run reports at the final time is the mesh it wrote out. + + The count comes from the solver, the comparison from the HDF5 file: a + sparsity index computed on ghosts, or on the mesh of one level, would pass + every arithmetic check and fail this one. + """ + metrics, output = run_case_with_metrics( + "euler_2d", tmp_path, "lax_liu", riemann_config=3, min_level=3, max_level=6, Tf=0.05 + ) + _, volume, _ = read(output) + + assert metrics["final_cells"] == volume.size + assert volume.size < metrics["uniform_cells"], "the mesh never coarsened, this exercises nothing" + assert metrics["final_sparsity"] == pytest.approx(100.0 * volume.size / metrics["uniform_cells"]) + assert metrics["mcu_per_second"] == pytest.approx(1e-6 * metrics["cell_updates"] / metrics["run_time"]) + + +def test_the_time_to_solution_holds_the_adaptation_and_not_the_output(tmp_path): + """Adapting is part of solving; writing files is not. + + The second half is what makes two runs comparable: --nfiles is a choice of + whoever runs the solver, and a table of times to solution that moved with it + would compare nothing. + """ + metrics, _ = run_case_with_metrics( + "euler_2d", tmp_path, "lax_liu", riemann_config=3, min_level=3, max_level=6, Tf=0.05, nfiles=5 + ) + + assert 0.0 < metrics["adapt_time"] < metrics["run_time"] + assert metrics["adapt_fraction"] == pytest.approx(100.0 * metrics["adapt_time"] / metrics["run_time"]) + assert metrics["output_time"] > 0.0, "five files were asked for and none was timed" diff --git a/tests/test_validation.py b/tests/test_validation.py index 462c241..d77c4df 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -17,7 +17,7 @@ import numpy as np import pytest -from util import ROOT, level_count, read, run_case, sedov_blast_energy +from util import ROOT, level_count, read, run_case, run_case_with_metrics, sedov_blast_energy sys.path.insert(0, str(ROOT / "python")) from error_analysis import errors, exact_vortex # noqa: E402 @@ -414,3 +414,44 @@ def test_limited_reconstruction_stays_close_to_second_order(tmp_path): orders = [np.log2(a / b) for a, b in zip(l1, l1[1:])] assert orders[-1] > 1.8, f"L1 errors {l1}, orders {orders}" + + +# --------------------------------------------------------------------------- +# What the adaptation saves, in the units the article reports +# --------------------------------------------------------------------------- +SPARSITY_LEVELS = [6, 7, 8] + + +def test_the_sparsity_index_falls_with_the_resolution(tmp_path): + """The finer the mesh, the smaller the share of it the solution needs. + + This is the reading the performance table of the article turns on, and the + one that makes a single sparsity index meaningless on its own: what the + multiresolution keeps is a neighbourhood of the discontinuities, which are + curves in a plane, so their share of a square mesh falls as the resolution + rises. A sparsity index quoted without the max-level it was measured at + compares nothing. + + Run on the reference case of the article, configuration 3 of Lax & Liu, to + its final time. + """ + sparsity = [] + for level in SPARSITY_LEVELS: + metrics, _ = run_case_with_metrics( + "euler_2d", + tmp_path / f"level{level}", + "lax_liu", + riemann_config=3, + min_level=3, + max_level=level, + Tf=0.8, + order=2, + ) + assert metrics["uniform_cells"] == 4 ** level + sparsity.append(metrics["final_sparsity"]) + + assert all(b < a for a, b in zip(sparsity, sparsity[1:])), f"sparsity indices {sparsity}" + # Measured 92%, 83%, 69% at levels 6, 7 and 8: the threshold below is loose + # enough that a better or a worse adaptation still passes, and only a + # sparsity that has stopped following the resolution fails. + assert sparsity[-1] < 0.9 * sparsity[0], f"sparsity indices {sparsity}" diff --git a/tests/util.py b/tests/util.py index 44fbbe3..a036f37 100644 --- a/tests/util.py +++ b/tests/util.py @@ -9,6 +9,7 @@ other. """ +import json import os import re import subprocess @@ -74,6 +75,20 @@ def run_case(binary, workdir, case, scheme="hllc", label=None, **options): return out, stem +METRICS = "metrics.json" + + +def run_case_with_metrics(binary, workdir, case, **options): + """Run one case with --metrics-file, and read the metrics back. + + Returns (metrics, output file stem): what the run reports about itself, and + what it computed, which is what makes the two comparable. + """ + out, stem = run_case(binary, workdir, case, metrics_file=METRICS, **options) + with open(Path(workdir) / METRICS) as handle: + return json.load(handle), out / stem + + def read(h5file): """Return cell centers, cell volumes and the primitive fields of one output.