diff --git a/make/tests b/make/tests
index 0fbe233945b..aa425b70fe9 100644
--- a/make/tests
+++ b/make/tests
@@ -38,7 +38,7 @@ endif
 
 ##
 # Adding a test for multiple translation units. If this fails,
-# a new function is probably missing an inline
+# a public umbrella header likely exposes a new non-inline external definition.
 ##
 
 ifneq ($(OS),Windows_NT)
diff --git a/runChecks.py b/runChecks.py
index 1cb3636a579..23898e24a9d 100755
--- a/runChecks.py
+++ b/runChecks.py
@@ -142,6 +142,47 @@ def check_non_test_files_in_test():
     return errors
 
 
+def check_rev_test_fixtures():
+    test_files = [
+        x for x in files_in_folder("test/unit/math/rev")
+        if os.path.isfile(x) and x.endswith(testsfx)
+    ]
+    errors = []
+    for filepath in test_files:
+        line_num = 0
+        multi_line_comment = False
+        old_state_multi_line_comment = False
+        with open(filepath, "r") as f:
+            for line in f:
+                line_num += 1
+                if multi_line_comment:
+                    if re.search("\*/", line):
+                        multi_line_comment = False
+                else:
+                    if re.search("/\*", line):
+                        multi_line_comment = True
+                if not multi_line_comment or (
+                    multi_line_comment and not old_state_multi_line_comment
+                ):
+                    if (
+                        not re.search(r".*\bTEST\(.*\*/.*", line)
+                        and not re.search(r".*/\*.*\bTEST\(", line)
+                        and not re.search(r".*//.*\bTEST\(", line)
+                        and re.search(r"\bTEST\(", line)
+                    ):
+                        errors.append(
+                            filepath
+                            + " at line "
+                            + str(line_num)
+                            + ":\n\t[rev-tests] Reverse-mode tests in "
+                            + "test/unit/math/rev must use a cleanup fixture. "
+                            + "Replace raw TEST(...) with TEST_F(AgradRev, ...) "
+                            + "or another approved fixture-based form."
+                        )
+                old_state_multi_line_comment = multi_line_comment
+    return errors
+
+
 def main():
     errors = []
     # Check for files inside stan/math/prim that contain stan/math/rev or stan/math/fwd
@@ -242,6 +283,7 @@ def main():
     )
 
     errors.extend(check_non_test_files_in_test())
+    errors.extend(check_rev_test_fixtures())
 
     errors.extend(check_non_unique_test_names())
 
diff --git a/test/unit/math/rev/core/precomputed_gradients_test.cpp b/test/unit/math/rev/core/precomputed_gradients_test.cpp
index d7d1007407a..a0f8fd1cabe 100644
--- a/test/unit/math/rev/core/precomputed_gradients_test.cpp
+++ b/test/unit/math/rev/core/precomputed_gradients_test.cpp
@@ -131,8 +131,7 @@ TEST_F(AgradRev, StanAgradRevInternal_precomputed_gradients_containers) {
   stan::math::recover_memory();
 }
 
-TEST(StanAgradRevInternal,
-     precomputed_gradients_containers_direct_construction) {
+TEST_F(AgradRev, precomputed_gradients_containers_direct_construction) {
   double value = 1;
   std::vector<stan::math::var> vars;
   std::vector<double> gradients;
diff --git a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp
index 39a9bfe145f..6d7453dab33 100644
--- a/test/unit/math/rev/functor/integrate_1d_impl_test.cpp
+++ b/test/unit/math/rev/functor/integrate_1d_impl_test.cpp
@@ -237,7 +237,7 @@ inline void test_derivatives(const F &f, double a, double b,
   }
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_test_integer_arguments) {
+TEST_F(AgradRev, integrate_1d_impl_rev_test_integer_arguments) {
   stan::math::var v;
   std::vector<stan::math::var> theta = {0.5};
   std::vector<double> x_r;
@@ -250,7 +250,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_test_integer_arguments) {
                                                     theta, x_r, x_i));
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_easy) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_easy) {
   // Easy integrals
   using stan::math::var;
   test_derivatives<double, double, var>(f1{}, 0.2, 0.7, {0.75}, {}, {},
@@ -267,7 +267,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_easy) {
                                         {0.0});
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_zero_crossing) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_zero_crossing) {
   // Zero crossing integral + test x_r + vars at endpoints
   using stan::math::var;
   test_derivatives<var, var, var>(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9}, {2.5, 3.0},
@@ -278,8 +278,8 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_zero_crossing) {
                                   -19.06340613646808, 21.41380852375568);
 }
 
