twostep_time_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 
20 #include "libmesh/diff_system.h"
21 #include "libmesh/euler_solver.h"
22 #include "libmesh/numeric_vector.h"
23 
24 namespace libMesh
25 {
26 
27 
28 
31 
32 {
33  // We start with a reasonable time solver: implicit Euler
34  core_time_solver.reset(new EulerSolver(s));
35 }
36 
37 
38 
40 {
41 }
42 
43 
44 
46 {
47  libmesh_assert(core_time_solver.get());
48 
49  // The core_time_solver will handle any first_solve actions
50  first_solve = false;
51 
52  // We may have to repeat timesteps entirely if our error is bad
53  // enough
54  bool max_tolerance_met = false;
55 
56  // Calculating error values each time
57  Real single_norm(0.), double_norm(0.), error_norm(0.),
58  relative_error(0.);
59 
60  while (!max_tolerance_met)
61  {
62  // If we've been asked to reduce deltat if necessary, make sure
63  // the core timesolver does so
64  core_time_solver->reduce_deltat_on_diffsolver_failure =
66 
67  if (!quiet)
68  {
69  libMesh::out << "\n === Computing adaptive timestep === "
70  << std::endl;
71  }
72 
73  // Use the double-length timestep first (so the
74  // old_nonlinear_solution won't have to change)
75  core_time_solver->solve();
76 
77  // Save a copy of the double-length nonlinear solution
78  // and the old nonlinear solution
79  std::unique_ptr<NumericVector<Number>> double_solution =
80  _system.solution->clone();
81  std::unique_ptr<NumericVector<Number>> old_solution =
82  _system.get_vector("_old_nonlinear_solution").clone();
83 
84  double_norm = calculate_norm(_system, *double_solution);
85  if (!quiet)
86  {
87  libMesh::out << "Double norm = " << double_norm << std::endl;
88  }
89 
90  // Then reset the initial guess for our single-length calcs
91  *(_system.solution) = _system.get_vector("_old_nonlinear_solution");
92 
93  // Call two single-length timesteps
94  // Be sure that the core_time_solver does not change the
95  // timestep here. (This is unlikely because it just succeeded
96  // with a timestep twice as large!)
97  // FIXME: even if diffsolver failure is unlikely, we ought to
98  // do *something* if it happens
99  core_time_solver->reduce_deltat_on_diffsolver_failure = 0;
100 
101  Real old_time = _system.time;
102  Real old_deltat = _system.deltat;
103  _system.deltat *= 0.5;
104  core_time_solver->solve();
105  core_time_solver->advance_timestep();
106  core_time_solver->solve();
107 
108  single_norm = calculate_norm(_system, *_system.solution);
109  if (!quiet)
110  {
111  libMesh::out << "Single norm = " << single_norm << std::endl;
112  }
113 
114  // Reset the core_time_solver's reduce_deltat... value.
115  core_time_solver->reduce_deltat_on_diffsolver_failure =
117 
118  // But then back off just in case our advance_timestep() isn't
119  // called.
120  // FIXME: this probably doesn't work with multistep methods
121  _system.get_vector("_old_nonlinear_solution") = *old_solution;
122  _system.time = old_time;
123  _system.deltat = old_deltat;
124 
125  // Find the relative error
126  *double_solution -= *(_system.solution);
127  error_norm = calculate_norm(_system, *double_solution);
128  relative_error = error_norm / _system.deltat /
129  std::max(double_norm, single_norm);
130 
131  // If the relative error makes no sense, we're done
132  if (!double_norm && !single_norm)
133  return;
134 
135  if (!quiet)
136  {
137  libMesh::out << "Error norm = " << error_norm << std::endl;
138  libMesh::out << "Local relative error = "
139  << (error_norm /
140  std::max(double_norm, single_norm))
141  << std::endl;
142  libMesh::out << "Global relative error = "
143  << (error_norm / _system.deltat /
144  std::max(double_norm, single_norm))
145  << std::endl;
146  libMesh::out << "old delta t = " << _system.deltat << std::endl;
147  }
148 
149  // If our upper tolerance is negative, that means we want to set
150  // it based on the first successful time step
151  if (this->upper_tolerance < 0)
152  this->upper_tolerance = -this->upper_tolerance * relative_error;
153 
154  // If we haven't met our upper error tolerance, we'll have to
155  // repeat this timestep entirely
156  if (this->upper_tolerance && relative_error > this->upper_tolerance)
157  {
158  // Reset the initial guess for our next try
159  *(_system.solution) =
160  _system.get_vector("_old_nonlinear_solution");
161 
162  // Chop delta t in half
163  _system.deltat /= 2.;
164 
165  if (!quiet)
166  {
167  libMesh::out << "Failed to meet upper error tolerance"
168  << std::endl;
169  libMesh::out << "Retrying with delta t = "
170  << _system.deltat << std::endl;
171  }
172  }
173  else
174  max_tolerance_met = true;
175  }
176 
177 
178  // Otherwise, compare the relative error to the tolerance
179  // and adjust deltat
181 
182  // If our target tolerance is negative, that means we want to set
183  // it based on the first successful time step
184  if (this->target_tolerance < 0)
185  this->target_tolerance = -this->target_tolerance * relative_error;
186 
187  const Real global_shrink_or_growth_factor =
188  std::pow(this->target_tolerance / relative_error,
189  static_cast<Real>(1. / core_time_solver->error_order()));
190 
191  const Real local_shrink_or_growth_factor =
192  std::pow(this->target_tolerance /
193  (error_norm/std::max(double_norm, single_norm)),
194  static_cast<Real>(1. / (core_time_solver->error_order()+1.)));
195 
196  if (!quiet)
197  {
198  libMesh::out << "The global growth/shrink factor is: "
199  << global_shrink_or_growth_factor << std::endl;
200  libMesh::out << "The local growth/shrink factor is: "
201  << local_shrink_or_growth_factor << std::endl;
202  }
203 
204  // The local s.o.g. factor is based on the expected **local**
205  // truncation error for the timestepping method, the global
206  // s.o.g. factor is based on the method's **global** truncation
207  // error. You can shrink/grow the timestep to attempt to satisfy
208  // either a global or local time-discretization error tolerance.
209 
210  Real shrink_or_growth_factor =
211  this->global_tolerance ? global_shrink_or_growth_factor :
212  local_shrink_or_growth_factor;
213 
214  if (this->max_growth && this->max_growth < shrink_or_growth_factor)
215  {
216  if (!quiet && this->global_tolerance)
217  {
218  libMesh::out << "delta t is constrained by max_growth" << std::endl;
219  }
220  shrink_or_growth_factor = this->max_growth;
221  }
222 
223  _system.deltat *= shrink_or_growth_factor;
224 
225  // Restrict deltat to max-allowable value if necessary
226  if ((this->max_deltat != 0.0) && (_system.deltat > this->max_deltat))
227  {
228  if (!quiet)
229  {
230  libMesh::out << "delta t is constrained by maximum-allowable delta t."
231  << std::endl;
232  }
233  _system.deltat = this->max_deltat;
234  }
235 
236  // Restrict deltat to min-allowable value if necessary
237  if ((this->min_deltat != 0.0) && (_system.deltat < this->min_deltat))
238  {
239  if (!quiet)
240  {
241  libMesh::out << "delta t is constrained by minimum-allowable delta t."
242  << std::endl;
243  }
244  _system.deltat = this->min_deltat;
245  }
246 
247  if (!quiet)
248  {
249  libMesh::out << "new delta t = " << _system.deltat << std::endl;
250  }
251 }
252 
253 } // namespace libMesh
virtual Real calculate_norm(System &, NumericVector< Number > &)
virtual std::unique_ptr< NumericVector< T > > clone() const =0
long double max(long double a, double b)
sys_type & _system
Definition: time_solver.h:258
const NumericVector< Number > & get_vector(const std::string &vec_name) const
Definition: system.C:774
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
unsigned int reduce_deltat_on_diffsolver_failure
Definition: time_solver.h:221
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void solve() override
OStreamProxy out(std::cout)
std::unique_ptr< UnsteadySolver > core_time_solver