diff --git a/src/ipc/candidates/candidates.cpp b/src/ipc/candidates/candidates.cpp
index 060edea16..a048df6dc 100644
--- a/src/ipc/candidates/candidates.cpp
+++ b/src/ipc/candidates/candidates.cpp
@@ -9,6 +9,7 @@
 #include <ipc/distance/point_plane.hpp>
 #include <ipc/distance/point_point.hpp>
 #include <ipc/utils/eigen_ext.hpp>
+#include <ipc/utils/profiler.hpp>
 #include <ipc/utils/save_obj.hpp>
 
 #include <igl/remove_unreferenced.h>
@@ -45,6 +46,8 @@ void Candidates::build(
     const double inflation_radius,
     BroadPhase* broad_phase)
 {
+    IPC_TOOLKIT_PROFILE_BLOCK("Candidates::build(static)");
+
     std::unique_ptr<BroadPhase> default_broad_phase;
     if (broad_phase == nullptr) {
         default_broad_phase = make_default_broad_phase();
@@ -130,6 +133,8 @@ void Candidates::build(
     const double inflation_radius,
     BroadPhase* broad_phase)
 {
+    IPC_TOOLKIT_PROFILE_BLOCK("Candidates::build(dynamic)");
+
     std::unique_ptr<BroadPhase> default_broad_phase;
     if (broad_phase == nullptr) {
         default_broad_phase = make_default_broad_phase();
@@ -255,6 +260,7 @@ double Candidates::compute_collision_free_stepsize(
 {
     assert(vertices_t0.rows() == mesh.num_vertices());
     assert(vertices_t1.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("Candidates::compute_collision_free_stepsize");
 
     if (empty()) {
         return 1; // No possible collisions, so can take full step.
diff --git a/src/ipc/collisions/normal/normal_collisions.cpp b/src/ipc/collisions/normal/normal_collisions.cpp
index da63d4605..d24630e06 100644
--- a/src/ipc/collisions/normal/normal_collisions.cpp
+++ b/src/ipc/collisions/normal/normal_collisions.cpp
@@ -2,6 +2,7 @@
 
 #include <ipc/collisions/normal/normal_collisions_builder.hpp>
 #include <ipc/utils/local_to_global.hpp>
+#include <ipc/utils/profiler.hpp>
 
 #include <tbb/blocked_range.h>
 #include <tbb/enumerable_thread_specific.h>
@@ -24,6 +25,7 @@ void NormalCollisions::build(
     BroadPhase* broad_phase)
 {
     assert(vertices.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("NormalCollisions::build");
 
     const double inflation_radius = 0.5 * (dhat + dmin);
 
@@ -41,6 +43,7 @@ void NormalCollisions::build(
     const double dmin)
 {
     assert(vertices.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("NormalCollisions::build(candidates)");
 
     clear();
 
diff --git a/src/ipc/collisions/tangential/tangential_collisions.cpp b/src/ipc/collisions/tangential/tangential_collisions.cpp
index 85702e546..1946afcc3 100644
--- a/src/ipc/collisions/tangential/tangential_collisions.cpp
+++ b/src/ipc/collisions/tangential/tangential_collisions.cpp
@@ -3,6 +3,7 @@
 #include <ipc/distance/edge_edge_mollifier.hpp>
 #include <ipc/friction/smooth_mu.hpp>
 #include <ipc/utils/local_to_global.hpp>
+#include <ipc/utils/profiler.hpp>
 
 #include <tbb/blocked_range.h>
 #include <tbb/enumerable_thread_specific.h>
@@ -69,6 +70,7 @@ void TangentialCollisions::build(
 {
     assert(mu_s.size() == vertices.rows());
     assert(mu_k.size() == vertices.rows());
+    IPC_TOOLKIT_PROFILE_BLOCK("TangentialCollisions::build");
 
     const Eigen::MatrixXi& edges = mesh.edges();
     const Eigen::MatrixXi& faces = mesh.faces();
@@ -180,6 +182,7 @@ void TangentialCollisions::build(
 {
     assert(mu_k.size() == vertices.rows());
     assert(mu_s.size() == vertices.rows());
+    IPC_TOOLKIT_PROFILE_BLOCK("TangentialCollisions::build(smooth)");
 
     const Eigen::MatrixXi& edges = mesh.edges();
     const Eigen::MatrixXi& faces = mesh.faces();
diff --git a/src/ipc/ipc.cpp b/src/ipc/ipc.cpp
index b34c09238..b838f33db 100644
--- a/src/ipc/ipc.cpp
+++ b/src/ipc/ipc.cpp
@@ -11,6 +11,8 @@
 #include <scalable_ccd/cuda/ipc_ccd_strategy.hpp>
 #endif
 
+#include <ipc/utils/profiler.hpp>
+
 #include <igl/predicates/segment_segment_intersect.h>
 
 namespace ipc {
@@ -25,6 +27,7 @@ bool is_step_collision_free(
 {
     assert(vertices_t0.rows() == mesh.num_vertices());
     assert(vertices_t1.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("is_step_collision_free");
 
     // Broad phase
     Candidates candidates;
@@ -49,6 +52,7 @@ double compute_collision_free_stepsize(
 {
     assert(vertices_t0.rows() == mesh.num_vertices());
     assert(vertices_t1.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("compute_collision_free_stepsize");
 
     std::unique_ptr<BroadPhase> default_broad_phase;
     if (broad_phase == nullptr) {
@@ -104,6 +108,7 @@ bool has_intersections(
     BroadPhase* broad_phase)
 {
     assert(vertices.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK("has_intersections");
 
     std::unique_ptr<BroadPhase> default_broad_phase;
     if (broad_phase == nullptr) {
diff --git a/src/ipc/potentials/barrier_potential.hpp b/src/ipc/potentials/barrier_potential.hpp
index 06611355f..bf5e23a70 100644
--- a/src/ipc/potentials/barrier_potential.hpp
+++ b/src/ipc/potentials/barrier_potential.hpp
@@ -74,6 +74,8 @@ class BarrierPotential : public NormalPotential {
         m_barrier = barrier;
     }
 
+    std::string name() const override { return "BarrierPotential"; }
+
     using Super::operator();
     using Super::gradient;
     using Super::hessian;
diff --git a/src/ipc/potentials/friction_potential.hpp b/src/ipc/potentials/friction_potential.hpp
index cf0f1afea..54f473827 100644
--- a/src/ipc/potentials/friction_potential.hpp
+++ b/src/ipc/potentials/friction_potential.hpp
@@ -24,6 +24,8 @@ class FrictionPotential : public TangentialPotential {
         m_eps_v = eps_v;
     }
 
+    std::string name() const override { return "FrictionPotential"; }
+
 protected:
     /// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth mollifier.
     /// @param x The tangential relative speed.
diff --git a/src/ipc/potentials/normal_adhesion_potential.hpp b/src/ipc/potentials/normal_adhesion_potential.hpp
index bb55aa649..92cdf4872 100644
--- a/src/ipc/potentials/normal_adhesion_potential.hpp
+++ b/src/ipc/potentials/normal_adhesion_potential.hpp
@@ -15,6 +15,8 @@ class NormalAdhesionPotential : public NormalPotential {
         const double _Y,
         const double _eps_c);
 
+    std::string name() const override { return "NormalAdhesionPotential"; }
+
     using Super::operator();
     using Super::gradient;
     using Super::hessian;
diff --git a/src/ipc/potentials/potential.cpp b/src/ipc/potentials/potential.cpp
index 548fc5af9..839ae594c 100644
--- a/src/ipc/potentials/potential.cpp
+++ b/src/ipc/potentials/potential.cpp
@@ -3,6 +3,7 @@
 #include <ipc/collisions/normal/normal_collisions.hpp>
 #include <ipc/collisions/tangential/tangential_collisions.hpp>
 #include <ipc/utils/local_to_global.hpp>
+#include <ipc/utils/profiler.hpp>
 
 #include <tbb/blocked_range.h>
 #include <tbb/combinable.h>
@@ -39,6 +40,7 @@ double Potential<TCollisions>::operator()(
     Eigen::ConstRef<Eigen::MatrixXd> X) const
 {
     assert(X.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK(this->name() + "::operator()");
 
     return tbb::parallel_reduce(
         tbb::blocked_range<size_t>(size_t(0), collisions.size()), 0.0,
@@ -61,6 +63,7 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
     Eigen::ConstRef<Eigen::MatrixXd> X) const
 {
     assert(X.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK(this->name() + "::gradient()");
 
     if (collisions.empty()) {
         return Eigen::VectorXd::Zero(X.size());
@@ -70,20 +73,25 @@ Eigen::VectorXd Potential<TCollisions>::gradient(
 
     tbb::combinable<Eigen::VectorXd> grad(Eigen::VectorXd::Zero(X.size()));
 
-    tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
-        const TCollision& collision = collisions[i];
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("compute local gradients");
+        tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
+            const TCollision& collision = collisions[i];
 
-        const VectorMaxNd local_grad = this->gradient(
-            collision, collision.dof(X, mesh.edges(), mesh.faces()));
+            const VectorMaxNd local_grad = this->gradient(
+                collision, collision.dof(X, mesh.edges(), mesh.faces()));
 
-        local_gradient_to_global_gradient(
-            local_grad, collision.vertex_ids(mesh.edges(), mesh.faces()), dim,
-            grad.local());
-    });
+            local_gradient_to_global_gradient(
+                local_grad, collision.vertex_ids(mesh.edges(), mesh.faces()),
+                dim, grad.local());
+        });
+    }
 
-    return grad.combine([](const Eigen::VectorXd& a, const Eigen::VectorXd& b) {
-        return a + b;
-    });
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("combine local gradients");
+        return grad.combine([](const Eigen::VectorXd& a,
+                               const Eigen::VectorXd& b) { return a + b; });
+    }
 }
 
 template <class TCollisions>
@@ -94,6 +102,7 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     const PSDProjectionMethod project_hessian_to_psd) const
 {
     assert(X.rows() == mesh.num_vertices());
+    IPC_TOOLKIT_PROFILE_BLOCK(this->name() + "::hessian()");
 
     if (collisions.empty()) {
         return Eigen::SparseMatrix<double>(X.size(), X.size());
@@ -111,20 +120,30 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     tbb::enumerable_thread_specific<LocalThreadMatStorage> storage(
         LocalThreadMatStorage(buffer_size, ndof, ndof));
 
-    tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
-        auto& hess_triplets = storage.local();
-
-        const TCollision& collision = collisions[i];
-
-        const MatrixMaxNd local_hess = this->hessian(
-            collisions[i], collisions[i].dof(X, edges, faces),
-            project_hessian_to_psd);
-
-        local_hessian_to_global_triplets(
-            local_hess, collision.vertex_ids(edges, faces), dim,
-            *(hess_triplets.cache), mesh.num_vertices());
-    });
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("compute local hessians and triplets");
+        tbb::parallel_for(size_t(0), collisions.size(), [&](size_t i) {
+            auto& hess_triplets = storage.local();
+
+            const TCollision& collision = collisions[i];
+
+            MatrixMaxNd local_hess;
+            {
+                IPC_TOOLKIT_PROFILE_BLOCK("compute local hessian");
+                local_hess = this->hessian(
+                    collision, collision.dof(X, edges, faces),
+                    project_hessian_to_psd);
+            }
 
+            {
+                IPC_TOOLKIT_PROFILE_BLOCK(
+                    "map local hessian to global triplets");
+                local_hessian_to_global_triplets(
+                    local_hess, collision.vertex_ids(edges, faces), dim,
+                    *(hess_triplets.cache), mesh.num_vertices());
+            }
+        });
+    }
     if (storage.empty()) {
         return Eigen::SparseMatrix<double>();
     }
@@ -132,9 +151,12 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
     // Assemble the stiffness matrix by concatenating the tuples in each local
     // storage
 
-    tbb::parallel_for_each(
-        storage.begin(), storage.end(),
-        [](const auto& local_storage) { local_storage.cache->prune(); });
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("prune local storages");
+        tbb::parallel_for_each(
+            storage.begin(), storage.end(),
+            [](const auto& local_storage) { local_storage.cache->prune(); });
+    }
 
     // Prepares for parallel concatenation
     std::vector<size_t> offsets(storage.size());
@@ -164,26 +186,37 @@ Eigen::SparseMatrix<double> Potential<TCollisions>::hessian(
         return hess;
     }
 
-    triplets.resize(triplet_count);
+    // Allocate triplets
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("allocate triplets");
+        triplets.resize(triplet_count);
+    }
 
     // Parallel copy into triplets
-    tbb::parallel_for(size_t(0), storage.size(), [&](size_t i) {
-        const SparseMatrixCache& cache = dynamic_cast<const SparseMatrixCache&>(
-            *((storage.begin() + i)->cache));
-        size_t offset = offsets[i];
-
-        std::copy(
-            cache.entries().begin(), cache.entries().end(),
-            triplets.begin() + offset);
-        offset += cache.entries().size();
-
-        if (cache.mat().nonZeros() > 0) {
-            set_triplets(cache.mat(), triplets, offset);
-        }
-    });
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("parallel copy into triplets");
+        tbb::parallel_for(size_t(0), storage.size(), [&](size_t i) {
+            const SparseMatrixCache& cache =
+                dynamic_cast<const SparseMatrixCache&>(
+                    *((storage.begin() + i)->cache));
+            size_t offset = offsets[i];
+
+            std::copy(
+                cache.entries().begin(), cache.entries().end(),
+                triplets.begin() + offset);
+            offset += cache.entries().size();
+
+            if (cache.mat().nonZeros() > 0) {
+                set_triplets(cache.mat(), triplets, offset);
+            }
+        });
+    }
 
     // Sort and assemble
-    hess.setFromTriplets(triplets.begin(), triplets.end());
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("assemble hessian from triplets");
+        hess.setFromTriplets(triplets.begin(), triplets.end());
+    }
 
     return hess;
 }
diff --git a/src/ipc/potentials/potential.hpp b/src/ipc/potentials/potential.hpp
index cedde81bd..2c6c42785 100644
--- a/src/ipc/potentials/potential.hpp
+++ b/src/ipc/potentials/potential.hpp
@@ -3,6 +3,8 @@
 #include <ipc/collision_mesh.hpp>
 #include <ipc/utils/eigen_ext.hpp>
 
+#include <string>
+
 namespace ipc {
 
 /// @brief Base class for potentials.
@@ -19,6 +21,9 @@ template <class TCollisions> class Potential {
     Potential() = default;
     virtual ~Potential() = default;
 
+    /// @brief The name of this potential (used for profiling).
+    virtual std::string name() const = 0;
+
     // -- Cumulative methods ---------------------------------------------------
 
     /// @brief Compute the potential for a set of collisions.
diff --git a/src/ipc/potentials/tangential_adhesion_potential.hpp b/src/ipc/potentials/tangential_adhesion_potential.hpp
index c3290e4d4..ee5250f8f 100644
--- a/src/ipc/potentials/tangential_adhesion_potential.hpp
+++ b/src/ipc/potentials/tangential_adhesion_potential.hpp
@@ -24,6 +24,8 @@ class TangentialAdhesionPotential : public TangentialPotential {
         m_eps_a = eps_a;
     }
 
+    std::string name() const override { return "TangentialAdhesionPotential"; }
+
 protected:
     /// @brief Compute the value of the ∫ μ(y) f₁(y) dy, where f₁ is the first derivative of the smooth mollifier.
     /// @param x The tangential relative speed.
diff --git a/src/ipc/utils/profiler.cpp b/src/ipc/utils/profiler.cpp
index 12bbfaeb6..a5e61d7b2 100644
--- a/src/ipc/utils/profiler.cpp
+++ b/src/ipc/utils/profiler.cpp
@@ -7,7 +7,7 @@
 
 namespace ipc {
 
-Profiler::Profiler() { }
+Profiler::Profiler() : m_main_thread_id(std::this_thread::get_id()) { }
 
 Profiler& profiler()
 {
@@ -15,13 +15,22 @@ Profiler& profiler()
     return instance;
 }
 
-void Profiler::clear() { m_data.clear(); }
+void Profiler::clear()
+{
+    m_data = nlohmann::json::object();
+    m_current_scope = nlohmann::json::json_pointer(); // root
+}
 
 void Profiler::start(const std::string& name)
 {
-    current_scope.push_back(name);
-    if (!m_data.contains(current_scope)) {
-        m_data[current_scope] = {
+    if (std::this_thread::get_id() != m_main_thread_id) {
+        return;
+    }
+
+    m_current_scope.push_back(name);
+
+    if (!m_data.contains(m_current_scope)) {
+        m_data[m_current_scope] = {
             { "time_ms", 0 },
             { "count", 0 },
         };
@@ -30,28 +39,34 @@ void Profiler::start(const std::string& name)
 
 void Profiler::stop(const double time_ms)
 {
+    if (std::this_thread::get_id() != m_main_thread_id) {
+        return;
+    }
+
     const static std::string log_fmt_text = fmt::format(
         "[{}] {{}} {{:.6f}} ms",
         fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing"));
 
     logger().trace(
-        fmt::runtime(log_fmt_text), current_scope.to_string(), time_ms);
-
-    assert(m_data.contains(current_scope));
-    assert(m_data.at(current_scope).contains("time_ms"));
-    assert(m_data.at(current_scope).contains("count"));
-    m_data[current_scope]["time_ms"] =
-        m_data[current_scope]["time_ms"].get<double>() + time_ms;
-    m_data[current_scope]["count"] =
-        m_data[current_scope]["count"].get<size_t>() + 1;
-    current_scope.pop_back();
+        fmt::runtime(log_fmt_text), m_current_scope.to_string(), time_ms);
+
+    assert(m_data.contains(m_current_scope));
+    assert(m_data.at(m_current_scope).contains("time_ms"));
+    assert(m_data.at(m_current_scope).contains("count"));
+    m_data[m_current_scope]["time_ms"] =
+        m_data[m_current_scope]["time_ms"].get<double>() + time_ms;
+    m_data[m_current_scope]["count"] =
+        m_data[m_current_scope]["count"].get<size_t>() + 1;
+
+    m_current_scope.pop_back();
 }
 
 void Profiler::reset()
 {
+    m_main_thread_id = std::this_thread::get_id();
     m_data.clear();
     // reset the calling thread's scope
-    current_scope = nlohmann::json::json_pointer(); // root
+    m_current_scope = nlohmann::json::json_pointer(); // root
 }
 
 void Profiler::print() const
@@ -118,4 +133,4 @@ void Profiler::write_csv(std::ostream& os) const
 
 } // namespace ipc
 
-#endif
\ No newline at end of file
+#endif
diff --git a/src/ipc/utils/profiler.hpp b/src/ipc/utils/profiler.hpp
index 2f560247e..766fd2097 100644
--- a/src/ipc/utils/profiler.hpp
+++ b/src/ipc/utils/profiler.hpp
@@ -13,6 +13,7 @@
 
 #include <chrono>
 #include <iostream>
+#include <thread>
 
 // Helper macro to stringify/paste after expansion
 #define IPC_TOOLKIT_PROFILE_BLOCK_CONCAT_IMPL(a, b) a##b
@@ -63,12 +64,21 @@ class Profiler {
     /// @brief Access the profiling data as a JSON object.
     nlohmann::json& data() { return m_data; }
 
+    bool is_recording_thread() const
+    {
+        return std::this_thread::get_id() == m_main_thread_id;
+    }
+
 protected:
     /// @brief The profiling data stored as a JSON object.
     nlohmann::json m_data;
 
     /// @brief The global scope pointer into the JSON data.
-    nlohmann::json::json_pointer current_scope;
+    nlohmann::json::json_pointer m_current_scope;
+
+    /// @brief The thread that records data; calls from all other threads are
+    ///        silently ignored, giving a single-thread estimate of block costs.
+    std::thread::id m_main_thread_id;
 };
 
 Profiler& profiler();
@@ -96,20 +106,26 @@ template <class Timer = ChronoTimer> class ProfilePoint {
 
     ProfilePoint(Profiler& p_profiler, const std::string& name)
         : m_profiler(p_profiler)
+        , m_active(p_profiler.is_recording_thread())
     {
-        m_profiler.start(name);
-        timer.start();
+        if (m_active) {
+            m_profiler.start(name);
+            timer.start();
+        }
     }
 
     ~ProfilePoint()
     {
-        timer.stop();
-        m_profiler.stop(timer.getElapsedTimeInMilliSec());
+        if (m_active) {
+            timer.stop();
+            m_profiler.stop(timer.getElapsedTimeInMilliSec());
+        }
     }
 
 protected:
     Profiler& m_profiler;
     Timer timer;
+    bool m_active;
 };
 
 } // namespace ipc
diff --git a/tests/src/tests/utils/CMakeLists.txt b/tests/src/tests/utils/CMakeLists.txt
index e04072a47..8acd17887 100644
--- a/tests/src/tests/utils/CMakeLists.txt
+++ b/tests/src/tests/utils/CMakeLists.txt
@@ -1,8 +1,9 @@
 set(SOURCES
   # Tests
-  test_utils.cpp
-  test_matrixcache.cpp
   test_local_to_global.cpp
+  test_matrixcache.cpp
+  test_profiler.cpp
+  test_utils.cpp
 
   # Benchmarks
 
diff --git a/tests/src/tests/utils/test_profiler.cpp b/tests/src/tests/utils/test_profiler.cpp
new file mode 100644
index 000000000..c24cb20b7
--- /dev/null
+++ b/tests/src/tests/utils/test_profiler.cpp
@@ -0,0 +1,192 @@
+#include <tests/config.hpp>
+
+#ifdef IPC_TOOLKIT_WITH_PROFILER
+
+#include <ipc/ipc.hpp>
+#include <ipc/ccd/additive_ccd.hpp>
+#include <ipc/collisions/normal/normal_collisions.hpp>
+#include <ipc/collisions/tangential/tangential_collisions.hpp>
+#include <ipc/potentials/barrier_potential.hpp>
+#include <ipc/potentials/friction_potential.hpp>
+#include <ipc/utils/profiler.hpp>
+
+#include <igl/edges.h>
+#include <igl/read_triangle_mesh.h>
+
+#include <catch2/catch_all.hpp>
+#include <tbb/parallel_for.h>
+#include <tbb/task_arena.h>
+
+using namespace ipc;
+
+TEST_CASE("Profiler", "[profiler]")
+{
+    constexpr int sleep_time_ms = 100;
+    constexpr int num_threads = 10;
+
+    {
+        IPC_TOOLKIT_PROFILE_BLOCK("Block 1");
+        std::this_thread::sleep_for(std::chrono::milliseconds(sleep_time_ms));
+    }
+
+    nlohmann::json data = profiler().data();
+
+    CHECK(data.size() == 1);
+
+    REQUIRE(data.contains("Block 1"));
+    nlohmann::json block1 = data.at("Block 1");
+
+    CHECK(block1.size() == 2); // count, time_ms
+    CHECK(
+        block1["time_ms"].get<double>()
+        == Catch::Approx(sleep_time_ms).margin(10));
+    CHECK(block1["count"].get<int>() == 1);
+
+    // ---------------------------------------------------------------------
+
+    tbb::parallel_for(0, num_threads, [&](int) {
+        // for (int i = 0; i < num_threads; ++i) {
+        IPC_TOOLKIT_PROFILE_BLOCK("Block 2");
+        std::this_thread::sleep_for(
+            std::chrono::milliseconds(sleep_time_ms / num_threads));
+        // }
+    });
+
+    data = profiler().data();
+    profiler().print();
+
+    CHECK(data.size() == 2);
+
+    REQUIRE(data.contains("Block 2"));
+    nlohmann::json block2 = data.at("Block 2");
+    CHECK(block2.size() == 2); // count, time_ms
+    // Only the coordinator thread's iterations are recorded
+    CHECK(block2["count"].get<int>() >= 1);
+
+    // ---------------------------------------------------------------------
+
+    auto foo = []() {
+        IPC_TOOLKIT_PROFILE_BLOCK("Block 3");
+        {
+            IPC_TOOLKIT_PROFILE_BLOCK("Block 4");
+            tbb::parallel_for(0, num_threads, [&](int) {
+                // for (int i = 0; i < num_threads; ++i) {
+                {
+                    IPC_TOOLKIT_PROFILE_BLOCK("Block 5");
+                    std::this_thread::sleep_for(
+                        std::chrono::milliseconds(sleep_time_ms / num_threads));
+                }
+
+                {
+                    IPC_TOOLKIT_PROFILE_BLOCK("Block 6");
+                    std::this_thread::sleep_for(
+                        std::chrono::milliseconds(sleep_time_ms / num_threads));
+                }
+                // }
+            });
+        }
+    };
+
+    foo();
+    foo(); // Run it twice to check that counts are aggregated correctly
+
+    data = profiler().data();
+    profiler().print();
+
+    CHECK(data.size() == 3);
+
+    REQUIRE(data.contains("Block 3"));
+    CHECK_FALSE(data.contains("Block 4"));
+    nlohmann::json block3 = data.at("Block 3");
+    CHECK(block3.size() == 3); // count, time_ms, Block 4
+    CHECK(block3["count"].get<int>() == 2);
+
+    REQUIRE(block3.contains("Block 4"));
+    nlohmann::json block4 = block3.at("Block 4");
+
+    REQUIRE(block4.contains("Block 5"));
+    nlohmann::json block5 = block4.at("Block 5");
+    CHECK(block5.size() == 2); // count, time_ms
+    // foo() runs twice; coordinator records at least once per call
+    CHECK(block5["count"].get<int>() >= 2);
+    CHECK(block5["time_ms"].get<double>() < block4["time_ms"].get<double>());
+
+    REQUIRE(block4.contains("Block 6"));
+    nlohmann::json block6 = block4.at("Block 6");
+    CHECK(block6.size() == 2); // count, time_ms
+    CHECK(block6["count"].get<int>() >= 2);
+    CHECK(block6["time_ms"].get<double>() < block4["time_ms"].get<double>());
+
+    CHECK(
+        block5["time_ms"].get<double>() + block6["time_ms"].get<double>()
+        < block4["time_ms"].get<double>());
+}
+
+TEST_CASE("Profile full pipeline", "[!benchmark][profiler]")
+{
+    // const std::string mesh_t0 = "cloth_ball92.ply";
+    // const std::string mesh_t1 = "cloth_ball93.ply";
+    const std::string mesh_t0 = "rod-twist/3036.ply";
+    const std::string mesh_t1 = "rod-twist/3037.ply";
+
+    Eigen::MatrixXd V0_full, V1_full;
+    Eigen::MatrixXi F0, F1;
+
+    const bool loaded_t0 = igl::read_triangle_mesh(
+        (ipc::tests::DATA_DIR / mesh_t0).string(), V0_full, F0);
+    const bool loaded_t1 = igl::read_triangle_mesh(
+        (ipc::tests::DATA_DIR / mesh_t1).string(), V1_full, F1);
+
+    if (!loaded_t0 || !loaded_t1) {
+        WARN("Skipping profiler test: rod-twist data not found");
+        return;
+    }
+
+    REQUIRE(F0.rows() == F1.rows());
+    REQUIRE(V0_full.rows() == V1_full.rows());
+
+    Eigen::MatrixXi E;
+    igl::edges(F0, E);
+
+    CollisionMesh mesh = CollisionMesh::build_from_full_mesh(V0_full, E, F0);
+    Eigen::MatrixXd V0 = mesh.vertices(V0_full);
+    Eigen::MatrixXd V1 = mesh.vertices(V1_full);
+
+    profiler().reset();
+
+    const double dhat = 1e-3;
+    const double mu = 0.3;
+    const double eps_v = 1e-3;
+
+    // CCD step size
+    const double step = compute_collision_free_stepsize(
+        mesh, V0, V1, 0.0, nullptr, AdditiveCCD());
+
+    // Normal collisions and barrier potential
+    NormalCollisions collisions;
+    collisions.build(mesh, V0, dhat);
+
+    BarrierPotential bp(dhat, /*stiffness=*/1e4);
+    bp(collisions, mesh, V0);
+    bp.gradient(collisions, mesh, V0);
+    bp.hessian(collisions, mesh, V0);
+
+    // Tangential (friction) collisions and friction potential
+    TangentialCollisions tangential;
+    tangential.build(mesh, V0, collisions, bp, mu);
+
+    FrictionPotential fp(eps_v);
+    fp(tangential, mesh, V0);
+    fp.gradient(tangential, mesh, V0);
+    fp.hessian(tangential, mesh, V0);
+
+    const std::string output = "profiler_output.csv";
+    profiler().write_csv(output);
+
+    logger().info(
+        "Profiler output written to: {}\n"
+        "  vertices: {}, faces: {}, collisions: {}, step_size: {:.6g}\n",
+        output, V0.rows(), F0.rows(), collisions.size(), step);
+}
+
+#endif
\ No newline at end of file