libMesh
newmark_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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  UniquePtr<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  UniquePtr<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 
138  _system.project_vector( old_solution_accel, f, g );
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 bool side_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on side of elem to elem_residual.
Definition: diff_physics.h:190
Real _gamma
The value for to employ.
void swap(DenseMatrix< T > &other_matrix)
STL-like swap method.
Definition: dense_matrix.h:746
const DenseVector< Number > & get_elem_solution() const
Accessor for element solution.
Definition: diff_context.h:110
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseVector< T3 > &vec)
Adds factor times vec to this vector.
Definition: dense_vector.h:431
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
bool _initial_accel_set
This method requires an initial acceleration.
virtual Real error_order() const libmesh_override
Error convergence order: 2 for , 1 otherwise.
virtual bool side_residual(bool request_jacobian, DiffContext &) libmesh_override
This method uses the DifferentiablePhysics&#39; side_time_derivative() and side_constraint() to build a f...
virtual bool side_damping_residual(bool request_jacobian, DiffContext &)
Subtracts a damping vector contribution on side of elem from elem_residual.
Definition: diff_physics.h:390
const DenseVector< Number > & get_elem_solution_accel() const
Accessor for element solution accel of change w.r.t.
Definition: diff_context.h:179
virtual bool nonlocal_residual(bool request_jacobian, DiffContext &) libmesh_override
This method uses the DifferentiablePhysics&#39; nonlocal_time_derivative() and nonlocal_constraint() to b...
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const DenseVector< Number > & get_elem_fixed_solution() const
Accessor for element fixed solution.
Definition: diff_context.h:214
void project_initial_accel(FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr)
Specify non-zero initial acceleration.
virtual bool nonlocal_time_derivative(bool request_jacobian, DiffContext &)
Adds any nonlocal time derivative contributions (e.g.
Definition: diff_physics.h:208
virtual bool damping_residual(bool request_jacobian, DiffContext &)
Subtracts a damping vector contribution on elem from elem_residual.
Definition: diff_physics.h:373
virtual void solve() libmesh_override
This method solves for the solution at the next timestep.
The libMesh namespace provides an interface to certain functionality in the library.
void(DiffContext::* ReinitFuncType)(Real)
Definition: time_solver.h:272
Number old_solution_accel(const dof_id_type global_dof_number) const
virtual void advance_timestep() libmesh_override
This method advances the solution to the next timestep, after a solve() has been performed.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual bool nonlocal_constraint(bool request_jacobian, DiffContext &)
Adds any nonlocal constraint contributions (e.g.
Definition: diff_physics.h:226
virtual void elem_reinit(Real)
Gives derived classes the opportunity to reinitialize data (FE objects in FEMSystem, for example) needed for an interior integration at a new point within a timestep.
Definition: diff_context.h:75
bool _eulerian_time_deriv(bool request_jacobian, DiffContext &)
This method simply combines element_time_derivative() and eulerian_residual(), which makes its addres...
Definition: diff_physics.C:102
virtual void compute_initial_accel()
This method uses the specified initial displacement and velocity to compute the initial acceleration ...
virtual ~NewmarkSolver()
Destructor.
virtual bool mass_residual(bool request_jacobian, DiffContext &)
Subtracts a mass vector contribution on elem from elem_residual.
Definition: diff_physics.h:317
NewmarkSolver(sys_type &s)
Constructor.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class provides a specific system class.
Definition: diff_system.h:53
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
sys_type & _system
A reference to the system we are solving.
Definition: time_solver.h:257
virtual bool _general_residual(bool request_jacobian, DiffContext &, ResFuncType mass, ResFuncType damping, ResFuncType time_deriv, ResFuncType constraint, ReinitFuncType reinit)
This method is the underlying implementation of the public residual methods.
Number old_nonlinear_solution(const dof_id_type global_dof_number) const
Real elem_solution_rate_derivative
The derivative of elem_solution_rate with respect to the current nonlinear solution, for use by systems with non default mass_residual terms.
Definition: diff_context.h:490
const DofMap & get_dof_map() const
Definition: system.h:2030
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
bool use_fixed_solution
A boolean to be set to true by systems using elem_fixed_solution, for optional use by e...
Definition: system.h:1493
virtual void adjoint_advance_timestep() libmesh_override
This method advances the adjoint solution to the previous timestep, after an adjoint_solve() has been...
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
void set_initial_accel_avail(bool initial_accel_set)
Allow the user to (re)set whether the initial acceleration is available.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
bool(DifferentiablePhysics::* ResFuncType)(bool, DiffContext &)
Definitions of argument types for use in refactoring subclasses.
Definition: time_solver.h:270
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:794
UniquePtr< NumericVector< Number > > _old_local_solution_rate
Serial vector of previous time step velocity .
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:169
Generic class from which second order UnsteadySolvers should subclass.
const DenseVector< Number > & get_elem_solution_rate() const
Accessor for element solution rate of change w.r.t.
Definition: diff_context.h:144
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
UniquePtr< NumericVector< Number > > _old_local_solution_accel
Serial vector of previous time step acceleration .
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void elem_side_reinit(Real)
Gives derived classes the opportunity to reinitialize data needed for a side integration at a new poi...
Definition: diff_context.h:81
virtual bool nonlocal_damping_residual(bool request_jacobian, DiffContext &)
Subtracts any nonlocal damping vector contributions (e.g.
Definition: diff_physics.h:405
virtual bool element_residual(bool request_jacobian, DiffContext &) libmesh_override
This method uses the DifferentiablePhysics&#39; element_time_derivative() and element_constraint() to bui...
bool _is_accel_solve
Need to be able to indicate to _general_residual if we are doing an acceleration solve or not...
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=libmesh_nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
UniquePtr< DiffSolver > _diff_solver
An implicit linear or nonlinear solver to use at each timestep.
Definition: time_solver.h:247
Number old_solution_rate(const dof_id_type global_dof_number) const
virtual bool side_mass_residual(bool request_jacobian, DiffContext &)
Subtracts a mass vector contribution on side of elem from elem_residual.
Definition: diff_physics.h:334
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
virtual bool side_time_derivative(bool request_jacobian, DiffContext &)
Adds the time derivative contribution on side of elem to elem_residual.
Definition: diff_physics.h:170
virtual void advance_timestep() libmesh_override
This method advances the solution to the next timestep, after a solve() has been performed.
virtual bool element_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on elem to elem_residual.
Definition: diff_physics.h:141
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:394
Real elem_solution_accel_derivative
The derivative of elem_solution_accel with respect to the current nonlinear solution, for use by systems with non default mass_residual terms.
Definition: diff_context.h:497
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
virtual bool nonlocal_mass_residual(bool request_jacobian, DiffContext &c)
Subtracts any nonlocal mass vector contributions (e.g.
Definition: diff_physics.C:63
virtual void solve() libmesh_override
This method solves for the solution at the next timestep.
Real _beta
The value for the to employ.
virtual UniquePtr< NumericVector< T > > clone() const =0
bool first_solve
A bool that will be true the first time solve() is called, and false thereafter.
virtual void nonlocal_reinit(Real)
Gives derived classes the opportunity to reinitialize data needed for nonlocal calculations at a new ...
Definition: diff_context.h:93