-TEST(StanMath_integrate_1d_impl_rev,
-     TestDerivatives_var_right_endpoint_var_params) {
+TEST_F(AgradRev,
+       integrate_1d_impl_rev_TestDerivatives_var_right_endpoint_var_params) {
   // Zero crossing integral + test x_r + vars at right endpoint
   using stan::math::var;
   test_derivatives<double, var, var>(
@@ -289,8 +289,8 @@ TEST(StanMath_integrate_1d_impl_rev,
       {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, 0.0, 21.41380852375568);
 }
 
-TEST(StanMath_integrate_1d_impl_rev,
-     TestDerivatives_var_left_endpoint_var_params) {
+TEST_F(AgradRev,
+       integrate_1d_impl_rev_TestDerivatives_var_left_endpoint_var_params) {
   // Zero crossing integral + test x_r + var at left endpoint
   using stan::math::var;
   test_derivatives<var, double, var>(
@@ -300,7 +300,7 @@ TEST(StanMath_integrate_1d_impl_rev,
       {5 * pow(0.5, 1.5), 12 * 1.75 * 1.75, 4.0}, -19.06340613646808, 0.0);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_no_param_vars) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_no_param_vars) {
   // No param vars
   using stan::math::var;
   test_derivatives<var, var, double>(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9},
@@ -310,8 +310,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_no_param_vars) {
                                      {}, -19.06340613646808, 21.41380852375568);
 }
 
-TEST_F(AgradRev,
-       StanMath_integrate_1d_impl_rev_TestDerivatives_left_limit_var) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_left_limit_var) {
   // No param vars, only left limit var
   using stan::math::var;
   test_derivatives<var, double, double>(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9},
@@ -321,8 +320,7 @@ TEST_F(AgradRev,
                                         {}, -19.06340613646808, 0.0);
 }
 
-TEST_F(AgradRev,
-       StanMath_integrate_1d_impl_rev_TestDerivatives_right_limit_var) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_right_limit_var) {
   // No param vars, only right limit var
   using stan::math::var;
   test_derivatives<double, var, double>(f3{}, -1.0, 1.0, {0.5, 1.75, 3.9},
@@ -332,7 +330,7 @@ TEST_F(AgradRev,
                                         {}, 0.0, 21.41380852375568);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_tricky1) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_tricky1) {
   // Tricky integral from Boost docs + limit at infinity + no gradients
   using stan::math::var;
   test_derivatives<double, double, var>(f4{}, 0.0,
@@ -340,7 +338,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_tricky1) {
                                         {}, {}, {}, 1.772453850905516, {});
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_tricky2) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_tricky2) {
   // Tricky integral from Boost docs + limit at infinity with gradients
   using stan::math::var;
   test_derivatives<double, double, var>(
@@ -350,15 +348,14 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_tricky2) {
        -1.772453850905516 * 0.5 / (2 * pow(0.5 * 3.0, 1.5))});
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_tricky3) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_tricky3) {
   // Tricky integral from Boost docs
   using stan::math::var;
   test_derivatives<double, double, var>(
       f6{}, 0.0, 1.0, {0.75}, {}, {}, 0.851926727945904, {0.4814066053874294});
 }
 
