FEMSystem TimeSolver Fixes

include/solvers/newmark_solver.h

@@ -141,6 +141,19 @@ class NewmarkSolver : public SecondOrderUnsteadySolver
                                   DiffContext &) libmesh_override;
 
 
+  /**
+   * Setter for \f$ \beta \f$
+   */
+  void set_beta ( Real beta )
+  { _beta = beta; }
+
+  /**
+   * Setter for \f$ \gamma \f$. Note method is second order
+   * only for \f$ \gamma = 1/2 \f$
+   */
+  void set_gamma ( Real gamma )
+  { _gamma = gamma; }
+
 protected:
 
   /**
@@ -179,7 +192,8 @@ class NewmarkSolver : public SecondOrderUnsteadySolver
                                   ResFuncType mass,
                                   ResFuncType damping,
                                   ResFuncType time_deriv,
-                                  ResFuncType constraint);
+                                  ResFuncType constraint,
+                                  ReinitFuncType reinit);
 };
 
 } // namespace libMesh

src/solvers/euler2_solver.C

@@ -134,6 +134,9 @@ bool Euler2Solver::_general_residual (bool request_jacobian,
   context.get_elem_solution_rate() *=
     context.elem_solution_rate_derivative;
 
+  // Move the mesh into place first if necessary, set t = t_{n+1}
+  (context.*reinit_func)(1.);
+
   // First, evaluate time derivative at the new timestep.
   // The element should already be in the proper place
   // even for a moving mesh problem.
@@ -167,7 +170,7 @@ bool Euler2Solver::_general_residual (bool request_jacobian,
   context.get_elem_solution().swap(old_elem_solution);
   context.elem_solution_derivative = 0.0;
 
-  // Move the mesh into place first if necessary
+  // Move the mesh into place if necessary, set t = t_{n}
   (context.*reinit_func)(0.);
 
   jacobian_computed =
@@ -194,7 +197,7 @@ bool Euler2Solver::_general_residual (bool request_jacobian,
   context.get_elem_solution().swap(old_elem_solution);
   context.elem_solution_derivative = 1;
 
-  // Restore the elem position if necessary
+  // Restore the elem position if necessary, set t = t_{n+1}
   (context.*reinit_func)(1.);
 
   // Add back (or restore) the old residual/jacobian

src/solvers/euler_solver.C

@@ -129,7 +129,7 @@ bool EulerSolver::_general_residual (bool request_jacobian,
   // Move theta_->elem_, elem_->theta_
   context.get_elem_solution().swap(theta_solution);
 
-  // Move the mesh into place first if necessary
+  // Move the mesh into place first if necessary, set t = t_{\theta}
   (context.*reinit_func)(theta);
 
   // Get the time derivative at t_theta
@@ -139,7 +139,7 @@ bool EulerSolver::_general_residual (bool request_jacobian,
   jacobian_computed = (_system.*mass)(jacobian_computed, context) &&
     jacobian_computed;
 
-  // Restore the elem position if necessary
+  // Restore the elem position if necessary, set t = t_{n+1}
   (context.*reinit_func)(1);
 
   // Move elem_->elem_, theta_->theta_

src/solvers/newmark_solver.C

@@ -168,7 +168,8 @@ bool NewmarkSolver::element_residual (bool request_jacobian,
                                  &DifferentiablePhysics::mass_residual,
                                  &DifferentiablePhysics::damping_residual,
                                  &DifferentiablePhysics::_eulerian_time_deriv,
-                                 &DifferentiablePhysics::element_constraint);
+                                 &DifferentiablePhysics::element_constraint,
+                                 &DiffContext::elem_reinit);
 }
 
 
@@ -181,7 +182,8 @@ bool NewmarkSolver::side_residual (bool request_jacobian,
                                  &DifferentiablePhysics::side_mass_residual,
                                  &DifferentiablePhysics::side_damping_residual,
                                  &DifferentiablePhysics::side_time_derivative,
-                                 &DifferentiablePhysics::side_constraint);
+                                 &DifferentiablePhysics::side_constraint,
+                                 &DiffContext::elem_side_reinit);
 }
 
 
@@ -194,7 +196,8 @@ bool NewmarkSolver::nonlocal_residual (bool request_jacobian,
                                  &DifferentiablePhysics::nonlocal_mass_residual,
                                  &DifferentiablePhysics::nonlocal_damping_residual,
                                  &DifferentiablePhysics::nonlocal_time_derivative,
-                                 &DifferentiablePhysics::nonlocal_constraint);
+                                 &DifferentiablePhysics::nonlocal_constraint,
+                                 &DiffContext::nonlocal_reinit);
 }
 
 
@@ -204,7 +207,8 @@ bool NewmarkSolver::_general_residual (bool request_jacobian,
                                        ResFuncType mass,
                                        ResFuncType damping,
                                        ResFuncType time_deriv,
-                                       ResFuncType constraint)
+                                       ResFuncType constraint,
+                                       ReinitFuncType reinit_func)
 {
   unsigned int n_dofs = context.get_elem_solution().size();
 
@@ -303,6 +307,9 @@ bool NewmarkSolver::_general_residual (bool request_jacobian,
       context.get_elem_solution_accel() *= context.elem_solution_accel_derivative;
       context.get_elem_solution_accel().add(-1.0/(_beta*dt), old_elem_solution_rate);
       context.get_elem_solution_accel().add(-(1.0/(2.0*_beta)-1.0), old_elem_solution_accel);
+
+      // Move the mesh into place first if necessary, set t = t_{n+1}
+      (context.*reinit_func)(1.);
     }
 
   // If a fixed solution is requested, we'll use x_{n+1}

src/systems/fem_context.C

@@ -1305,18 +1305,20 @@ void FEMContext::elem_reinit(Real theta)
 
   // Handle a moving element if necessary.
   if (_mesh_sys)
-    // We assume that the ``default'' state
-    // of the mesh is its final, theta=1.0
-    // position, so we don't bother with
-    // mesh motion in that case.
-    if (theta != 1.0)
-      {
-        // FIXME - ALE is not threadsafe yet!
-        libmesh_assert_equal_to (libMesh::n_threads(), 1);
+    {
+      // We assume that the ``default'' state
+      // of the mesh is its final, theta=1.0
+      // position, so we don't bother with
+      // mesh motion in that case.
+      if (theta != 1.0)
+        {
+          // FIXME - ALE is not threadsafe yet!
+          libmesh_assert_equal_to (libMesh::n_threads(), 1);
 
-        elem_position_set(theta);
-      }
-  elem_fe_reinit();
+          elem_position_set(theta);
+        }
+      elem_fe_reinit();
+   }
 }
 
 

tests/solvers/first_order_unsteady_solver_test.C

@@ -125,6 +125,7 @@ public:
   CPPUNIT_TEST_SUITE( Euler2SolverTest );
 
   CPPUNIT_TEST( testEuler2SolverConstantFirstOrderODE );
+  CPPUNIT_TEST( testEuler2SolverLinearTimeFirstOrderODE );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -138,6 +139,12 @@ public:
     this->run_test_with_exact_soln<ConstantFirstOrderODE>(0.5,10);
   }
 
+  void testEuler2SolverLinearTimeFirstOrderODE()
+  {
+    // Need \theta = 0.5 since this has t in F.
+    this->set_theta(0.5);
+    this->run_test_with_exact_soln<LinearTimeFirstOrderODE>(0.5,10);
+  }
 };
 
 CPPUNIT_TEST_SUITE_REGISTRATION( EulerSolverTest );

tests/solvers/second_order_unsteady_solver_test.C

@@ -49,13 +49,57 @@ public:
   { return 2.71/3.14*0.5*t*t; }
 };
 
+//! Implements ODE: 1.0\ddot{u} = 6.0*t+2.0, u(0) = 0, \dot{u}(0) = 0
+class LinearTimeSecondOrderODE : public SecondOrderScalarSystemBase
+{
+public:
+  LinearTimeSecondOrderODE(EquationSystems & es,
+                           const std::string & name_in,
+                           const unsigned int number_in)
+    : SecondOrderScalarSystemBase(es, name_in, number_in)
+  {}
+
+  virtual Number F( FEMContext & context, unsigned int /*qp*/ )
+  { return -6.0*context.get_time()-2.0; }
+
+  virtual Number C( FEMContext & /*context*/, unsigned int /*qp*/ )
+  { return 0.0; }
+
+  virtual Number M( FEMContext & /*context*/, unsigned int /*qp*/ )
+  { return 1.0; }
+
+  virtual Number u( Real t )
+  { return t*t*t+t*t; }
+};
+
+class NewmarkSolverTestBase : public TimeSolverTestImplementation<NewmarkSolver>
+{
+public:
+  NewmarkSolverTestBase()
+    : TimeSolverTestImplementation<NewmarkSolver>(),
+    _beta(0.25)
+  {}
+
+protected:
+
+  virtual void aux_time_solver_init( NewmarkSolver & time_solver )
+  { time_solver.set_beta(_beta);
+    time_solver.compute_initial_accel(); }
+
+  void set_beta( Real beta )
+  { _beta = beta; }
+
+  Real _beta;
+};
+
 class NewmarkSolverTest : public CppUnit::TestCase,
-                          public TimeSolverTestImplementation<NewmarkSolver>
+                          public NewmarkSolverTestBase
 {
 public:
   CPPUNIT_TEST_SUITE( NewmarkSolverTest );
 
   CPPUNIT_TEST( testNewmarkSolverConstantSecondOrderODE );
+  CPPUNIT_TEST( testNewmarkSolverLinearTimeSecondOrderODE );
 
   CPPUNIT_TEST_SUITE_END();
 
@@ -66,10 +110,14 @@ public:
     this->run_test_with_exact_soln<ConstantSecondOrderODE>(0.5,10);
   }
 
-protected:
-
-  virtual void aux_time_solver_init( NewmarkSolver & time_solver )
-  { time_solver.compute_initial_accel(); }
+  void testNewmarkSolverLinearTimeSecondOrderODE()
+  {
+    // For \beta = 1/6, we have the "linear acceleration method" for which
+    // we should be able to exactly integrate linear (in time) acceleration
+    // functions.
+    this->set_beta(1.0/6.0);
+    this->run_test_with_exact_soln<LinearTimeSecondOrderODE>(0.5,10);
+  }
 
 };