cl_IWG_Timestep.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 *
 * See the top-level LICENSE file for the complete license and disclaimer.
 */
 
#ifndef BELFEM_CL_IWG_TIMESTEP_HPP
#define BELFEM_CL_IWG_TIMESTEP_HPP
 
 
#include "cl_IWG.hpp"
#include "en_SolverEnums.hpp"
#include "cl_TimestepMatrices.hpp"
 
namespace belfem
{
    namespace fem
    {
        class Element ;

        class IWG_Timestep: public IWG
        {
            // Unit-test seam, no runtime cost and no public API widened.
            //
            // The BDF coefficients ( mAlpha, mBeta ) and the step-size history
            // that determines them ( mH, mStepCount ) are private on purpose.
            // The production writers are shift_fields, reset_fields,
            // restore_savepoint — each needing an initialized DofManager —
            // and, since 2026-08-15, the validated memdump seam
            // restore_history_state ( which needs none ). Without this
            // declaration the
            // variable-step coefficient formulas below are reachable by no
            // test at all: they are pure arithmetic, but their inputs cannot
            // be posed and their outputs cannot be read.
            //
            // The probe is defined in tests/fem/test_BdfTimestepMethod.cpp and
            // exists only there. Nothing in src/ may use it.
            friend class BdfCoefficientProbe ;
 
            EulerMethod mMethod = EulerMethod::BackwardDifference1 ;
 
 
 
            Cell< Cell< Vector< real > * > > mFieldData ;
 
            real mAlpha = 1.0 ;
            Vector< real > mBeta ;
            Vector< real > mH ;
 
            // true if BDF coefficients must be recomputed before the next assembly
            bool mCoeffsDirty = true ;
 
            // savepoint storage: deep copies of the dof fields and the step-size
            // history, so a rejected coupled timestep can restore the state
            // across any number of shift_fields calls ( see make_savepoint )
            Cell< Cell< Vector< real > > > mFieldSnapshot ;
            Vector< real > mHSnapshot ;
            real mDeltaTimeSnapshot = BELFEM_QUIET_NAN ;
            real mHDroppedSnapshot = BELFEM_QUIET_NAN ;
            uint mStepCountSnapshot = 0 ;
            bool mHaveSavepoint = false ;
 
            // the step size that falls off the end of mH on a shift: the
            // rotation destroys it, but reset_fields must restore it for a
            // BDF5 retry, which reads mH( 3 )
            real mHDropped = BELFEM_QUIET_NAN ;
 
            // number of timesteps started ( shifted ), drives the startup order ramp
            uint mStepCount = 0 ;
 
            // order of the scheme actually running this step ( startup ramp
            // may run below the configured mOrder ); consumed by collect_qhist
            uint mOrderActive = 1 ;
 
            // remember whether the scheme was configured with a stiffness matrix
            bool mHaveStiffness = true ;
 
            // timestepping pointer as configured by set_timestepping_method
            void
            ( IWG_Timestep:: * mTimestep )(
                    Matrix< real > & aJ,
                    Vector< real > & aRHS );
 
            // scheme that actually runs this step; differs from mTimestep only
            // during the startup ramp of a BDF-p run ( step 1: BDF1, step 2: BDF2, ... )
            void
            ( IWG_Timestep:: * mTimestepActive )(
                    Matrix< real > & aJ,
                    Vector< real > & aRHS ) = nullptr ;
 
//------------------------------------------------------------------------------
        protected:
//------------------------------------------------------------------------------
            //real mDeltaTime = 1.0 ;
            uint mOrder = 1 ;
 
            TimestepMatrices * mTimeStepMatrices ;
 
//------------------------------------------------------------------------------
        public:
//------------------------------------------------------------------------------
 
            IWG_Timestep( const IwgType aType,
                          const ModelDimensionality aDimensionality,
                          const IwgMode aMode=IwgMode::Iterative,
                          const SymmetryMode aSymmetryMode=SymmetryMode::Unsymmetric,
                          const DofMode      aDofMode=DofMode::AllBlocksEqual,
                          const SideSetDofLinkMode aSideSetDofLinkMode=SideSetDofLinkMode::FacetOnly );
 
 
            ~IWG_Timestep() override ;
 
            // set the timestepping method
            void
            set_timestepping_method( const EulerMethod aMethod, const bool aHaveStiffness = true ) override;
 
            // return the timestepping method
            EulerMethod
            method() const override;
 
            // BDF order actually running this step ( startup ramp )
            uint
            order_active() const;
 
            // BDF order the UPCOMING assembly will run; ask this, not
            // order_active(), before the step has been assembled
            uint
            order_pending() const;
 
//------------------------------------------------------------------------------
 
            void
            compute_jacobian_and_rhs(
                    Element        * aElement,
                    Matrix< real > & aJacobian,
                    Vector< real > & aRHS ) override;
 
//------------------------------------------------------------------------------
 
            void
            shift_fields() override;
 
//------------------------------------------------------------------------------
 
//------------------------------------------------------------------------------

            void
            compute_bdf_coefficients();
 
//------------------------------------------------------------------------------

            const Vector< real > &
            collect_qhist();
 
//------------------------------------------------------------------------------
 