-TEST_F(AgradRev,
-       StanMath_integrate_1d_impl_rev_TestDerivatives_zero_crossing2) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_zero_crossing2) {
   // Zero crossing integral + limit at infinity + var at left limit
   using stan::math::var;
   test_derivatives<var, double, var>(
@@ -367,8 +364,7 @@ TEST_F(AgradRev,
       std::numeric_limits<double>::quiet_NaN());
 }
 
-TEST_F(AgradRev,
-       StanMath_integrate_1d_impl_rev_TestDerivatives_zero_crossing3) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_zero_crossing3) {
   // Zero crossing integral + limit at negative infinity + var at right limit
   using stan::math::var;
   test_derivatives<double, var, var>(
@@ -377,7 +373,7 @@ TEST_F(AgradRev,
       std::numeric_limits<double>::quiet_NaN(), 1808.042414456063);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_indefinite) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_indefinite) {
   // Both limits at infinity + test x_r/x_i + no gradients
   using stan::math::var;
   test_derivatives<double, double, var>(
@@ -386,8 +382,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_indefinite) {
       2.536571480364399, {});
 }
 
-TEST_F(AgradRev,
-       StanMath_integrate_1d_impl_rev_TestDerivatives_endpoint_precision) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_endpoint_precision) {
   // Various integrals of beta function
   using stan::math::var;
   test_derivatives<double, double, var>(f11{}, 0.0, 1.0, {0.1, 0.1}, {}, {},
@@ -407,7 +402,7 @@ TEST_F(AgradRev,
       {-0.01040816326530613, -0.004852607709750566});
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_gaussian) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDerivatives_gaussian) {
   // Check Gaussian integrates to 1.0 always
   using stan::math::var;
   test_derivatives<double, double, var>(
@@ -416,8 +411,8 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDerivatives_gaussian) {
       {0.0, 0.0});
 }
 
-TEST(StanMath_integrate_1d_impl_rev,
-     TestDerivativesSameVarAtEndpointAndInParams) {
+TEST_F(AgradRev,
+       integrate_1d_impl_rev_TestDerivativesSameVarAtEndpointAndInParams) {
   using stan::math::var;
 
   var a = 2.0;
@@ -435,7 +430,7 @@ TEST(StanMath_integrate_1d_impl_rev,
   EXPECT_LT(std::abs(12.0 - b.adj()), 1e-8);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestBeta) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestBeta) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -459,7 +454,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestBeta) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestCauchy) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestCauchy) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -485,7 +480,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestCauchy) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestChiSquare) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestChiSquare) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -508,7 +503,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestChiSquare) {
   EXPECT_FLOAT_EQ(1, 1 + g[0]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDoubleExponential) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestDoubleExponential) {
   using stan::math::exp;
   using stan::math::integrate_1d;
   using stan::math::var;
@@ -536,7 +531,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestDoubleExponential) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestExponential) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestExponential) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -559,7 +554,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestExponential) {
   EXPECT_FLOAT_EQ(1, 1 + g[0]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestFrechet) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestFrechet) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -585,7 +580,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestFrechet) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestGamma) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestGamma) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -611,7 +606,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestGamma) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestGumbel) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestGumbel) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -637,7 +632,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestGumbel) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestInvChiSquared) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestInvChiSquared) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -661,7 +656,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestInvChiSquared) {
   EXPECT_FLOAT_EQ(1, 1 + g[0]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestLogistic) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestLogistic) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -687,7 +682,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestLogistic) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestLogNormal) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestLogNormal) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -713,7 +708,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestLogNormal) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestNormal) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestNormal) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -739,7 +734,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestNormal) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestPareto) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestPareto) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -765,7 +760,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestPareto) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestPareto2) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestPareto2) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -793,7 +788,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestPareto2) {
   EXPECT_FLOAT_EQ(1, 1 + g[2]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestRayleigh) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestRayleigh) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -816,7 +811,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestRayleigh) {
   EXPECT_FLOAT_EQ(1, 1 + g[0]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestScaledInvChiSquare) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestScaledInvChiSquare) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -842,7 +837,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestScaledInvChiSquare) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestStudentT) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestStudentT) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -870,7 +865,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestStudentT) {
   EXPECT_FLOAT_EQ(1, 1 + g[2]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestUniform) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestUniform) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -894,7 +889,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestUniform) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestVonMises) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestVonMises) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
