Line data    Source code

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