diff --git a/integration/ROS2S/Make.package b/integration/ROS2S/Make.package
new file mode 100644
index 0000000000..51fdcba9cd
--- /dev/null
+++ b/integration/ROS2S/Make.package
@@ -0,0 +1,7 @@
+CEXE_headers += ros2s_type.H
+CEXE_headers += ros2s_integrator.H
+ifeq ($(USE_ALL_SDC), TRUE)
+  CEXE_headers += actual_integrator_sdc.H
+else
+  CEXE_headers += actual_integrator.H
+endif
diff --git a/integration/ROS2S/_parameters b/integration/ROS2S/_parameters
new file mode 100644
index 0000000000..d9df057822
--- /dev/null
+++ b/integration/ROS2S/_parameters
@@ -0,0 +1,21 @@
+@namespace: integrator
+
+# for the step rejection logic on mass fractions, we only consider
+# species that are > X_reject_buffer * atol_spec
+X_reject_buffer              real         1.0
+
+# Rosenbrock tableau selector:
+#   0 = existing 3-stage ROS2S tableau
+#   1 = Verwer et al. (1999) 2-stage ROS2 tableau
+#   2 = Sandu (2000) positivity-favoring method A
+#   3 = Sandu (2000) positivity-favoring method B
+#   4 = Sandu (2000) positivity-favoring method C
+#   5 = Sandu (2000) positivity-favoring method D
+rosenbrock_tableau           int          0
+
+# H211b timestep controller parameters
+h211b_b                      real         4.0
+h211b_k                      real         2.5
+h211b_fac_min                real         0.2
+h211b_fac_max                real         6.0
+h211b_reduction_fac          real         0.5
diff --git a/integration/ROS2S/actual_integrator.H b/integration/ROS2S/actual_integrator.H
new file mode 100644
index 0000000000..561e91d88e
--- /dev/null
+++ b/integration/ROS2S/actual_integrator.H
@@ -0,0 +1,28 @@
+#ifndef actual_integrator_H
+#define actual_integrator_H
+
+#include <AMReX_REAL.H>
+
+#include <network.H>
+#include <burn_type.H>
+
+#include <integrator_data.H>
+#include <integrator_setup_strang.H>
+#include <ros2s_integrator.H>
+#include <ros2s_type.H>
+
+template <typename BurnT>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false)
+{
+    constexpr int int_neqs = integrator_neqs<BurnT>();
+
+    auto ros2s_state = integrator_setup<BurnT, ros2s_t<int_neqs>>(state, dt, is_retry);
+    auto state_save = integrator_backup(state);
+
+    int istate = ros2s_integrator(state, ros2s_state);
+
+    integrator_cleanup(ros2s_state, state, istate, state_save, dt);
+}
+
+#endif
diff --git a/integration/ROS2S/actual_integrator_sdc.H b/integration/ROS2S/actual_integrator_sdc.H
new file mode 100644
index 0000000000..5e187acd96
--- /dev/null
+++ b/integration/ROS2S/actual_integrator_sdc.H
@@ -0,0 +1,28 @@
+#ifndef actual_integrator_H
+#define actual_integrator_H
+
+#include <AMReX_REAL.H>
+
+#include <network.H>
+#include <burn_type.H>
+
+#include <integrator_data.H>
+#include <integrator_setup_sdc.H>
+#include <ros2s_integrator.H>
+#include <ros2s_type.H>
+
+template <typename BurnT>
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void actual_integrator (BurnT& state, amrex::Real dt, bool is_retry=false)
+{
+    constexpr int int_neqs = integrator_neqs<BurnT>();
+
+    auto ros2s_state = integrator_setup<BurnT, ros2s_t<int_neqs>>(state, dt, is_retry);
+    auto state_save = integrator_backup(state);
+
+    int istate = ros2s_integrator(state, ros2s_state);
+
+    integrator_cleanup(ros2s_state, state, istate, state_save, dt);
+}
+
+#endif
diff --git a/integration/ROS2S/ros2s_integrator.H b/integration/ROS2S/ros2s_integrator.H
new file mode 100644
index 0000000000..e75b931656
--- /dev/null
+++ b/integration/ROS2S/ros2s_integrator.H
@@ -0,0 +1,805 @@
+#ifndef ROS2S_INTEGRATOR_H
+#define ROS2S_INTEGRATOR_H
+
+#include <AMReX_Algorithm.H>
+
+#include <cmath>
+#include <limits>
+
+#include <network.H>
+#include <burn_type.H>
+#include <linpack.H>
+
+#ifdef STRANG
+#include <integrator_rhs_strang.H>
+#endif
+#ifdef SDC
+#include <integrator_rhs_sdc.H>
+#endif
+
+#include <integrator_data.H>
+#include <ros2s_type.H>
+
+namespace ros2s {
+
+struct ros2s_tableau {
+    static constexpr int stages = 3;
+    static constexpr amrex::Real gamma = 0.292893218813452_rt;
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        return i == 2 ? 0.585786437626905_rt : 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return 2.0000000000000036_rt;
+        }
+        if (i == 3 && j == 1) {
+            return 6.828427124746214_rt;
+        }
+        if (i == 3 && j == 2) {
+            return 3.4142135623731007_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return -6.828427124746214_rt;
+        }
+        if (i == 3 && j == 1) {
+            return -10.949747468305889_rt;
+        }
+        if (i == 3 && j == 2) {
+            return -7.535533905932761_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return 6.828427124746214_rt;
+        }
+        if (i == 2) {
+            return 3.414213562373101_rt;
+        }
+        return 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        if (i == 1) {
+            return -0.23570226039551292_rt;
+        }
+        if (i == 2) {
+            return -0.23570226039551567_rt;
+        }
+        return -0.13807118745769906_rt;
+    }
+};
+
+struct ros2_tableau {
+    static constexpr int stages = 2;
+    static constexpr amrex::Real gamma = 1.7071067811865475_rt;
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        (void) i;
+        return 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return 0.585786437626905_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return -1.1715728752538102_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return 0.8786796564403575_rt;
+        }
+        return 0.2928932188134525_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        (void) i;
+        return 0.2928932188134525_rt;
+    }
+};
+
+struct sandu_a_tableau {
+    static constexpr int stages = 3;
+    static constexpr amrex::Real gamma = 0.7886751345948129_rt; // (3 + sqrt(3)) / 6
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        return i == 1 ? 0.0_rt : 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if ((i == 2 || i == 3) && j == 1) {
+            return 1.2679491924311228_rt; // 3 - sqrt(3)
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return 0.9282032302755092_rt; // -6 + 4 * sqrt(3)
+        }
+        if (i == 3 && (j == 1 || j == 2)) {
+            return -0.4641016151377544_rt; // 3 - 2 * sqrt(3)
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return 1.2679491924311228_rt; // 3 - sqrt(3)
+        }
+        if (i == 3) {
+            return 1.0_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        return i == 3 ? 1.0_rt : 0.0_rt;
+    }
+};
+
+struct sandu_b_tableau {
+    static constexpr int stages = 3;
+    static constexpr amrex::Real gamma = 0.7886751345948129_rt; // (3 + sqrt(3)) / 6
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        return i == 1 ? 0.0_rt : 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if ((i == 2 || i == 3) && j == 1) {
+            return 1.2679491924311228_rt; // 3 - sqrt(3)
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 3 && j == 1) {
+            return -0.5358983848622456_rt; // -4 + 2 * sqrt(3)
+        }
+        if (i == 3 && j == 2) {
+            return -0.7320508075688772_rt; // 1 - sqrt(3)
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return 1.2679491924311228_rt; // 3 - sqrt(3)
+        }
+        if (i == 3) {
+            return 1.0_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        return i == 3 ? 1.0_rt : 0.0_rt;
+    }
+};
+
+struct sandu_c_tableau {
+    static constexpr int stages = 3;
+    static constexpr amrex::Real gamma = 0.78867513459481275_rt; // (3 + sqrt(3)) / 6
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        if (i == 3) {
+            return 0.6666666666666666_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if (i == 3 && j == 1) {
+            return -0.4335531486700170_rt; // (1272 - 823 * sqrt(3)) / 354
+        }
+        if (i == 3 && j == 2) {
+            return 1.7222282832648295_rt; // 3 * (-51 + 49 * sqrt(3)) / 59
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return -0.3264246859454366_rt; // -1 + 7 * sqrt(3) / 18
+        }
+        if (i == 3 && j == 1) {
+            return 0.4138899169340995_rt; // (25 - 13 * sqrt(3)) / 6
+        }
+        if (i == 3 && j == 2) {
+            return 1.6076951545867360_rt; // 12 - 6 * sqrt(3)
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return -7.1285171234651825_rt; // (-12089 + 5037 * sqrt(3)) / 472
+        }
+        if (i == 2) {
+            return 8.4030361763035106_rt; // 9 * (344 - 135 * sqrt(3)) / 118
+        }
+        return 0.9509618943233421_rt; // 3 * (3 - sqrt(3)) / 4
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        if (i == 1) {
+            return -8.396466315896305_rt;
+        }
+        if (i == 2) {
+            return 8.4030361763035106_rt;
+        }
+        return 0.9509618943233421_rt;
+    }
+};
+
+struct sandu_d_tableau {
+    static constexpr int stages = 4;
+    static constexpr amrex::Real gamma = 0.5_rt;
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real ct (const int i) {
+        return i == 1 ? 0.0_rt : 1.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real a (const int i, const int j) {
+        if (i > 1 && j == 1) {
+            return 2.0_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real c (const int i, const int j) {
+        if (i == 2 && j == 1) {
+            return -1.3333333333333333_rt; // -4 / 3
+        }
+        if (i == 3 && j == 1) {
+            return -3.3333333333333333_rt; // -10 / 3
+        }
+        if (i == 3 && j == 2) {
+            return -2.0_rt;
+        }
+        if (i == 4 && j == 1) {
+            return -0.5_rt;
+        }
+        if (i == 4 && j == 3) {
+            return 1.5_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real b (const int i) {
+        if (i == 1) {
+            return 2.0_rt;
+        }
+        if (i == 4) {
+            return 1.0_rt;
+        }
+        return 0.0_rt;
+    }
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static constexpr amrex::Real e (const int i) {
+        if (i == 1) {
+            return -0.6666666666666666_rt; // 2 - 8 / 3
+        }
+        if (i == 2 || i == 3) {
+            return -1.0_rt;
+        }
+        return 1.3333333333333333_rt; // 1 - (-1 / 3)
+    }
+};
+
+using coefficients = ros2s_tableau;
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real& stage (ros2s_t<int_neqs>& ros2s, const int i, const int n)
+{
+    if (i == 1) {
+        return ros2s.ak1(n);
+    }
+    if (i == 2) {
+        return ros2s.ak2(n);
+    }
+    if (i == 3) {
+        return ros2s.ak3(n);
+    }
+    if (i == 4) {
+        return ros2s.ak4(n);
+    }
+    return ros2s.work(n);
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real stage (const ros2s_t<int_neqs>& ros2s, const int i, const int n)
+{
+    if (i == 1) {
+        return ros2s.ak1(n);
+    }
+    if (i == 2) {
+        return ros2s.ak2(n);
+    }
+    if (i == 3) {
+        return ros2s.ak3(n);
+    }
+    if (i == 4) {
+        return ros2s.ak4(n);
+    }
+    return ros2s.work(n);
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+RArray1D& stage_vector (ros2s_t<int_neqs>& ros2s, const int i)
+{
+    if (i == 1) {
+        return ros2s.ak1;
+    }
+    if (i == 2) {
+        return ros2s.ak2;
+    }
+    if (i == 3) {
+        return ros2s.ak3;
+    }
+    if (i == 4) {
+        return ros2s.ak4;
+    }
+    return ros2s.work;
+}
+
+template <typename Ros2sT>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real rtol_for (const Ros2sT& ros2s, const int n)
+{
+    return n == net_ienuc ? ros2s.rtol_enuc : ros2s.rtol_spec;
+}
+
+template <typename Ros2sT>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real atol_for (const Ros2sT& ros2s, const int n)
+{
+    return n == net_ienuc ? ros2s.atol_enuc : ros2s.atol_spec;
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+amrex::Real error_norm (const ros2s_t<int_neqs>& ros2s)
+{
+    amrex::Real err = 0.0_rt;
+    for (int n = 1; n <= int_neqs; ++n) {
+        const amrex::Real sk = atol_for(ros2s, n) +
+            rtol_for(ros2s, n) * amrex::max(std::abs(ros2s.y(n)), std::abs(ros2s.ynew(n)));
+        const amrex::Real term = ros2s.work(n) / sk;
+        err += term * term;
+    }
+    return std::sqrt(err / static_cast<amrex::Real>(int_neqs));
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+bool species_change_valid (const ros2s_t<int_neqs>& ros2s)
+{
+#ifdef STRANG
+    constexpr amrex::Real increase_change_factor = 4.0_rt;
+    constexpr amrex::Real decrease_change_factor = 0.25_rt;
+
+    for (int i = 1; i <= NumSpec; ++i) {
+        if (std::abs(ros2s.y(i)) > integrator_rp::X_reject_buffer * ros2s.atol_spec &&
+            std::abs(ros2s.ynew(i)) > integrator_rp::X_reject_buffer * ros2s.atol_spec &&
+            (std::abs(ros2s.ynew(i)) > increase_change_factor * std::abs(ros2s.y(i)) ||
+             std::abs(ros2s.ynew(i)) < decrease_change_factor * std::abs(ros2s.y(i)))) {
+            return false;
+        }
+    }
+#endif
+
+    return true;
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+bool valid_integrator_state (const amrex::Array1D<amrex::Real, 1, int_neqs>& y)
+{
+    for (int n = 1; n <= int_neqs; ++n) {
+        const amrex::Real yn = y(n);
+        if (yn != yn || std::abs(yn) == std::numeric_limits<amrex::Real>::infinity()) {
+            return false;
+        }
+    }
+
+#ifdef STRANG
+    for (int n = 1; n <= NumSpec; ++n) {
+        if (y(n) < 0.0_rt) {
+            return false;
+        }
+
+        if (! integrator_rp::use_number_densities && y(n) > 1.0_rt) {
+            return false;
+        }
+    }
+
+    if (y(net_ienuc) <= 0.0_rt) {
+        return false;
+    }
+#endif
+
+    return true;
+}
+
+template <typename BurnT, int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void evaluate_jacobian (BurnT& state, ros2s_t<int_neqs>& ros2s, const amrex::Real time)
+{
+    constexpr bool in_jacobian = true;
+
+    rhs(time, state, ros2s, ros2s.ak1, in_jacobian);
+    ros2s.n_rhs += 1;
+
+    if (ros2s.jacobian_type == 1) {
+        jac(time, state, ros2s, ros2s.jac);
+        ros2s.n_jac += 1;
+        return;
+    }
+
+    amrex::Real fac = 0.0_rt;
+    for (int n = 1; n <= int_neqs; ++n) {
+        const amrex::Real weight = 1.0_rt /
+            (rtol_for(ros2s, n) * std::abs(ros2s.y(n)) + atol_for(ros2s, n));
+        fac += (ros2s.ak1(n) * weight) * (ros2s.ak1(n) * weight);
+    }
+    fac = std::sqrt(fac / static_cast<amrex::Real>(int_neqs));
+
+    amrex::Real r0 = 1000.0_rt * std::abs(ros2s.dt) *
+        std::numeric_limits<amrex::Real>::epsilon() *
+        static_cast<amrex::Real>(int_neqs) * fac;
+    if (r0 == 0.0_rt) {
+        r0 = 1.0_rt;
+    }
+
+    for (int j = 1; j <= int_neqs; ++j) {
+        const amrex::Real yj = ros2s.y(j);
+        const amrex::Real weight = rtol_for(ros2s, j) * std::abs(yj) + atol_for(ros2s, j);
+        const amrex::Real dy = amrex::max(std::sqrt(std::numeric_limits<amrex::Real>::epsilon()) *
+                                          std::abs(yj), r0 * weight);
+        ros2s.y(j) += dy;
+
+        rhs(time, state, ros2s, ros2s.work, in_jacobian);
+        const amrex::Real dyinv = 1.0_rt / dy;
+        for (int i = 1; i <= int_neqs; ++i) {
+            ros2s.jac(i, j) = (ros2s.work(i) - ros2s.ak1(i)) * dyinv;
+        }
+
+        ros2s.y(j) = yj;
+    }
+
+    ros2s.n_rhs += int_neqs;
+    ros2s.n_jac += 1;
+}
+
+template <typename BurnT, int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void rhs_at (const amrex::Real time, BurnT& state, ros2s_t<int_neqs>& ros2s,
+             const amrex::Array1D<amrex::Real, 1, int_neqs>& y,
+             amrex::Array1D<amrex::Real, 1, int_neqs>& ydot)
+{
+    amrex::Array1D<amrex::Real, 1, int_neqs> ysave;
+
+    for (int n = 1; n <= int_neqs; ++n) {
+        ysave(n) = ros2s.y(n);
+        ros2s.y(n) = y(n);
+    }
+
+    rhs(time, state, ros2s, ydot);
+
+    for (int n = 1; n <= int_neqs; ++n) {
+        ros2s.y(n) = ysave(n);
+    }
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int decompose (ros2s_t<int_neqs>& ros2s, const amrex::Real fac)
+{
+    for (int j = 1; j <= int_neqs; ++j) {
+        for (int i = 1; i <= int_neqs; ++i) {
+            ros2s.jac(i, j) = -ros2s.jac(i, j);
+        }
+        ros2s.jac(j, j) += fac;
+    }
+
+    int ierr_linpack = 0;
+    if (integrator_rp::linalg_do_pivoting == 1) {
+        constexpr bool allow_pivot{true};
+        dgefa<int_neqs, allow_pivot>(ros2s.jac, ros2s.pivot, ierr_linpack);
+    } else {
+        constexpr bool allow_pivot{false};
+        dgefa<int_neqs, allow_pivot>(ros2s.jac, ros2s.pivot, ierr_linpack);
+    }
+
+    return ierr_linpack;
+}
+
+template <int int_neqs>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void solve (ros2s_t<int_neqs>& ros2s, RArray1D& b)
+{
+    if (integrator_rp::linalg_do_pivoting == 1) {
+        constexpr bool allow_pivot{true};
+        dgesl<int_neqs, allow_pivot>(ros2s.jac, ros2s.pivot, b);
+    } else {
+        constexpr bool allow_pivot{false};
+        dgesl<int_neqs, allow_pivot>(ros2s.jac, ros2s.pivot, b);
+    }
+}
+
+} // namespace ros2s
+
+template <typename BurnT, typename Tableau>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int rosenbrock_integrator (BurnT& state, ros2s_t<integrator_neqs<BurnT>()>& ros2s)
+{
+    using C = Tableau;
+    constexpr int int_neqs = integrator_neqs<BurnT>();
+
+    if (ros2s.tout == ros2s.t) {
+        return IERR_SUCCESS;
+    }
+
+    for (int n = 1; n <= int_neqs; ++n) {
+        if (ros2s::atol_for(ros2s, n) <= 0.0_rt ||
+            ros2s::rtol_for(ros2s, n) <= 10.0_rt * std::numeric_limits<amrex::Real>::epsilon()) {
+            return IERR_TOO_MUCH_ACCURACY_REQUESTED;
+        }
+    }
+
+    ros2s.n_rhs = 0;
+    ros2s.n_jac = 0;
+    ros2s.n_step = 0;
+
+    ros2s.dt = ros2s.tout - ros2s.t;
+
+    const amrex::Real posneg = ros2s.tout >= ros2s.t ? 1.0_rt : -1.0_rt;
+    const amrex::Real hmax = amrex::min(integrator_rp::ode_max_dt,
+                                        std::abs(ros2s.tout - ros2s.t));
+    amrex::Real h = amrex::min(std::abs(ros2s.dt), hmax) * posneg;
+    if (std::abs(h) <= 10.0_rt * std::numeric_limits<amrex::Real>::epsilon()) {
+        h = 1.e-6_rt * posneg;
+    }
+
+    bool reject = false;
+    bool last = false;
+    int nsing = 0;
+    int n_reject = 0;
+    amrex::Real facold = 1.0_rt;
+    amrex::Real errold = 1.0_rt;
+    amrex::Real hopt = h;
+    amrex::Real x = ros2s.t;
+    amrex::Array1D<amrex::Real, 1, int_neqs> rhs_tmp;
+
+    while (true) {
+        if (ros2s.n_step > integrator_rp::ode_max_steps) {
+            ros2s.t = x;
+            ros2s.dt = h;
+            return IERR_TOO_MANY_STEPS;
+        }
+
+        if (0.1_rt * std::abs(h) <= std::abs(x) * std::numeric_limits<amrex::Real>::epsilon()) {
+            ros2s.t = x;
+            ros2s.dt = h;
+            return IERR_DT_UNDERFLOW;
+        }
+
+        if (last) {
+            ros2s.t = x;
+            ros2s.dt = hopt;
+            return IERR_SUCCESS;
+        }
+
+        hopt = h;
+        if ((x + h * (1.0_rt + timestep_safety_factor) - ros2s.tout) * posneg >= 0.0_rt) {
+            h = ros2s.tout - x;
+            last = true;
+        }
+
+        while (true) {
+            ros2s.dt = h;
+
+            if (0.1_rt * std::abs(h) <= std::abs(x) * std::numeric_limits<amrex::Real>::epsilon()) {
+                ros2s.t = x;
+                ros2s.dt = h;
+                return IERR_DT_UNDERFLOW;
+            }
+
+            const amrex::Real fac = 1.0_rt / (h * C::gamma);
+            ros2s::evaluate_jacobian(state, ros2s, x);
+            int ierr_linpack = ros2s::decompose(ros2s, fac);
+
+            if (ierr_linpack != 0) {
+                nsing += 1;
+                if (nsing >= 5) {
+                    ros2s.t = x;
+                    ros2s.dt = h;
+                    return IERR_LU_DECOMPOSITION_ERROR;
+                }
+                h *= 0.5_rt;
+                reject = true;
+                last = false;
+                continue;
+            }
+
+            ros2s::solve(ros2s, ros2s.ak1);
+
+            bool valid_stage_state = true;
+            for (int istage = 2; istage <= C::stages; ++istage) {
+                for (int n = 1; n <= int_neqs; ++n) {
+                    ros2s.ynew(n) = ros2s.y(n);
+                    ros2s::stage(ros2s, istage, n) = 0.0_rt;
+
+                    for (int jstage = 1; jstage < istage; ++jstage) {
+                        const amrex::Real akj = ros2s::stage(ros2s, jstage, n);
+                        ros2s.ynew(n) += C::a(istage, jstage) * akj;
+                        ros2s::stage(ros2s, istage, n) +=
+                            (C::c(istage, jstage) / h) * akj;
+                    }
+                }
+
+                if (! ros2s::valid_integrator_state(ros2s.ynew)) {
+                    h *= 0.25_rt;
+                    reject = true;
+                    n_reject += 1;
+                    last = false;
+                    valid_stage_state = false;
+                    break;
+                }
+
+                ros2s::rhs_at(x + C::ct(istage) * h, state, ros2s, ros2s.ynew, rhs_tmp);
+                ros2s.n_rhs += 1;
+
+                for (int n = 1; n <= int_neqs; ++n) {
+                    ros2s::stage(ros2s, istage, n) += rhs_tmp(n);
+                }
+                ros2s::solve(ros2s, ros2s::stage_vector(ros2s, istage));
+            }
+
+            if (! valid_stage_state) {
+                continue;
+            }
+
+            for (int n = 1; n <= int_neqs; ++n) {
+                amrex::Real y_update = 0.0_rt;
+                amrex::Real err = 0.0_rt;
+
+                for (int istage = 1; istage <= C::stages; ++istage) {
+                    const amrex::Real aki = ros2s::stage(ros2s, istage, n);
+                    y_update += C::b(istage) * aki;
+                    err += C::e(istage) * aki;
+                }
+
+                ros2s.ynew(n) = ros2s.y(n) + y_update;
+                ros2s.work(n) = err;
+            }
+
+            ros2s.n_step += 1;
+
+            const bool valid_state = ros2s::valid_integrator_state(ros2s.ynew);
+            if (! valid_state) {
+                reject = true;
+                n_reject += 1;
+                last = false;
+                h *= 0.25_rt;
+                continue;
+            }
+
+            const amrex::Real err = ros2s::error_norm(ros2s);
+            const bool valid_update = ros2s::species_change_valid(ros2s);
+            constexpr amrex::Real err_min = 1.e-10_rt;
+            const amrex::Real controller_fac = amrex::Clamp(
+                std::pow(1.0_rt / amrex::max(err, err_min),
+                         1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) *
+                std::pow(1.0_rt / amrex::max(errold, err_min),
+                         1.0_rt / (integrator_rp::h211b_b * integrator_rp::h211b_k)) *
+                std::pow(facold, -1.0_rt / integrator_rp::h211b_b),
+                integrator_rp::h211b_fac_min, integrator_rp::h211b_fac_max);
+            facold = controller_fac;
+            errold = err;
+
+            amrex::Real hnew = h * controller_fac;
+
+            if (err <= 1.0_rt && valid_update) {
+                for (int n = 1; n <= int_neqs; ++n) {
+                    ros2s.y(n) = ros2s.ynew(n);
+                }
+                x += h;
+
+                if (std::abs(hnew) > hmax) {
+                    hnew = posneg * hmax;
+                }
+                if (reject) {
+                    hnew = posneg * amrex::min(std::abs(hnew), std::abs(h));
+                }
+                reject = false;
+                n_reject = 0;
+                h = hnew;
+                break;
+            }
+
+            reject = true;
+            n_reject += 1;
+            last = false;
+            if (n_reject >= 2) {
+                hnew *= integrator_rp::h211b_reduction_fac;
+            }
+            h = valid_update ? hnew : 0.25_rt * h;
+        }
+    }
+}
+
+template <typename BurnT>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+int ros2s_integrator (BurnT& state, ros2s_t<integrator_neqs<BurnT>()>& ros2s)
+{
+    if (integrator_rp::rosenbrock_tableau == 1) {
+        return rosenbrock_integrator<BurnT, ros2s::ros2_tableau>(state, ros2s);
+    }
+    if (integrator_rp::rosenbrock_tableau == 2) {
+        return rosenbrock_integrator<BurnT, ros2s::sandu_a_tableau>(state, ros2s);
+    }
+    if (integrator_rp::rosenbrock_tableau == 3) {
+        return rosenbrock_integrator<BurnT, ros2s::sandu_b_tableau>(state, ros2s);
+    }
+    if (integrator_rp::rosenbrock_tableau == 4) {
+        return rosenbrock_integrator<BurnT, ros2s::sandu_c_tableau>(state, ros2s);
+    }
+    if (integrator_rp::rosenbrock_tableau == 5) {
+        return rosenbrock_integrator<BurnT, ros2s::sandu_d_tableau>(state, ros2s);
+    }
+
+    return rosenbrock_integrator<BurnT, ros2s::ros2s_tableau>(state, ros2s);
+}
+
+#endif
diff --git a/integration/ROS2S/ros2s_type.H b/integration/ROS2S/ros2s_type.H
new file mode 100644
index 0000000000..b33266d873
--- /dev/null
+++ b/integration/ROS2S/ros2s_type.H
@@ -0,0 +1,54 @@
+#ifndef ROS2S_TYPE_H
+#define ROS2S_TYPE_H
+
+#include <AMReX_Array.H>
+#include <AMReX_REAL.H>
+
+#include <ArrayUtilities.H>
+#include <network.H>
+#include <integrator_data.H>
+
+#ifdef STRANG
+#include <integrator_type_strang.H>
+#endif
+#ifdef SDC
+#include <integrator_type_sdc.H>
+#endif
+
+// When checking the integration time to see if we're done,
+// be careful with roundoff issues.
+
+const amrex::Real timestep_safety_factor = 1.0e-4_rt;
+
+template <int int_neqs>
+struct ros2s_t {
+
+    amrex::Real t;      // the starting time
+    amrex::Real tout;   // the stopping time
+    amrex::Real dt;     // next internal timestep
+
+    int n_step;
+    int n_rhs;
+    int n_jac;
+
+    amrex::Real atol_spec;
+    amrex::Real rtol_spec;
+
+    amrex::Real atol_enuc;
+    amrex::Real rtol_enuc;
+
+    amrex::Array1D<amrex::Real, 1, int_neqs> y;
+    amrex::Array1D<amrex::Real, 1, int_neqs> ynew;
+    amrex::Array1D<amrex::Real, 1, int_neqs> ak1;
+    amrex::Array1D<amrex::Real, 1, int_neqs> ak2;
+    amrex::Array1D<amrex::Real, 1, int_neqs> ak3;
+    amrex::Array1D<amrex::Real, 1, int_neqs> ak4;
+    amrex::Array1D<amrex::Real, 1, int_neqs> work;
+
+    ArrayUtil::MathArray2D<jac_t, 1, int_neqs, 1, int_neqs> jac;
+    IArray1D pivot;
+
+    short jacobian_type;
+};
+
+#endif
diff --git a/unit_test/burn_cell_primordial_chem/GNUmakefile b/unit_test/burn_cell_primordial_chem/GNUmakefile
index 669e9839c7..9578116d3f 100644
--- a/unit_test/burn_cell_primordial_chem/GNUmakefile
+++ b/unit_test/burn_cell_primordial_chem/GNUmakefile
@@ -2,11 +2,11 @@ PRECISION  = DOUBLE
 PROFILE    = FALSE
 
 # Set DEBUG to TRUE if debugging
-DEBUG      = TRUE
+DEBUG      = FALSE
 
 DIM        = 1
 
-COMP       = gnu
+COMP       = llvm
 
 USE_MPI    = FALSE
 USE_OMP    = FALSE
@@ -15,7 +15,7 @@ USE_CUDA   = FALSE
 USE_REACT = TRUE
 
 # Set USE_MICROPHYSICS_DEBUG to TRUE if debugging
-USE_MICROPHYSICS_DEBUG = TRUE
+USE_MICROPHYSICS_DEBUG = FALSE
 
 EBASE = main
 
@@ -30,7 +30,8 @@ NETWORK_DIR := primordial_chem
 
 CONDUCTIVITY_DIR := stellar
 
-INTEGRATOR_DIR =  VODE
+INTEGRATOR_DIR = VODE
+USERSuffix := .$(INTEGRATOR_DIR)
 
 EXTERN_SEARCH += .
 
diff --git a/unit_test/burn_cell_primordial_chem/burn_cell.H b/unit_test/burn_cell_primordial_chem/burn_cell.H
index f5ced4cb5e..ba06060ac4 100644
--- a/unit_test/burn_cell_primordial_chem/burn_cell.H
+++ b/unit_test/burn_cell_primordial_chem/burn_cell.H
@@ -9,6 +9,7 @@
 #include <extern_parameters.H>
 #include <network.H>
 #include <actual_network.H>
+#include <AMReX_ParallelDescriptor.H>
 
 amrex::Real grav_constant = 6.674e-8;
 
@@ -147,6 +148,7 @@ auto burn_cell_c() -> int {
   // the loop below is similar to that used in krome and pynucastro
   amrex::Real dd = rhotot;
   amrex::Real dd1 = 0.0_rt;
+  amrex::Real collapse_loop_walltime = amrex::ParallelDescriptor::second();
 
   for (int n = 0; n < unit_test_rp::nsteps; n++) {
 
@@ -237,6 +239,14 @@ auto burn_cell_c() -> int {
     }
     state_over_time << std::endl;
   }
+  collapse_loop_walltime =
+      amrex::ParallelDescriptor::second() - collapse_loop_walltime;
+  const int IOProc = amrex::ParallelDescriptor::IOProcessorNumber();
+  amrex::ParallelDescriptor::ReduceRealMax(collapse_loop_walltime, IOProc);
+  if (amrex::ParallelDescriptor::IOProcessor()) {
+    std::cout << "collapse loop walltime: " << collapse_loop_walltime << " s"
+              << std::endl;
+  }
   state_over_time.close();
 
   // output diagnostics to the terminal
diff --git a/unit_test/burn_cell_primordial_chem/inputs_primordial_chem b/unit_test/burn_cell_primordial_chem/inputs_primordial_chem
index 27305ba634..b1f86cdde4 100644
--- a/unit_test/burn_cell_primordial_chem/inputs_primordial_chem
+++ b/unit_test/burn_cell_primordial_chem/inputs_primordial_chem
@@ -46,9 +46,17 @@ integrator.X_reject_buffer = 1e100
 integrator.jacobian = 1
 # do you want to normalize abundances within VODE? (you don't!)
 integrator.renormalize_abundances = 0
-# tolerances
+
+# species tolerances
 integrator.rtol_spec = 1.0e-4
-integrator.atol_spec = 1.0e-4
+integrator.atol_spec = 1.0e-8
+
+# energy tolerances
+integrator.rtol_enuc = 1.0e-6
+integrator.atol_enuc = 1.0
+
+# tableau (0 = ROS2S, 1 = ROS2)
+integrator.rosenbrock_tableau = 1
 
 #assumed redshift for Pop III star formation
 network.redshift = 30.0
diff --git a/unit_test/burn_cell_primordial_chem/pareto_vode_ros2s.py b/unit_test/burn_cell_primordial_chem/pareto_vode_ros2s.py
new file mode 100755
index 0000000000..8991e63c92
--- /dev/null
+++ b/unit_test/burn_cell_primordial_chem/pareto_vode_ros2s.py
@@ -0,0 +1,531 @@
+#!/usr/bin/env python3
+"""Run a VODE/ROS2S tolerance sweep using the internal collapse-loop timer."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import math
+import re
+import statistics
+import subprocess
+import sys
+from datetime import datetime
+from pathlib import Path
+
+
+FLOAT_RE = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eEdD][-+]?\d+)?"
+WALLTIME_RE = re.compile(rf"collapse loop walltime:\s*({FLOAT_RE})\s*s")
+SUCCESS_RE = re.compile(r"successful\?\s+([01])")
+FINAL_RE = {
+    "final_T": re.compile(rf"T final\s*=\s*({FLOAT_RE})"),
+    "final_Eint": re.compile(rf"Eint final\s*=\s*({FLOAT_RE})"),
+    "final_rho": re.compile(rf"rho final\s*=\s*({FLOAT_RE})"),
+}
+
+TOLERANCE_KEYS = (
+    "integrator.rtol_spec",
+    "integrator.atol_spec",
+    "integrator.rtol_enuc",
+    "integrator.atol_enuc",
+)
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description=(
+            "Sweep rtol_spec for the VODE and ROS2S burn_cell executables. "
+            "Other tolerances are scaled from their values in the inputs file."
+        )
+    )
+    parser.add_argument("--input", default="inputs_primordial_chem")
+    parser.add_argument("--vode-exe", default="./main1d.llvm.VODE.ex")
+    parser.add_argument("--ros-exe", default="./main1d.llvm.ROS2S.ex")
+    parser.add_argument("--rtol-min", type=float, default=1.0e-4)
+    parser.add_argument("--rtol-max", type=float, default=1.0e-2)
+    parser.add_argument("--num-rtols", type=int, default=9)
+    parser.add_argument(
+        "--rtols",
+        help="Comma-separated rtol_spec values. Overrides --rtol-min/--rtol-max/--num-rtols.",
+    )
+    parser.add_argument("--repeats", type=int, default=3)
+    parser.add_argument("--output", default="pareto_vode_ros2s.csv")
+    parser.add_argument("--runs-output", default="pareto_vode_ros2s_runs.csv")
+    parser.add_argument("--log-root", default="pareto_vode_ros2s_logs")
+    parser.add_argument(
+        "--reference-output",
+        help=(
+            "Existing burn_cell stdout to use as the accuracy reference. "
+            "If omitted, the script tries reference_solution_vode_rtol_1e-7.out, "
+            "then reference_solution.out, then the strictest successful VODE run."
+        ),
+    )
+    parser.add_argument(
+        "--no-auto-reference",
+        action="store_true",
+        help="Do not automatically use local reference_solution*.out files.",
+    )
+    parser.add_argument("--timeout", type=float, default=300.0)
+    parser.add_argument(
+        "--dry-run",
+        action="store_true",
+        help="Print the commands that would run without executing them.",
+    )
+    return parser.parse_args()
+
+
+def parse_float(value: str) -> float:
+    return float(value.replace("D", "E").replace("d", "e"))
+
+
+def format_float(value: float) -> str:
+    return f"{value:.16e}"
+
+
+def slug_float(value: float) -> str:
+    return f"{value:.3e}".replace("+", "").replace("-", "m").replace(".", "p")
+
+
+def logspace(start: float, stop: float, count: int) -> list[float]:
+    if count <= 0:
+        raise ValueError("--num-rtols must be positive")
+    if start <= 0.0 or stop <= 0.0:
+        raise ValueError("--rtol-min and --rtol-max must be positive")
+    if count == 1:
+        return [start]
+    log_start = math.log10(start)
+    log_stop = math.log10(stop)
+    return [
+        10.0 ** (log_start + i * (log_stop - log_start) / (count - 1))
+        for i in range(count)
+    ]
+
+
+def requested_rtols(args: argparse.Namespace) -> list[float]:
+    if args.rtols:
+        values = [parse_float(value.strip()) for value in args.rtols.split(",")]
+        if not values:
+            raise ValueError("--rtols did not contain any values")
+        return values
+    return logspace(args.rtol_min, args.rtol_max, args.num_rtols)
+
+
+def read_base_tolerances(inputs: Path) -> dict[str, float]:
+    text = inputs.read_text()
+    values: dict[str, float] = {}
+    for key in TOLERANCE_KEYS:
+        pattern = re.compile(rf"^\s*{re.escape(key)}\s*=\s*({FLOAT_RE})", re.MULTILINE)
+        match = pattern.search(text)
+        if match is None:
+            raise ValueError(f"{inputs} is missing {key}")
+        values[key] = parse_float(match.group(1))
+    return values
+
+
+def parse_output(text: str) -> dict[str, object]:
+    result: dict[str, object] = {}
+
+    walltime_match = WALLTIME_RE.search(text)
+    result["walltime_s"] = (
+        parse_float(walltime_match.group(1)) if walltime_match is not None else math.nan
+    )
+
+    success_matches = list(SUCCESS_RE.finditer(text))
+    result["successful"] = (
+        int(success_matches[-1].group(1)) if success_matches else 0
+    )
+
+    for key, pattern in FINAL_RE.items():
+        match = pattern.search(text)
+        result[key] = parse_float(match.group(1)) if match is not None else math.nan
+
+    result["number_densities"] = parse_number_densities(text)
+    return result
+
+
+def parse_number_densities(text: str) -> list[float]:
+    values: list[float] = []
+    in_block = False
+    for line in text.splitlines():
+        if "New number densities" in line:
+            in_block = True
+            continue
+        if not in_block:
+            continue
+        stripped = line.strip()
+        if not stripped:
+            continue
+        try:
+            values.append(parse_float(stripped.split()[0]))
+        except ValueError:
+            if values:
+                break
+    return values
+
+
+def final_vector(parsed: dict[str, object]) -> list[float]:
+    vector = [
+        float(parsed.get("final_T", math.nan)),
+        float(parsed.get("final_Eint", math.nan)),
+        float(parsed.get("final_rho", math.nan)),
+    ]
+    vector.extend(float(value) for value in parsed.get("number_densities", []))
+    return vector
+
+
+def symmetric_linf_relative_error(values: list[float], reference: list[float]) -> float:
+    if not values or len(values) != len(reference):
+        return math.nan
+    error = 0.0
+    for value, ref in zip(values, reference):
+        if not math.isfinite(value) or not math.isfinite(ref):
+            return math.nan
+        denom = abs(value) + abs(ref) + 1.0e-300
+        error = max(error, 2.0 * abs(value - ref) / denom)
+    return error
+
+
+def load_reference(args: argparse.Namespace, cwd: Path) -> tuple[list[float] | None, str]:
+    candidates: list[Path] = []
+    if args.reference_output:
+        candidates.append(Path(args.reference_output))
+    elif not args.no_auto_reference:
+        candidates.extend(
+            [
+                cwd / "reference_solution_vode_rtol_1e-7.out",
+                cwd / "reference_solution.out",
+            ]
+        )
+
+    for path in candidates:
+        if not path.exists():
+            continue
+        parsed = parse_output(path.read_text(errors="replace"))
+        vector = final_vector(parsed)
+        if vector and all(math.isfinite(value) for value in vector):
+            return vector, str(path)
+    return None, ""
+
+
+def command_for_run(
+    executable: Path,
+    inputs: Path,
+    tolerances: dict[str, float],
+) -> list[str]:
+    return [
+        str(executable),
+        str(inputs),
+        *[f"{key}={format_float(value)}" for key, value in tolerances.items()],
+    ]
+
+
+def run_case(
+    *,
+    integrator: str,
+    executable: Path,
+    inputs: Path,
+    tolerances: dict[str, float],
+    rtol_spec: float,
+    repeat: int,
+    run_dir: Path,
+    timeout: float,
+) -> dict[str, object]:
+    run_dir.mkdir(parents=True, exist_ok=True)
+    cmd = command_for_run(executable, inputs, tolerances)
+    (run_dir / "command.txt").write_text(" ".join(cmd) + "\n")
+
+    completed = subprocess.run(
+        cmd,
+        cwd=run_dir,
+        text=True,
+        capture_output=True,
+        timeout=timeout,
+        check=False,
+    )
+    (run_dir / "stdout.txt").write_text(completed.stdout)
+    (run_dir / "stderr.txt").write_text(completed.stderr)
+
+    parsed = parse_output(completed.stdout + "\n" + completed.stderr)
+    row: dict[str, object] = {
+        "integrator": integrator,
+        "repeat": repeat,
+        "rtol_spec": rtol_spec,
+        "atol_spec": tolerances["integrator.atol_spec"],
+        "rtol_enuc": tolerances["integrator.rtol_enuc"],
+        "atol_enuc": tolerances["integrator.atol_enuc"],
+        "returncode": completed.returncode,
+        "successful": int(parsed["successful"]) if completed.returncode == 0 else 0,
+        "walltime_s": parsed["walltime_s"],
+        "final_T": parsed["final_T"],
+        "final_Eint": parsed["final_Eint"],
+        "final_rho": parsed["final_rho"],
+        "stdout": str(run_dir / "stdout.txt"),
+        "stderr": str(run_dir / "stderr.txt"),
+        "state_over_time": str(run_dir / "state_over_time.txt"),
+    }
+    for index, value in enumerate(parsed["number_densities"], start=1):
+        row[f"number_density_{index}"] = value
+    row["_final_vector"] = final_vector(parsed)
+    return row
+
+
+def finite(value: object) -> bool:
+    try:
+        return math.isfinite(float(value))
+    except (TypeError, ValueError):
+        return False
+
+
+def summarize(
+    rows: list[dict[str, object]],
+    repeats: int,
+    reference_vector: list[float] | None,
+) -> list[dict[str, object]]:
+    grouped: dict[tuple[str, float], list[dict[str, object]]] = {}
+    for row in rows:
+        grouped.setdefault((str(row["integrator"]), float(row["rtol_spec"])), []).append(row)
+
+    if reference_vector is None:
+        reference_vector = best_sweep_reference(grouped)
+
+    summary_rows: list[dict[str, object]] = []
+    for (integrator, rtol_spec), group in sorted(grouped.items()):
+        successful = [
+            row
+            for row in group
+            if int(row["successful"]) == 1 and finite(row["walltime_s"])
+        ]
+        walltimes = [float(row["walltime_s"]) for row in successful]
+        representative = successful[0] if successful else group[0]
+        final_vector_values = representative.get("_final_vector", [])
+        error = (
+            symmetric_linf_relative_error(final_vector_values, reference_vector)
+            if reference_vector is not None
+            else math.nan
+        )
+        summary_rows.append(
+            {
+                "integrator": integrator,
+                "rtol_spec": rtol_spec,
+                "atol_spec": representative["atol_spec"],
+                "rtol_enuc": representative["rtol_enuc"],
+                "atol_enuc": representative["atol_enuc"],
+                "repeats": repeats,
+                "successful_runs": len(successful),
+                "successful": int(len(successful) == repeats),
+                "walltime_median_s": statistics.median(walltimes)
+                if walltimes
+                else math.nan,
+                "walltime_mean_s": statistics.mean(walltimes)
+                if walltimes
+                else math.nan,
+                "walltime_min_s": min(walltimes) if walltimes else math.nan,
+                "walltime_max_s": max(walltimes) if walltimes else math.nan,
+                "walltime_stdev_s": statistics.stdev(walltimes)
+                if len(walltimes) > 1
+                else 0.0,
+                "error_linf_rel": error,
+                "pareto_integrator": 0,
+                "pareto_global": 0,
+            }
+        )
+
+    mark_pareto(summary_rows, "pareto_global")
+    for integrator in sorted({row["integrator"] for row in summary_rows}):
+        mark_pareto(
+            [row for row in summary_rows if row["integrator"] == integrator],
+            "pareto_integrator",
+        )
+    return summary_rows
+
+
+def best_sweep_reference(
+    grouped: dict[tuple[str, float], list[dict[str, object]]]
+) -> list[float] | None:
+    candidates: list[tuple[int, float, list[float]]] = []
+    for (integrator, rtol_spec), group in grouped.items():
+        for row in group:
+            if int(row["successful"]) != 1:
+                continue
+            priority = 0 if integrator == "VODE" else 1
+            vector = row.get("_final_vector", [])
+            if vector:
+                candidates.append((priority, rtol_spec, vector))
+                break
+    if not candidates:
+        return None
+    candidates.sort(key=lambda item: (item[0], item[1]))
+    return candidates[0][2]
+
+
+def mark_pareto(rows: list[dict[str, object]], column: str) -> None:
+    eligible = [
+        row
+        for row in rows
+        if int(row["successful"]) == 1
+        and finite(row["walltime_median_s"])
+        and finite(row["error_linf_rel"])
+    ]
+    for row in eligible:
+        time = float(row["walltime_median_s"])
+        error = float(row["error_linf_rel"])
+        dominated = False
+        for other in eligible:
+            if other is row:
+                continue
+            other_time = float(other["walltime_median_s"])
+            other_error = float(other["error_linf_rel"])
+            if (
+                other_time <= time
+                and other_error <= error
+                and (other_time < time or other_error < error)
+            ):
+                dominated = True
+                break
+        row[column] = int(not dominated)
+
+
+def write_csv(path: Path, rows: list[dict[str, object]], fieldnames: list[str]) -> None:
+    with path.open("w", newline="") as csvfile:
+        writer = csv.DictWriter(csvfile, fieldnames=fieldnames, extrasaction="ignore")
+        writer.writeheader()
+        for row in rows:
+            writer.writerow(
+                {
+                    key: format_csv_value(row.get(key, ""))
+                    for key in fieldnames
+                }
+            )
+
+
+def format_csv_value(value: object) -> object:
+    if isinstance(value, float):
+        return format_float(value) if math.isfinite(value) else "nan"
+    return value
+
+
+def main() -> int:
+    args = parse_args()
+    cwd = Path.cwd()
+    inputs = Path(args.input).resolve()
+    executables = {
+        "VODE": Path(args.vode_exe).resolve(),
+        "ROS2S": Path(args.ros_exe).resolve(),
+    }
+
+    if not inputs.exists():
+        raise FileNotFoundError(inputs)
+    for executable in executables.values():
+        if not executable.exists():
+            raise FileNotFoundError(executable)
+    if args.repeats <= 0:
+        raise ValueError("--repeats must be positive")
+
+    base_tolerances = read_base_tolerances(inputs)
+    rtols = requested_rtols(args)
+    base_rtol = base_tolerances["integrator.rtol_spec"]
+
+    reference_vector, reference_source = load_reference(args, cwd)
+    if reference_source:
+        print(f"using accuracy reference: {reference_source}")
+    else:
+        print("using strictest successful sweep run as accuracy reference")
+
+    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
+    log_root = Path(args.log_root).resolve() / timestamp
+
+    dry_commands: list[list[str]] = []
+    run_rows: list[dict[str, object]] = []
+
+    for rtol_spec in rtols:
+        scale = rtol_spec / base_rtol
+        tolerances = {
+            "integrator.rtol_spec": rtol_spec,
+            "integrator.atol_spec": base_tolerances["integrator.atol_spec"] * scale,
+            "integrator.rtol_enuc": base_tolerances["integrator.rtol_enuc"] * scale,
+            "integrator.atol_enuc": base_tolerances["integrator.atol_enuc"] * scale,
+        }
+        for integrator, executable in executables.items():
+            for repeat in range(1, args.repeats + 1):
+                run_dir = (
+                    log_root
+                    / f"{integrator}_rtol_{slug_float(rtol_spec)}_run_{repeat:02d}"
+                )
+                if args.dry_run:
+                    dry_commands.append(command_for_run(executable, inputs, tolerances))
+                    continue
+                print(
+                    f"running {integrator} rtol_spec={rtol_spec:.3e} "
+                    f"repeat {repeat}/{args.repeats}"
+                )
+                run_rows.append(
+                    run_case(
+                        integrator=integrator,
+                        executable=executable,
+                        inputs=inputs,
+                        tolerances=tolerances,
+                        rtol_spec=rtol_spec,
+                        repeat=repeat,
+                        run_dir=run_dir,
+                        timeout=args.timeout,
+                    )
+                )
+
+    if args.dry_run:
+        for cmd in dry_commands:
+            print(" ".join(cmd))
+        return 0
+
+    run_fieldnames = [
+        "integrator",
+        "repeat",
+        "rtol_spec",
+        "atol_spec",
+        "rtol_enuc",
+        "atol_enuc",
+        "returncode",
+        "successful",
+        "walltime_s",
+        "final_T",
+        "final_Eint",
+        "final_rho",
+        *[f"number_density_{index}" for index in range(1, 15)],
+        "stdout",
+        "stderr",
+        "state_over_time",
+    ]
+    summary_fieldnames = [
+        "integrator",
+        "rtol_spec",
+        "atol_spec",
+        "rtol_enuc",
+        "atol_enuc",
+        "repeats",
+        "successful_runs",
+        "successful",
+        "walltime_median_s",
+        "walltime_mean_s",
+        "walltime_min_s",
+        "walltime_max_s",
+        "walltime_stdev_s",
+        "error_linf_rel",
+        "pareto_integrator",
+        "pareto_global",
+    ]
+
+    summary_rows = summarize(run_rows, args.repeats, reference_vector)
+    write_csv(Path(args.runs_output), run_rows, run_fieldnames)
+    write_csv(Path(args.output), summary_rows, summary_fieldnames)
+
+    print(f"wrote {args.output}")
+    print(f"wrote {args.runs_output}")
+    print(f"logs are under {log_root}")
+    return 0
+
+
+if __name__ == "__main__":
+    try:
+        raise SystemExit(main())
+    except Exception as exc:
+        print(f"error: {exc}", file=sys.stderr)
+        raise SystemExit(1)
diff --git a/unit_test/burn_cell_primordial_chem/plot_pareto_vode_ros2s.py b/unit_test/burn_cell_primordial_chem/plot_pareto_vode_ros2s.py
new file mode 100755
index 0000000000..3c6eafb809
--- /dev/null
+++ b/unit_test/burn_cell_primordial_chem/plot_pareto_vode_ros2s.py
@@ -0,0 +1,90 @@
+#!/usr/bin/env python3
+"""Plot VODE/ROS2S walltime versus rtol_spec on log-log axes."""
+
+from __future__ import annotations
+
+import argparse
+import csv
+from collections import defaultdict
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Plot median walltime vs rtol_spec from pareto_vode_ros2s.py output."
+    )
+    parser.add_argument("csv", nargs="?", default="pareto_vode_ros2s.csv")
+    parser.add_argument("-o", "--output", default="pareto_vode_ros2s.png")
+    parser.add_argument(
+        "--title",
+        default="VODE vs ROS2S tolerance sweep",
+    )
+    return parser.parse_args()
+
+
+def read_rows(path: Path) -> list[dict[str, str]]:
+    with path.open(newline="") as csvfile:
+        reader = csv.DictReader(csvfile)
+        required = {
+            "integrator",
+            "rtol_spec",
+            "walltime_median_s",
+            "error_linf_rel",
+            "pareto_global",
+        }
+        missing = required - set(reader.fieldnames or [])
+        if missing:
+            raise ValueError(f"{path} is missing column(s): {', '.join(sorted(missing))}")
+        return list(reader)
+
+
+def main() -> int:
+    args = parse_args()
+    rows = read_rows(Path(args.csv))
+
+    by_integrator: dict[str, list[dict[str, str]]] = defaultdict(list)
+    for row in rows:
+        if int(row["successful"]) != 1:
+            continue
+        rtol = float(row["rtol_spec"])
+        walltime = float(row["walltime_median_s"])
+        if rtol <= 0.0 or walltime <= 0.0:
+            continue
+        by_integrator[row["integrator"]].append(row)
+
+    fig, ax = plt.subplots(figsize=(7.0, 5.0), constrained_layout=True)
+
+    markers = {"VODE": "o", "ROS2S": "s"}
+    for integrator, group in sorted(by_integrator.items()):
+        group = sorted(group, key=lambda row: float(row["rtol_spec"]))
+        x = [float(row["rtol_spec"]) for row in group]
+        y = [float(row["walltime_median_s"]) for row in group]
+        ax.plot(
+            x,
+            y,
+            marker=markers.get(integrator, "o"),
+            linewidth=1.5,
+            markersize=5,
+            label=integrator,
+        )
+
+    ax.set_xscale("log")
+    ax.set_yscale("log")
+    ax.set_xlabel("rtol_spec")
+    ax.set_ylabel("median walltime [s]")
+    ax.set_title(args.title)
+    ax.grid(True, which="both", linestyle=":", linewidth=0.7)
+    ax.legend()
+
+    fig.savefig(args.output, dpi=200)
+    print(f"wrote {args.output}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
diff --git a/unit_test/burn_cell_primordial_chem/run_perf_comparison.py b/unit_test/burn_cell_primordial_chem/run_perf_comparison.py
new file mode 100755
index 0000000000..7fdc801934
--- /dev/null
+++ b/unit_test/burn_cell_primordial_chem/run_perf_comparison.py
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+./pareto_vode_ros2s.py
+python3 plot_pareto_vode_ros2s.py