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_NEWMARK_SOLVER_H
      21             : #define LIBMESH_NEWMARK_SOLVER_H
      22             : 
      23             : // Local includes
      24             : #include "libmesh/second_order_unsteady_solver.h"
      25             : 
      26             : namespace libMesh
      27             : {
      28             : /**
      29             :  * This class defines a Newmark time integrator for
      30             :  * second order (in time) DifferentiableSystems.
      31             :  * There are two parameters
      32             :  * \f$\gamma\f$ and \f$\beta\f$ (defaulting to 0.5 and 0.25,
      33             :  * respectively).
      34             :  *
      35             :  * \note Projectors are included for the initial velocity and
      36             :  * acceleration; the initial displacement can be set by calling
      37             :  * System::project_solution().
      38             :  *
      39             :  * This class is part of the new DifferentiableSystem framework,
      40             :  * which is still experimental.  Users of this framework should
      41             :  * beware of bugs and future API changes.
      42             :  *
      43             :  * \author Paul T. Bauman
      44             :  * \date 2015
      45             :  */
      46          48 : class NewmarkSolver : public SecondOrderUnsteadySolver
      47             : {
      48             : public:
      49             :   /**
      50             :    * The parent class
      51             :    */
      52             :   typedef UnsteadySolver Parent;
      53             : 
      54             :   /**
      55             :    * Constructor. Requires a reference to the system
      56             :    * to be solved.
      57             :    */
      58             :   explicit
      59             :   NewmarkSolver (sys_type & s);
      60             : 
      61             :   /**
      62             :    * Destructor.
      63             :    */
      64             :   virtual ~NewmarkSolver ();
      65             : 
      66             :   /**
      67             :    * This method advances the solution to the next timestep, after a
      68             :    * solve() has been performed.  Often this will be done after every
      69             :    * UnsteadySolver::solve(), but adaptive mesh refinement and/or adaptive
      70             :    * time step selection may require some solve() steps to be repeated.
      71             :    */
      72             :   virtual void advance_timestep () override;
      73             : 
      74             :   /**
      75             :    * This method advances the adjoint solution to the previous
      76             :    * timestep, after an adjoint_solve() has been performed.  This will
      77             :    * be done before every UnsteadySolver::adjoint_solve().
      78             :    */
      79             :   virtual void adjoint_advance_timestep () override;
      80             : 
      81             :   /**
      82             :    * This method uses the specified initial displacement and velocity
      83             :    * to compute the initial acceleration \f$a_0\f$.
      84             :    */
      85             :   virtual void compute_initial_accel();
      86             : 
      87             :   /**
      88             :    * Specify non-zero initial acceleration. Should be called before solve().
      89             :    * This is an alternative to compute_initial_acceleration() if the
      90             :    * initial acceleration is actually known.
      91             :    * The function value f and its gradient g are user-provided cloneable functors.
      92             :    * A gradient g is only required/used for projecting onto finite element spaces
      93             :    * with continuous derivatives.
      94             :    */
      95             :   void project_initial_accel (FunctionBase<Number> * f,
      96             :                               FunctionBase<Gradient> * g = nullptr);
      97             : 
      98             :   /**
      99             :    * Allow the user to (re)set whether the initial acceleration is available.
     100             :    * This is not needed if either compute_initial_accel() or project_initial_accel()
     101             :    * are called. This is useful is the user is restarting a calculation and the acceleration
     102             :    * is available from the restart.
     103             :    */
     104             :   void set_initial_accel_avail (bool initial_accel_set);
     105             : 
     106             :   /**
     107             :    * Error convergence order: 2 for \f$\gamma=0.5\f$, 1 otherwise
     108             :    */
     109             :   virtual Real error_order() const override;
     110             : 
     111             :   /**
     112             :    * This method solves for the solution at the next timestep.
     113             :    * Usually we will only need to solve one (non)linear system per timestep,
     114             :    * but more complex subclasses may override this.
     115             :    */
     116             :   virtual void solve () override;
     117             : 
     118             :   /**
     119             :    * This method uses the DifferentiablePhysics'
     120             :    * element_time_derivative() and element_constraint()
     121             :    * to build a full residual on an element.  What combination
     122             :    * it uses will depend on theta.
     123             :    */
     124             :   virtual bool element_residual (bool request_jacobian,
     125             :                                  DiffContext &) override;
     126             : 
     127             :   /**
     128             :    * This method uses the DifferentiablePhysics'
     129             :    * side_time_derivative() and side_constraint()
     130             :    * to build a full residual on an element's side.
     131             :    * What combination it uses will depend on theta.
     132             :    */
     133             :   virtual bool side_residual (bool request_jacobian,
     134             :                               DiffContext &) override;
     135             : 
     136             :   /**
     137             :    * This method uses the DifferentiablePhysics'
     138             :    * nonlocal_time_derivative() and nonlocal_constraint()
     139             :    * to build a full residual for non-local terms.
     140             :    * What combination it uses will depend on theta.
     141             :    */
     142             :   virtual bool nonlocal_residual (bool request_jacobian,
     143             :                                   DiffContext &) override;
     144             : 
     145             : 
     146             :   /**
     147             :    * Setter for \f$ \beta \f$
     148             :    */
     149             :   void set_beta ( Real beta )
     150             :   { _beta = beta; }
     151             : 
     152             :   /**
     153             :    * Setter for \f$ \gamma \f$.
     154             :    *
     155             :    * \note The method is second order only for \f$ \gamma = 1/2 \f$.
     156             :    */
     157             :   void set_gamma ( Real gamma )
     158             :   { _gamma = gamma; }
     159             : 
     160             : protected:
     161             : 
     162             :   /**
     163             :    * The value for the \f$\beta\f$ to employ. Method is
     164             :    * unconditionally stable for
     165             :    * \f$ \beta \ge \frac{1}{4} \left( \gamma + \frac{1}{2}\right)^2 \f$
     166             :    */
     167             :   Real _beta;
     168             : 
     169             :   /**
     170             :    * The value for \f$\gamma\f$ to employ. Newmark
     171             :    * is 2nd order iff \f$\gamma=0.5\f$.
     172             :    */
     173             :   Real _gamma;
     174             : 
     175             :   /**
     176             :    * Need to be able to indicate to _general_residual if we
     177             :    * are doing an acceleration solve or not.
     178             :    */
     179             :   bool _is_accel_solve;
     180             : 
     181             : 
     182             :   /**
     183             :    * This method requires an initial acceleration. So, we force the
     184             :    * user to call either compute_initial_accel or project_initial_accel
     185             :    * to set the initial acceleration.
     186             :    */
     187             :   bool _initial_accel_set;
     188             : 
     189             :   /**
     190             :    * This method is the underlying implementation of the public
     191             :    * residual methods.
     192             :    */
     193             :   virtual bool _general_residual (bool request_jacobian,
     194             :                                   DiffContext &,
     195             :                                   ResFuncType mass,
     196             :                                   ResFuncType damping,
     197             :                                   ResFuncType time_deriv,
     198             :                                   ResFuncType constraint,
     199             :                                   ReinitFuncType reinit);
     200             : };
     201             : 
     202             : } // namespace libMesh
     203             : 
     204             : #endif // LIBMESH_NEWMARK_SOLVER_H