euler_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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  : UnsteadySolver(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 {
53  unsigned int n_dofs = context.get_elem_solution().size();
54 
55  // Local nonlinear solution at old timestep
56  DenseVector<Number> old_elem_solution(n_dofs);
57  for (unsigned int i=0; i != n_dofs; ++i)
58  old_elem_solution(i) =
60 
61  // Local nonlinear solution at time t_theta
62  DenseVector<Number> theta_solution(context.get_elem_solution());
63  theta_solution *= theta;
64  theta_solution.add(1. - theta, old_elem_solution);
65 
66  // Technically the elem_solution_derivative is either theta
67  // or -1.0 in this implementation, but we scale the former part
68  // ourselves
69  context.elem_solution_derivative = 1.0;
70 
71 // Technically the fixed_solution_derivative is always theta,
72 // but we're scaling a whole jacobian by theta after these first
73 // evaluations
74  context.fixed_solution_derivative = 1.0;
75 
76  // If a fixed solution is requested, we'll use theta_solution
78  context.get_elem_fixed_solution() = theta_solution;
79 
80  // Move theta_->elem_, elem_->theta_
81  context.get_elem_solution().swap(theta_solution);
82 
83  // Move the mesh into place first if necessary
84  context.elem_reinit(theta);
85 
86  // We're going to compute just the change in elem_residual
87  // (and possibly elem_jacobian), then add back the old values
88  DenseVector<Number> old_elem_residual(context.get_elem_residual());
89  DenseMatrix<Number> old_elem_jacobian;
90  if (request_jacobian)
91  {
92  old_elem_jacobian = context.get_elem_jacobian();
93  context.get_elem_jacobian().zero();
94  }
95  context.get_elem_residual().zero();
96 
97  // Get the time derivative at t_theta
98  bool jacobian_computed =
99  _system.element_time_derivative(request_jacobian, context);
100 
101  // For a moving mesh problem we may need the pseudoconvection term too
102  jacobian_computed =
103  _system.eulerian_residual(jacobian_computed, context) && jacobian_computed;
104 
105  // Scale the time-dependent residual and jacobian correctly
106  context.get_elem_residual() *= _system.deltat;
107  if (jacobian_computed)
108  context.get_elem_jacobian() *= (theta * _system.deltat);
109 
110  // The fixed_solution_derivative is always theta,
111  // and now we're done scaling jacobians
113 
114  // We evaluate mass_residual with the change in solution
115  // to get the mass matrix, reusing old_elem_solution to hold that
116  // delta_solution. We're solving dt*F(u) - du = 0, so our
117  // delta_solution is old_solution - new_solution.
118  // We're still keeping elem_solution in theta_solution for now
119  old_elem_solution -= theta_solution;
120 
121  // Move old_->elem_, theta_->old_
122  context.get_elem_solution().swap(old_elem_solution);
123 
124  // We do a trick here to avoid using a non-1
125  // elem_solution_derivative:
126  context.get_elem_jacobian() *= -1.0;
127  context.fixed_solution_derivative *= -1.0;
128  jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
129  jacobian_computed;
130  context.get_elem_jacobian() *= -1.0;
131  context.fixed_solution_derivative *= -1.0;
132 
133  // Move elem_->elem_, old_->theta_
134  context.get_elem_solution().swap(theta_solution);
135 
136  // Restore the elem position if necessary
137  context.elem_reinit(1.);
138 
139  // Add the constraint term
140  jacobian_computed = _system.element_constraint(jacobian_computed, context) &&
141  jacobian_computed;
142 
143  // Add back the old residual and jacobian
144  context.get_elem_residual() += old_elem_residual;
145  if (request_jacobian)
146  {
147  if (jacobian_computed)
148  context.get_elem_jacobian() += old_elem_jacobian;
149  else
150  context.get_elem_jacobian().swap(old_elem_jacobian);
151  }
152 
153  return jacobian_computed;
154 }
155 
156 
157 
158 bool EulerSolver::side_residual (bool request_jacobian,
159  DiffContext &context)
160 {
161  unsigned int n_dofs = context.get_elem_solution().size();
162 
163  // Local nonlinear solution at old timestep
164  DenseVector<Number> old_elem_solution(n_dofs);
165  for (unsigned int i=0; i != n_dofs; ++i)
166  old_elem_solution(i) =
168 
169  // Local nonlinear solution at time t_theta
170  DenseVector<Number> theta_solution(context.get_elem_solution());
171  theta_solution *= theta;
172  theta_solution.add(1. - theta, old_elem_solution);
173 
174  // Technically the elem_solution_derivative is either theta
175  // or 1.0 in this implementation, but we scale the former part
176  // ourselves
177  context.elem_solution_derivative = 1.0;
178 
179 // Technically the fixed_solution_derivative is always theta,
180 // but we're scaling a whole jacobian by theta after these first
181 // evaluations
182  context.fixed_solution_derivative = 1.0;
183 
184  // If a fixed solution is requested, we'll use theta_solution
186  context.get_elem_fixed_solution() = theta_solution;
187 
188  // Move theta_->elem_, elem_->theta_
189  context.get_elem_solution().swap(theta_solution);
190 
191  // Move the mesh into place first if necessary
192  context.elem_side_reinit(theta);
193 
194  // We're going to compute just the change in elem_residual
195  // (and possibly elem_jacobian), then add back the old values
196  DenseVector<Number> old_elem_residual(context.get_elem_residual());
197  DenseMatrix<Number> old_elem_jacobian;
198  if (request_jacobian)
199  {
200  old_elem_jacobian = context.get_elem_jacobian();
201  context.get_elem_jacobian().zero();
202  }
203  context.get_elem_residual().zero();
204 
205  // Get the time derivative at t_theta
206  bool jacobian_computed =
207  _system.side_time_derivative(request_jacobian, context);
208 
209  // Scale the time-dependent residual and jacobian correctly
210  context.get_elem_residual() *= _system.deltat;
211  if (jacobian_computed)
212  context.get_elem_jacobian() *= (theta * _system.deltat);
213 
214  // The fixed_solution_derivative is always theta,
215  // and now we're done scaling jacobians
217 
218  // We evaluate side_mass_residual with the change in solution
219  // to get the mass matrix, reusing old_elem_solution to hold that
220  // delta_solution. We're solving dt*F(u) - du = 0, so our
221  // delta_solution is old_solution - new_solution.
222  // We're still keeping elem_solution in theta_solution for now
223  old_elem_solution -= theta_solution;
224 
225  // Move old_->elem_, theta_->old_
226  context.get_elem_solution().swap(old_elem_solution);
227 
228  // We do a trick here to avoid using a non-1
229  // elem_solution_derivative:
230  context.get_elem_jacobian() *= -1.0;
231  jacobian_computed = _system.side_mass_residual(jacobian_computed, context) &&
232  jacobian_computed;
233  context.get_elem_jacobian() *= -1.0;
234 
235  // Move elem_->elem_, old_->theta_
236  context.get_elem_solution().swap(theta_solution);
237 
238  // Restore the elem position if necessary
239  context.elem_side_reinit(1.);
240 
241  // Add the constraint term
242  jacobian_computed = _system.side_constraint(jacobian_computed, context) &&
243  jacobian_computed;
244 
245  // Add back the old residual and jacobian
246  context.get_elem_residual() += old_elem_residual;
247  if (request_jacobian)
248  {
249  if (jacobian_computed)
250  context.get_elem_jacobian() += old_elem_jacobian;
251  else
252  context.get_elem_jacobian().swap(old_elem_jacobian);
253  }
254 
255  return jacobian_computed;
256 }
257 
258 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo