Implement under-relaxation for the tight-coupling iterative solver

docs/source/user/api_change.rst

@@ -9,6 +9,20 @@ The changes are tabulated according to the module input file, line number, and f
 The line number corresponds to the resulting line number after all changes are implemented.
 Thus, be sure to implement each in order so that subsequent line numbers are correct.
 
+OpenFAST v5.0.x to OpenFAST  v5.1.0
+-----------------------------------
+
+Under-relaxation is introduced for the tight-coupling iterative solver to improve numerical stability, requiring two new inputs in the main OpenFAST input file.
+
+============================================= ======== ==================== ==========================================================================================================================================================================================================================================================================================================
+Added in OpenFAST `5.1.0`
+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+Module                                        Line     Flag Name            Example Value
+============================================= ======== ==================== ==========================================================================================================================================================================================================================================================================================================
+OpenFAST                                      14       AutoRelax            default  AutoRelax   - Adaptive under-relaxation for the tight-coupling iterative solver (flag) [default=true]
+OpenFAST                                      15       RelaxFactor          default  RelaxFactor - Constant or initial (if AutoRelax) under-relaxation factor for the tight-coupling iterative solver (-) [>0 and <=1; default=0.7 if AutoRelax=false; default=0.3 if AutoRelax=true]
+============================================= ======== ==================== ==========================================================================================================================================================================================================================================================================================================
+
 OpenFAST v4.2.x to OpenFAST  v5.0.0
 -----------------------------------
 

docs/source/user/glue-code/solver.rst

@@ -50,19 +50,38 @@ All solver parameters are set in the main OpenFAST input file
        maximum damping (first-order accurate).  Typical value: **0.9**.
        Reducing ``RhoInf`` below 1 damps high-frequency numerical noise at the
        cost of slightly reduced accuracy.
+   * - ``ConvTol``
+     - real
+     - Convergence tolerance.  The iteration stops when the average
+       `L2`-norm of the Newton update vector falls below this value.
+       Typical value: ``1.0e-4``.  Tighter tolerances increase
+       computational cost but may be needed for stiff problems.
    * - ``MaxConvIter``
      - integer
      - Maximum number of Newton convergence iterations per time step before the
        solver declares convergence failure.  Typical value: **20**.
        With ``ModCoupling=2`` or ``1``, a fatal error is issued on failure;
        with ``ModCoupling=3`` the Jacobian is rebuilt first and the step is
        retried before a warning is emitted.
-   * - ``ConvTol``
+   * - ``AutoRelax``
+     - bool
+     - Adaptive under-relaxation for the tight-coupling iterative solver.
+       If set to ``true``, OpenFAST will only use the user provided
+       ``RelaxFactor`` in the first iteration of each predictor or corrector
+       step. In subsequent iterations, OpenFAST will increase the relaxation
+       factor from the previous iteration by a factor of 1.2 if converging
+       (residual decreasing) or halve the relaxation factor if diverging
+       (residual increasing). The relaxation factor is bounded between a
+       minimum of 0.01 and a maximum of 0.8. If set to ``false``, OpenFAST
+       will use a constant relaxation factor given by ``RelaxFactor``.
+       Default is true.
+   * - ``RelaxFactor``
      - real
-     - Convergence tolerance.  The iteration stops when the average
-       `L2`-norm of the Newton update vector falls below this value.
-       Typical value: ``1.0e-4``.  Tighter tolerances increase
-       computational cost but may be needed for stiff problems.
+     - Constant (if ``AutoRelax`` is ``False``) or initial (if ``AutoRelax``
+       is ``True``) under-relaxation factor for the tight-coupling iterative
+       solver. Must be a positive number less than or equal to **1.0**. If
+       ``AutoRelax`` is ``False``, default value is **0.7**. If ``AutoRelax``
+       is ``True``, default value is **0.3**.
    * - ``DT_UJac``
      - real
      - Time interval (seconds) between Jacobian rebuilds when

modules/nwtc-library/src/ModVar.f90

@@ -787,15 +787,22 @@ subroutine MV_ExtrapInterp(VarAry, y, tin, y_out, tin_out, ErrStat, ErrMsg)
 
 end subroutine
 
