cl_FEM_Controller.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 CL_FEM_CONTROLLER_HPP
#define CL_FEM_CONTROLLER_HPP
#include "cl_FEM_DofMgr_SolverData.hpp"
#include "typedefs.hpp"
#include "cl_Timer.hpp"
#include "cl_Input_Section.hpp"
#include "cl_Circuit.hpp"
#include "cl_FEM_PhysicalBoundaryCondition.hpp"
#include "cl_FEM_Kernel.hpp"
#include "cl_ShiftRegister.hpp"
 
namespace belfem
{
    class Mesh ;
 
    namespace fem
    {
        class Kernel ;
        class IWG_Timestep ;

        class Controller
        {
            const proc_t   mCommRank ;
            const proc_t   mCommSize ;
            Kernel       * mKernel = nullptr ;
            IWG_Timestep * mEquation = nullptr ;
            Mesh         * mMesh = nullptr ;
            Kernel       * mKernel2   = nullptr ;
            IWG_Timestep * mEquation2 = nullptr ;
            Mesh         * mMesh2     = nullptr ;
            Circuit      * mCircuit   = nullptr ;
 
            real mTime = 0.0 ;
            real mDeltaTime = BELFEM_QUIET_NAN ;
            real mDeltaTimeTemporary = BELFEM_QUIET_NAN ; //Temporary time step to reach the time steps to be saved
            real mDeltaTime0 = BELFEM_QUIET_NAN ;

            real mDeltaTimeInitial = BELFEM_QUIET_NAN ;
            uint mRunningTimeStep = 0 ;
            uint mMeshTimeStep = 1 ;
 
            real  mTime2 =  0.0 ;
            real  mDeltaTime2 = BELFEM_QUIET_NAN ;
            real mDeltaTimeTemporary2 = BELFEM_QUIET_NAN ; //Temporary time step to reach the time steps to be saved
 
            // - - - - - - begin user settings - - - - - -
 
            real mAlpha = 0.5 ;
            real mBeta  = 1.1 ;
            real mGamma = 0.4 ;
 
            uint mMinNumIterations = 2 ;
            uint mMaxNumIterations = 100 ;
            uint mIterationTarget = 20 ;

            uint mDivergenceStrikes = 0 ;
 
            // stagnation guard: if the mean absolute deviation of the last
            // mStallWindow residuals (in dB) stays below mStallBand while the
            // residual is still above tolerance, the timestep is cut instead of
            // grinding on to mMaxNumIterations (see iterate_coupled).
            // The band must sit ABOVE the noise of a floored residual, or the
            // Newton->Picard demotion this guard exists for never fires: the
            // sidecoatings floor wandered +-0.05..0.2 dB against the old
            // 0.001 dB band and burned the full watchdog window instead.
            // Genuine convergence moves >= O(1) dB per iterate, so 0.2 dB
            // separates the two regimes. Deck override: "stall tolerance"
            uint mStallWindow = 5 ;
            real mStallBand   = 0.2 ; // dB
 
            // progress watchdog: cut the timestep when the residual has not
            // reached a NEW MINIMUM for this many iterations while still above
            // 10x tolerance. Catches what the flat-band guard cannot see:
            // omega-sawtooth limit cycles ( locally "converging" on any short
            // window ) and slow monotone creep of the thermal residual at the
            // relaxation floor. 0 disables. Input key: "watchdog window".
            uint mWatchdogWindow  = 30 ;
            uint mWatchdogWindow2 = 30 ;

            real mWatchdogOmegaPrev  = BELFEM_QUIET_NAN ;
            real mWatchdogOmegaPrev2 = BELFEM_QUIET_NAN ;
 
            uint mMinNumIterations2 = 2 ;
            uint mMaxNumIterations2 = 100 ;
            uint mIterationTarget2 = 20 ;
 
            uint mMaxNumIterationsDiv = 10 ; //Number of diverging iterations before resetting the relaxation
            uint mNumIterationsDiv = 0; //Number of diverging iterations
            // conditioning diagnostic, per field. It describes a property of
            // ONE solve, and the two fields configure their solvers
            // separately ( "linear magnetic" / "linear thermal", falling back
            // to "linear" ), so the flag lives with them.
            // "compute conditioning" in the solver section still sets both.
            //
            // These gate the EIGENVALUE estimate and nothing else. The MUMPS
            // Arioli/Demmel/Duff numbers have their own flag below -- the two
            // diagnostics are independent, and either may be asked for alone
            bool mComputeConditioning  = false ; // magnetic; always costs ARPACK
            bool mComputeConditioning2 = false ; // thermal
 
