1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 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 20 #ifndef LIBMESH_ADAPTIVE_TIME_SOLVER_H 21 #define LIBMESH_ADAPTIVE_TIME_SOLVER_H 22 23 // Local includes 24 #include "libmesh/system_norm.h" 25 #include "libmesh/first_order_unsteady_solver.h" 26 27 // C++ includes 28 29 namespace libMesh 30 { 31 32 // Forward declarations 33 class System; 34 35 /** 36 * This class wraps another UnsteadySolver derived class, and compares 37 * the results of timestepping with deltat and timestepping with 38 * 2*deltat to adjust future timestep lengths. 39 * 40 * Currently this class only works on fully coupled Systems 41 * 42 * This class is part of the new DifferentiableSystem framework, 43 * which is still experimental. Users of this framework should 44 * beware of bugs and future API changes. 45 * 46 * \author Roy H. Stogner 47 * \date 2007 48 */ 49 class AdaptiveTimeSolver : public FirstOrderUnsteadySolver 50 { 51 public: 52 /** 53 * The parent class 54 */ 55 typedef FirstOrderUnsteadySolver Parent; 56 57 /** 58 * Constructor. Requires a reference to the system 59 * to be solved. 60 */ 61 explicit 62 AdaptiveTimeSolver (sys_type & s); 63 64 /** 65 * Destructor. 66 */ 67 virtual ~AdaptiveTimeSolver (); 68 69 virtual void init() override; 70 71 virtual void reinit() override; 72 73 virtual void solve() override = 0; 74 75 virtual std::pair<unsigned int, Real> adjoint_solve (const QoISet & qoi_indices) override = 0; 76 77 virtual void advance_timestep() override; 78 79 virtual void adjoint_advance_timestep() override; 80 81 virtual void retrieve_timestep() override; 82 83 virtual void integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities) override = 0; 84 last_complete_deltat()85 virtual Real last_complete_deltat() override { return completedtimestep_deltat; } 86 87 /** 88 * This method is passed on to the core_time_solver 89 */ 90 virtual Real error_order () const override; 91 92 /** 93 * This method is passed on to the core_time_solver 94 */ 95 virtual bool element_residual (bool get_jacobian, 96 DiffContext &) override; 97 98 /** 99 * This method is passed on to the core_time_solver 100 */ 101 virtual bool side_residual (bool get_jacobian, 102 DiffContext &) override; 103 104 /** 105 * This method is passed on to the core_time_solver 106 */ 107 virtual bool nonlocal_residual (bool get_jacobian, 108 DiffContext &) override; 109 110 /** 111 * An implicit linear or nonlinear solver to use at each timestep. 112 */ 113 virtual std::unique_ptr<DiffSolver> & diff_solver() override; 114 115 /** 116 * An implicit linear solver to use for adjoint and sensitivity 117 * problems. 118 */ 119 virtual std::unique_ptr<LinearSolver<Number>> & linear_solver() override; 120 121 /** 122 * This object is used to take timesteps 123 */ 124 std::unique_ptr<UnsteadySolver> core_time_solver; 125 126 /** 127 * Error calculations are done in this norm, DISCRETE_L2 by default. 128 */ 129 SystemNorm component_norm; 130 131 /** 132 * If component_norms is non-empty, each variable's contribution to the error 133 * of a system will also be scaled by component_scale[var], unless 134 * component_scale is empty in which case all variables will be weighted 135 * equally. 136 */ 137 std::vector<float> component_scale; 138 139 /** 140 * This tolerance is the target relative error between an exact 141 * time integration and a single time step output, scaled by deltat. 142 * integrator, scaled by deltat. If the estimated error exceeds or 143 * undershoots the target error tolerance, future timesteps will be 144 * run with deltat shrunk or grown to compensate. 145 * 146 * The default value is 1.0e-2; obviously users should select their 147 * own tolerance. 148 * 149 * If a *negative* target_tolerance is specified, then its absolute 150 * value is used to scale the estimated error from the first 151 * simulation time step, and this becomes the target tolerance of 152 * all future time steps. 153 */ 154 Real target_tolerance; 155 156 /** 157 * This tolerance is the maximum relative error between an exact 158 * time integration and a single time step output, scaled by deltat. 159 * If this error tolerance is exceeded by the estimated error of the 160 * current time step, that time step will be repeated with a smaller 161 * deltat. 162 * 163 * If you use the default upper_tolerance=0.0, then the current 164 * time step will not be repeated regardless of estimated error. 165 * 166 * If a *negative* upper_tolerance is specified, then its absolute 167 * value is used to scale the estimated error from the first 168 * simulation time step, and this becomes the upper tolerance of 169 * all future time steps. 170 */ 171 Real upper_tolerance; 172 173 /** 174 * Do not allow the adaptive time solver to select deltat > max_deltat. 175 * If you use the default max_deltat=0.0, then deltat is unlimited. 176 */ 177 Real max_deltat; 178 179 /** 180 * Do not allow the adaptive time solver to select deltat < min_deltat. 181 * The default value is 0.0. 182 */ 183 Real min_deltat; 184 185 /** 186 * Do not allow the adaptive time solver to select 187 * a new deltat greater than max_growth times the old deltat. 188 * If you use the default max_growth=0.0, then the deltat growth is 189 * unlimited. 190 */ 191 Real max_growth; 192 193 /** 194 * The adaptive time solver's have two notions of deltat. The deltat the solver 195 * ended up using for the completed timestep. And the deltat the solver determined 196 * would be workable for the coming timestep. The latter gets set as system.deltat. 197 * We need a variable to save the deltat used for the completed timesteo. 198 * */ 199 Real completedtimestep_deltat; 200 201 /** 202 * This flag, which is true by default, grows (shrinks) the timestep 203 * based on the expected global accuracy of the timestepping scheme. 204 * Global in this sense means the cumulative final-time accuracy of 205 * the scheme. For example, the backward Euler scheme's truncation 206 * error is locally of order 2, so that after N timesteps of size 207 * deltat, the result is first-order accurate. If you set this to 208 * false, you can grow (shrink) your timestep based on the local 209 * accuracy rather than the global accuracy of the core TimeSolver. 210 * 211 * \note By setting this value to false you may fail to achieve 212 * the predicted convergence in time of the underlying method, however 213 * it may be possible to get more fine-grained control over step sizes 214 * as well. 215 */ 216 bool global_tolerance; 217 218 protected: 219 220 /** 221 * We need to store the value of the last deltat used, so 222 * that advance_timestep() will increment the system time 223 * correctly. 224 */ 225 Real last_deltat; 226 227 /** 228 * A helper function to calculate error norms 229 */ 230 virtual Real calculate_norm(System &, NumericVector<Number> &); 231 }; 232 233 234 } // namespace libMesh 235 236 237 #endif // LIBMESH_ADAPTIVE_TIME_SOLVER_H 238