Merge pull request #1230 from LLNL/tupek/initial_acceleration

Easier initial guess for acceleration in dynamics.

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!

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.");

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;
 };

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...);