diff --git a/autotest/out_baseline.dat b/autotest/out_baseline.dat
index 0928712..8cdebe5 100644
--- a/autotest/out_baseline.dat
+++ b/autotest/out_baseline.dat
@@ -185,19 +185,19 @@ Final mass u:  0.8143001155
 Max value u:   0.9995762608
 
 --- Product remap 2D (FCT)
-mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps
-Final mass us: 0.1796082878
-Mass loss us:  4.15262e-07
+mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13
+Final mass us: 0.1767475452
+Mass loss us:  0.00286033
 
 --- Product remap 2D (ClipScale)
-mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 1
-Final mass us: 0.1814773886
-Mass loss us:  0.00186952
+mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13
+Final mass us: 0.1782170448
+Mass loss us:  0.00139083
 
 --- Product remap 2D (FCTProject)
-mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 1
-Final mass us: 0.1814920761
-Mass loss us:  0.0018842
+mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13
+Final mass us: 0.1776386793
+Mass loss us:  0.00196919
 
 --- BLAST sharpening test - Pacman remap auto-dt
 mpirun -np 2 ./remhos -no-vis --verify-bounds -m ./data/inline-quad.mesh -p 14 -rs 1 -dt -1 -tf 0.75 -ho 3 -lo 5 -fct 4 -bt 1 -dtc 1
diff --git a/autotest/test.sh b/autotest/test.sh
index 70cd450..79029fc 100755
--- a/autotest/test.sh
+++ b/autotest/test.sh
@@ -66,17 +66,17 @@ for method in "${methods[@]}"; do
 done
 
 echo -e '\n'"--- Product remap 2D (FCT)" >> $file
-run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps"
+run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 1 -ps -s 13"
 echo -e $run_line >> $file
 $run_line | grep -e 'mass us' -e 'loss us'>> $file
 
 echo -e '\n'"--- Product remap 2D (ClipScale)" >> $file
-run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 1"
+run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 2 -ps -s 13"
 echo -e $run_line >> $file
 $run_line | grep -e 'mass us' -e 'loss us'>> $file
 
 echo -e '\n'"--- Product remap 2D (FCTProject)" >> $file
-run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 1"
+run_line=$command" -m ./data/inline-quad.mesh -p 14 -rs 2 -dt 0.005 -tf 0.75 -ho 3 -lo 1 -fct 4 -ps -s 13"
 echo -e $run_line >> $file
 $run_line | grep -e 'mass us' -e 'loss us'>> $file
 
diff --git a/remhos.cpp b/remhos.cpp
index 93ea0f3..054ae73 100644
--- a/remhos.cpp
+++ b/remhos.cpp
@@ -47,16 +47,20 @@
 using namespace std;
 using namespace mfem;
 
-enum class HOSolverType {None, Neumann, CG, LocalInverse};
-enum class FCTSolverType {None, FluxBased, ClipScale,
-                          NonlinearPenalty, FCTProject
-                         };
-enum class LOSolverType {None,    DiscrUpwind,    DiscrUpwindPrec,
-                         ResDist, ResDistSubcell, MassBased
-                        };
-enum class MonolithicSolverType {None, ResDistMono, ResDistMonoSubcell};
+enum class HOSolverType
+{ None, Neumann, CG, LocalInverse };
 
-enum class TimeStepControl {FixedTimeStep, LOBoundsError};
+enum class LOSolverType
+{ None, DiscrUpwind, DiscrUpwindPrec, ResDist, ResDistSubcell, MassBased };
+
+enum class FCTSolverType
+{ None, FluxBased, ClipScale, NonlinearPenalty, FCTProject };
+
+enum class MonolithicSolverType
+{ None, ResDistMono, ResDistMonoSubcell };
+
+enum class TimeStepControl
+{ FixedTimeStep, LOBoundsError };
 
 // Choice for the problem setup. The fluid velocity, initial condition and
 // inflow boundary condition are chosen based on this parameter.
@@ -306,18 +310,25 @@ int main(int argc, char *argv[])
 
    // Define the ODE solver used for time integration. Several explicit
    // Runge-Kutta methods are available.
