From fe32f8bf60518afa886504a3a6d0c258c7c1f434 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 09:50:26 -0400 Subject: [PATCH 1/7] Only call elem_fe_reinit in elem_reinit if there's mesh motion Otherwise, it's unnecessary since elem_fe_reinit is called before assembling each element residual contributions. --- src/systems/fem_context.C | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/systems/fem_context.C b/src/systems/fem_context.C index 0579463a64b..4885d925085 100644 --- a/src/systems/fem_context.C +++ b/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(); + } } From 45744d01def9558fc4592126eeb877ab7a2c1481 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 09:58:13 -0400 Subject: [PATCH 2/7] Add initial reinit_func(1.) call in Euler2Solver We need the call to reinit_func to set the time, t, in the context to the correct value. Also added clarifying comments in EulerSolver and Euler2Solver that we're also setting the time in addition to possibly resetting the mesh if there's mesh motion. There is probably a way to make this more efficient such that we only call reinit_func twice in Euler2Solver, but I didn't put any thought into it. --- src/solvers/euler2_solver.C | 7 +++++-- src/solvers/euler_solver.C | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/solvers/euler2_solver.C b/src/solvers/euler2_solver.C index 5fc5a27f7d9..89fb818e596 100644 --- a/src/solvers/euler2_solver.C +++ b/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 diff --git a/src/solvers/euler_solver.C b/src/solvers/euler_solver.C index e378dbf438d..f96c52c0110 100644 --- a/src/solvers/euler_solver.C +++ b/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_ From 0b40def3ac27b8566cb61f856d5b73fa4a7613db Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 10:01:44 -0400 Subject: [PATCH 3/7] Add Euler2Solver unit test for second order This test exercises the reinit_func fix. Before, we didn't have the correct time in the context so we didn't integrate this ODE exactly as we should have. Now we do. --- tests/solvers/first_order_unsteady_solver_test.C | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/solvers/first_order_unsteady_solver_test.C b/tests/solvers/first_order_unsteady_solver_test.C index 35d89b18739..b9e4f0446eb 100644 --- a/tests/solvers/first_order_unsteady_solver_test.C +++ b/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(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(0.5,10); + } }; CPPUNIT_TEST_SUITE_REGISTRATION( EulerSolverTest ); From 13f6fcbcc5b204767f3575078942da440feef5e4 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 10:25:19 -0400 Subject: [PATCH 4/7] Add NewmarkSolver setters for beta, gamma --- include/solvers/newmark_solver.h | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/include/solvers/newmark_solver.h b/include/solvers/newmark_solver.h index 006856d5a0d..db4273c4e6d 100644 --- a/include/solvers/newmark_solver.h +++ b/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: /** From 48ecc2666792a66e05e7dd1fe3de67bb6801f607 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 10:26:11 -0400 Subject: [PATCH 5/7] Add reinit_func to NewmarkSolver, call in residual We need to call reinit_func in order to set t = t_{n+1} for evaluating the residual for Newmark. Of course, we'll also need it for the mesh motion if it's there. --- include/solvers/newmark_solver.h | 3 ++- src/solvers/newmark_solver.C | 15 +++++++++++---- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/include/solvers/newmark_solver.h b/include/solvers/newmark_solver.h index db4273c4e6d..3099a19a8f2 100644 --- a/include/solvers/newmark_solver.h +++ b/include/solvers/newmark_solver.h @@ -192,7 +192,8 @@ class NewmarkSolver : public SecondOrderUnsteadySolver ResFuncType mass, ResFuncType damping, ResFuncType time_deriv, - ResFuncType constraint); + ResFuncType constraint, + ReinitFuncType reinit); }; } // namespace libMesh diff --git a/src/solvers/newmark_solver.C b/src/solvers/newmark_solver.C index 0b6a7906e93..f7e7ab4d6b4 100644 --- a/src/solvers/newmark_solver.C +++ b/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} From 0657125114f48a7e8fcfba745a8f0c4833e0f282 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 10:28:11 -0400 Subject: [PATCH 6/7] Rework NewmarkSolverTest a bit This is to prepare for adding a new test where we reset beta. --- .../second_order_unsteady_solver_test.C | 26 +++++++++++++++---- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/tests/solvers/second_order_unsteady_solver_test.C b/tests/solvers/second_order_unsteady_solver_test.C index 374968a7359..4fbea74e90e 100644 --- a/tests/solvers/second_order_unsteady_solver_test.C +++ b/tests/solvers/second_order_unsteady_solver_test.C @@ -49,8 +49,28 @@ public: { return 2.71/3.14*0.5*t*t; } }; +class NewmarkSolverTestBase : public TimeSolverTestImplementation +{ +public: + NewmarkSolverTestBase() + : TimeSolverTestImplementation(), + _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 + public NewmarkSolverTestBase { public: CPPUNIT_TEST_SUITE( NewmarkSolverTest ); @@ -66,10 +86,6 @@ public: this->run_test_with_exact_soln(0.5,10); } -protected: - - virtual void aux_time_solver_init( NewmarkSolver & time_solver ) - { time_solver.compute_initial_accel(); } }; From 2dd3138add395e34791227b7c572cfd6cc8d8930 Mon Sep 17 00:00:00 2001 From: "Paul T. Bauman" Date: Wed, 24 Aug 2016 10:28:39 -0400 Subject: [PATCH 7/7] Add new NewmarkSolver unit test For this one, we should be able to integrate exactly a linear function (in time) of acceleration. This required both the beta setter (since you need beta = 1/6 for this case) and the reinit_func fix in the NewmarkSolver. --- .../second_order_unsteady_solver_test.C | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/tests/solvers/second_order_unsteady_solver_test.C b/tests/solvers/second_order_unsteady_solver_test.C index 4fbea74e90e..9a5a11658fa 100644 --- a/tests/solvers/second_order_unsteady_solver_test.C +++ b/tests/solvers/second_order_unsteady_solver_test.C @@ -49,6 +49,29 @@ 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 { public: @@ -76,6 +99,7 @@ public: CPPUNIT_TEST_SUITE( NewmarkSolverTest ); CPPUNIT_TEST( testNewmarkSolverConstantSecondOrderODE ); + CPPUNIT_TEST( testNewmarkSolverLinearTimeSecondOrderODE ); CPPUNIT_TEST_SUITE_END(); @@ -86,6 +110,14 @@ public: this->run_test_with_exact_soln(0.5,10); } + 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(0.5,10); + } };