euler2_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/euler2_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 Euler2Solver::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  // We evaluate mass_residual with the change in solution
62  // to get the mass matrix, reusing old_elem_solution to hold that
63  // delta_solution.
64  DenseVector<Number> delta_elem_solution(context.get_elem_solution());
65  delta_elem_solution -= old_elem_solution;
66 
67  // Our first evaluations are at the true elem_solution
68  context.elem_solution_derivative = 1.0;
69 
70  // If a fixed solution is requested, we'll use the elem_solution
71  // at the new timestep
73  context.get_elem_fixed_solution() = context.get_elem_solution();
74 
75  context.fixed_solution_derivative = 1.0;
76 
77  // We're going to compute just the change in elem_residual
78  // (and possibly elem_jacobian), then add back the old values
79  DenseVector<Number> total_elem_residual(context.get_elem_residual());
80  DenseMatrix<Number> old_elem_jacobian, total_elem_jacobian;
81  if (request_jacobian)
82  {
83  old_elem_jacobian = context.get_elem_jacobian();
84  total_elem_jacobian = context.get_elem_jacobian();
85  context.get_elem_jacobian().zero();
86  }
87  context.get_elem_residual().zero();
88 
89  // First, evaluate time derivative at the new timestep.
90  // The element should already be in the proper place
91  // even for a moving mesh problem.
92  bool jacobian_computed =
93  _system.element_time_derivative(request_jacobian, context);
94 
95  // For a moving mesh problem we may need the pseudoconvection term too
96  jacobian_computed =
97  _system.eulerian_residual(jacobian_computed, context) && jacobian_computed;
98 
99  // Scale the new time-dependent residual and jacobian correctly
100  context.get_elem_residual() *= (theta * _system.deltat);
101  total_elem_residual += context.get_elem_residual();
102  context.get_elem_residual().zero();
103 
104  if (jacobian_computed)
105  {
106  context.get_elem_jacobian() *= (theta * _system.deltat);
107  total_elem_jacobian += context.get_elem_jacobian();
108  context.get_elem_jacobian().zero();
109  }
110 
111  // Next, evaluate the mass residual at the new timestep,
112  // with the delta_solution.
113  // Evaluating the mass residual at both old and new timesteps will be
114  // redundant in most problems but may be necessary for time accuracy
115  // or stability in moving mesh problems or problems with user-overridden
116  // mass_residual functions
117 
118  // Move elem_->delta_, delta_->elem_
119  context.get_elem_solution().swap(delta_elem_solution);
120 
121  jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
122  jacobian_computed;
123 
124  context.get_elem_residual() *= -theta;
125  total_elem_residual += context.get_elem_residual();
126  context.get_elem_residual().zero();
127 
128  if (jacobian_computed)
129  {
130  // The minus sign trick here is to avoid using a non-1
131  // elem_solution_derivative:
132  context.get_elem_jacobian() *= -theta;
133  total_elem_jacobian += context.get_elem_jacobian();
134  context.get_elem_jacobian().zero();
135  }
136 
137  // Add the time-dependent term for the old solution
138 
139  // Make sure elem_solution is set up for elem_reinit to use
140  // Move delta_->old_, old_->elem_
141  context.get_elem_solution().swap(old_elem_solution);
142 
143  // Move the mesh into place first if necessary
144  context.elem_reinit(0.);
145 
147  {
148  context.elem_solution_derivative = 0.0;
149  jacobian_computed =
150  _system.element_time_derivative(jacobian_computed, context) &&
151  jacobian_computed;
152  jacobian_computed =
153  _system.eulerian_residual(jacobian_computed, context) &&
154  jacobian_computed;
155  context.elem_solution_derivative = 1.0;
156  context.get_elem_residual() *= ((1. - theta) * _system.deltat);
157  total_elem_residual += context.get_elem_residual();
158  if (jacobian_computed)
159  {
160  context.get_elem_jacobian() *= ((1. - theta) * _system.deltat);
161  total_elem_jacobian += context.get_elem_jacobian();
162  context.get_elem_jacobian().zero();
163  }
164  }
165  else
166  {
167  // FIXME - we should detect if element_time_derivative() edits
168  // elem_jacobian and lies about it!
169  _system.element_time_derivative(false, context);
170  _system.eulerian_residual(false, context);
171  context.get_elem_residual() *= ((1. - theta) * _system.deltat);
172  total_elem_residual += context.get_elem_residual();
173  }
174 
175  context.get_elem_residual().zero();
176 
177  // Add the mass residual term for the old solution
178 
179  // Move old_->old_, delta_->elem_
180  context.get_elem_solution().swap(old_elem_solution);
181 
182  jacobian_computed = _system.mass_residual(jacobian_computed, context) &&
183  jacobian_computed;
184 
185  context.get_elem_residual() *= -(1. - theta);
186  total_elem_residual += context.get_elem_residual();
187  context.get_elem_residual().zero();
188 
189  if (jacobian_computed)
190  {
191  // The minus sign trick here is to avoid using a non-1
192  // *_solution_derivative:
193  context.get_elem_jacobian() *= -(1. - theta);
194  total_elem_jacobian += context.get_elem_jacobian();
195  context.get_elem_jacobian().zero();
196  }
197 
198  // Restore the elem_solution
199  // Move elem_->elem_, delta_->delta_
200  context.get_elem_solution().swap(delta_elem_solution);
201 
202  // Restore the elem position if necessary
203  context.elem_reinit(1.);
204 
205  // Add the constraint term
206  jacobian_computed = _system.element_constraint(jacobian_computed, context) &&
207  jacobian_computed;
208 
209  // Add back the previously accumulated residual and jacobian
210  context.get_elem_residual() += total_elem_residual;
211  if (request_jacobian)
212  {
213  if (jacobian_computed)
214  context.get_elem_jacobian() += total_elem_jacobian;
215  else
216  context.get_elem_jacobian().swap(old_elem_jacobian);
217  }
218 
219  return jacobian_computed;
220 }
221 
222 
223 
224 bool Euler2Solver::side_residual (bool request_jacobian,
225  DiffContext &context)
226 {
227  unsigned int n_dofs = context.get_elem_solution().size();
228 
229  // Local nonlinear solution at old timestep
230  DenseVector<Number> old_elem_solution(n_dofs);
231  for (unsigned int i=0; i != n_dofs; ++i)
232  old_elem_solution(i) =
234 
235  // We evaluate mass_residual with the change in solution
236  // to get the mass matrix, reusing old_elem_solution to hold that
237  // delta_solution.
238  DenseVector<Number> delta_elem_solution(context.get_elem_solution());
239  delta_elem_solution -= old_elem_solution;
240 
241  // Our first evaluations are at the true elem_solution
242  context.elem_solution_derivative = 1.0;
243 
244  // If a fixed solution is requested, we'll use the elem_solution
245  // at the new timestep
247  context.get_elem_fixed_solution() = context.get_elem_solution();
248 
249  context.fixed_solution_derivative = 1.0;
250 
251  // We're going to compute just the change in elem_residual
252  // (and possibly elem_jacobian), then add back the old values
253  DenseVector<Number> total_elem_residual(context.get_elem_residual());
254  DenseMatrix<Number> old_elem_jacobian, total_elem_jacobian;
255  if (request_jacobian)
256  {
257  old_elem_jacobian = context.get_elem_jacobian();
258  total_elem_jacobian = context.get_elem_jacobian();
259  context.get_elem_jacobian().zero();
260  }
261  context.get_elem_residual().zero();
262 
263  // First, evaluate time derivative at the new timestep.
264  // The element should already be in the proper place
265  // even for a moving mesh problem.
266  bool jacobian_computed =
267  _system.side_time_derivative(request_jacobian, context);
268 
269  // Scale the time-dependent residual and jacobian correctly
270  context.get_elem_residual() *= (theta * _system.deltat);
271  total_elem_residual += context.get_elem_residual();
272  context.get_elem_residual().zero();
273 
274  if (jacobian_computed)
275  {
276  context.get_elem_jacobian() *= (theta * _system.deltat);
277  total_elem_jacobian += context.get_elem_jacobian();
278  context.get_elem_jacobian().zero();
279  }
280 
281  // Next, evaluate the mass residual at the new timestep,
282  // with the delta_solution.
283  // Evaluating the mass residual at both old and new timesteps will be
284  // redundant in most problems but may be necessary for time accuracy
285  // or stability in moving mesh problems or problems with user-overridden
286  // mass_residual functions
287 
288  // Move elem_->delta_, delta_->elem_
289  context.get_elem_solution().swap(delta_elem_solution);
290 
291  jacobian_computed = _system.side_mass_residual(jacobian_computed, context) &&
292  jacobian_computed;
293 
294  context.get_elem_residual() *= -theta;
295  total_elem_residual += context.get_elem_residual();
296  context.get_elem_residual().zero();
297 
298  if (jacobian_computed)
299  {
300  // The minus sign trick here is to avoid using a non-1
301  // elem_solution_derivative:
302  context.get_elem_jacobian() *= -theta;
303  total_elem_jacobian += context.get_elem_jacobian();
304  context.get_elem_jacobian().zero();
305  }
306 
307  // Add the time-dependent term for the old solution
308 
309  // Make sure elem_solution is set up for elem_reinit to use
310  // Move delta_->old_, old_->elem_
311  context.get_elem_solution().swap(old_elem_solution);
312 
313  // Move the mesh into place first if necessary
314  context.elem_side_reinit(0.);
315 
317  {
318  context.elem_solution_derivative = 0.0;
319  jacobian_computed =
320  _system.side_time_derivative(jacobian_computed, context) &&
321  jacobian_computed;
322  context.elem_solution_derivative = 1.0;
323  context.get_elem_residual() *= ((1. - theta) * _system.deltat);
324  total_elem_residual += context.get_elem_residual();
325  if (jacobian_computed)
326  {
327  context.get_elem_jacobian() *= ((1. - theta) * _system.deltat);
328  total_elem_jacobian += context.get_elem_jacobian();
329  context.get_elem_jacobian().zero();
330  }
331  }
332  else
333  {
334  // FIXME - we should detect if side_time_derivative() edits
335  // elem_jacobian and lies about it!
336  _system.side_time_derivative(false, context);
337  context.get_elem_residual() *= ((1. - theta) * _system.deltat);
338  total_elem_residual += context.get_elem_residual();
339  }
340 
341  context.get_elem_residual().zero();
342 
343  // Add the mass residual term for the old solution
344 
345  // Move old_->old_, delta_->elem_
346  context.get_elem_solution().swap(old_elem_solution);
347 
348  jacobian_computed =
349  _system.side_mass_residual(jacobian_computed, context) &&
350  jacobian_computed;
351 
352  context.get_elem_residual() *= -(1. - theta);
353  total_elem_residual += context.get_elem_residual();
354  context.get_elem_residual().zero();
355 
356  if (jacobian_computed)
357  {
358  // The minus sign trick here is to avoid using a non-1
359  // *_solution_derivative:
360  context.get_elem_jacobian() *= -(1. - theta);
361  total_elem_jacobian += context.get_elem_jacobian();
362  context.get_elem_jacobian().zero();
363  }
364 
365  // Restore the elem_solution
366  // Move elem_->elem_, delta_->delta_
367  context.get_elem_solution().swap(delta_elem_solution);
368 
369  // Restore the elem position if necessary
370  context.elem_side_reinit(1.);
371 
372  // Add the constraint term
373  jacobian_computed = _system.side_constraint(jacobian_computed, context) &&
374  jacobian_computed;
375 
376  // Add back the previously accumulated residual and jacobian
377  context.get_elem_residual() += total_elem_residual;
378  if (request_jacobian)
379  {
380  if (jacobian_computed)
381  context.get_elem_jacobian() += total_elem_jacobian;
382  else
383  context.get_elem_jacobian().swap(old_elem_jacobian);
384  }
385 
386  return jacobian_computed;
387 }
388 
389 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo