euler2_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 
19 #include "libmesh/diff_system.h"
20 #include "libmesh/euler2_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 Euler2Solver::element_residual (bool request_jacobian,
51  DiffContext & context)
52 {
54 
55  return this->_general_residual(request_jacobian,
56  context,
63 }
64 
65 
66 
67 bool Euler2Solver::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 Euler2Solver::nonlocal_residual (bool request_jacobian,
83  DiffContext & context)
84 {
86 
87  return this->_general_residual(request_jacobian,
88  context,
95 }
96 
97 
98 
99 bool Euler2Solver::_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,
106  bool compute_second_order_eqns)
107 {
108  unsigned int n_dofs = context.get_elem_solution().size();
109 
110  // Local nonlinear solution at old timestep
111  DenseVector<Number> old_elem_solution(n_dofs);
112  for (unsigned int i=0; i != n_dofs; ++i)
113  old_elem_solution(i) =
115 
116  // Local time derivative of solution
117  context.get_elem_solution_rate() = context.get_elem_solution();
118  context.get_elem_solution_rate() -= old_elem_solution;
120  context.get_elem_solution_rate() *=
122 
123  // Our first evaluations are at the final elem_solution
124  context.elem_solution_derivative = 1.0;
125 
126  // If a fixed solution is requested, we'll use the elem_solution
127  // at the new timestep
128  // FIXME - should this be the theta solution instead?
130  context.get_elem_fixed_solution() = context.get_elem_solution();
131 
132  context.fixed_solution_derivative = 1.0;
133 
134  // We need to save the old jacobian and old residual since we'll be
135  // multiplying some of the new contributions by theta or 1-theta
136  DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
137  DenseVector<Number> old_elem_residual(n_dofs);
138  old_elem_residual.swap(context.get_elem_residual());
139  if (request_jacobian)
140  old_elem_jacobian.swap(context.get_elem_jacobian());
141 
142  // Local time derivative of solution
143  context.get_elem_solution_rate() = context.get_elem_solution();
144  context.get_elem_solution_rate() -= old_elem_solution;
146  context.get_elem_solution_rate() *=
148 
149  // If we are asked to compute residuals for second order variables,
150  // we also populate the acceleration part so the user can use that.
152  this->prepare_accel(context);
153 
154  // Move the mesh into place first if necessary, set t = t_{n+1}
155  (context.*reinit_func)(1.);
156 
157  // First, evaluate time derivative at the new timestep.
158  // The element should already be in the proper place
159  // even for a moving mesh problem.
160  bool jacobian_computed =
161  (_system.get_physics()->*time_deriv)(request_jacobian, context);
162 
163  // Next, evaluate the mass residual at the new timestep
164 
165  jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
166  jacobian_computed;
167 
168  // If we have second-order variables, we need to get damping terms
169  // and the velocity equations
171  {
172  jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
173  jacobian_computed;
174 
175  jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
176  jacobian_computed;
177  }
178 
179  // Add the constraint term
180  jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
181  jacobian_computed;
182 
183  // The new solution's contribution is scaled by theta
184  context.get_elem_residual() *= theta;
185  context.get_elem_jacobian() *= theta;
186 
187  // Save the new solution's term
188  DenseMatrix<Number> elem_jacobian_newterm(n_dofs, n_dofs);
189  DenseVector<Number> elem_residual_newterm(n_dofs);
190  elem_residual_newterm.swap(context.get_elem_residual());
191  if (request_jacobian)
192  elem_jacobian_newterm.swap(context.get_elem_jacobian());
193 
194  // Add the time-dependent term for the old solution
195 
196  // Make sure elem_solution is set up for elem_reinit to use
197  // Move elem_->old_, old_->elem_
198  context.get_elem_solution().swap(old_elem_solution);
199  context.elem_solution_derivative = 0.0;
200 
201  // Move the mesh into place if necessary, set t = t_{n}
202  (context.*reinit_func)(0.);
203 
204  jacobian_computed =
205  (_system.get_physics()->*time_deriv)(jacobian_computed, context) &&
206  jacobian_computed;
207 
208  // Add the mass residual term for the old solution
209 
210  // Evaluating the mass residual at both old and new timesteps will be
211  // redundant in most problems but may be necessary for time accuracy
212  // or stability in moving mesh problems or problems with user-overridden
213  // mass_residual functions
214 
215  jacobian_computed =
216  (_system.get_physics()->*mass)(jacobian_computed, context) &&
217  jacobian_computed;
218 
219  // If we have second-order variables, we need to get damping terms
220  // and the velocity equations
222  {
223  jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
224  jacobian_computed;
225 
226  jacobian_computed = this->compute_second_order_eqns(jacobian_computed, context) &&
227  jacobian_computed;
228  }
229 
230  // The old solution's contribution is scaled by (1-theta)
231  context.get_elem_residual() *= (1-theta);
232  context.get_elem_jacobian() *= (1-theta);
233 
234  // Restore the elem_solution
235  // Move elem_->elem_, old_->old_
236  context.get_elem_solution().swap(old_elem_solution);
237  context.elem_solution_derivative = 1;
238 
239  // Restore the elem position if necessary, set t = t_{n+1}
240  (context.*reinit_func)(1.);
241 
242  // Add back (or restore) the old residual/jacobian
243  context.get_elem_residual() += old_elem_residual;
244  if (request_jacobian)
245  {
246  if (jacobian_computed)
247  context.get_elem_jacobian() += old_elem_jacobian;
248  else
249  context.get_elem_jacobian().swap(old_elem_jacobian);
250  }
251 
252  // Add the saved new-solution terms
253  context.get_elem_residual() += elem_residual_newterm;
254  if (jacobian_computed)
255  context.get_elem_jacobian() += elem_jacobian_newterm;
256 
257  return jacobian_computed;
258 }
259 
260 
261 
262 } // namespace libMesh
virtual unsigned int size() const override
Definition: dense_vector.h:92
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
bool compute_second_order_eqns(bool compute_jacobian, DiffContext &c)
void swap(DenseMatrix< T > &other_matrix)
Definition: dense_matrix.h:762
virtual bool _general_residual(bool request_jacobian, DiffContext &, ResFuncType mass, ResFuncType damping, ResFuncType time_deriv, ResFuncType constraint, ReinitFuncType reinit, bool compute_second_order_eqns)
Definition: euler2_solver.C:99
Euler2Solver(sys_type &s)
Definition: euler2_solver.C:27
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 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
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 &)
void swap(DenseVector< T > &other_vector)
Definition: dense_vector.h:346
virtual bool mass_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:318
sys_type & _system
Definition: time_solver.h:258
const std::set< unsigned int > & get_second_order_vars() const
Definition: diff_physics.h:530
Real elem_solution_rate_derivative
Definition: diff_context.h:507
const DenseVector< Number > & get_elem_solution() const
Definition: diff_context.h:111
virtual bool element_residual(bool request_jacobian, DiffContext &) override
Definition: euler2_solver.C:50
bool use_fixed_solution
Definition: system.h:1493
virtual bool nonlocal_residual(bool request_jacobian, DiffContext &) override
Definition: euler2_solver.C:82
const std::vector< dof_id_type > & get_dof_indices() const
Definition: diff_context.h:367
virtual bool side_residual(bool request_jacobian, DiffContext &) override
Definition: euler2_solver.C:67
bool have_second_order_scalar_vars() const
Definition: diff_system.C:335
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:182
const DenseVector< Number > & get_elem_residual() const
Definition: diff_context.h:249
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)
virtual bool side_mass_residual(bool request_jacobian, DiffContext &)
Definition: diff_physics.h:335
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
virtual Real error_order() const override
Definition: euler2_solver.C:40
virtual void nonlocal_reinit(Real)
Definition: diff_context.h:94