            // "mumps error analysis", per field, read from the same three
            // places as the flag above. Gates MUMPS ICNTL(11) ONLY -- the
            // COND1 / COND2 / omega2 triple that feeds the "MUMPS ADD" footer
            // rows. A no-op on any other library, which is why the key names
            // MUMPS: see the warnings raised for a field that asks for it
            // without one
            bool mMumpsErrorAnalysis  = false ; // magnetic
            bool mMumpsErrorAnalysis2 = false ; // thermal
 
            // guards check_thermal_diagnostics() against warning twice: both
            // attach paths call it, and a caller can use both
            bool mThermalDiagnosticsChecked = false ;
 
            // whether set_params has run. The thermal diagnostics must not be
            // judged against the DEFAULT flags, so their check requires both
            // prerequisites -- parsed params AND an attached kernel -- and
            // runs from whichever call arrives second, in any order
            bool mParamsSet = false ;

            real mOmegaNoiseBand = 0.05 ;
 
            // MIT-3a: consecutive first-trial Newton steps that at least
            // halved the residual. The geometric omega recovery
            // ( tGrowth = 2.0 ) requires TWO of them, because a single
            // qualifying step is not evidence of contraction near a noise
            // floor -- on tapestack3d the doubling fired once ( omega
            // 0.483 -> 0.966 ) and regressed on the very next iterate.
            // SCOPE: this addresses the x2 rule ONLY; a later loss in the
            // same trace came from ordinary arctan growth and is NOT fixed
            // here. A healthily
            // contracting Newton earns the streak on its second step and
            // keeps the acceleration; a lucky step never does.
            //
            // Reset on: timestep start ( both drivers ), reject,
            // promotion, escalation, and every demotion ON THE COUPLED
            // PATH -- the evidence belongs to one uninterrupted Newton
            // run, not to the timestep. NOTE this is NOT simply
            // "wherever mFirstFlip is reset": iterate_magnetic's
            // promotion clears mFirstFlip without touching the streak.
            // That twin never increments and never consumes the streak
            // ( it has no x2 rule ), and initialize_magnetic zeroes at
            // every step start, so the omission is inert -- but do not
            // rely on the mFirstFlip locator when adding a path
            uint mNewtonTrustStreak = 0 ;
 
            real mEpsilonSwitch = 1e-4 ; // residual at which Picard hands off to Newton (Messe et al. 2023, Eq. 13)
            real mRelativeEpsilonTarget = 1e-6 ; // relative error criterion
            real mAbsoluteEpsilonTarget = 0.0 ;
 
            real mEpsilonSwitch2 = 1e-3 ; // thermal residual at which Picard hands off to Newton
            real mRelativeEpsilonTarget2 = 1e-6 ; // relative thermal error criterion
            real mAbsoluteEpsilonTarget2 = 0.0 ;  // absolute thermal error criterion ( see above: opt-in )

            real mThermalUpdateGate = BELFEM_REAL_MAX ;

            bool mThermalFrozen = false ;
 
            real mSimulationTime = BELFEM_REAL_MAX ;
            bool mAdaptTimestep = true ;
            real mSaveEvery = 0.0 ;
            bool mLastSave = false ;
            bool mSave = true ;
            
            // - - - - - - end user settings - - - - - -
 
            // todo: add flag if we want to reset omega at each timestep
 
            // end user settings
            real mOmegaNewton = 1.0 ;
            real mOmegaPicard = 1.0 ;
            real mOmegaNewton2 = 1.0 ;
            real mOmegaPicard2 = 1.0 ;
            bool mJustPicard = false ;
            bool mJustPicard2 = false ;

            uint mThermalFlipCount = 0 ;

            uint mAndersonDepth  = 0 ;
            uint mAndersonDepth2 = 0 ;

            uint mAndersonCommits  = 0 ;
            uint mAndersonCommits2 = 0 ;
 
            // Picard-breakdown rescue ( see try_escalate_to_newton ). When the magnetic
            // Picard iteration diverges at a residual far above mEpsilonSwitch -- so the
            // normal Picard->Newton promotion never fires -- the controller escalates to
            // the configured Newton tangent once before cutting the timestep. mForceNewton
            // pins the solver to Newton for the rest of the attempt; mNewtonEscalated caps
            // the escalation at one per timestep attempt so it cannot ping-pong.
            bool mForceNewton = false ;
            bool mNewtonEscalated = false ;
 