@@ -920,7 +915,7 @@ TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestVonMises) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
-TEST_F(AgradRev, StanMath_integrate_1d_impl_rev_TestWeibull) {
+TEST_F(AgradRev, integrate_1d_impl_rev_TestWeibull) {
   using stan::math::exp;
   using stan::math::integrate_1d_impl;
   using stan::math::var;
diff --git a/test/unit/math/rev/prob/categorical_logit_glm_lpmf_test.cpp b/test/unit/math/rev/prob/categorical_logit_glm_lpmf_test.cpp
index ac406b643fe..429533c1117 100644
--- a/test/unit/math/rev/prob/categorical_logit_glm_lpmf_test.cpp
+++ b/test/unit/math/rev/prob/categorical_logit_glm_lpmf_test.cpp
@@ -30,8 +30,7 @@ categorical_logit_glm_simple_lpmf(const std::vector<int>& y, const T_x& x,
   return lpmf;
 }
 
-TEST(ProbDistributionsCategoricalLogitGLM,
-     glm_matches_categorical_logit_doubles) {
+TEST_F(AgradRev, CategoricalLogitGLM_glm_matches_categorical_logit_doubles) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::MatrixXd;
diff --git a/test/unit/math/rev/prob/gamma_lccdf_test.cpp b/test/unit/math/rev/prob/gamma_lccdf_test.cpp
index 1066773e060..89347153337 100644
--- a/test/unit/math/rev/prob/gamma_lccdf_test.cpp
+++ b/test/unit/math/rev/prob/gamma_lccdf_test.cpp
@@ -1,9 +1,10 @@
 #include <stan/math/rev.hpp>
+#include <test/unit/math/rev/util.hpp>
 #include <gtest/gtest.h>
 #include <vector>
 #include <limits>
 
-TEST(ProbDistributionsGamma, lccdf_values) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_values) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -22,7 +23,7 @@ TEST(ProbDistributionsGamma, lccdf_values) {
   EXPECT_LE(lccdf_var.val(), 0.0);
 }
 
-TEST(ProbDistributionsGamma, lccdf_derivatives_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_derivatives_y) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -47,7 +48,7 @@ TEST(ProbDistributionsGamma, lccdf_derivatives_y) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lccdf_derivatives_alpha) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_derivatives_alpha) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -70,7 +71,7 @@ TEST(ProbDistributionsGamma, lccdf_derivatives_alpha) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lccdf_derivatives_beta) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_derivatives_beta) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -93,7 +94,7 @@ TEST(ProbDistributionsGamma, lccdf_derivatives_beta) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lccdf_derivatives_all_params) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_derivatives_all_params) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -123,7 +124,7 @@ TEST(ProbDistributionsGamma, lccdf_derivatives_all_params) {
   EXPECT_LT(grads[0], 0.0) << "d/dy should be negative";
 }
 
-TEST(ProbDistributionsGamma, lccdf_finite_diff_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_finite_diff_y) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -149,7 +150,7 @@ TEST(ProbDistributionsGamma, lccdf_finite_diff_y) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-5);
 }
 
-TEST(ProbDistributionsGamma, lccdf_finite_diff_alpha) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_finite_diff_alpha) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -176,7 +177,7 @@ TEST(ProbDistributionsGamma, lccdf_finite_diff_alpha) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-3);
 }
 
-TEST(ProbDistributionsGamma, lccdf_finite_diff_beta) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_finite_diff_beta) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -202,7 +203,7 @@ TEST(ProbDistributionsGamma, lccdf_finite_diff_beta) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-5);
 }
 
-TEST(ProbDistributionsGamma, lccdf_extreme_values_small) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_extreme_values_small) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -230,7 +231,7 @@ TEST(ProbDistributionsGamma, lccdf_extreme_values_small) {
   }
 }
 
-TEST(ProbDistributionsGamma, lccdf_extreme_values_large) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_extreme_values_large) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -258,7 +259,7 @@ TEST(ProbDistributionsGamma, lccdf_extreme_values_large) {
   }
 }
 
-TEST(ProbDistributionsGamma, lccdf_alpha_one_derivatives) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_alpha_one_derivatives) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -287,7 +288,7 @@ TEST(ProbDistributionsGamma, lccdf_alpha_one_derivatives) {
   EXPECT_NEAR(grads[2], -y_d, 1e-10);
 }
 