            void
            reset_fields() override;
 
//------------------------------------------------------------------------------

            void
            make_savepoint();
 
//------------------------------------------------------------------------------

            void
            restore_savepoint();
 
//------------------------------------------------------------------------------

            void
            save_history_state(
                Vector< real > & aH,
                uint           & aStepCount,
                real           & aLastDeltaTime ) const ;
 
//------------------------------------------------------------------------------

            void
            history_field_labels(
                Cell< string > & aParents,
                Cell< string > & aHistory,
                const uint       aDepth = 0 ) const ;
 
//------------------------------------------------------------------------------

            bool
            restore_history_state(
                const Vector< real > & aH,
                const uint             aStepCount,
                const real             aLastDeltaTime ) ;
 
 //------------------------------------------------------------------------------
 
            TimestepMatrices *
            matrices() override ;
 
//------------------------------------------------------------------------------
        protected:
//------------------------------------------------------------------------------

            void
            compute_mkf( Element * aElement) override;
 
//------------------------------------------------------------------------------
 
            void
            compute_timestep( Matrix< real > & aJ,
                              Vector< real > & aRHS );
 
//------------------------------------------------------------------------------
 
            void
            set_field( DofManagerBase * aField ) override;
 
//------------------------------------------------------------------------------
        private:
//------------------------------------------------------------------------------
 
            // The startup ramp, in one place. BDF-p needs p history states, so
            // the first steps run the highest order the history supports.
            // Const and side-effect free, so it can be asked before the lazy
            // recompute has fired; compute_bdf_coefficients assigns its result
            // to mOrderActive, and order_pending() reports it to the log.
            inline uint
            ramped_order() const
            {
                switch( mMethod )
                {
                    case( EulerMethod::BackwardDifference2 ) :
                    case( EulerMethod::BackwardDifference3 ) :
                    case( EulerMethod::BackwardDifference4 ) :
                    case( EulerMethod::BackwardDifference5 ) :
                    {
                        return ( mOrder > 1 && mStepCount < mOrder )
                                ? ( mStepCount > 0 ? mStepCount : 1 )
                                : mOrder ;
                    }
                    default :
                    {
                        // every other scheme is single-state
                        return 1 ;
                    }
                }
            }
 
//------------------------------------------------------------------------------
 
            void
            create_old_dof_fields();
 
//------------------------------------------------------------------------------
 
            void
            static_subfield( Matrix< real > & aM, Vector< real > & aRHS );
 
//------------------------------------------------------------------------------
            void
            explicit_euler( Matrix< real > & aM, Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            crank_nicolson( Matrix< real > & aM, Vector< real > & aRHS );
 
//------------------------------------------------------------------------------
 
            void
            galerkin( Matrix< real > & aM, Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            mass_matrix_only(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
//------------------------------------------------------------------------------
 
            void
            stiffness_matrix_only(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf1(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf2(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf3(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf4(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf5(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf1_nok(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf2_nok(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf3_nok(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf4_nok(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            bdf5_nok(
                Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            compute_bdf_coeffs_2();
 
//------------------------------------------------------------------------------
 
            void
            compute_bdf_coeffs_3();
 
//------------------------------------------------------------------------------
 
            void
            compute_bdf_coeffs_4();
 
//------------------------------------------------------------------------------
 
            void
            compute_bdf_coeffs_5();
 
//------------------------------------------------------------------------------
 
            void
            derivative( Matrix< real > & aJ,
                Vector< real > & aRHS  );
 
//------------------------------------------------------------------------------
 
            void
            euler_no_stiffness( Matrix< real > & aJ,
                                  Vector< real > & aRHS );
 
//------------------------------------------------------------------------------
 
            void
            derivative_no_stiffness( Matrix< real > & aJ,
                                  Vector< real > & aRHS );
 
 
//------------------------------------------------------------------------------
 
            uint
            timestepping_order() const override ;
 
//------------------------------------------------------------------------------
        };
 
//------------------------------------------------------------------------------
 
        // return the timestepping method
        inline EulerMethod
        IWG_Timestep::method() const
        {
            return mMethod ;
        }
 
//------------------------------------------------------------------------------
 
        // BDF order actually running this step ( startup ramp may sit below
        // the configured order ); the controller's step-growth clamp reads it
        inline uint
        IWG_Timestep::order_active() const
        {
            return mOrderActive ;
        }
 
//------------------------------------------------------------------------------
 
        // BDF order the UPCOMING assembly will run. The coefficients are
        // recomputed lazily, on the first element of a step, so between
        // shift_fields and that first assembly mOrderActive still holds the
        // PREVIOUS step's order. Anything reporting the order before assembly
        // -- the timestep header, above all -- must ask this instead of
        // order_active(), or it lags the ramp by one step.
        inline uint
        IWG_Timestep::order_pending() const
        {
            return mCoeffsDirty ? this->ramped_order() : mOrderActive ;
        }
 
//------------------------------------------------------------------------------
 
        inline uint
        IWG_Timestep::timestepping_order() const
        {
            return mOrder ;
        }
 
//------------------------------------------------------------------------------
    }
}
 
#endif //BELFEM_CL_IWG_TIMESTEP_HPP