            // Thermal twin of the pair above ( see try_escalate_thermal_to_newton,
            // ). The thermal Picard residual can floor bit-flat ABOVE
            // mEpsilonSwitch2, and since the coupled Picard->Newton promotion fires
            // only below the switch, Newton is gated behind progress that only
            // Newton can make. Same lifecycle as the magnetic pair: cleared in
            // initialize_timestep / initialize_thermal, at most one escalation
            // per timestep attempt.
            bool mForceNewton2 = false ;
            bool mNewtonEscalated2 = false ;
 
            SolverAlgorithm mAlgorithm = SolverAlgorithm::Picard ;
            SolverAlgorithm mAlgorithmThermal = SolverAlgorithm::Picard ;
            EulerMethod          mTimeStepping = EulerMethod::BackwardDifference1 ;
 
            real mOmega0 = BELFEM_QUIET_NAN ;
 
            real mOmegaMin = 0.001 ;
            real mOmegaMax = 1.0 ;
            real mOmegaMin2 = 0.1 ;
            real mOmegaMax2 = 1.0 ;
 
            real mTime0 = 0.0 ;
            real mTime02  = 0.0 ;
            Timer * mTimer = nullptr ;
 
            real mEpsilon  = BELFEM_REAL_MAX ;
            real mEpsilon0 = BELFEM_REAL_MAX ;
 
            real mEpsilon2 = BELFEM_REAL_MAX ;
            real mEpsilon20 = BELFEM_REAL_MAX ;

            real mEpsilonAbs2 = BELFEM_REAL_MAX ;

            real mEpsilonAbs = BELFEM_REAL_MAX ;

            uint mThermalFlatCount = 0 ;
            bool mThermalStalled = false ;

            bool mMagneticHitTarget = false ;

            bool mMagBodyRanLastTrip = true ;
 
            bool mTripExit = false ;

            real mBestEpsilon  = BELFEM_REAL_MAX ;
            real mBestEpsilon2 = BELFEM_REAL_MAX ;
            uint mBestEpsilonIteration  = 0 ;
            uint mBestEpsilonIteration2 = 0 ;

            uint mSolverFailCount = 0 ;

            bool mUseLegacyTimestepControl = false ;

            real mCtrlKp = 0.15 ;
            real mCtrlKi = 0.30 ;
            real mCtrlKd = 0.0 ;

            real mCtrlErr0 = 1.0 ;
            real mCtrlErr1 = 1.0 ;
            real mCtrlErr2 = 1.0 ;

            uint mPostFailureHold      = 0 ;
            uint mPostFailureHoldSteps = 2 ;

            uint mFloorRetries       = 0 ;
            uint mFloorEscalationCap = 4 ;

            uint mMaxFloorRetries = 20 ;

            uint mMaxNumIterationsDeck  = 100 ;
            uint mMaxNumIterations2Deck = 100 ;
            uint mWatchdogWindowDeck    = 30 ;
            uint mWatchdogWindow2Deck   = 30 ;

            bool mAllowRestart = true ;

            uint mIteration = 0 ;
            uint mIteration2 = 0 ;
            uint mIterationTime = 0 ;
            uint mIteration0 = 0 ;
            uint mIteration02 = 0 ;
            uint mEigenAnalysisTime = 0 ;
            uint mPostprocesingTime = 0 ;


            // slot 0 : the eigen estimate ( spectral ratio, or kappa_2 when
            //          the symmetric driver earned the name )
            // slot 1 : MUMPS ADD COND1
            // slot 2 : MUMPS ADD COND2
            // slot 3 : MUMPS omega2 -- NOT reported, it is the discriminator
            //          that says whether slot 2 is a measurement at all
            Vector< real > mConditionNumbers1  = { BELFEM_QUIET_NAN, BELFEM_QUIET_NAN, BELFEM_QUIET_NAN, BELFEM_QUIET_NAN };
            Vector< real > mConditionNumbers2  = { BELFEM_QUIET_NAN, BELFEM_QUIET_NAN, BELFEM_QUIET_NAN, BELFEM_QUIET_NAN };
 