-subroutine MV_AddDelta(VarAry, DeltaAry, DataAry)
+subroutine MV_AddDelta(VarAry, DeltaAry, DataAry, RelaxIn)
    type(ModVarType), intent(in)  :: VarAry(:)      ! Array of variables
    real(R8Ki), intent(in)        :: DeltaAry(:)    ! Array of delta values
    real(R8Ki), intent(inout)     :: DataAry(:)     ! Array to be modified
+   real(R8Ki), intent(in),optional :: RelaxIn
    integer(IntKi)                :: i, j, k
-   real(R8Ki)                    :: quat_base(3), quat_delta(3), rvec(3), dcm(3, 3)
+   real(R8Ki)                    :: quat_base(3), quat_delta(3), rvec(3), dcm(3, 3), Relax
    integer(IntKi)                :: ErrStat2
    character(ErrMsgLen)          :: ErrMsg2
 
+   if (present(RelaxIn)) then
+      Relax = RelaxIn
+   else
+      Relax = 1.0_R8Ki
+   endif
+
    ! Loop through variables
    do i = 1, size(VarAry)
       associate (iLoc => VarAry(i)%iLoc)
@@ -812,7 +819,7 @@ subroutine MV_AddDelta(VarAry, DeltaAry, DataAry)
                quat_base = DataAry(k:k + 2)
 
                ! Get rotation vector delta
-               rvec = DeltaAry(k:k + 2)
+               rvec = Relax*DeltaAry(k:k + 2)
 
                if (UseSmallRotAngles) then
                   call SmllRotTrans('linearization perturbation', rvec(1), rvec(2), rvec(3), dcm, ErrStat=ErrStat2, ErrMsg=ErrMsg2)
@@ -829,7 +836,7 @@ subroutine MV_AddDelta(VarAry, DeltaAry, DataAry)
             end do
 
          case default
-            DataAry(iLoc(1):iLoc(2)) = DataAry(iLoc(1):iLoc(2)) + DeltaAry(iLoc(1):iLoc(2))
+            DataAry(iLoc(1):iLoc(2)) = DataAry(iLoc(1):iLoc(2)) + Relax*DeltaAry(iLoc(1):iLoc(2))
          end select
       end associate
    end do

modules/openfast-library/src/FAST_Registry.txt

@@ -129,6 +129,8 @@ typedef	^	FAST_ParameterType	IntKi	ModCoupling	-	-	-	"Module coupling type {1=lo
 typedef	^	FAST_ParameterType	DbKi	RhoInf	-	-	-	"Numerical damping parameter for tight coupling generalized-alpha integrator (-) [0.0 to 1.0]"	-
 typedef	^	FAST_ParameterType	DbKi	ConvTol	-	-	-	"Convergence iteration error tolerance for tight coupling generalized alpha integrator (-)"	-
 typedef	^	FAST_ParameterType	IntKi	MaxConvIter	-	-	-	"Maximum number of convergence iterations for tight coupling generalized alpha integrator (-)"	-
+typedef ^	FAST_ParameterType	logical	AutoRelax	-	.true.	-	"Adaptive under-relaxation (flag)"	-
+typedef	^	FAST_ParameterType	DbKi	RelaxFactor	-	0.7_DbKi	-	"Constant or initial under-relaxation factor for the tight-coupling iterative solver (-) [>0 and <=1]"	-
 # Data for Jacobians:
 typedef	^	FAST_ParameterType	DbKi	DT_Ujac	-	-	-	"Time between when we need to re-calculate these Jacobians"	s
 typedef	^	FAST_ParameterType	Reki	UJacSclFact	-	-	-	"Scaling factor used to get similar magnitudes between accelerations, forces, and moments in Jacobians"	-

modules/openfast-library/src/FAST_Solver.f90

