newmark_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/diff_system.h"
19 #include "libmesh/newmark_solver.h"
20 #include "libmesh/dof_map.h"
21 #include "libmesh/diff_solver.h"
22 
23 namespace libMesh
24 {
27  _beta(0.25),
28  _gamma(0.5),
29  _is_accel_solve(false),
30  _initial_accel_set(false)
31 {}
32 
34 {}
35 
37 {
38  if (_gamma == 0.5)
39  return 2.;
40  return 1.;
41 }
42 
44 {
45  // We need to update velocity and acceleration before
46  // we update the nonlinear solution (displacement) and
47  // delta_t
48 
50  _system.get_vector("_old_solution_rate");
51 
53  _system.get_vector("_old_solution_accel");
54 
55  if (!first_solve)
56  {
57  NumericVector<Number> & old_nonlinear_soln =
58  _system.get_vector("_old_nonlinear_solution");
59 
60  NumericVector<Number> & nonlinear_solution =
61  *(_system.solution);
62 
63  // We need to cache the new solution_rate before updating the old_solution_rate
64  // so we can update acceleration with the proper old_solution_rate
65  // v_{n+1} = gamma/(beta*Delta t)*(x_{n+1}-x_n)
66  // - ((gamma/beta)-1)*v_n
67  // - (gamma/(2*beta)-1)*(Delta t)*a_n
68  std::unique_ptr<NumericVector<Number>> new_solution_rate = nonlinear_solution.clone();
69  (*new_solution_rate) -= old_nonlinear_soln;
70  (*new_solution_rate) *= (_gamma/(_beta*_system.deltat));
71  new_solution_rate->add( (1.0-_gamma/_beta), old_solution_rate );
72  new_solution_rate->add( (1.0-_gamma/(2.0*_beta))*_system.deltat, old_solution_accel );
73 
74  // a_{n+1} = (1/(beta*(Delta t)^2))*(x_{n+1}-x_n)
75  // - 1/(beta*Delta t)*v_n
76  // - (1-1/(2*beta))*a_n
77  std::unique_ptr<NumericVector<Number>> new_solution_accel = old_solution_accel.clone();
78  (*new_solution_accel) *= -(1.0/(2.0*_beta)-1.0);
79  new_solution_accel->add( -1.0/(_beta*_system.deltat), old_solution_rate );
80  new_solution_accel->add( 1.0/(_beta*_system.deltat*_system.deltat), nonlinear_solution );
81  new_solution_accel->add( -1.0/(_beta*_system.deltat*_system.deltat), old_nonlinear_soln );
82 
83  // Now update old_solution_rate
84  old_solution_rate = (*new_solution_rate);
85  old_solution_accel = (*new_solution_accel);
86  }
87 
88  // Localize updated vectors
89  old_solution_rate.localize
92 
93  old_solution_accel.localize
96 
97  // Now we can finish advancing the timestep
99 }
100 
102 {
103  libmesh_not_implemented();
104 }
105 
107 {
108  // We need to compute the initial acceleration based off of
109  // the initial position and velocity and, thus, acceleration
110  // is the unknown in diff_solver and not the displacement. So,
111  // We swap solution and acceleration. NewmarkSolver::_general_residual
112  // will check _is_accel_solve and provide the correct
113  // values to the FEMContext assuming this swap was made.
114  this->_is_accel_solve = true;
115 
116  //solution->accel, accel->solution
117  _system.solution->swap(_system.get_vector("_old_solution_accel"));
118  _system.update();
119 
120  this->_diff_solver->solve();
121 
122  // solution->solution, accel->accel
123  _system.solution->swap(_system.get_vector("_old_solution_accel"));
124  _system.update();
125 
126  // We're done, so no longer doing an acceleration solve
127  this->_is_accel_solve = false;
128 
129  this->set_initial_accel_avail(true);
130 }
131 
134 {
136  _system.get_vector("_old_solution_accel");
137 
139 
140  this->set_initial_accel_avail(true);
141 }
142 
143 void NewmarkSolver::set_initial_accel_avail( bool initial_accel_set )
144 {
145  _initial_accel_set = initial_accel_set;
146 }
147 
149 {
150  // First, check that the initial accel was set one way or another
151  if (!_initial_accel_set)
152  {
153  std::string error = "ERROR: Must first set initial acceleration using one of:\n";
154  error += "NewmarkSolver::compute_initial_accel()\n";
155  error += "NewmarkSolver::project_initial_accel()\n";
156  libmesh_error_msg(error);
157  }
158 
159  // That satisfied, now we can solve
161 }
162 
163 bool NewmarkSolver::element_residual (bool request_jacobian,
164  DiffContext & context)
165 {
166  return this->_general_residual(request_jacobian,
167  context,
173 }
174 
175 
176 
177 bool NewmarkSolver::side_residual (bool request_jacobian,
178  DiffContext & context)
179 {
180  return this->_general_residual(request_jacobian,
181  context,
187 }
188 
189 
190 
191 bool NewmarkSolver::nonlocal_residual (bool request_jacobian,
192  DiffContext & context)
193 {
194  return this->_general_residual(request_jacobian,
195  context,
201 }
202 
203 
204 
205 bool NewmarkSolver::_general_residual (bool request_jacobian,
206  DiffContext & context,
207  ResFuncType mass,
208  ResFuncType damping,
209  ResFuncType time_deriv,
210  ResFuncType constraint,
211  ReinitFuncType reinit_func)
212 {
213  unsigned int n_dofs = context.get_elem_solution().size();
214 
215  // We might need to save the old jacobian in case one of our physics
216  // terms later is unable to update it analytically.
217  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
218 
219  // Local velocity at old time step
220  DenseVector<Number> old_elem_solution_rate(n_dofs);
221  for (unsigned int i=0; i != n_dofs; ++i)
222  old_elem_solution_rate(i) =
223  old_solution_rate(context.get_dof_indices()[i]);
224 
225  // The user is computing the initial acceleration
226  // So upstream we've swapped _system.solution and _old_local_solution_accel
227  // So we need to give the context the correct entries since we're solving for
228  // acceleration here.
229  if (_is_accel_solve)
230  {
231  // System._solution is actually the acceleration right now so we need
232  // to reset the elem_solution to the right thing, which in this case
233  // is the initial guess for displacement, which got swapped into
234  // _old_solution_accel vector
235  DenseVector<Number> old_elem_solution(n_dofs);
236  for (unsigned int i=0; i != n_dofs; ++i)
237  old_elem_solution(i) =
238  old_solution_accel(context.get_dof_indices()[i]);
239 
240  context.elem_solution_derivative = 0.0;
241  context.elem_solution_rate_derivative = 0.0;
242  context.elem_solution_accel_derivative = 1.0;
243 
244  // Acceleration is currently the unknown so it's already sitting
245  // in elem_solution() thanks to FEMContext::pre_fe_reinit
246  context.get_elem_solution_accel() = context.get_elem_solution();
247 
248  // Now reset elem_solution() to what the user is expecting
249  context.get_elem_solution() = old_elem_solution;
250 
251  context.get_elem_solution_rate() = old_elem_solution_rate;
252 
253  // The user's Jacobians will be targeting derivatives w.r.t. u_{n+1}.
254  // Although the vast majority of cases will have the correct analytic
255  // Jacobians in this iteration, since we reset elem_solution_derivative*,
256  // if there are coupled/overlapping problems, there could be
257  // mismatches in the Jacobian. So we force finite differencing for
258  // the first iteration.
259  request_jacobian = false;
260  }
261  // Otherwise, the unknowns are the displacements and everything is straight
262  // forward and is what you think it is
263  else
264  {
265  if (request_jacobian)
266  old_elem_jacobian.swap(context.get_elem_jacobian());
267 
268  // Local displacement at old timestep
269  DenseVector<Number> old_elem_solution(n_dofs);
270  for (unsigned int i=0; i != n_dofs; ++i)
271  old_elem_solution(i) =
273 
274  // Local acceleration at old time step
275  DenseVector<Number> old_elem_solution_accel(n_dofs);
276  for (unsigned int i=0; i != n_dofs; ++i)
277  old_elem_solution_accel(i) =
278  old_solution_accel(context.get_dof_indices()[i]);
279 
280  // Convenience
282 
283  context.elem_solution_derivative = 1.0;
284 
285  // Local velocity at current time step
286  // v_{n+1} = gamma/(beta*Delta t)*(x_{n+1}-x_n)
287  // + (1-(gamma/beta))*v_n
288  // + (1-gamma/(2*beta))*(Delta t)*a_n
289  context.elem_solution_rate_derivative = (_gamma/(_beta*dt));
290 
291  context.get_elem_solution_rate() = context.get_elem_solution();
292  context.get_elem_solution_rate() -= old_elem_solution;
294  context.get_elem_solution_rate().add( (1.0-_gamma/_beta), old_elem_solution_rate);
295  context.get_elem_solution_rate().add( (1.0-_gamma/(2.0*_beta))*dt, old_elem_solution_accel);
296 
297 
298 
299  // Local acceleration at current time step
300  // a_{n+1} = (1/(beta*(Delta t)^2))*(x_{n+1}-x_n)
301  // - 1/(beta*Delta t)*v_n
302  // - (1/(2*beta)-1)*a_n
303  context.elem_solution_accel_derivative = 1.0/(_beta*dt*dt);
304 
305  context.get_elem_solution_accel() = context.get_elem_solution();
306  context.get_elem_solution_accel() -= old_elem_solution;
308  context.get_elem_solution_accel().add(-1.0/(_beta*dt), old_elem_solution_rate);
309  context.get_elem_solution_accel().add(-(1.0/(2.0*_beta)-1.0), old_elem_solution_accel);
310 
311  // Move the mesh into place first if necessary, set t = t_{n+1}
312  (context.*reinit_func)(1.);
313  }
314 
315  // If a fixed solution is requested, we'll use x_{n+1}
317  context.get_elem_fixed_solution() = context.get_elem_solution();
318 
319  // Get the time derivative at t_{n+1}, F(u_{n+1})
320  bool jacobian_computed = (_system.get_physics()->*time_deriv)(request_jacobian, context);
321 
322  // Damping at t_{n+1}, C(u_{n+1})
323  jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
324  jacobian_computed;
325 
326  // Mass at t_{n+1}, M(u_{n+1})
327  jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
328  jacobian_computed;
329 
330  // Add the constraint term
331  jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
332  jacobian_computed;
333 
334  // Add back (or restore) the old jacobian
335  if (request_jacobian)
336  {
337  if (jacobian_computed)
338  context.get_elem_jacobian() += old_elem_jacobian;
339  else
340  context.get_elem_jacobian().swap(old_elem_jacobian);
341  }
342 
343  return jacobian_computed;
344 }
345 
346 }
virtual unsigned int size() const override
Definition: dense_vector.h:92
virtual bool element_residual(bool request_jacobian, DiffContext &) override
virtual bool side_constraint(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:191
const DenseMatrix< Number > & get_elem_jacobian() const
Definition: diff_context.h:283
Number old_solution_accel(const dof_id_type global_dof_number) const
void swap(DenseMatrix< T > &other_matrix)
Definition: dense_matrix.h:762
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
Definition: dense_vector.h:436
std::unique_ptr< DiffSolver > _diff_solver
Definition: time_solver.h:248
const DenseVector< Number > & get_elem_fixed_solution() const
Definition: diff_context.h:215
const DenseVector< Number > & get_elem_solution_rate() const
Definition: diff_context.h:145
virtual bool side_damping_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:391
Number old_nonlinear_solution(const dof_id_type global_dof_number) const
virtual bool side_residual(bool request_jacobian, DiffContext &) override
virtual std::unique_ptr< NumericVector< T > > clone() const =0
virtual bool nonlocal_time_derivative(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:209
virtual bool damping_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:374
std::unique_ptr< NumericVector< Number > > _old_local_solution_accel
virtual bool nonlocal_constraint(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:227
virtual void elem_reinit(Real)
Definition: diff_context.h:76
bool _eulerian_time_deriv(bool request_jacobian, DiffContext &)
virtual void compute_initial_accel()
virtual bool mass_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:318
NewmarkSolver(sys_type &s)
sys_type & _system
Definition: time_solver.h:258
virtual bool _general_residual(bool request_jacobian, DiffContext &, ResFuncType mass, ResFuncType damping, ResFuncType time_deriv, ResFuncType constraint, ReinitFuncType reinit)
Real elem_solution_rate_derivative
Definition: diff_context.h:507
const DenseVector< Number > & get_elem_solution() const
Definition: diff_context.h:111
bool use_fixed_solution
Definition: system.h:1493
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:774
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
void set_initial_accel_avail(bool initial_accel_set)
const std::vector< dof_id_type > & get_dof_indices() const
Definition: diff_context.h:367
virtual void advance_timestep() override
Number old_solution_rate(const dof_id_type global_dof_number) const
virtual bool nonlocal_residual(bool request_jacobian, DiffContext &) override
virtual void solve() override
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:182
virtual void advance_timestep() override
virtual void update()
Definition: system.C:408
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void elem_side_reinit(Real)
Definition: diff_context.h:82
virtual bool nonlocal_damping_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:406
virtual bool nonlocal_mass_residual(bool request_jacobian, DiffContext &c)
const DenseVector< Number > & get_elem_solution_accel() const
Definition: diff_context.h:180
virtual Real error_order() const override
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr, int is_adjoint=-1) const
virtual void solve() override
virtual bool side_mass_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:335
std::unique_ptr< NumericVector< Number > > _old_local_solution_rate
virtual void add(const numeric_index_type i, const T value)=0
virtual bool side_time_derivative(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:171
virtual bool element_constraint(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:142
const DofMap & get_dof_map() const
Definition: system.h:2049
Real elem_solution_accel_derivative
Definition: diff_context.h:514
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:450
virtual void adjoint_advance_timestep() override
void project_initial_accel(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr)
virtual void nonlocal_reinit(Real)
Definition: diff_context.h:94