            real mDeltaTimeMax = BELFEM_REAL_MAX ;
            real mDeltaTimeMin = 1e-10 ;
            bool mReset = false ;
            bool mResetThermal = false ;
            bool mFirstFlip = false ;
            bool mFirstFlip2 = false ;
 
            bool mIsFullyCoupled = true ;
            real mCouplingFactor = 1 ; //Number of thermal time steps within each magnetic time steps
 
 
            Cell < PhysicalBoundaryCondition * > mCircuitCurrentBCs ;
            Cell < PhysicalBoundaryCondition * > mCircuitVoltageBCs ;
 
            Vector< real > mLHS ;
            Vector< real > mLHS0 ;
            bool mFirstIVSave = true ;

            Cell< string > mIVNamesI ;
            Cell< string > mIVNamesU ;
 
            // ring buffer of the last mStallWindow residuals (in dB)
            ShiftRegister< real > mResidualHistory ;
 
            Cell< Vector< real > > mBackupFields ;
            Vector< real > mBackupDofValues ;

            real mLastJJcMax = 0. ;

            bool mBoxSectionOpen = false ;

            bool mPostProcessed = false ;
 
        public:
 
            Controller( Kernel * aKernel,  Kernel * aKernel2 = nullptr ) ;
 
            ~Controller();
 
            void
            initialize_timestep();
 
            void
            initialize_magnetic();
 
            void
            initialize_thermal();
 
//------------------------------------------------------------------------------

            void
            solve_coupled();
 
            void
            solve_magnetic();
 
            void
            solve_thermal();
 
            bool
            solve_circuit();
 
            void
            finalize( const bool aPostProcess = true );
 
            const real &
            time() const;
 
            const real &
            time_thermal() const;
 
            real
            epsilon() const ;
 
            uint
            iteration() const ;
 
 
            bool
            is_fullycoupled() const ;
 
            bool
            reset() const ;
 
            void
            save( const std::string & aFilename );
 
            void
            save_IV( const std::string & aFilename );

            void
            create_iv_names();
 
            void
            set_params( const input::Section * aSection );
 
            void
            set_thermal_kernel( Kernel * aKernel );
 
            real
            simulation_time()  const ;
 
            bool
            save()  const ;
 
            void
            set_circuit( Circuit * aCircuit );
 
            Circuit *
            circuit() ;
 
            Kernel *
            kernel() ;
 
            Kernel *
            thermal_kernel() ;
 
            void
            save_memdump( const string & aPath );
 
            void
            load_memdump( const string & aPath );
 
            EulerMethod
            euler_method() const ;
 
 
        private:
 
            // imposes the voltage/current boundary conditions on the RHS ( rank 0 only,
            // after assembly and before the solve ); shared by the magnetic iterate paths
 
            void
            iterate_coupled();
 
            void
            iterate_magnetic();
 
            void
            iterate_thermal();
 
            void
            impose_voltage_bcs();
 
            // env-gated diagnostic dump of the assembled system, shared by both
            // iterate paths; off unless BELFEM_DUMP_SYSTEM is set.
            // aTag names the field in the file name, and aCount is that field's
            // OWN budget -- one shared counter let whichever field assembled
            // first spend all four dumps and starve the other
            void
            dump_system_if_requested( DofManager * aDofMgr,
                                      const char * aTag,
                                      int        & aCount );

            int mDumpCountMagnetic = 0 ;
            int mDumpCountThermal  = 0 ;
 
 
            // magnetic stagnation guard: if the residual window is flat ( deviation below
            // mStallBand ), Newton falls back to Picard internally; returns true only when
            // Picard itself has stalled, signalling the caller to reset_timestep()
            bool
            magnetic_stagnation_forces_reset();
 
            // progress watchdogs ( see mWatchdogWindow ): update the best-residual
            // tracker and return true when the attempt has made no new minimum for
            // the configured window while still far from tolerance -- UNLESS the
            // field's relaxation grew since the previous call, which marks a line
            // search recovering from an overshoot rather than a stall ( replayed
            // 2026-08-21: fires on the omega-pinned grinds, spares the AIMD
            // recovery whose cut used to cascade the timestep down ).
            //
            // aOmega must be the relaxation that PRODUCED the current residual --
            // the caller's post-line-search tOmega / tOmega2 -- never re-derived
            // from the live algorithm, which a same-call escalation or a
            // stagnation demotion may already have flipped or mutated.
            bool
            watchdog_magnetic( const real aOmega );
 