-TEST(ProbDistributionsGamma, lccdf_various_parameter_combinations) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_various_parameter_combinations) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -323,7 +324,7 @@ TEST(ProbDistributionsGamma, lccdf_various_parameter_combinations) {
   }
 }
 
-TEST(ProbDistributionsGamma, lccdf_second_derivative_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lccdf_second_derivative_y) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
@@ -345,7 +346,8 @@ TEST(ProbDistributionsGamma, lccdf_second_derivative_y) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lccdf_numerically_challenging_derivatives) {
+TEST_F(AgradRev,
+       ProbDistributionsGamma_lccdf_numerically_challenging_derivatives) {
   using stan::math::gamma_lccdf;
   using stan::math::var;
 
diff --git a/test/unit/math/rev/prob/gamma_lcdf_test.cpp b/test/unit/math/rev/prob/gamma_lcdf_test.cpp
index 218ef77a423..f3b5db8caa3 100644
--- a/test/unit/math/rev/prob/gamma_lcdf_test.cpp
+++ b/test/unit/math/rev/prob/gamma_lcdf_test.cpp
@@ -1,9 +1,10 @@
 #include <stan/math/rev.hpp>
+#include <test/unit/math/rev/util.hpp>
 #include <gtest/gtest.h>
 #include <vector>
 #include <limits>
 
-TEST(ProbDistributionsGamma, lcdf_values) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_values) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -22,7 +23,7 @@ TEST(ProbDistributionsGamma, lcdf_values) {
   EXPECT_LE(lcdf_var.val(), 0.0);  // log of probability
 }
 
-TEST(ProbDistributionsGamma, lcdf_derivatives_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_derivatives_y) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -47,7 +48,7 @@ TEST(ProbDistributionsGamma, lcdf_derivatives_y) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lcdf_derivatives_alpha) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_derivatives_alpha) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -70,7 +71,7 @@ TEST(ProbDistributionsGamma, lcdf_derivatives_alpha) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lcdf_derivatives_beta) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_derivatives_beta) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -93,7 +94,7 @@ TEST(ProbDistributionsGamma, lcdf_derivatives_beta) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lcdf_derivatives_all_params) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_derivatives_all_params) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -123,7 +124,7 @@ TEST(ProbDistributionsGamma, lcdf_derivatives_all_params) {
   EXPECT_GT(grads[0], 0.0) << "d/dy should be positive";
 }
 
-TEST(ProbDistributionsGamma, lcdf_finite_diff_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_finite_diff_y) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -149,7 +150,7 @@ TEST(ProbDistributionsGamma, lcdf_finite_diff_y) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-5);
 }
 
-TEST(ProbDistributionsGamma, lcdf_finite_diff_alpha) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_finite_diff_alpha) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -176,7 +177,7 @@ TEST(ProbDistributionsGamma, lcdf_finite_diff_alpha) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-3);
 }
 
-TEST(ProbDistributionsGamma, lcdf_finite_diff_beta) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_finite_diff_beta) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -202,7 +203,7 @@ TEST(ProbDistributionsGamma, lcdf_finite_diff_beta) {
   EXPECT_NEAR(grad_autodiff, grad_findiff, 1e-5);
 }
 
-TEST(ProbDistributionsGamma, lcdf_extreme_values_small) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_extreme_values_small) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -230,7 +231,7 @@ TEST(ProbDistributionsGamma, lcdf_extreme_values_small) {
   }
 }
 
-TEST(ProbDistributionsGamma, lcdf_extreme_values_large) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_extreme_values_large) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -258,7 +259,7 @@ TEST(ProbDistributionsGamma, lcdf_extreme_values_large) {
   }
 }
 
-TEST(ProbDistributionsGamma, lcdf_alpha_one_derivatives) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_alpha_one_derivatives) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -286,7 +287,7 @@ TEST(ProbDistributionsGamma, lcdf_alpha_one_derivatives) {
   EXPECT_NEAR(grads[0], expected_dy, 1e-10);
 }
 