-   IDPODESolver *ode_solver = NULL;
+   IDPODESolver *idp_ode_solver = nullptr;
+   ODESolver *ode_solver = nullptr;
    switch (ode_solver_type)
    {
-      case 1: ode_solver = new ForwardEulerIDPSolver(); break;
-      case 2: ode_solver = new RK2IDPSolver(); break;
-      case 3: ode_solver = new RK3IDPSolver(); break;
-      case 4: ode_solver = new RK4IDPSolver(); break;
-      case 6: ode_solver = new RK6IDPSolver(); break;
+      case 1: ode_solver = new ForwardEulerSolver(); break;
+      case 2: ode_solver = new RK2Solver(1.0); break;
+      case 3: ode_solver = new RK3SSPSolver(); break;
+      case 4: ode_solver = new RK4Solver(); break;
+      case 6: ode_solver = new RK6Solver(); break;
+      case 11: idp_ode_solver = new ForwardEulerIDPSolver(); break;
+      case 12: idp_ode_solver = new RK2IDPSolver(); break;
+      case 13: idp_ode_solver = new RK3IDPSolver(); break;
+      case 14: idp_ode_solver = new RK4IDPSolver(); break;
+      case 16: idp_ode_solver = new RK6IDPSolver(); break;
       default:
          cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
          return 3;
    }
+   if (ode_solver_type > 10) { ode_solver = idp_ode_solver; }
 
    // Check if the input mesh is periodic.
    const bool periodic = pmesh.GetNodes() != NULL &&
@@ -1452,6 +1463,8 @@ void AdvectionOperator::ComputeMask(const Vector &x, Array<bool> &mask) const
    {
       for (int d = 0; d < ndim; d++)
       {
+         // Note that this puts a mask on the first field as well, meaning
+         // that propagation in new elements will be through Forward Euler.
          mask[i+ndofs*d] = bool_dofs[i];
       }
    }
diff --git a/remhos_solvers.cpp b/remhos_solvers.cpp
index 9c9f4cf..43a40d1 100644
--- a/remhos_solvers.cpp
+++ b/remhos_solvers.cpp
@@ -212,13 +212,19 @@ void RKIDPSolver::Step(Vector &x, double &t, double &dt)
       // Update mask with the HO update
       UpdateMask(x, dxs[i], dct, mask);
 
-      // Convert the HO update to Forward Euler
+      // Form the unlimited update for the stage.
+      // Note that it converts eq. (2.16) in JLG's paper into an update using
+      // the previous limited updates.
       if (d_i[i] != 1.)
       {
+         // for mask = 0, we get dxs (nothing happens).
+         //               the loop below won't change it -> Forward Euler.
+         // for mask = 1, we scale dxs by d_i[i].
          AddMasked(mask, d_i[i]-1., dxs[i], dxs[i]);
       }
       for (int j = 0; j < i; j++)
       {
+         // Use all previous limited updates.
          AddMasked(mask, d_i[j], dxs[j], dxs[i]);
       }
 
diff --git a/remhos_solvers.hpp b/remhos_solvers.hpp
index b37a96f..d3a983d 100644
--- a/remhos_solvers.hpp
+++ b/remhos_solvers.hpp
@@ -69,7 +69,12 @@ class IDPODESolver : public ODESolver
    LimitedTimeDependentOperator *f;  // f(.,t) : R^n --> R^n
 
    void Init(TimeDependentOperator &f_) override
-   { MFEM_ABORT("Limited time-dependent operator must be assigned!"); }
+   {
+      auto lo = dynamic_cast<LimitedTimeDependentOperator *>(&f_);
+      if (lo) { Init(*lo); }
+      else    { MFEM_ABORT("LimitedTimeDependentOperator must be assigned!"); }
+   }
+
 public:
    IDPODESolver() : ODESolver(), f(NULL) { }
    virtual ~IDPODESolver() { }
@@ -94,6 +99,9 @@ class RKIDPSolver : public IDPODESolver
    Vector *dxs;
    Array<bool> mask;
 
+   // This function constructs coefficients that transform eq. (2.16) from
+   // JLG's paper to an update that only uses the previous limited updates.
+   // This function does not depend on the Operator f in any way.
    void ConstructD();
 
    /// Adds only DOFs that have mask = true.