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_TIME_SOLVER_H 21 #define LIBMESH_TIME_SOLVER_H 22 23 // Local includes 24 #include "libmesh/libmesh_common.h" 25 #include "libmesh/reference_counted_object.h" 26 27 // C++ includes 28 #include <memory> 29 30 namespace libMesh 31 { 32 33 // Forward Declarations 34 class DiffContext; 35 class DiffSolver; 36 class DifferentiablePhysics; 37 class DifferentiableSystem; 38 class ParameterVector; 39 class SensitivityData; 40 class SolutionHistory; 41 class SystemNorm; 42 class QoISet; 43 44 template <typename T> 45 class LinearSolver; 46 47 /** 48 * This is a generic class that defines a solver to handle 49 * time integration of DifferentiableSystems. 50 * 51 * A user can define a solver by deriving from this class and 52 * implementing certain functions. 53 * 54 * This class is part of the new DifferentiableSystem framework, 55 * which is still experimental. Users of this framework should 56 * beware of bugs and future API changes. 57 * 58 * \author Roy H. Stogner 59 * \date 2006 60 */ 61 class TimeSolver : public ReferenceCountedObject<TimeSolver> 62 { 63 public: 64 /** 65 * The type of system 66 */ 67 typedef DifferentiableSystem sys_type; 68 69 /** 70 * Constructor. Requires a reference to the system 71 * to be solved. 72 */ 73 explicit 74 TimeSolver (sys_type & s); 75 76 /** 77 * Destructor. 78 */ 79 virtual ~TimeSolver (); 80 81 /** 82 * The initialization function. This method is used to 83 * initialize internal data structures before a simulation begins. 84 */ 85 virtual void init (); 86 87 /** 88 * The data initialization function. This method is used to 89 * initialize internal data structures after the underlying System 90 * has been initialized 91 */ 92 virtual void init_data (); 93 94 /** 95 * The reinitialization function. This method is used after 96 * changes in the mesh 97 */ 98 virtual void reinit (); 99 100 /** 101 * This method solves for the solution at the next timestep (or 102 * solves for a steady-state solution). Usually we will only need 103 * to solve one (non)linear system per timestep, but more complex 104 * subclasses may override this. 105 */ 106 virtual void solve (); 107 108 /** 109 * This method advances the solution to the next timestep, after a 110 * solve() has been performed. Often this will be done after every 111 * UnsteadySolver::solve(), but adaptive mesh refinement and/or adaptive 112 * time step selection may require some solve() steps to be repeated. 113 */ 114 virtual void advance_timestep (); 115 116 /** 117 * This method solves for the adjoint solution at the next adjoint timestep 118 * (or a steady state adjoint solve) 119 */ 120 virtual std::pair<unsigned int, Real> adjoint_solve (const QoISet & qoi_indices); 121 122 /** 123 * This method advances the adjoint solution to the previous 124 * timestep, after an adjoint_solve() has been performed. This will 125 * be done before every UnsteadySolver::adjoint_solve(). 126 */ 127 virtual void adjoint_advance_timestep (); 128 129 /** 130 * This method retrieves all the stored solutions at the current 131 * system.time 132 */ 133 virtual void retrieve_timestep(); 134 135 /** 136 * A method to integrate the adjoint sensitivity w.r.t a given parameter 137 * vector. int_{tstep_start}^{tstep_end} dQ/dp dt = int_{tstep_start}^{tstep_end} (\partialQ / \partial p) - ( \partial R (u,z) / \partial p ) dt 138 */ 139 virtual void integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities); 140 141 /** 142 * This method uses the DifferentiablePhysics 143 * element_time_derivative(), element_constraint(), and 144 * mass_residual() to build a full residual on an element. What 145 * combination 146 * 147 * it uses will depend on the type of solver. See 148 * the subclasses for more details. 149 */ 150 virtual bool element_residual (bool request_jacobian, 151 DiffContext &) = 0; 152 153 /** 154 * This method uses the DifferentiablePhysics 155 * side_time_derivative(), side_constraint(), and 156 * side_mass_residual() to build a full residual on an element's 157 * side. 158 * 159 * What combination it uses will depend on the type 160 * of solver. See the subclasses for more details. 161 */ 162 virtual bool side_residual (bool request_jacobian, 163 DiffContext &) = 0; 164 165 /** 166 * This method uses the DifferentiablePhysics 167 * nonlocal_time_derivative(), nonlocal_constraint(), and 168 * nonlocal_mass_residual() to build a full residual of non-local 169 * terms. 170 * 171 * What combination it uses will depend on the type 172 * of solver. See the subclasses for more details. 173 */ 174 virtual bool nonlocal_residual (bool request_jacobian, 175 DiffContext &) = 0; 176 177 /** 178 * This method is for subclasses or users to override 179 * to do arbitrary processing between timesteps 180 */ before_timestep()181 virtual void before_timestep () {} 182 183 /** 184 * \returns A constant reference to the system we are solving. 185 */ system()186 const sys_type & system () const { return _system; } 187 188 /** 189 * \returns A writable reference to the system we are solving. 190 */ system()191 sys_type & system () { return _system; } 192 193 /** 194 * An implicit linear or nonlinear solver to use at each timestep. 195 */ diff_solver()196 virtual std::unique_ptr<DiffSolver> & diff_solver() { return _diff_solver; } 197 198 /** 199 * An implicit linear solver to use for adjoint and sensitivity problems. 200 */ linear_solver()201 virtual std::unique_ptr<LinearSolver<Number>> & linear_solver() { return _linear_solver; } 202 203 /** 204 * Print extra debugging information if quiet == false. 205 */ 206 bool quiet; 207 208 /** 209 * Computes the size of ||u^{n+1} - u^{n}|| in some norm. 210 * 211 * \note While you can always call this function, its 212 * result may or may not be very meaningful. For example, if 213 * you call this function right after calling advance_timestep() 214 * then you'll get a result of zero since old_nonlinear_solution 215 * is set equal to nonlinear_solution in this function. 216 */ 217 virtual Real du(const SystemNorm & norm) const = 0; 218 219 /** 220 * Is this effectively a steady-state solver? 221 */ 222 virtual bool is_steady() const = 0; 223 224 /** 225 * This value (which defaults to zero) is the number of times the 226 * TimeSolver is allowed to halve deltat and let the DiffSolver 227 * repeat the latest failed solve with a reduced timestep. 228 * 229 * \note This has no effect for SteadySolvers. 230 * \note You must set at least one of the DiffSolver flags 231 * "continue_after_max_iterations" or 232 * "continue_after_backtrack_failure" to allow the TimeSolver to 233 * retry the solve. 234 */ 235 unsigned int reduce_deltat_on_diffsolver_failure; 236 237 /** 238 * A setter function users will employ if they need to do something 239 * other than save no solution history 240 */ 241 void set_solution_history(const SolutionHistory & _solution_history); 242 243 /** 244  * A getter function that returns a reference to the solution history 245  * object owned by TimeSolver 246  * */ 247 SolutionHistory & get_solution_history(); 248 249 /** 250 * Accessor for querying whether we need to do a primal 251 * or adjoint solve 252 */ is_adjoint()253 bool is_adjoint() const 254 { return _is_adjoint; } 255 256 /** 257 * Accessor for setting whether we need to do a primal 258 * or adjoint solve 259 */ set_is_adjoint(bool _is_adjoint_value)260 void set_is_adjoint(bool _is_adjoint_value) 261 { _is_adjoint = _is_adjoint_value; } 262 263 /** 264 * Function to return 'last_deltat()', returns system.deltat if 265 * fixed timestep solver is used, completed_deltat if the adaptive 266 * time solver is used. 267 */ 268 virtual Real last_complete_deltat(); 269 270 271 protected: 272 273 /** 274 * An implicit linear or nonlinear solver to use at each timestep. 275 */ 276 std::unique_ptr<DiffSolver> _diff_solver; 277 278 /** 279 * An implicit linear solver to use for adjoint problems. 280 */ 281 std::unique_ptr<LinearSolver<Number>> _linear_solver; 282 283 /** 284 * A reference to the system we are solving. 285 */ 286 sys_type & _system; 287 288 /** 289 * A std::unique_ptr to a SolutionHistory object. Default is 290 * NoSolutionHistory, which the user can override by declaring a 291 * different kind of SolutionHistory in the application 292 */ 293 std::unique_ptr<SolutionHistory> solution_history; 294 295 /** 296 * Definitions of argument types for use in refactoring subclasses. 297 */ 298 299 typedef bool (DifferentiablePhysics::*ResFuncType) (bool, DiffContext &); 300 301 typedef void (DiffContext::*ReinitFuncType) (Real); 302 303 private: 304 305 /** 306 * This boolean tells the TimeSolver whether we are solving a primal or 307 * adjoint problem 308 */ 309 bool _is_adjoint; 310 311 }; 312 313 314 } // namespace libMesh 315 316 317 #endif // LIBMESH_TIME_SOLVER_H 318