-TEST(ProbDistributionsGamma, lcdf_various_parameter_combinations) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_various_parameter_combinations) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -322,7 +323,7 @@ TEST(ProbDistributionsGamma, lcdf_various_parameter_combinations) {
   }
 }
 
-TEST(ProbDistributionsGamma, lcdf_consistency_with_lccdf) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_consistency_with_lccdf) {
   using stan::math::gamma_lccdf;
   using stan::math::gamma_lcdf;
   using stan::math::var;
@@ -357,7 +358,7 @@ TEST(ProbDistributionsGamma, lcdf_consistency_with_lccdf) {
   EXPECT_LT(grads_lccdf[0], 0.0);
 }
 
-TEST(ProbDistributionsGamma, lcdf_second_derivative_y) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_second_derivative_y) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -379,7 +380,8 @@ TEST(ProbDistributionsGamma, lcdf_second_derivative_y) {
   EXPECT_TRUE(std::isfinite(grads[0]));
 }
 
-TEST(ProbDistributionsGamma, lcdf_numerically_challenging_derivatives) {
+TEST_F(AgradRev,
+       ProbDistributionsGamma_lcdf_numerically_challenging_derivatives) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
@@ -416,7 +418,7 @@ TEST(ProbDistributionsGamma, lcdf_numerically_challenging_derivatives) {
   }
 }
 