            bool
            watchdog_thermal( const real aOmega );
 
            // Anderson bookkeeping ( plan §3.2 ): SolverData stages the pair,
            // the controller decides its fate. Helpers keep the commit
            // counters and the per-kernel SolverData in sync; all are cheap
            // no-ops while the respective depth is 0.
            void anderson_commit_magnetic();
            void anderson_discard_magnetic();
            void anderson_clear_magnetic();
            void anderson_commit_thermal();
            void anderson_clear_thermal();
 
            // Picard-breakdown rescue: if the configured terminal algorithm is Newton and
            // the magnetic Picard iteration has broken down ( residual stuck far above
            // mEpsilonSwitch, so the normal promotion never fires ), switch to the Newton
            // tangent for the rest of this timestep attempt instead of cutting Δt. Returns
            // true if it escalated ( caller keeps iterating ), false otherwise ( caller
            // resets ). Escalates at most once per attempt, guarded by mNewtonEscalated.
            bool
            try_escalate_to_newton();
 
            // Thermal twin: called from the coupled stagnation
            // detector when the thermal Picard residual sits bit-flat above
            // mEpsilonSwitch2 with Newton configured but not running. Returns
            // true if it escalated. At most once per attempt ( mNewtonEscalated2 ).
            bool
            try_escalate_thermal_to_newton();
 
            void
            reset_timestep();
 
            void
            reset_thermal();
 
            void
            adjust_timestep();
 
            void
            compute_circuit_current() ;
 
            void
            compute_conditioning();
 
            // MUMPS computes its condition estimate during the SOLVE, so the
            // error analysis is armed for exactly ONE solve per timestep --
            // the first iterate -- and disarmed as soon as the value is read.
            // Leaving it armed ( the pre-2026-08-10 behaviour ) paid for the
            // analysis on every solve of every iteration and read one value.
            // The first iterate is also the better sample: it is assembled at
            // the converged previous step, so kappa is measured at a
            // comparable state in every timestep and the series is a trend.
            // All four are no-ops unless that FIELD's error-analysis flag is
            // set ( mMumpsErrorAnalysis / mMumpsErrorAnalysis2 -- NOT the
            // eigen flags, which gate a different diagnostic ) and its own
            // solver is MUMPS. They run on EVERY rank, unguarded:
            // the ICNTL setting must match across the communicator, and the
            // branch is taken identically everywhere
            void
            arm_conditioning_magnetic();
 
            void
            arm_conditioning_thermal();
 
            void
            capture_conditioning_magnetic();
 
            void
            capture_conditioning_thermal();

            void
            setup_thermal_eigen();

            void
            check_magnetic_diagnostics();

            void
            check_thermal_diagnostics();
 
//-----------------------------------------------------------------------------

            void
            check_iterative_solver_headroom(
                    DofManager * aDofMgr,
                    const real   aNonlinTol,
                    const char * aFieldName ) const ;
 
//-----------------------------------------------------------------------------

            void
            check_compression_headroom(
                    DofManager * aDofMgr,
                    const real   aNonlinTol,
                    const char * aFieldName ) const ;
 
//-----------------------------------------------------------------------------

            void
            synchronize_history_fields(
                    IWG_Timestep * aEquation,
                    DofManager   * aDofMgr,
                    const uint     aStepCount,
                    const string & aPath,
                    const char   * aEquationName );
 
            void
            print_header();
 
            void
            print_line( const real aOmega, const real aOmega2 );
 
            void
            print_line_magnetic( const real aOmega );
 
            void
            print_line_thermal( const real aOmega );
 
            void
            print_footer();

            void
            print_thermal_stall_warning();
 
            real
            get_Tmax();
 
            real
            get_Imax();

            string
            excitation_cell();
 
            real
            get_JJCmax();
 
            void
            print_physics_stats();
 
            void
            reset_dotQ();
 
            void
            collect_dotQ();
 
        };
 
        inline
        Circuit *
        Controller::circuit()
        {
            return mCircuit;
        }
 
        inline
        Kernel *
        Controller::kernel()
        {
            return mKernel;
        }
 
        inline
        Kernel *
        Controller::thermal_kernel()
        {
            return mKernel2;
        }
 
        inline
        EulerMethod
        Controller::euler_method() const
        {
            return mTimeStepping;
        }
 
    }
}
 
#endif //CL_FEM_CONTROLLER_HPP