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_UNSTEADY_SOLVER_H
      21             : #define LIBMESH_UNSTEADY_SOLVER_H
      22             : 
      23             : // Local includes
      24             : #include "libmesh/libmesh_common.h"
      25             : #include "libmesh/time_solver.h"
      26             : 
      27             : // C++ includes
      28             : #include <memory>
      29             : 
      30             : namespace libMesh
      31             : {
      32             : 
      33             : // Forward declarations
      34             : template <typename T> class NumericVector;
      35             : 
      36             : /**
      37             :  * This is a generic class that defines a solver to handle
      38             :  * time integration of DifferentiableSystems.
      39             :  *
      40             :  * A user can define a solver for unsteady problems by deriving
      41             :  * from this class and implementing certain functions.
      42             :  *
      43             :  * This class is part of the new DifferentiableSystem framework,
      44             :  * which is still experimental.  Users of this framework should
      45             :  * beware of bugs and future API changes.
      46             :  *
      47             :  * \author Roy H. Stogner
      48             :  * \date 2008
      49             :  */
      50         264 : class UnsteadySolver : public TimeSolver
      51             : {
      52             : public:
      53             :   /**
      54             :    * Constructor. Requires a reference to the system
      55             :    * to be solved.
      56             :    */
      57             :   explicit
      58             :   UnsteadySolver (sys_type & s);
      59             : 
      60             :   /**
      61             :    * Destructor.
      62             :    */
      63             :   virtual ~UnsteadySolver ();
      64             : 
      65             :   /**
      66             :    * The initialization function.  This method is used to
      67             :    * initialize internal data structures before a simulation begins.
      68             :    */
      69             :   virtual void init () override;
      70             : 
      71             :   /**
      72             :    * @brief Add adjoint vectors and old_adjoint_vectors
      73             :    * as per the indices of QoISet
      74             :    */
      75             :   virtual void init_adjoints () override;
      76             : 
      77             :   /**
      78             :    * The data initialization function.  This method is used to
      79             :    * initialize internal data structures after the underlying System
      80             :    * has been initialized
      81             :    */
      82             :   virtual void init_data () override;
      83             : 
      84             :   /**
      85             :    * The reinitialization function.  This method is used to
      86             :    * resize internal data vectors after a mesh change.
      87             :    */
      88             :   virtual void reinit () override;
      89             : 
      90             :   /**
      91             :    * This method solves for the solution at the next timestep.
      92             :    * Usually we will only need to solve one (non)linear system per timestep,
      93             :    * but more complex subclasses may override this.
      94             :    */
      95             :   virtual void solve () override;
      96             : 
      97             :   /**
      98             :    * This method advances the solution to the next timestep, after a
      99             :    * solve() has been performed.  Often this will be done after every
     100             :    * UnsteadySolver::solve(), but adaptive mesh refinement and/or adaptive
     101             :    * time step selection may require some solve() steps to be repeated.
     102             :    */
     103             :   virtual void advance_timestep () override;
     104             : 
     105             :   void update();
     106             : 
     107             :   /**
     108             :    * This method solves for the adjoint solution at the next adjoint timestep
     109             :    * (or a steady state adjoint solve)
     110             :    */
     111             :   virtual std::pair<unsigned int, Real> adjoint_solve (const QoISet & qoi_indices) override;
     112             : 
     113             :   /**
     114             :    * This method advances the adjoint solution to the previous
     115             :    * timestep, after an adjoint_solve() has been performed.  This will
     116             :    * be done before every UnsteadySolver::adjoint_solve().
     117             :    */
     118             :   virtual void adjoint_advance_timestep () override;
     119             : 
     120             :   /**
     121             :    * This method retrieves all the stored solutions at the current
     122             :    * system.time
     123             :    */
     124             :   virtual void retrieve_timestep () override;
     125             : 
     126             :   /**
     127             :    * A method to integrate the system::QoI functionals.
     128             :    */
     129             :   virtual void integrate_qoi_timestep() override;
     130             : 
     131             :   /**
     132             :    * A method to integrate the adjoint sensitivity w.r.t a given parameter
     133             :    * vector. int_{tstep_start}^{tstep_end} dQ/dp dt = int_{tstep_start}^{tstep_end} (\partialQ / \partial p) - ( \partial R (u,z) / \partial p ) dt
     134             :    * The trapezoidal rule is used to numerically integrate the timestep.
     135             :    */
     136             :   virtual void integrate_adjoint_sensitivity(const QoISet & qois, const ParameterVector & parameter_vector, SensitivityData & sensitivities) override;
     137             : 
     138             : #ifdef LIBMESH_ENABLE_AMR
     139             :   /**
     140             :    * A method to compute the adjoint refinement error estimate at the current timestep.
     141             :    * int_{tstep_start}^{tstep_end} R(u^h,z) dt
     142             :    * The user provides an initialized ARefEE object.
     143             :    * Fills in an ErrorVector that contains the weighted sum of errors from all the QoIs and can be used to guide AMR.
     144             :    * CURRENTLY ONLY SUPPORTED for Backward Euler.
     145             :    */
     146             :   virtual void integrate_adjoint_refinement_error_estimate(AdjointRefinementEstimator & /*adjoint_refinement_error_estimator*/, ErrorVector & /*QoI_elementwise_error*/) override;
     147             : #endif // LIBMESH_ENABLE_AMR
     148             : 
     149             :   /**
     150             :    * This method should return the expected convergence order of the
     151             :    * (non-local) error of the time discretization scheme - e.g. 2 for the
     152             :    * O(deltat^2) Crank-Nicholson, or 1 for the O(deltat) Backward Euler.
     153             :    *
     154             :    * Useful for adaptive timestepping schemes.
     155             :    */
     156             :   virtual Real error_order () const = 0;
     157             : 
     158             :   /**
     159             :    * \returns The maximum order of time derivatives for which the
     160             :    * UnsteadySolver subclass is capable of handling.
     161             :    *
     162             :    * For example, EulerSolver will have \p time_order() = 1 and
     163             :    * NewmarkSolver will have \p time_order() = 2.
     164             :    */
     165             :   virtual unsigned int time_order () const = 0;
     166             : 
     167             :   /**
     168             :    * \returns The old nonlinear solution for the specified global
     169             :    * DOF.
     170             :    */
     171             :   Number old_nonlinear_solution (const dof_id_type global_dof_number) const;
     172             : 
     173             :   /**
     174             :    * Serial vector of _system.get_vector("_old_nonlinear_solution")
     175             :    * This is a shared_ptr so that it can be shared between different
     176             :    * derived class instances, as in e.g. AdaptiveTimeSolver.
     177             :    */
     178             :   std::shared_ptr<NumericVector<Number>> old_local_nonlinear_solution;
     179             : 
     180             :   /**
     181             :    * Computes the size of ||u^{n+1} - u^{n}|| in some norm.
     182             :    *
     183             :    * \note While you can always call this function, its
     184             :    * result may or may not be very meaningful.  For example, if
     185             :    * you call this function right after calling advance_timestep()
     186             :    * then you'll get a result of zero since old_nonlinear_solution
     187             :    * is set equal to nonlinear_solution in this function.
     188             :    */
     189             :   virtual Real du(const SystemNorm & norm) const override;
     190             : 
     191             :   /**
     192             :    * This is not a steady-state solver.
     193             :    */
     194    93335107 :   virtual bool is_steady() const override { return false; }
     195             : 
     196             :   /**
     197             :    * A setter for the first_adjoint_step boolean. Needed for nested time solvers.
     198             :    */
     199         120 :   void set_first_adjoint_step(bool first_adjoint_step_setting)
     200             :   {
     201        4200 :     first_adjoint_step = first_adjoint_step_setting;
     202         120 :   }
     203             : 
     204             :   /*
     205             :    * A setter for the first_solver boolean. Useful for example if we are using
     206             :    * a nested time solver, and the outer solver wants to tell the inner one that
     207             :    * the initial conditions have already been handled.
     208             :    */
     209          12 :   void set_first_solve(bool first_solve_setting)
     210             :   {
     211         420 :     first_solve = first_solve_setting;
     212         408 :   }
     213             : 
     214             : protected:
     215             : 
     216             :   /**
     217             :    * A bool that will be true the first time solve() is called,
     218             :    * and false thereafter
     219             :    */
     220             :   bool first_solve;
     221             : 
     222             :   /**
     223             :    * A bool that will be true the first time adjoint_advance_timestep() is called,
     224             :    * (when the primal solution is to be used to set adjoint boundary conditions) and false thereafter
     225             :    */
     226             :   bool first_adjoint_step;
     227             : 
     228             :   /**
     229             :    * A vector of pointers to vectors holding the adjoint solution at the last time step
     230             :    */
     231             :   std::vector< std::unique_ptr<NumericVector<Number>> > old_adjoints;
     232             : 
     233             :   /**
     234             :    * We will need to move the system.time around to ensure that residuals
     235             :    * are built with the right deltat and the right time.
     236             :    */
     237             :   Real last_step_deltat;
     238             :   Real next_step_deltat;
     239             : };
     240             : 
     241             : 
     242             : 
     243             : } // namespace libMesh
     244             : 
     245             : 
     246             : #endif // LIBMESH_UNSTEADY_SOLVER_H