libMesh
euler_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 
19 #include "libmesh/diff_system.h"
20 #include "libmesh/euler_solver.h"
21 
22 namespace libMesh
23 {
24 
25 
26 
28  : FirstOrderUnsteadySolver(s), theta(1.)
29 {
30 }
31 
32 
33 
35 {
36 }
37 
38 
39 
41 {
42  if (theta == 0.5)
43  return 2.;
44  return 1.;
45 }
46 
47 
48 
49 
50 bool EulerSolver::element_residual (bool request_jacobian,
51  DiffContext & context)
52 {
54 
55  return this->_general_residual(request_jacobian,
56  context,
62  compute_second_order_eqns);
63 }
64 
65 
66 
67 bool EulerSolver::side_residual (bool request_jacobian,
68  DiffContext & context)
69 {
70  return this->_general_residual(request_jacobian,
71  context,
77  false);
78 }
79 
80 
81 
82 bool EulerSolver::nonlocal_residual (bool request_jacobian,
83  DiffContext & context)
84 {
86 
87  return this->_general_residual(request_jacobian,
88  context,
94  compute_second_order_eqns);
95 }
96 
97 
98 
99 bool EulerSolver::_general_residual (bool request_jacobian,
100  DiffContext & context,
101  ResFuncType mass,
102  ResFuncType damping,
103  ResFuncType time_deriv,
104  ResFuncType constraint,
105  ReinitFuncType reinit_func,
107 {
108  unsigned int n_dofs = context.get_elem_solution().size();
109 
110  // We might need to save the old jacobian in case one of our physics
111  // terms later is unable to update it analytically.
112  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
113  if (request_jacobian)
114  old_elem_jacobian.swap(context.get_elem_jacobian());
115 
116  // Local nonlinear solution at old timestep
117  DenseVector<Number> old_elem_solution(n_dofs);
118  for (unsigned int i=0; i != n_dofs; ++i)
119  old_elem_solution(i) =
121 
122  // Local time derivative of solution
123  context.get_elem_solution_rate() = context.get_elem_solution();
124  context.get_elem_solution_rate() -= old_elem_solution;
126  context.get_elem_solution_rate() *=
128 
129  // If we are asked to compute residuals for second order variables,
130  // we also populate the acceleration part so the user can use that.
131  if (compute_second_order_eqns)
132  this->prepare_accel(context);
133 
134  // Local nonlinear solution at time t_theta
135  DenseVector<Number> theta_solution(context.get_elem_solution());
136  theta_solution *= theta;
137  theta_solution.add(1. - theta, old_elem_solution);
138 
141 
142  // If a fixed solution is requested, we'll use theta_solution
144  context.get_elem_fixed_solution() = theta_solution;
145 
146  // Move theta_->elem_, elem_->theta_
147  context.get_elem_solution().swap(theta_solution);
148 
149  // Move the mesh into place first if necessary, set t = t_{\theta}
150  (context.*reinit_func)(theta);
151 
152  // Get the time derivative at t_theta
153  bool jacobian_computed =
154  (_system.get_physics()->*time_deriv)(request_jacobian, context);
155 
156  jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
157  jacobian_computed;
158 
159  // If we have second-order variables, we need to get damping terms
160  // and the velocity equations
161  if (compute_second_order_eqns)
162  {
163  jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
164  jacobian_computed;
165 
166  jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
167  jacobian_computed;
168  }
169 
170  // Restore the elem position if necessary, set t = t_{n+1}
171  (context.*reinit_func)(1);
172 
173  // Move elem_->elem_, theta_->theta_
174  context.get_elem_solution().swap(theta_solution);
175  context.elem_solution_derivative = 1;
176 
177  // Add the constraint term
178  jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
179  jacobian_computed;
180 
181  // Add back (or restore) the old jacobian
182  if (request_jacobian)
183  {
184  if (jacobian_computed)
185  context.get_elem_jacobian() += old_elem_jacobian;
186  else
187  context.get_elem_jacobian().swap(old_elem_jacobian);
188  }
189 
190  return jacobian_computed;
191 }
192 
193 
194 } // namespace libMesh
virtual bool side_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on side of elem to elem_residual.
Definition: diff_physics.h:190
virtual ~EulerSolver()
Destructor.
Definition: euler_solver.C:34
bool compute_second_order_eqns(bool compute_jacobian, DiffContext &c)
If there are second order variables, then we need to compute their residual equations and correspondi...
Real fixed_solution_derivative
The derivative of elem_fixed_solution with respect to the nonlinear solution, for use by systems cons...
Definition: diff_context.h:504
void swap(DenseMatrix< T > &other_matrix)
STL-like swap method.
Definition: dense_matrix.h:746
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...
Definition: euler_solver.C:67
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
void prepare_accel(DiffContext &context)
If there are second order variables in the system, then we also prepare the accel for those variables...
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
Real theta
The value for the theta method to employ: 1.0 corresponds to backwards Euler, 0.0 corresponds to forw...
Definition: euler_solver.h:100
virtual bool nonlocal_residual(bool request_jacobian, DiffContext &) libmesh_override
This method uses the DifferentiablePhysics&#39; nonlocal_time_derivative() and nonlocal_constraint() to b...
Definition: euler_solver.C:82
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
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 bool _general_residual(bool request_jacobian, DiffContext &, ResFuncType mass, ResFuncType damping, ResFuncType time_deriv, ResFuncType constraint, ReinitFuncType reinit, bool compute_second_order_eqns)
This method is the underlying implementation of the public residual methods.
Definition: euler_solver.C:99
The libMesh namespace provides an interface to certain functionality in the library.
void(DiffContext::* ReinitFuncType)(Real)
Definition: time_solver.h:272
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
bool have_second_order_scalar_vars() const
Check for any second order vars that are also belong to FEFamily::SCALAR.
Definition: diff_system.C:348
virtual bool nonlocal_constraint(bool request_jacobian, DiffContext &)
Adds any nonlocal constraint contributions (e.g.
Definition: diff_physics.h:226
bool have_second_order_vars() const
Definition: diff_physics.h:523
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
void swap(DenseVector< T > &other_vector)
STL-like swap method.
Definition: dense_vector.h:341
virtual bool mass_residual(bool request_jacobian, DiffContext &)
Subtracts a mass vector contribution on elem from elem_residual.
Definition: diff_physics.h:317
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
Number old_nonlinear_solution(const dof_id_type global_dof_number) const
Generic class from which first order UnsteadySolvers should subclass.
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
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
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
bool(DifferentiablePhysics::* ResFuncType)(bool, DiffContext &)
Definitions of argument types for use in refactoring subclasses.
Definition: time_solver.h:270
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:169
const DenseVector< Number > & get_elem_solution_rate() const
Accessor for element solution rate of change w.r.t.
Definition: diff_context.h:144
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 Real error_order() const libmesh_override
Error convergence order: 2 for Crank-Nicolson, 1 otherwise.
Definition: euler_solver.C:40
virtual bool element_residual(bool request_jacobian, DiffContext &) libmesh_override
This method uses the DifferentiablePhysics&#39; element_time_derivative() and element_constraint() to bui...
Definition: euler_solver.C:50
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 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 bool element_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on elem to elem_residual.
Definition: diff_physics.h:141
virtual bool nonlocal_mass_residual(bool request_jacobian, DiffContext &c)
Subtracts any nonlocal mass vector contributions (e.g.
Definition: diff_physics.C:63
EulerSolver(sys_type &s)
Constructor.
Definition: euler_solver.C:27
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