-TEST(ProbDistributionsGamma, lcdf_monotonic_derivative) {
+TEST_F(AgradRev, ProbDistributionsGamma_lcdf_monotonic_derivative) {
   using stan::math::gamma_lcdf;
   using stan::math::var;
 
diff --git a/test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp b/test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp
index fb8cf1d02d4..142de0147ca 100644
--- a/test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp
+++ b/test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp
@@ -7,8 +7,7 @@
 
 //  We check that the values of the new regression match those of one built
 //  from existing primitives.
-TEST(ProbDistributionsNegBinomial2LogGLM,
-     glm_matches_neg_binomial_2_log_doubles) {
+TEST_F(AgradRev, NegBinomial2LogGLM_glm_matches_neg_binomial_2_log_doubles) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::var;
@@ -32,8 +31,8 @@ TEST(ProbDistributionsNegBinomial2LogGLM,
 }
 //  We check that the values of the new regression match those of one built
 //  from existing primitives.
-TEST(ProbDistributionsNegBinomial2LogGLM,
-     glm_matches_neg_binomial_2_log_doubles_rand) {
+TEST_F(AgradRev,
+       NegBinomial2LogGLM_glm_matches_neg_binomial_2_log_doubles_rand) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::var;
@@ -541,8 +540,8 @@ TYPED_TEST(ProbDistributionsNegBinomial2LogGLM,
 }
 
 //  We check that the right errors are thrown.
-TEST(ProbDistributionsNegBinomial2LogGLM,
-     glm_matches_neg_binomial_2_log_error_checking) {
+TEST_F(AgradRev,
+       NegBinomial2LogGLM_glm_matches_neg_binomial_2_log_error_checking) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::var;
diff --git a/test/unit/math/rev/prob/ordered_logistic_glm_lpmf_test.cpp b/test/unit/math/rev/prob/ordered_logistic_glm_lpmf_test.cpp
index 64a7c3b4991..89d5889aedd 100644
--- a/test/unit/math/rev/prob/ordered_logistic_glm_lpmf_test.cpp
+++ b/test/unit/math/rev/prob/ordered_logistic_glm_lpmf_test.cpp
@@ -14,8 +14,7 @@ ordered_logistic_glm_simple_lpmf(const std::vector<int>& y, T_x&& x,
   return stan::math::ordered_logistic_lpmf<propto>(y, location, cuts);
 }
 
-TEST(ProbDistributionsOrderedLogisticGLM,
-     glm_matches_ordered_logistic_doubles) {
+TEST_F(AgradRev, OrderedLogisticGLM_glm_matches_ordered_logistic_doubles) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using std::vector;
@@ -33,8 +32,8 @@ TEST(ProbDistributionsOrderedLogisticGLM,
                   ordered_logistic_glm_simple_lpmf<false>(y, x, beta, cuts));
 }
 
-TEST(ProbDistributionsOrderedLogisticGLM,
-     glm_matches_ordered_logistic_doubles_broadcast_y) {
+TEST_F(AgradRev,
+       OrderedLogisticGLM_glm_matches_ordered_logistic_doubles_broadcast_y) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using std::vector;
diff --git a/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp b/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
index 3dc41c72845..64641e76981 100644
--- a/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
+++ b/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
@@ -26,8 +26,7 @@ inline auto skew_de_ccdf_test(const T1& y, const T2& mu, const T3& sigma,
   }
 }
 
-TEST(RevProbDistributionsSkewedDoubleExponential,
-     lccdf_computes_correct_gradients) {
+TEST_F(AgradRev, SkewedDoubleExponential_lccdf_computes_correct_gradients) {
   using stan::math::skew_double_exponential_lccdf;
 
   for (double ys : {-1.7, 0.2, 0.5, 0.9, 1.1, 3.2, 8.3}) {
diff --git a/test/unit/math/rev/util.hpp b/test/unit/math/rev/util.hpp
index b7a21185def..a25aff6561f 100644
--- a/test/unit/math/rev/util.hpp
+++ b/test/unit/math/rev/util.hpp
@@ -40,6 +40,7 @@ struct var_matrix_types {
   using matrix_v = stan::math::test::cond_var_matrix_t<UseVarMat>;
   using row_vector_v = stan::math::test::cond_var_row_vector_t<UseVarMat>;
   using vector_v = stan::math::test::cond_var_vector_t<UseVarMat>;
+  virtual ~var_matrix_types() { stan::math::recover_memory(); }
 };
 
 template <class T>
diff --git a/test/unit/multiple_translation_units1.cpp b/test/unit/multiple_translation_units1.cpp
index c1f4f92125b..3e9a69ef42d 100644
--- a/test/unit/multiple_translation_units1.cpp
+++ b/test/unit/multiple_translation_units1.cpp
@@ -1,3 +1,3 @@
-#include <stan/math.hpp>
+#include <test/unit/multiple_translation_units_includes.hpp>
 
 inline stan::math::var function1() { return 0; }
diff --git a/test/unit/multiple_translation_units2.cpp b/test/unit/multiple_translation_units2.cpp
index cb11d730165..8958b90aaa0 100644
--- a/test/unit/multiple_translation_units2.cpp
+++ b/test/unit/multiple_translation_units2.cpp
@@ -1,3 +1,3 @@
-#include <stan/math.hpp>
+#include <test/unit/multiple_translation_units_includes.hpp>
 
 inline stan::math::var function2() { return 0; }
diff --git a/test/unit/multiple_translation_units_includes.hpp b/test/unit/multiple_translation_units_includes.hpp
new file mode 100644
index 00000000000..20b4db63e2f
--- /dev/null
+++ b/test/unit/multiple_translation_units_includes.hpp
@@ -0,0 +1,14 @@
+#ifndef TEST_UNIT_MULTIPLE_TRANSLATION_UNITS_INCLUDES_HPP
+#define TEST_UNIT_MULTIPLE_TRANSLATION_UNITS_INCLUDES_HPP
+
+// Keep this list aligned with the public Stan Math umbrella headers. Each
+// translation unit test source includes this header so linking catches ODR
+// violations from any externally linked header-defined function reachable from
+// the public entry points.
+#include <stan/math.hpp>
+#include <stan/math/fwd.hpp>
+#include <stan/math/mix.hpp>
+#include <stan/math/prim.hpp>
+#include <stan/math/rev.hpp>
+
+#endif
diff --git a/test/unit/multiple_translation_units_test.cpp b/test/unit/multiple_translation_units_test.cpp
index 787c4031333..8894851aba2 100644
--- a/test/unit/multiple_translation_units_test.cpp
+++ b/test/unit/multiple_translation_units_test.cpp
@@ -2,5 +2,5 @@
 
 TEST(multiple_translation_units, compile) {
   SUCCEED() << "this test compiling indicates that compiling the math library "
-            << "with multiple translation units is ok.";
+            << "public umbrella headers with multiple translation units is ok.";
 }