diff --git a/src/serac/numerics/equation_solver.cpp b/src/serac/numerics/equation_solver.cpp
index 14744cd05..b407ba27b 100644
--- a/src/serac/numerics/equation_solver.cpp
+++ b/src/serac/numerics/equation_solver.cpp
@@ -505,7 +505,7 @@ class TrustRegion : public mfem::NewtonSolver {
     settings.min_cg_iterations                          = static_cast<size_t>(nonlinear_options.min_iterations);
     settings.max_cg_iterations                          = static_cast<size_t>(linear_options.max_iterations);
     settings.cg_tol                                     = 0.5 * norm_goal;
-    double tr_size                                      = 0.1 * std::sqrt(X.Size());
+    double tr_size                                      = nonlinear_options.trust_region_scaling * std::sqrt(X.Size());
     size_t cumulative_cg_iters_from_last_precond_update = 0;
 
     auto& d  = trResults.d;   // reuse, maybe dangerous!
diff --git a/src/serac/numerics/odes.cpp b/src/serac/numerics/odes.cpp
index 2590a1bc4..55a01916e 100644
--- a/src/serac/numerics/odes.cpp
+++ b/src/serac/numerics/odes.cpp
@@ -235,11 +235,11 @@ void SecondOrderODE::Solve(const double time, const double c0, const double c1,
   dU_dt_.SetSubVectorComplement(constrained_dofs, 0.0);
   state_.du_dt += dU_dt_;
 
-  // use the previous solution as our starting guess
+  // use 0 as our starting guess for unconstrained dofs
   d2u_dt2 = state_.d2u_dt2;
   d2u_dt2.SetSubVector(constrained_dofs, 0.0);
-  d2U_dt2_.SetSubVectorComplement(constrained_dofs, 0.0);
   d2u_dt2 += d2U_dt2_;
+  d2u_dt2.SetSubVectorComplement(constrained_dofs, 0.0);
 
   solver_.solve(d2u_dt2);
   SLIC_WARNING_ROOT_IF(!solver_.nonlinearSolver().GetConverged(), "Newton Solver did not converge.");
diff --git a/src/serac/numerics/solver_config.hpp b/src/serac/numerics/solver_config.hpp
index 7aaabf759..e895b2451 100644
--- a/src/serac/numerics/solver_config.hpp
+++ b/src/serac/numerics/solver_config.hpp
@@ -393,6 +393,9 @@ struct NonlinearSolverOptions {
   /// Debug print level
   int print_level = 0;
 
+  /// Scaling for the initial trust region size
+  double trust_region_scaling = 0.1;
+
   /// Should the gradient be converted to a monolithic matrix
   bool force_monolithic = false;
 };
diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp
index 8feaed7ee..3ba8818ab 100644
--- a/src/serac/physics/solid_mechanics.hpp
+++ b/src/serac/physics/solid_mechanics.hpp
@@ -1740,7 +1740,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
       du_[j] -= displacement_(j);
     }
 
-    if (use_warm_start_ && is_quasistatic_) {
+    if (use_warm_start_) {
       // Update the linearized Jacobian matrix
       auto r = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_,
                             *parameters_[parameter_indices].state...);
diff --git a/tests b/tests
index ba34e958e..662c093dd 160000
--- a/tests
+++ b/tests
@@ -1 +1 @@
-Subproject commit ba34e958ec41f4a8d5d7197056cdb4acb38dd734
+Subproject commit 662c093dddcff184983547ea0212adce3632a751