@@ -64,6 +64,12 @@ subroutine FAST_SolverInit(p_FAST, p, m, GlueModData, GlueModMaps, Turbine, ErrS
    ! Convergence tolerance
    p%ConvTol = p_FAST%ConvTol
 
+   ! Adaptive under-relaxation
+   p%AutoRelax = p_FAST%AutoRelax
+
+   ! Under-relaxation factor
+   p%RelaxFactor = p_FAST%RelaxFactor
+
    ! Solver time step
    p%h = p_FAST%DT
 
@@ -1209,14 +1215,15 @@ subroutine FAST_SolverStep(n_t_global, t_initial, p, m, GlueModData, GlueModMaps
    logical, parameter         :: IsSolve = .true.
    integer(IntKi)             :: ConvIter, CorrIter, TotalIter
    integer(IntKi)             :: NumUJac, NumCorrections
-   real(R8Ki)                 :: ConvError
+   real(R8Ki)                 :: ConvError, ConvErrorLast
    real(DbKi)                 :: t_global_next     ! next simulation time (m_FAST%t_global + p_FAST%dt)
    integer(IntKi)             :: n_t_global_next   ! n_t_global + 1
    integer(IntKi)             :: i, j, k
    integer(IntKi)             :: iMod
    integer(IntKi)             :: ConvUJac          ! Jacobian updated for convergence
    integer(IntKi)             :: MaxConvUJac       ! Max times Jacobian can be updated for convergence
    real(R8Ki)                 :: RotDiff(3, 3)
+   real(R8Ki)                 :: RelaxFactor
 
    ErrStat = ErrID_None
    ErrMsg = ''
@@ -1398,6 +1405,7 @@ subroutine FAST_SolverStep(n_t_global, t_initial, p, m, GlueModData, GlueModMaps
       !-------------------------------------------------------------------------
 
       ! Loop through convergence iterations
+      RelaxFactor = p%RelaxFactor
       do ConvIter = 0, p%MaxConvIter
 
          ! Increment total number of convergence iterations in step
@@ -1546,7 +1554,19 @@ subroutine FAST_SolverStep(n_t_global, t_initial, p, m, GlueModData, GlueModMaps
 
          ! If at least one convergence iteration has been done and
          ! the RHS norm is less than convergence tolerance, exit loop
-         if ((ConvIter > 0) .and. (ConvError < p%ConvTol)) exit
+         if (ConvIter == 0) then
+            ConvErrorLast = ConvError
+         else
+            if (ConvError < p%ConvTol) exit
+            if (p%AutoRelax) then
+               if (ConvError<ConvErrorLast) then
+                  RelaxFactor = min(1.2_R8Ki*RelaxFactor,0.80_R8Ki)
+               else ! Residual is stagnant or diverging
+                  RelaxFactor = max(0.5_R8Ki*RelaxFactor,0.01_R8Ki)
+               endif
+            endif
+            ConvErrorLast = ConvError
+         endif
 
          ! Remove load condition conditioning on input changes
          if (p%iJL(1) > 0) m%XB(p%iJL(1):p%iJL(2), 1) = m%XB(p%iJL(1):p%iJL(2), 1)*p%Scale_UJac
@@ -1562,7 +1582,7 @@ subroutine FAST_SolverStep(n_t_global, t_initial, p, m, GlueModData, GlueModMaps
          !----------------------------------------------------------------------
 
          ! Add change in inputs
-         if (p%iJU(1) > 0) call MV_AddDelta(m%Mod%Vars%u, m%XB(p%iJU(1):p%iJU(2), 1), m%Mod%Lin%u)
+         if (p%iJU(1) > 0) call MV_AddDelta(m%Mod%Vars%u, m%XB(p%iJU(1):p%iJU(2), 1), m%Mod%Lin%u, RelaxIn=RelaxFactor)
 
          !----------------------------------------------------------------------
          ! TC and Option 1: Transfer updated states and inputs to modules
@@ -1639,6 +1659,7 @@ subroutine FAST_CalcOutputsAndSolveForInputs(p, m, GlueModData, GlueModMaps, Thi
    integer(IntKi)                         :: ErrStat2
    character(ErrMsgLen)                   :: ErrMsg2
    integer(IntKi)                         :: i
+   real(R8Ki)                             :: RelaxFactor, ConvErrorLast
 
    ErrStat = ErrID_None
    ErrMsg = ''
@@ -1705,7 +1726,7 @@ subroutine FAST_CalcOutputsAndSolveForInputs(p, m, GlueModData, GlueModMaps, Thi
    !----------------------------------------------------------------------------
    ! Convergence Iterations
    !----------------------------------------------------------------------------
-
+   RelaxFactor = p%RelaxFactor
    ! Loop through convergence iterations
    do ConvIter = 0, p%MaxConvIter
 
@@ -1823,6 +1844,19 @@ subroutine FAST_CalcOutputsAndSolveForInputs(p, m, GlueModData, GlueModMaps, Thi
          exit
       end if
 
+      if (ConvIter == 0) then
+         ConvErrorLast = ConvError
+      else
+         if (p%AutoRelax) then
+            if (ConvError<ConvErrorLast) then
+               RelaxFactor = min(1.2_R8Ki*RelaxFactor,0.80_R8Ki)
+            else ! Residual is stagnant or diverging
+               RelaxFactor = max(0.5_R8Ki*RelaxFactor,0.01_R8Ki)
+            endif
+         endif
+         ConvErrorLast = ConvError
+      endif
+
       ! Remove load condition conditioning on input changes
       if (p%iUL(1) > 0) m%IO_X(p%iUL(1):p%iUL(2), 1) = m%IO_X(p%iUL(1):p%iUL(2), 1)*p%Scale_UJac
 
@@ -1831,7 +1865,7 @@ subroutine FAST_CalcOutputsAndSolveForInputs(p, m, GlueModData, GlueModMaps, Thi
       !-------------------------------------------------------------------------
 
       ! Add change in inputs
-      call MV_AddDelta(m%Mod%Vars%u, m%IO_X(:, 1), m%Mod%Lin%u)
+      call MV_AddDelta(m%Mod%Vars%u, m%IO_X(:, 1), m%Mod%Lin%u, RelaxIn=RelaxFactor)
 
       ! Transfer updated TC and Option 1 inputs to modules
       do i = 1, size(m%Mod%ModData)

modules/openfast-library/src/FAST_Subs.f90

@@ -1982,6 +1982,10 @@ SUBROUTINE ValidateInputData(p, m_FAST, ErrStat, ErrMsg)
       CALL SetErrStat( ErrID_Fatal, 'MaxIter must be at least 1.', ErrStat, ErrMsg, RoutineName )
    END IF
 
+   IF ( (p%RelaxFactor <= 0.0_DbKi) .or. (p%RelaxFactor > 1.0_DbKi) ) THEN
+      CALL SetErrStat( ErrID_Fatal, 'RelaxFactor must be positive and less than or equal to 1.', ErrStat, ErrMsg, RoutineName )
+   END IF
+
       ! Check that InputFileData%OutFmt is a valid format specifier and will fit over the column headings
    CALL ChkRealFmtStr( p%OutFmt, 'OutFmt', p%FmtWidth, ErrStat2, ErrMsg2 )
       call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
@@ -2837,7 +2841,20 @@ SUBROUTINE FAST_ReadPrimaryFile( InputFile, p, m_FAST, OverrideAbortErrLev, ErrS
    CALL ReadVar( UnIn, InputFile, p%MaxConvIter, "MaxConvIter", "Maximum number of convergence iterations "//&
                      "for tight coupling generalized alpha integrator (-)", ErrStat2, ErrMsg2, UnEc)
    if (Failed()) return
-
+
+      ! AutoRelax - Adaptive under-relaxation (flag)
+   CALL ReadVarWDefault( UnIn, InputFile, p%AutoRelax, "AutoRelax", "Adaptive under-relaxation (flag)", .true., ErrStat2, ErrMsg2, UnEc)
+   if (Failed()) return
+
+      ! RelaxFactor - Constant or initial under-relaxation factor for the iterative solver (-) [>0 and <=1]
+   if (p%AutoRelax) then
+      CALL ReadVarWDefault( UnIn, InputFile, p%RelaxFactor, "RelaxFactor", "Constant or initial under-relaxation factor for the iterative solver (-) [>0 and <=1]", 0.3_R8Ki, ErrStat2, ErrMsg2, UnEc)
+      if (Failed()) return
+   else
+      CALL ReadVarWDefault( UnIn, InputFile, p%RelaxFactor, "RelaxFactor", "Constant or initial under-relaxation factor for the iterative solver (-) [>0 and <=1]", 0.7_R8Ki, ErrStat2, ErrMsg2, UnEc)
+      if (Failed()) return
+   endif
+
       ! DT_UJac - Time between calls to get Jacobians (s)
    CALL ReadVar( UnIn, InputFile, p%DT_UJac, "DT_UJac", "Time between calls to get Jacobians (s)", ErrStat2, ErrMsg2, UnEc)
    if (Failed()) return

modules/openfast-library/src/FAST_Types.f90

@@ -145,6 +145,8 @@ MODULE FAST_Types
     REAL(DbKi)  :: RhoInf = 0.0_R8Ki      !< Numerical damping parameter for tight coupling generalized-alpha integrator (-) [0.0 to 1.0] [-]
     REAL(DbKi)  :: ConvTol = 0.0_R8Ki      !< Convergence iteration error tolerance for tight coupling generalized alpha integrator (-) [-]
     INTEGER(IntKi)  :: MaxConvIter = 0_IntKi      !< Maximum number of convergence iterations for tight coupling generalized alpha integrator (-) [-]
+    LOGICAL  :: AutoRelax = .true.      !< Adaptive under-relaxation (flag) [-]
+    REAL(DbKi)  :: RelaxFactor = 0.7_DbKi      !< Constant or initial under-relaxation factor for the tight-coupling iterative solver (-) [>0 and <=1] [-]
     REAL(DbKi)  :: DT_Ujac = 0.0_R8Ki      !< Time between when we need to re-calculate these Jacobians [s]
     REAL(ReKi)  :: UJacSclFact = 0.0_ReKi      !< Scaling factor used to get similar magnitudes between accelerations, forces, and moments in Jacobians [-]
     INTEGER(IntKi) , DIMENSION(1:9)  :: SizeJac_Opt1 = 0_IntKi      !< (1)=size of matrix; (2)=size of ED portion; (3)=size of SD portion [2 meshes]; (4)=size of HD portion; (5)=size of BD portion blade 1; (6)=size of BD portion blade 2; (7)=size of BD portion blade 3; (8)=size of Orca portion; (9)=size of ExtPtfm portion; [-]
@@ -1162,6 +1164,8 @@ subroutine FAST_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg)
    DstParamData%RhoInf = SrcParamData%RhoInf
    DstParamData%ConvTol = SrcParamData%ConvTol
    DstParamData%MaxConvIter = SrcParamData%MaxConvIter
+   DstParamData%AutoRelax = SrcParamData%AutoRelax
+   DstParamData%RelaxFactor = SrcParamData%RelaxFactor
    DstParamData%DT_Ujac = SrcParamData%DT_Ujac
    DstParamData%UJacSclFact = SrcParamData%UJacSclFact
    DstParamData%SizeJac_Opt1 = SrcParamData%SizeJac_Opt1
@@ -1391,6 +1395,8 @@ subroutine FAST_PackParam(RF, Indata)
    call RegPack(RF, InData%RhoInf)
    call RegPack(RF, InData%ConvTol)
    call RegPack(RF, InData%MaxConvIter)
+   call RegPack(RF, InData%AutoRelax)
+   call RegPack(RF, InData%RelaxFactor)
    call RegPack(RF, InData%DT_Ujac)
    call RegPack(RF, InData%UJacSclFact)
    call RegPack(RF, InData%SizeJac_Opt1)
@@ -1513,6 +1519,8 @@ subroutine FAST_UnPackParam(RF, OutData)
    call RegUnpack(RF, OutData%RhoInf); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%ConvTol); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%MaxConvIter); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%AutoRelax); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%RelaxFactor); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%DT_Ujac); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%UJacSclFact); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%SizeJac_Opt1); if (RegCheckErr(RF, RoutineName)) return

modules/openfast-library/src/Glue_Registry.txt

@@ -74,6 +74,8 @@ typedef ^ ^                     R8Ki            ConvTol                 -  -  -
 typedef ^ ^                     IntKi           ModCoupling             -  -  -     "Module coupling method {1=loose; 2=tight with fixed Jacobian updates (DT_UJac); 3=tight with automatic Jacobian updates}" -
 typedef ^ ^                     IntKi           NumCrctn                -  -  -     "" -
 typedef ^ ^                     IntKi           MaxConvIter             -  -  -     "" -
+typedef ^ ^                     logical         AutoRelax               -  .true.    -     "Adaptive under-relaxation (flag)" -
+typedef ^ ^                     R8Ki            RelaxFactor             -  0.7_R8Ki  -     "Constant or initial under-relaxation factor for the tight-coupling iterative solver" -
 typedef ^ ^                     IntKi           NIter_UJac              -  -  -     "Number of solution iterations between updating the Jacobian" -
 typedef ^ ^                     IntKi           NStep_UJac              -  -  -     "Number of global time steps between updating the Jacobian" -
 typedef ^ ^                     R8Ki            Scale_UJac              -  -  -     "Jacobian load scaling factor" -