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.