catch already reversible solvers

diffrax/_adjoint.py

@@ -26,6 +26,9 @@
     AbstractSolver,
     AbstractStratonovichSolver,
     AbstractWrappedSolver,
+    LeapfrogMidpoint,
+    ReversibleHeun,
+    SemiImplicitEuler,
 )
 from ._term import AbstractTerm, AdjointTerm
 
@@ -1126,6 +1129,12 @@ def loop(
                 "`diffrax.ReversibleAdjoint` is not compatible with events."
             )
 
+        if isinstance(solver, (SemiImplicitEuler, ReversibleHeun, LeapfrogMidpoint)):
+            raise ValueError(
+                "`diffrax.ReversibleAdjoint` is not compatible with solvers that are "
+                f"intrinsically algebraically reversible, such as {solver}."
+            )
+
         solver = _Reversible(solver, self.l)
         tprev = init_state.tprev
         tnext = init_state.tnext

test/test_reversible.py

@@ -354,3 +354,29 @@ def test_incorrect_saveat():
             stepsize_controller=diffrax.ConstantStepSize(),
             pytree_state=False,
         )
+
+
+def test_incorrect_solver():
+    y0 = (jnp.array([0.9, 5.4]), jnp.array([0.9, 5.4]))
+    args = (0.1, -1)
+    terms = (
+        diffrax.ODETerm(_VectorField(nondiff_arg=1, diff_arg=-0.1)),
+        diffrax.ODETerm(_VectorField(nondiff_arg=1, diff_arg=-0.1)),
+    )
+    incompatible_solvers = (
+        diffrax.SemiImplicitEuler(),
+        diffrax.ReversibleHeun(),
+        diffrax.LeapfrogMidpoint(),
+    )
+    for solver in incompatible_solvers:
+        with pytest.raises(ValueError):
+            diffrax.diffeqsolve(
+                terms,
+                solver,
+                t0=0,
+                t1=5,
+                dt0=0.01,
+                y0=y0,
+                args=args,
+                adjoint=diffrax.ReversibleAdjoint(),
+            )