From a9760727fda0713cf0984b6b40f647719b5860f6 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 10 Apr 2026 14:59:57 +0000 Subject: [PATCH 01/14] Defined new variables of Syrope model --- modules/moordyn/src/MoorDyn_Registry.txt | 13 ++++ modules/moordyn/src/MoorDyn_Types.f90 | 80 ++++++++++++++++++++++++ 2 files changed, 93 insertions(+) diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 210b98cf0b..5bafb79e61 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -83,6 +83,9 @@ typedef ^ ^ DbKi dampYs {MD_MaxNCoef} typedef ^ ^ IntKi nEIpoints - - - "number of values in bending stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi bstiffXs {MD_MaxNCoef} - - "x array for stress-strain lookup table (up to MD_MaxNCoef)" typedef ^ ^ DbKi bstiffYs {MD_MaxNCoef} - - "y array for stress-strain lookup table" +typedef ^ ^ IntKi syropeWCForm - - - "Syrope working curve formula (1 = LINEAR, 2 = QUADRATIC, 3 = EXP)" +typedef ^ ^ DbKi syropeWCK1 - - - "Coefficient for the Syrope working curve formula" +typedef ^ ^ DbKi syropeWCK2 - - - "Coefficient for the Syrope working curve formula" # rod properties from rod dictionary input typedef ^ MD_RodProp IntKi IdNum - - - "integer identifier of this set of rod properties" @@ -258,6 +261,10 @@ typedef ^ ^ DbKi cF - typedef ^ ^ IntKi nEApoints - - - "number of values in stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi stiffXs {MD_MaxNCoef} - - "x array for stress-strain lookup table (up to MD_MaxNCoef)" typedef ^ ^ DbKi stiffYs {MD_MaxNCoef} - - "y array for stress-strain lookup table" +typedef ^ ^ DbKi stiffZs {MD_MaxNCoef} - - "z array for stress-strain original working curve (Syrope model only)" +typedef ^ ^ DbKi stiffxs_ {MD_MaxNCoef} - - "auxiliary x array for Syrope working curve" +typedef ^ ^ DbKi stiffys_ {MD_MaxNCoef} - - "auxiliary y array for Syrope working curve" +typedef ^ ^ DbKi stiffzs_ {MD_MaxNCoef} - - "auxiliary z array for Syrope working curve" typedef ^ ^ IntKi nBApoints - - - "number of values in stress-strainrate lookup table (0 means using constant c)" typedef ^ ^ DbKi dampXs {MD_MaxNCoef} - - "x array for stress-strainrate lookup table (up to MD_MaxNCoef)" typedef ^ ^ DbKi dampYs {MD_MaxNCoef} - - "y array for stress-strainrate lookup table" @@ -283,6 +290,11 @@ typedef ^ ^ DbKi zeta {:} typedef ^ ^ DbKi PDyn {:} - - "water dynamic pressure at node" "[Pa]" typedef ^ ^ DbKi T {:}{:} - - "segment tension vectors" "[N]" typedef ^ ^ DbKi Td {:}{:} - - "segment internal damping force vectors" "[N]" +typedef ^ ^ IntKi syropeWCForm - - - "Syrope working curve formula (1 = LINEAR, 2 = QUADRATIC, 3 = EXP)" +typedef ^ ^ DbKi syropeWCK1 - - - "Coefficient for the Syrope working curve formula" +typedef ^ ^ DbKi syropeWCK2 - - - "Coefficient for the Syrope working curve formula" +typedef ^ ^ DbKi Tmax {:} - - "segment preceding highest mean tensions" "[N]" +typedef ^ ^ DbKi Tmean {:} - - "segment mean tensions" "[N]" typedef ^ ^ DbKi W {:}{:} - - "weight/buoyancy vectors" "[N]" typedef ^ ^ DbKi Dp {:}{:} - - "node drag (transverse)" "[N]" typedef ^ ^ DbKi Dq {:}{:} - - "node drag (axial)" "[N]" @@ -366,6 +378,7 @@ typedef ^ ^ IntKi nPointsExtra - typedef ^ ^ IntKi nBodies - 0 - "number of Body objects" "" typedef ^ ^ IntKi nRods - 0 - "number of Rod objects" "" typedef ^ ^ IntKi nLines - 0 - "number of Line objects" "" +typedef ^ ^ IntKi nSyropeLineICs - 0 - "number of Syrope Line initial conditions" "" typedef ^ ^ IntKi nExtLds - 0 - "number of external loads or damping" "" typedef ^ ^ IntKi nCtrlChans - 0 - "number of distinct control channels specified for use as inputs" "" typedef ^ ^ IntKi nFails - 0 - "number of failure conditions" "" diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index 86703d7aa7..5dc8e34859 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -96,6 +96,9 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: nEIpoints = 0_IntKi !< number of values in bending stress-strain lookup table (0 means using constant E) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: bstiffXs = 0.0_R8Ki !< x array for stress-strain lookup table (up to MD_MaxNCoef) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: bstiffYs = 0.0_R8Ki !< y array for stress-strain lookup table [-] + INTEGER(IntKi) :: syropeWCForm = 0_IntKi !< Syrope working curve formula (1 = LINEAR, 2 = QUADRATIC, 3 = EXP) [-] + REAL(DbKi) :: syropeWCK1 = 0.0_R8Ki !< Coefficient for the Syrope working curve formula [-] + REAL(DbKi) :: syropeWCK2 = 0.0_R8Ki !< Coefficient for the Syrope working curve formula [-] END TYPE MD_LineProp ! ======================= ! ========= MD_RodProp ======= @@ -279,6 +282,10 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: nEApoints = 0_IntKi !< number of values in stress-strain lookup table (0 means using constant E) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffXs = 0.0_R8Ki !< x array for stress-strain lookup table (up to MD_MaxNCoef) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffYs = 0.0_R8Ki !< y array for stress-strain lookup table [-] + REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffZs = 0.0_R8Ki !< z array for stress-strain original working curve (Syrope model only) [-] + REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffxs_ = 0.0_R8Ki !< auxiliary x array for Syrope working curve [-] + REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffys_ = 0.0_R8Ki !< auxiliary y array for Syrope working curve [-] + REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: stiffzs_ = 0.0_R8Ki !< auxiliary z array for Syrope working curve [-] INTEGER(IntKi) :: nBApoints = 0_IntKi !< number of values in stress-strainrate lookup table (0 means using constant c) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: dampXs = 0.0_R8Ki !< x array for stress-strainrate lookup table (up to MD_MaxNCoef) [-] REAL(DbKi) , DIMENSION(1:MD_MaxNCoef) :: dampYs = 0.0_R8Ki !< y array for stress-strainrate lookup table [-] @@ -304,6 +311,11 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: PDyn !< water dynamic pressure at node [[Pa]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: T !< segment tension vectors [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Td !< segment internal damping force vectors [[N]] + INTEGER(IntKi) :: syropeWCForm = 0_IntKi !< Syrope working curve formula (1 = LINEAR, 2 = QUADRATIC, 3 = EXP) [-] + REAL(DbKi) :: syropeWCK1 = 0.0_R8Ki !< Coefficient for the Syrope working curve formula [-] + REAL(DbKi) :: syropeWCK2 = 0.0_R8Ki !< Coefficient for the Syrope working curve formula [-] + REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: Tmax !< segment preceding highest mean tensions [[N]] + REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: Tmean !< segment mean tensions [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: W !< weight/buoyancy vectors [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Dp !< node drag (transverse) [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Dq !< node drag (axial) [[N]] @@ -403,6 +415,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: nBodies = 0 !< number of Body objects [] INTEGER(IntKi) :: nRods = 0 !< number of Rod objects [] INTEGER(IntKi) :: nLines = 0 !< number of Line objects [] + INTEGER(IntKi) :: nSyropeLineICs = 0 !< number of Syrope Line initial conditions [] INTEGER(IntKi) :: nExtLds = 0 !< number of external loads or damping [] INTEGER(IntKi) :: nCtrlChans = 0 !< number of distinct control channels specified for use as inputs [] INTEGER(IntKi) :: nFails = 0 !< number of failure conditions [] @@ -796,6 +809,9 @@ subroutine MD_CopyLineProp(SrcLinePropData, DstLinePropData, CtrlCode, ErrStat, DstLinePropData%nEIpoints = SrcLinePropData%nEIpoints DstLinePropData%bstiffXs = SrcLinePropData%bstiffXs DstLinePropData%bstiffYs = SrcLinePropData%bstiffYs + DstLinePropData%syropeWCForm = SrcLinePropData%syropeWCForm + DstLinePropData%syropeWCK1 = SrcLinePropData%syropeWCK1 + DstLinePropData%syropeWCK2 = SrcLinePropData%syropeWCK2 end subroutine subroutine MD_DestroyLineProp(LinePropData, ErrStat, ErrMsg) @@ -840,6 +856,9 @@ subroutine MD_PackLineProp(RF, Indata) call RegPack(RF, InData%nEIpoints) call RegPack(RF, InData%bstiffXs) call RegPack(RF, InData%bstiffYs) + call RegPack(RF, InData%syropeWCForm) + call RegPack(RF, InData%syropeWCK1) + call RegPack(RF, InData%syropeWCK2) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -876,6 +895,9 @@ subroutine MD_UnPackLineProp(RF, OutData) call RegUnpack(RF, OutData%nEIpoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%bstiffXs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%bstiffYs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCForm); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCK1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCK2); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine MD_CopyRodProp(SrcRodPropData, DstRodPropData, CtrlCode, ErrStat, ErrMsg) @@ -1790,6 +1812,10 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) DstLineData%nEApoints = SrcLineData%nEApoints DstLineData%stiffXs = SrcLineData%stiffXs DstLineData%stiffYs = SrcLineData%stiffYs + DstLineData%stiffZs = SrcLineData%stiffZs + DstLineData%stiffxs_ = SrcLineData%stiffxs_ + DstLineData%stiffys_ = SrcLineData%stiffys_ + DstLineData%stiffzs_ = SrcLineData%stiffzs_ DstLineData%nBApoints = SrcLineData%nBApoints DstLineData%dampXs = SrcLineData%dampXs DstLineData%dampYs = SrcLineData%dampYs @@ -2013,6 +2039,33 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) end if DstLineData%Td = SrcLineData%Td end if + DstLineData%syropeWCForm = SrcLineData%syropeWCForm + DstLineData%syropeWCK1 = SrcLineData%syropeWCK1 + DstLineData%syropeWCK2 = SrcLineData%syropeWCK2 + if (allocated(SrcLineData%Tmax)) then + LB(1:1) = lbound(SrcLineData%Tmax) + UB(1:1) = ubound(SrcLineData%Tmax) + if (.not. allocated(DstLineData%Tmax)) then + allocate(DstLineData%Tmax(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%Tmax.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%Tmax = SrcLineData%Tmax + end if + if (allocated(SrcLineData%Tmean)) then + LB(1:1) = lbound(SrcLineData%Tmean) + UB(1:1) = ubound(SrcLineData%Tmean) + if (.not. allocated(DstLineData%Tmean)) then + allocate(DstLineData%Tmean(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%Tmean.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%Tmean = SrcLineData%Tmean + end if if (allocated(SrcLineData%W)) then LB(1:2) = lbound(SrcLineData%W) UB(1:2) = ubound(SrcLineData%W) @@ -2274,6 +2327,12 @@ subroutine MD_DestroyLine(LineData, ErrStat, ErrMsg) if (allocated(LineData%Td)) then deallocate(LineData%Td) end if + if (allocated(LineData%Tmax)) then + deallocate(LineData%Tmax) + end if + if (allocated(LineData%Tmean)) then + deallocate(LineData%Tmean) + end if if (allocated(LineData%W)) then deallocate(LineData%W) end if @@ -2359,6 +2418,10 @@ subroutine MD_PackLine(RF, Indata) call RegPack(RF, InData%nEApoints) call RegPack(RF, InData%stiffXs) call RegPack(RF, InData%stiffYs) + call RegPack(RF, InData%stiffZs) + call RegPack(RF, InData%stiffxs_) + call RegPack(RF, InData%stiffys_) + call RegPack(RF, InData%stiffzs_) call RegPack(RF, InData%nBApoints) call RegPack(RF, InData%dampXs) call RegPack(RF, InData%dampYs) @@ -2384,6 +2447,11 @@ subroutine MD_PackLine(RF, Indata) call RegPackAlloc(RF, InData%PDyn) call RegPackAlloc(RF, InData%T) call RegPackAlloc(RF, InData%Td) + call RegPack(RF, InData%syropeWCForm) + call RegPack(RF, InData%syropeWCK1) + call RegPack(RF, InData%syropeWCK2) + call RegPackAlloc(RF, InData%Tmax) + call RegPackAlloc(RF, InData%Tmean) call RegPackAlloc(RF, InData%W) call RegPackAlloc(RF, InData%Dp) call RegPackAlloc(RF, InData%Dq) @@ -2447,6 +2515,10 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpack(RF, OutData%nEApoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%stiffXs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%stiffYs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stiffZs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stiffxs_); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stiffys_); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%stiffzs_); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nBApoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%dampXs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%dampYs); if (RegCheckErr(RF, RoutineName)) return @@ -2472,6 +2544,11 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpackAlloc(RF, OutData%PDyn); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%T); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Td); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCForm); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCK1); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%syropeWCK2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Tmax); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Tmean); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%W); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Dp); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Dq); if (RegCheckErr(RF, RoutineName)) return @@ -3017,6 +3094,7 @@ subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%nBodies = SrcParamData%nBodies DstParamData%nRods = SrcParamData%nRods DstParamData%nLines = SrcParamData%nLines + DstParamData%nSyropeLineICs = SrcParamData%nSyropeLineICs DstParamData%nExtLds = SrcParamData%nExtLds DstParamData%nCtrlChans = SrcParamData%nCtrlChans DstParamData%nFails = SrcParamData%nFails @@ -3413,6 +3491,7 @@ subroutine MD_PackParam(RF, Indata) call RegPack(RF, InData%nBodies) call RegPack(RF, InData%nRods) call RegPack(RF, InData%nLines) + call RegPack(RF, InData%nSyropeLineICs) call RegPack(RF, InData%nExtLds) call RegPack(RF, InData%nCtrlChans) call RegPack(RF, InData%nFails) @@ -3521,6 +3600,7 @@ subroutine MD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%nBodies); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nRods); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nLines); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%nSyropeLineICs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nExtLds); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nCtrlChans); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nFails); if (RegCheckErr(RF, RoutineName)) return From 4244f3339599d9f42a3ec24ce6db37cd82734b34 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 10 Apr 2026 15:00:14 +0000 Subject: [PATCH 02/14] 1D interpolation from x and y arrays --- modules/moordyn/src/MoorDyn_Misc.f90 | 70 ++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index 3c661f1766..3f39a1a7a7 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -53,6 +53,7 @@ MODULE MoorDyn_Misc PUBLIC :: calculate4Dinterpolation PUBLIC :: calculate3Dinterpolation PUBLIC :: calculate2Dinterpolation + PUBLIC :: calculate1DinterpolationXY PUBLIC :: getDepthFromBathymetry @@ -824,6 +825,75 @@ SUBROUTINE calculate1Dinterpolation(f, ix0, fx, c) c = c0*(1.0-fx) + c1*fx END SUBROUTINE calculate1Dinterpolation + SUBROUTINE calculate1DinterpolationXY(x, y, xi, yi, ErrStat, ErrMsg) + REAL(DbKi), INTENT(IN) :: x(:) + REAL(DbKi), INTENT(IN) :: y(:) + REAL(DbKi), INTENT(IN) :: xi + REAL(DbKi), INTENT(OUT) :: yi + INTEGER(IntKi), INTENT(OUT) :: ErrStat + CHARACTER(*), INTENT(OUT) :: ErrMsg + CHARACTER(*), PARAMETER :: RoutineName = 'calculate1DinterpolationXY' + + INTEGER(IntKi) :: n + INTEGER(IntKi) :: lo, hi, mid + INTEGER(IntKi) :: ix0 + REAL(DbKi) :: fx, dx + + ErrStat = ErrID_None + ErrMsg = '' + yi = 0.0_DbKi + + n = SIZE(x) + + IF (n /= SIZE(y)) THEN + CALL SetErrStat(ErrID_Fatal, 'x and y must have the same size.', ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + IF (n < 2) THEN + CALL SetErrStat(ErrID_Fatal, 'At least two points are required for interpolation.', ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + IF (ANY(x(2:n) <= x(1:n-1))) THEN + CALL SetErrStat(ErrID_Fatal, 'x array must be strictly increasing.', ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + ! Clamp to endpoints + IF (xi <= x(1)) THEN + yi = y(1) + RETURN + ELSE IF (xi >= x(n)) THEN + yi = y(n) + RETURN + END IF + + ! Binary search for interval + lo = 1 + hi = n + DO WHILE (hi - lo > 1) + mid = lo + (hi - lo) / 2 + IF (x(mid) <= xi) THEN + lo = mid + ELSE + hi = mid + END IF + END DO + + ix0 = lo + dx = x(ix0+1) - x(ix0) + + IF (ABS(dx) <= TINY(1.0_DbKi)) THEN + CALL SetErrStat(ErrID_Fatal, 'Zero or too-small x interval encountered.', ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + fx = (xi - x(ix0)) / dx + + CALL calculate1Dinterpolation(y, ix0, fx, yi) + + END SUBROUTINE calculate1DinterpolationXY From fc7b7387c74320489d28a366c97ec2203bdc5ea8 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 10 Apr 2026 15:00:30 +0000 Subject: [PATCH 03/14] Read Syrope settings --- modules/moordyn/src/MoorDyn.f90 | 329 ++++++++++++++++++++++++++++++-- 1 file changed, 313 insertions(+), 16 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 67191d4609..7d5a468553 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -127,8 +127,11 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CHARACTER(40) :: TempString3 ! CHARACTER(40) :: TempString4 ! CHARACTER(40) :: TempString5 ! + CHARACTER(40) :: TempString6 ! CHARACTER(40) :: TempStrings(6) ! Array of 6 strings used when parsing comma-separated items ! CHARACTER(1024) :: FileName ! + LOGICAL :: hasSyrope = .FALSE. ! flag for if a line type has a Syrope model specified + LOGICAL :: hasSyropeICs = .FALSE. ! flag for whether any line type uses Syrope model REAL(DbKi) :: depth ! local water depth interpolated from bathymetry grid [m] @@ -155,6 +158,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er INTEGER :: QuoteCh ! Character position. CHARACTER(*), PARAMETER :: RoutineName = 'MD_Init' + ! for Syrope working curves + CHARACTER(40) :: owcPath ! + CHARACTER(40) :: wcFormula ! + REAL(DbKi) :: wcK1 ! + REAL(DbKi) :: wcK2 ! + + ! for Syrope line initial conditions + REAL(DbKi) :: Tmax0 + REAL(DbKi) :: Tmean0 + ! Initialize Err stat ErrStat = ErrID_None @@ -371,6 +384,21 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er p%nLines = p%nLines + 1 Line = NextLine(i) END DO + + else if ((INDEX(Line, "SYROPE IC") > 0) ) then ! if Syrope Line initial conditions + + IF (wordy > 1) print *, " Reading Syrope line initial conditions: "; + + ! skip following two lines (label line and unit line) + Line = NextLine(i) + Line = NextLine(i) + + ! find how many elements of this type there are + Line = NextLine(i) + DO while (INDEX(Line, "---") == 0) ! while we DON'T find another header line + p%nSyropeLineICs = p%nSyropeLineICs + 1 + Line = NextLine(i) + END DO else if (INDEX(Line, "EXTERNAL LOADS") > 0) then ! if external load and damping header @@ -698,25 +726,76 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! process stiffness coefficients CALL SplitByBars(tempString1, N, tempStrings) - if (N > 3) then - CALL SetErrStat( ErrID_Fatal, 'A line type EA entry can have at most 3 (bar-separated) values.', ErrStat, ErrMsg, RoutineName ) - CALL CleanUp() - RETURN - else if (N==3) then ! visco-elastic case, load dependent dynamic stiffness! - m%LineTypeList(l)%ElasticMod = 3 + + ! check if the first entry follows the pattern "SYROPE:" for a Syrope model + ! if so, enable the Syrope model for this line type + tempString5 = adjustl(tempStrings(1)) + ! check length of tempString5 to out-of-bounds error + hasSyrope = (len(tempString5) > 7) .and. (tempString5(1:7) == 'SYROPE:') + + if (hasSyrope) then + if (N /= 3) then + CALL SetErrStat( ErrID_Fatal, 'A line type EA entry for a SYROPE model must have exactly 3 (bar-separated) values.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + end if + + m%LineTypeList(l)%ElasticMod = 4 read(tempStrings(2), *) m%LineTypeList(l)%alphaMBL read(tempStrings(3), *) m%LineTypeList(l)%vbeta - else if (N==2) then ! visco-elastic case, constant dynamic stiffness! - m%LineTypeList(l)%ElasticMod = 2 - read(tempStrings(2), *) m%LineTypeList(l)%EA_D - else - m%LineTypeList(l)%ElasticMod = 1 ! normal case - end if - ! get the regular/static coefficient or relation in all cases (can be from a lookup table) - CALL getCoefficientOrCurve(tempStrings(1), m%LineTypeList(l)%EA, & + + tempString6 = adjustl(tempString5(8:)) ! get the filepath of the SYROPE curves + + CALL MD_ReadSyropeWorkingCurves(tempString6, owcPath, wcFormula, wcK1, wcK2, ErrStat2, ErrMsg2) + + if (ErrStat2 /= ErrID_None ) then + CALL SetErrStat( ErrID_Fatal, 'Failed to read SYROPE working curves for line type '//trim(Num2LStr(l))//'. Check error message for details.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + end if + + if (wordy > 0) print *, " Read Syrope working curve for Line type ", trim(Num2LStr(l)) + + if (trim(wcFormula) == 'LINEAR') then + m%LineTypeList(l)%syropeWCForm = 1 + else if (trim(wcFormula) == 'QUADRATIC') then + m%LineTypeList(l)%syropeWCForm = 2 + else if (trim(wcFormula) == 'EXP') then + m%LineTypeList(l)%syropeWCForm = 3 + end if + m%LineTypeList(l)%syropeWCK1 = wcK1 + m%LineTypeList(l)%syropeWCK2 = wcK2 + + ! get the original working curve from the loopup table owcPath + CALL getCoefficientOrCurve(owcPath, m%LineTypeList(l)%EA, & m%LineTypeList(l)%nEApoints, & m%LineTypeList(l)%stiffXs, & m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2) + + hasSyrope = .false. ! reset hasSyrope flag for later use when processing + + else ! process other visco-elastic or normal cases + + if (N > 3) then + CALL SetErrStat( ErrID_Fatal, 'A line type EA entry can have at most 3 (bar-separated) values.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + else if (N==3) then ! visco-elastic case, load dependent dynamic stiffness! + m%LineTypeList(l)%ElasticMod = 3 + read(tempStrings(2), *) m%LineTypeList(l)%alphaMBL + read(tempStrings(3), *) m%LineTypeList(l)%vbeta + else if (N==2) then ! visco-elastic case, constant dynamic stiffness! + m%LineTypeList(l)%ElasticMod = 2 + read(tempStrings(2), *) m%LineTypeList(l)%EA_D + else + m%LineTypeList(l)%ElasticMod = 1 ! normal case + end if + ! get the regular/static coefficient or relation in all cases (can be from a lookup table) + CALL getCoefficientOrCurve(tempStrings(1), m%LineTypeList(l)%EA, & + m%LineTypeList(l)%nEApoints, & + m%LineTypeList(l)%stiffXs, & + m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2) + end if ! process damping coefficients @@ -1688,6 +1767,60 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF (wordy > 0) print *, "Set up Line", l, "of type", m%LineList(l)%PropsIdNum END DO ! l = 1,p%nLines + + !------------------------------------------------------------------------------------------- + else if ((INDEX(Line, "SYROPE IC") > 0) ) then ! if Syrope initial conditions header + hasSyropeICs = .false. + do l = 1, p%nLineTypes + if (m%LineTypeList(l)%ElasticMod == 4) then + hasSyropeICs = .true. + exit + end if + end do + + if (.not. hasSyropeICs) then + CALL SetErrStat( ErrID_Warn, 'No line type with Syrope model is defined, but SYROPE IC section is present.' , ErrStat, ErrMsg, RoutineName ) + end if + + if (wordy > 0) print *, " Reading Syrope line initial inputs" + + ! skip following two lines (label line and unit line) + Line = NextLine(i) + Line = NextLine(i) + + DO l = 1, p%nSyropeLineICs + + !read into a line + Line = NextLine(i) + + ! count commas to determine how many line IDs specified for this IC + N = count(transfer(Line, 'a', len(Line)) == ",") + 1 ! number of line IDs given + + ! check for correct number of columns in current line + IF ( CountWords( Line ) /= N+2 ) THEN + CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be 5 columns.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + ! parse out entries: SyropeLineID Tmax Tmean + IF (ErrStat2 == 0) THEN + READ(Line,*,IOSTAT=ErrStat2) TempIDnums(1:N), Tmax0, Tmean0 + END IF + + DO il = 1, N + if (TempIDnums(il) <= p%nLines) then ! ensure line ID is in range + DO j = 1, m%LineList(TempIDnums(il))%N + m%LineList(TempIDnums(il))%Tmax(j) = Tmax0 + m%LineList(TempIDnums(il))%Tmean(j) = Tmean0 + END DO + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse Syrope line initial conditions '//trim(Num2LStr(l))//'. Line number '//TRIM(Int2LStr(TempIDnums(il)))//' out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + END DO + END DO !------------------------------------------------------------------------------------------- else if (INDEX(Line, "EXTERNAL LOADS") > 0) then ! if external load header @@ -2668,8 +2801,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! print *, m%LineList(l)%r(:,I) END DO - ! if using viscoelastic model, initialize the internal states - if (m%LineList(l)%ElasticMod == 2) then + ! if using viscoelastic model or Syrope model, initialize the internal states + if (m%LineList(l)%ElasticMod == 2 .or. m%LineList(l)%ElasticMod == 4) then do I = 1,N x%states(m%LineStateIs1(l) + 6*N-6 + I-1) = m%LineList(l)%dl_1(I) ! should be zero end do @@ -4914,4 +5047,168 @@ SUBROUTINE MD_JacobianPConstrState( t, u, p, x, xd, z, OtherState, y, m, ErrStat END IF END SUBROUTINE MD_JacobianPConstrState +SUBROUTINE MD_ReadSyropeWorkingCurves(inputString, owcPath, wcFormula, k1, k2, ErrStat3, ErrMsg3) + + CHARACTER(*), INTENT(IN) :: inputString + CHARACTER(*), INTENT(OUT) :: owcPath + CHARACTER(*), INTENT(OUT) :: wcFormula + REAL(DbKi), INTENT(OUT) :: k1 + REAL(DbKi), INTENT(OUT) :: k2 + + INTEGER(IntKi), INTENT(OUT) :: ErrStat3 + CHARACTER(*), INTENT(OUT) :: ErrMsg3 + + ! local variables + TYPE(FileInfoType) :: FileInfo + CHARACTER(1024) :: NextLine + CHARACTER(64) :: OptString + CHARACTER(256) :: OptValue + LOGICAL :: FoundOWC = .false. + LOGICAL :: FoundWCFormula = .false. + LOGICAL :: FoundK1 = .false. + LOGICAL :: FoundK2 = .false. + CHARACTER(*), PARAMETER :: RoutineName = 'MD_ReadSyropeWorkingCurves' + + INTEGER(IntKi) :: i, ios + INTEGER(IntKi) :: ErrStat4 + CHARACTER(1024) :: ErrMsg4 + + ! initialize outputs + owcPath = '' + wcFormula = '' + k1 = 0.0_DbKi + k2 = 0.0_DbKi + ErrStat3 = ErrID_None + ErrMsg3 = '' + + ! read file + + CALL ProcessComFile(inputString, FileInfo, ErrStat4, ErrMsg4) + IF (ErrStat4 /= ErrID_None) THEN + CALL SetErrStat(ErrID_Fatal, & + 'Error processing Syrope working curve file ' // TRIM(inputString) // ': ' // TRIM(ErrMsg4), & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + CALL WrScr(' Loading SYROPE working curves from ' // TRIM(inputString)) + + ! parse lines + DO i = 1, FileInfo%NumLines + + NextLine = TRIM(FileInfo%Lines(i)) + IF (LEN_TRIM(NextLine) == 0) CYCLE + + READ(NextLine, *, IOSTAT=ios) OptValue, OptString + IF (ios /= 0) THEN + CALL SetErrStat(ErrID_Fatal, & + 'Failed to read options in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + CALL Conv2UC(OptString) + + IF (OptString == 'OWC') THEN + owcPath = TRIM(OptValue) + FoundOWC = .true. + ELSE IF (OptString == 'WCTYPE') THEN + wcFormula = TRIM(OptValue) + FoundWCFormula = .true. + ELSE IF (OptString == 'K1') THEN + READ(OptValue, *, IOSTAT=ios) k1 + IF (ios /= 0) THEN + CALL SetErrStat(ErrID_Fatal, & + 'Invalid K1 value in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + FoundK1 = .true. + ELSE IF (OptString == 'K2') THEN + READ(OptValue, *, IOSTAT=ios) k2 + IF (ios /= 0) THEN + CALL SetErrStat(ErrID_Fatal, & + 'Invalid K2 value in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + FoundK2 = .true. + ELSE + CALL SetErrStat(ErrID_Warn, & + 'Unknown keyword "' // TRIM(OptString) // '" in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + END IF + + END DO + + ! required checks + IF (.NOT. FoundOWC) THEN + CALL SetErrStat(ErrID_Fatal, & + 'OWC option not found in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + IF (.NOT. FoundWCFormula) THEN + CALL SetErrStat(ErrID_Fatal, & + 'WCTYPE option not found in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + IF (.NOT. FoundK1) THEN + CALL SetErrStat(ErrID_Fatal, & + 'K1 option not found in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + IF (.NOT. FoundK2) THEN + CALL SetErrStat(ErrID_Fatal, & + 'K2 option not found in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + ! validate formula + CALL Conv2UC(wcFormula) + + IF (wcFormula == 'LINEAR') THEN + + IF (k1 < 0.0_DbKi) THEN + CALL SetErrStat(ErrID_Fatal, & + 'LINEAR working curve requires k1 >= 0.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + ELSE IF (wcFormula == 'QUADRATIC') THEN + + IF (k1 < 0.0_DbKi .OR. k1 >= 1.0_DbKi .OR. k2 <= 0.0_DbKi .OR. k2 > 1.0_DbKi) THEN + CALL SetErrStat(ErrID_Fatal, & + 'QUADRATIC working curve requires 0 <= k1 < 1 and 0 < k2 <= 1.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + ELSE IF (wcFormula == 'EXP') THEN + + IF (k1 < 0.0_DbKi .OR. k1 >= 1.0_DbKi .OR. k2 <= 0.0_DbKi) THEN + CALL SetErrStat(ErrID_Fatal, & + 'EXP working curve requires 0 <= k1 < 1 and k2 > 0.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + END IF + + ELSE + + CALL SetErrStat(ErrID_Fatal, & + 'Invalid WCTYPE value "' // TRIM(wcFormula) // '" in Syrope working curve file ' // TRIM(inputString) // '.', & + ErrStat3, ErrMsg3, RoutineName) + RETURN + + END IF + +END SUBROUTINE MD_ReadSyropeWorkingCurves + END MODULE MoorDyn From 7fe99b98fff15a7b388b55b74f97f1cfc2d44c7f Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Fri, 10 Apr 2026 15:00:47 +0000 Subject: [PATCH 04/14] Implementation for the Syrope model --- modules/moordyn/src/MoorDyn_Line.f90 | 222 ++++++++++++++++++++++++++- 1 file changed, 217 insertions(+), 5 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 50e93350bf..3e4c567ed6 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -90,9 +90,9 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! copy over elasticity data Line%ElasticMod = LineProp%ElasticMod - if (Line%ElasticMod > 3) then + if (Line%ElasticMod > 4) then ErrStat = ErrID_Fatal - ErrMsg = "Line ElasticMod > 3. This is not possible." + ErrMsg = "Line ElasticMod > 4. This is not possible." RETURN endif @@ -101,6 +101,20 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%stiffXs(I) = LineProp%stiffXs(I) Line%stiffYs(I) = LineProp%stiffYs(I) ! note: this does not convert to E (not EA) like done in C version END DO + + ! find static strain of the slow spring on the original working curve if Syrope model is used + if (Line%ElasticMod == 4) then + DO I = 1,Line%nEApoints + Line%stiffZs(I) = Line%stiffXs(I) - 1.0 / Line%vbeta * log(1.0 + Line%vbeta / Line%alphaMBL * Line%stiffYs(I)) + END DO + end if + + ! set working curve formula + if (Line%ElasticMod == 4) then + Line%syropeWCForm = LineProp%syropeWCForm + line%syropeWCK1 = LineProp%syropeWCK1 + line%syropeWCK2 = LineProp%syropeWCK2 + endif Line%nBApoints = LineProp%nBApoints DO I = 1,Line%nBApoints @@ -241,6 +255,16 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) RETURN END IF + ! allocate Syrope state variables if Syrope model is being used + if (Line%ElasticMod == 4) then + ALLOCATE ( Line%Tmax(N), Line%Tmean(N) ) + IF ( ErrStat /= ErrID_None ) THEN + ErrMsg = ' Error allocating Syrope arrays, Tmax and Tmean.' + !CALL CleanUp() + RETURN + END IF + end if + ! allocate node force vectors ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), & Line%Ap(3, 0:N), Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), & @@ -301,6 +325,9 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) REAL(DbKi), ALLOCATABLE :: LNodesZ(:) INTEGER(IntKi) :: N + REAL(DbKi) :: WCK1 ! Syrope working curve parameter 1 + REAL(DbKi) :: WCK2 ! Syrope working curve parameter 2 + N = Line%N ! for convenience @@ -412,10 +439,10 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) ENDIF - ! If using the viscoelastic model, initalize the deltaL_1 to the delta L of the segment. + ! If using the viscoelastic model (not Syrope), initalize the deltaL_1 to the delta L of the segment. ! This is required here to initalize the state as non-zero, which avoids an initial ! transient where the segment goes from unstretched to stretched in one time step. - IF (Line%ElasticMod > 1) THEN + IF (Line%ElasticMod > 1 .and. Line%ElasticMod <= 4) THEN DO I = 1, N ! calculate current (Stretched) segment lengths and unit tangent vectors (qs) for each segment (this is used for bending calculations) CALL UnitVector(Line%r(:,I-1), Line%r(:,I), Line%qs(:,I), Line%lstr(I)) @@ -423,6 +450,26 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) END DO ENDIF + ! If using Syrope model, initialize dl_1 to the static stretch based on the working curve and mean tension + IF (Line%ElasticMod == 4) THEN + DO I = 1, N + CALL UnitVector(Line%r(:,I-1), Line%r(:,I), Line%qs(:,I), Line%lstr(I)) + CALL SetupWorkingCurve(Line, Line%Tmax(I), ErrStat2, ErrMsg2) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + CALL calculate1DinterpolationXY( Line%stiffys_(1:Line%nEApoints), & + Line%stiffzs_(1:Line%nEApoints), & + Line%Tmean(I), & + Line%dl_1(I), ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + END DO + END IF + ! initialize phi to unique values on the range 0-2pi if VIV model is active IF (Line%Cl > 0) THEN DO I=0, Line%N @@ -1172,6 +1219,10 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(120) :: RoutineName = 'Line_GetStateDeriv' + Real(DbKi) :: WCK1 ! coefficient for the Syrope working curve formula + Real(DbKi) :: WCK2 ! coefficient for the Syrope working curve formula + Real(DbKi) :: dl_s ! static stretch (both slow and fast components) + ErrStat = ErrID_None ErrMsg = "" @@ -1366,7 +1417,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai MagTd = Line%BA* ( Line%lstrd(I) - Line%lstr(I)*Line%ld(I)/Line%l(I) )/Line%l(I) ! viscoelastic model from https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018 - else if (Line%ElasticMod > 1) then + else if (Line%ElasticMod > 1 .and. Line%ElasticMod < 4) then if (Line%ElasticMod == 3) then if (Line%dl_1(I) > 0.0) then @@ -1420,6 +1471,76 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! update state derivative for static stiffness stretch (last N entries in the line state vector if no VIV model, otherwise N entries past the 6*N-6 entries in this vector) Xd( 6*N-6 + I) = ld_1 + else if (Line%ElasticMod == 4) then + + dl = Line%lstr(I) - Line%l(I) ! delta l of this segment + + if (EqualRealNos(Line%Tmax(I), Line%Tmean(I))) then + ! on the original working curve + CALL calculate1DinterpolationXY( Line%stiffZs(1:Line%nEApoints), & + Line%stiffYs(1:Line%nEApoints), & + Line%dl_1(I) / Line%l(I), & + Line%Tmean(I), ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + + if (dl >= 0.0) then + MagT = Line%Tmean(I) + else + MagT = 0.0_DbKi ! cable can't "push" + endif + + CALL calculate1DinterpolationXY( Line%stiffYs(1:Line%nEApoints), & + Line%stiffXs(1:Line%nEApoints), & + Line%Tmean(I), & + dl_s, ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + else + ! on the active working curve + CALL calculate1DinterpolationXY( Line%stiffzs_(1:Line%nEApoints), & + Line%stiffys_(1:Line%nEApoints), & + Line%dl_1(I) / Line%l(I), & + Line%Tmean(I), ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + + if (dl >= 0.0) then + MagT = Line%Tmean(I) + else + MagT = 0.0_DbKi ! cable can't "push" + endif + + CALL calculate1DinterpolationXY( Line%stiffys_(1:Line%nEApoints), & + Line%stiffxs_(1:Line%nEApoints), & + Line%Tmean(I), & + dl_s, ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + endif + + ld_1 = ((dl - dl_s * Line%l(I)) * (Line%alphaMBL + Line%vbeta * Line%Tmean(I)) + Line%BA_D * Line%lstrd(I)) / (Line%BA_D + Line%BA) + MagTd = Line%BA * ld_1 / Line%l(I) + Xd( 6*N-6 + I) = ld_1 + + ! update Tmax and working curve if Tmean > Tmax + if (Line%Tmean(I) > Line%Tmax(I)) then + Line%Tmax(I) = Line%Tmean(I) + CALL SetupWorkingCurve(Line, Line%Tmax(I), ErrStat2, ErrMsg2) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + RETURN + END IF + end if + end if @@ -2053,6 +2174,97 @@ subroutine VisLinesMesh_Update(p,m,y,ErrStat,ErrMsg) enddo end subroutine VisLinesMesh_Update + SUBROUTINE SetupWorkingCurve(Line, Tmax, ErrStat, ErrMsg) + + TYPE(MD_Line), INTENT(INOUT) :: Line + REAL(DbKi), INTENT(IN ) :: Tmax ! Maximum tension used to build the working curve + INTEGER(IntKi), INTENT( OUT) :: ErrStat + CHARACTER(*), INTENT( OUT) :: ErrMsg + CHARACTER(*), PARAMETER :: RoutineName = 'SetupWorkingCurve' + + INTEGER(IntKi) :: I + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + REAL(DbKi) :: eps_max, eps_min ! maximum and minimum strains for the working curve + REAL(DbKi) :: xi ! normalized strain on the working curve + REAL(DbKi) :: denom + REAL(DbKi) :: log_arg + REAL(DbKi) :: exp_denom + + ErrStat = ErrID_None + ErrMsg = '' + ErrStat2 = ErrID_None + ErrMsg2 = '' + + ! input validation + IF (Tmax < Line%StiffYs(1) .OR. Tmax > Line%StiffYs(Line%nEApoints)) THEN + CALL SetErrStat(ErrID_Fatal, "Tmax is outside the working curve range.", ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + ! eps_max is computed from Tmax and the original working curve + CALL calculate1DinterpolationXY( Line%StiffYs(1:Line%nEApoints), & + Line%StiffXs(1:Line%nEApoints), & + Tmax, eps_max, ErrStat2, ErrMsg2 ) + IF (ErrStat2 /= ErrID_None) THEN + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + ! eps_min is the strain at zero tension + eps_min = Line%StiffXs(1) + Line%syropeWCK1 * (eps_max - Line%StiffXs(1)) + + ! for linear working curve formula, if Line%syropeWCK1 > 1, treat it as slope/stiffness of the working curve + ! in general, this value is much larger than 1 + IF (Line%syropeWCForm == 1 .AND. Line%syropeWCK1 >= 1.0_DbKi) THEN + eps_min = eps_max - Tmax / Line%syropeWCK1 + END IF + + ! make sure eps_min is less than eps_max and larger than zero + IF (eps_min < 0.0_DbKi .OR. eps_min >= eps_max) THEN + CALL SetErrStat(ErrID_Fatal, "Invalid working curve parameters: the zero tension strain on the working curve is out of range.", ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + denom = eps_max - eps_min + IF (ABS(denom) <= TINY(1.0_DbKi)) THEN + CALL SetErrStat(ErrID_Fatal, & + "Vertical Syrope working curve detected.", & + ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + DO I = 1, Line%nEApoints + Line%Stiffxs_(I) = eps_min + denom * REAL(I - 1, DbKi) / REAL(Line%nEApoints - 1, DbKi) + xi = (Line%Stiffxs_(I) - eps_min) / denom ! normalized strain + + IF (Line%syropeWCForm == 1) THEN ! linear working curve + Line%Stiffys_(I) = Tmax * xi + + ELSE IF (Line%syropeWCForm == 2) THEN ! quadratic working curve + Line%Stiffys_(I) = Tmax * xi * (Line%syropeWCK2 * xi + (1.0_DbKi - Line%syropeWCK2)) + + ELSE IF (Line%syropeWCForm == 3) THEN ! exponential working curve + exp_denom = 1.0_DbKi - EXP(Line%syropeWCK2) + Line%Stiffys_(I) = Tmax * (1.0_DbKi - EXP(Line%syropeWCK2 * xi)) / exp_denom + + ELSE + CALL SetErrStat(ErrID_Fatal, "Invalid working curve formula specified.", ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + log_arg = 1.0_DbKi + (Line%vbeta / Line%alphaMBL) * Line%Stiffys_(I) + IF (log_arg <= 0.0_DbKi) THEN + CALL SetErrStat(ErrID_Fatal, & + "Non-positive argument in log() when setting working curve in Syrope model.", & + ErrStat, ErrMsg, RoutineName) + RETURN + END IF + + Line%Stiffzs_(I) = Line%Stiffxs_(I) - (1.0_DbKi / Line%vbeta) * LOG(log_arg) + END DO + + END SUBROUTINE SetupWorkingCurve END MODULE MoorDyn_Line From 38b3dc69b9cd36e373736c43bc12148ab976b667 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 20 Apr 2026 13:34:00 +0000 Subject: [PATCH 05/14] Allocate memory with stat --- modules/moordyn/src/MoorDyn_Line.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 3e4c567ed6..813a9d6cc6 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -257,7 +257,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! allocate Syrope state variables if Syrope model is being used if (Line%ElasticMod == 4) then - ALLOCATE ( Line%Tmax(N), Line%Tmean(N) ) + ALLOCATE ( Line%Tmax(N), Line%Tmean(N), STAT = ErrStat ) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating Syrope arrays, Tmax and Tmean.' !CALL CleanUp() @@ -442,7 +442,7 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) ! If using the viscoelastic model (not Syrope), initalize the deltaL_1 to the delta L of the segment. ! This is required here to initalize the state as non-zero, which avoids an initial ! transient where the segment goes from unstretched to stretched in one time step. - IF (Line%ElasticMod > 1 .and. Line%ElasticMod <= 4) THEN + IF (Line%ElasticMod > 1 .and. Line%ElasticMod < 4) THEN DO I = 1, N ! calculate current (Stretched) segment lengths and unit tangent vectors (qs) for each segment (this is used for bending calculations) CALL UnitVector(Line%r(:,I-1), Line%r(:,I), Line%qs(:,I), Line%lstr(I)) From 25cca6dc533bde9a7cd838f6645c7387151cbd1b Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 6 May 2026 07:23:47 +0000 Subject: [PATCH 06/14] Updated the error message in the SYROPE IC section --- modules/moordyn/src/MoorDyn.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 7d5a468553..cc4379fc40 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -1798,7 +1798,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! check for correct number of columns in current line IF ( CountWords( Line ) /= N+2 ) THEN - CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be 5 columns.', ErrStat, ErrMsg, RoutineName ) + CALL SetErrStat( ErrID_Fatal, ' Unable to parse Syrope line initial conditions '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be '//trim(Num2LStr(N+2))//' columns (, Tmax, Tmean).', ErrStat, ErrMsg, RoutineName ) CALL CleanUp() RETURN END IF From 3c146e49776b85d1e2bbfef47e422f575f95eb4f Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 6 May 2026 07:39:27 +0000 Subject: [PATCH 07/14] Line ID must be >=1 and used a robust way to read the ICs --- modules/moordyn/src/MoorDyn.f90 | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index cc4379fc40..a096cc2ca7 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -1804,12 +1804,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END IF ! parse out entries: SyropeLineID Tmax Tmean - IF (ErrStat2 == 0) THEN - READ(Line,*,IOSTAT=ErrStat2) TempIDnums(1:N), Tmax0, Tmean0 + ErrStat2 = 0 + READ(Line,*,IOSTAT=ErrStat2) TempIDnums(1:N), Tmax0, Tmean0 + IF (ErrStat2 /= 0) THEN + CALL SetErrStat( ErrID_Fatal, ' Unable to parse Syrope line initial conditions '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN END IF DO il = 1, N - if (TempIDnums(il) <= p%nLines) then ! ensure line ID is in range + if (TempIDnums(il) >= 1 .and. TempIDnums(il) <= p%nLines) then ! ensure line ID is in range DO j = 1, m%LineList(TempIDnums(il))%N m%LineList(TempIDnums(il))%Tmax(j) = Tmax0 m%LineList(TempIDnums(il))%Tmean(j) = Tmean0 From 94e001cd133b46ad12310a4debdd4fbc7d35d843 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 6 May 2026 08:16:28 +0000 Subject: [PATCH 08/14] Removed unused varaibles WCK1 and WCK2 --- modules/moordyn/src/MoorDyn_Line.f90 | 6 ------ 1 file changed, 6 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 813a9d6cc6..750a1a45c8 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -325,10 +325,6 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) REAL(DbKi), ALLOCATABLE :: LNodesZ(:) INTEGER(IntKi) :: N - REAL(DbKi) :: WCK1 ! Syrope working curve parameter 1 - REAL(DbKi) :: WCK2 ! Syrope working curve parameter 2 - - N = Line%N ! for convenience ! try to calculate initial line profile using catenary routine (from FAST v.7) @@ -1219,8 +1215,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(120) :: RoutineName = 'Line_GetStateDeriv' - Real(DbKi) :: WCK1 ! coefficient for the Syrope working curve formula - Real(DbKi) :: WCK2 ! coefficient for the Syrope working curve formula Real(DbKi) :: dl_s ! static stretch (both slow and fast components) ErrStat = ErrID_None From 4dda4d2403b0484710766d9dd6b563464336baef Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 6 May 2026 08:33:04 +0000 Subject: [PATCH 09/14] Handle Syrope working-curve warnings without aborting MD_Init --- modules/moordyn/src/MoorDyn.f90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index a096cc2ca7..3d6892479c 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -747,12 +747,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er tempString6 = adjustl(tempString5(8:)) ! get the filepath of the SYROPE curves CALL MD_ReadSyropeWorkingCurves(tempString6, owcPath, wcFormula, wcK1, wcK2, ErrStat2, ErrMsg2) - - if (ErrStat2 /= ErrID_None ) then - CALL SetErrStat( ErrID_Fatal, 'Failed to read SYROPE working curves for line type '//trim(Num2LStr(l))//'. Check error message for details.', ErrStat, ErrMsg, RoutineName ) - CALL CleanUp() - RETURN - end if + IF (ErrStat2 >= AbortErrLev) THEN + ErrMsg2 = 'Failed to read SYROPE working curves for line type '//trim(Num2LStr(l))//'.'//NewLine//TRIM(ErrMsg2) + END IF + CALL CheckError( ErrStat2, ErrMsg2 ) + IF (ErrStat >= AbortErrLev) RETURN + ErrStat2 = ErrID_None + ErrMsg2 = '' if (wordy > 0) print *, " Read Syrope working curve for Line type ", trim(Num2LStr(l)) From d78fa63cf373566064f5f0b03586e284c879bffa Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Wed, 6 May 2026 15:32:29 +0000 Subject: [PATCH 10/14] New tests passed --- reg_tests/CTestList.cmake | 2 ++ reg_tests/r-test | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index efe0587cab..8676ffb955 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -571,6 +571,8 @@ md_regression("md_waterkin3" "moordyn") py_md_regression("py_md_5MW_OC4Semi" "moordyn;python") # the following tests are excessively slow in double precision, so skip these in normal testing #md_regression("md_Single_Line_Quasi_Static_Test" "moordyn") +md_regression("md_viscoelastic" "moordyn") +md_regression("md_syrope" "moordyn") # OpenFAST IO Library regression tests py_openfast_io_library_pytest("openfast_io_library" "openfast_io;python") diff --git a/reg_tests/r-test b/reg_tests/r-test index ff9b577119..818aa5a855 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit ff9b577119f5d056168d711b291d0156ee033955 +Subproject commit 818aa5a855a9adfbf6d9913775453b48dd4a6b2c From aca2b2f1841594b86c429bbd4a1219753897957b Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Thu, 7 May 2026 11:16:23 +0000 Subject: [PATCH 11/14] dl_1 denotes stretch instead of strain --- modules/moordyn/src/MoorDyn_Line.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 750a1a45c8..a491e67596 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -459,6 +459,7 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) Line%stiffzs_(1:Line%nEApoints), & Line%Tmean(I), & Line%dl_1(I), ErrStat2, ErrMsg2 ) + Line%dl_1(I) = Line%dl_1(I) * Line%l(I) IF (ErrStat2 /= ErrID_None) THEN CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') RETURN From 494104cd488f813ab000d6a908d2a16d1d9aa3eb Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 11 May 2026 07:31:45 +0000 Subject: [PATCH 12/14] Initialize Tmax and Tmean to zero --- modules/moordyn/src/MoorDyn_Line.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index a491e67596..055d0d9190 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -263,6 +263,8 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) !CALL CleanUp() RETURN END IF + Line%Tmax = 0.0_DbKi + Line%Tmean = 0.0_DbKi end if ! allocate node force vectors From 6974791cb5f136be4507413026b59d6818b9c93e Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 11 May 2026 08:01:20 +0000 Subject: [PATCH 13/14] Check Syrope input validity --- modules/moordyn/src/MoorDyn_Line.f90 | 36 +++++++++++++++++++++++++--- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 055d0d9190..8daa5c34d2 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -463,7 +463,12 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) Line%dl_1(I), ErrStat2, ErrMsg2 ) Line%dl_1(I) = Line%dl_1(I) * Line%l(I) IF (ErrStat2 /= ErrID_None) THEN - CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + CALL SetErrStat(ErrStat2, & + 'The Syrope slow spring strain is not a single-valued function of mean tension for Line '// & + trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// & + 'The static stiffness on the working curve is smaller than the dynamic stiffness. '// & + trim(ErrMsg2), & + ErrStat, ErrMsg, '') RETURN END IF END DO @@ -1469,6 +1474,21 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Xd( 6*N-6 + I) = ld_1 else if (Line%ElasticMod == 4) then + + if (Line%dl_1(I) < 0.0) then + CALL SetErrStat(ErrID_Fatal,"Syrope model: Slow spring stretch cannot be less than zero",ErrStat,ErrMsg,RoutineName) + return + end if + + IF (Line%alphaMBL <= 0 .OR. Line%vbeta <= 0 .OR. Line%l(I) <= 0 .OR. Line%dl_1(I) <= 0) THEN + CALL SetErrStat(ErrID_Warn,"Syrope model: Assumptions violated",ErrStat,ErrMsg,RoutineName) + if (wordy > 2) then + print *, "Line%alphaMBL", Line%alphaMBL + print *, "Line%vbeta", Line%vbeta + print *, "Line%l(I)", Line%l(I) + print *, "Line%dl_1(I)", Line%dl_1(I) + endif + ENDIF dl = Line%lstr(I) - Line%l(I) ! delta l of this segment @@ -1479,7 +1499,12 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%dl_1(I) / Line%l(I), & Line%Tmean(I), ErrStat2, ErrMsg2 ) IF (ErrStat2 /= ErrID_None) THEN - CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + CALL SetErrStat(ErrStat2, & + 'The Syrope mean tension is not a single-valued function of slow spring strain for Line '// & + trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// & + 'The static stiffness on the original working curve is smaller than the dynamic stiffness. '// & + trim(ErrMsg2), & + ErrStat, ErrMsg, '') RETURN END IF @@ -1504,7 +1529,12 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%dl_1(I) / Line%l(I), & Line%Tmean(I), ErrStat2, ErrMsg2 ) IF (ErrStat2 /= ErrID_None) THEN - CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, '') + CALL SetErrStat(ErrStat2, & + 'The Syrope slow spring strain is not a single-valued function of mean tension for Line '// & + trim(Num2LStr(Line%IdNum))//', segment '//trim(Num2LStr(I))//'. '// & + 'The static stiffness on the working curve is smaller than the dynamic stiffness. '// & + trim(ErrMsg2), & + ErrStat, ErrMsg, '') RETURN END IF From 9951a306a94db60cde99daf600dedf4faa31f879 Mon Sep 17 00:00:00 2001 From: Zhilong Wei Date: Mon, 11 May 2026 14:11:05 +0000 Subject: [PATCH 14/14] Support relative path of the Syrope input files --- modules/moordyn/src/MoorDyn.f90 | 54 +++++++++++++++++++----------- modules/moordyn/src/MoorDyn_IO.f90 | 18 +++++++--- 2 files changed, 47 insertions(+), 25 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 3d6892479c..fd29d2d55f 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -122,13 +122,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CHARACTER(40) :: DepthValue ! Temporarily stores the optional WtrDpth setting for MD, which could be a number or a filename CHARACTER(40) :: WaterKinValue ! Temporarily stores the optional WaterKin setting for MD, which is typically a filename INTEGER(IntKi) :: nOpts ! number of options lines in input file - CHARACTER(40) :: TempString1 ! - CHARACTER(40) :: TempString2 ! - CHARACTER(40) :: TempString3 ! - CHARACTER(40) :: TempString4 ! - CHARACTER(40) :: TempString5 ! - CHARACTER(40) :: TempString6 ! - CHARACTER(40) :: TempStrings(6) ! Array of 6 strings used when parsing comma-separated items + CHARACTER(1024) :: TempString1 ! + CHARACTER(1024) :: TempString2 ! + CHARACTER(1024) :: TempString3 ! + CHARACTER(1024) :: TempString4 ! + CHARACTER(1024) :: TempString5 ! + CHARACTER(1024) :: TempString6 ! + CHARACTER(1024) :: TempStrings(6) ! Array of 6 strings used when parsing comma-separated items ! CHARACTER(1024) :: FileName ! LOGICAL :: hasSyrope = .FALSE. ! flag for if a line type has a Syrope model specified LOGICAL :: hasSyropeICs = .FALSE. ! flag for whether any line type uses Syrope model @@ -159,8 +159,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CHARACTER(*), PARAMETER :: RoutineName = 'MD_Init' ! for Syrope working curves - CHARACTER(40) :: owcPath ! - CHARACTER(40) :: wcFormula ! + CHARACTER(1024) :: owcPath ! + CHARACTER(1024) :: wcFormula ! REAL(DbKi) :: wcK1 ! REAL(DbKi) :: wcK2 ! @@ -745,6 +745,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er read(tempStrings(3), *) m%LineTypeList(l)%vbeta tempString6 = adjustl(tempString5(8:)) ! get the filepath of the SYROPE curves + IF ( PathIsRelative( tempString6 ) ) THEN + tempString6 = TRIM(p%PriPath)//TRIM(tempString6) + END IF CALL MD_ReadSyropeWorkingCurves(tempString6, owcPath, wcFormula, wcK1, wcK2, ErrStat2, ErrMsg2) IF (ErrStat2 >= AbortErrLev) THEN @@ -771,7 +774,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL getCoefficientOrCurve(owcPath, m%LineTypeList(l)%EA, & m%LineTypeList(l)%nEApoints, & m%LineTypeList(l)%stiffXs, & - m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2) + m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2, p%PriPath) hasSyrope = .false. ! reset hasSyrope flag for later use when processing @@ -795,7 +798,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL getCoefficientOrCurve(tempStrings(1), m%LineTypeList(l)%EA, & m%LineTypeList(l)%nEApoints, & m%LineTypeList(l)%stiffXs, & - m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2) + m%LineTypeList(l)%stiffYs, ErrStat2, ErrMsg2, p%PriPath) end if @@ -817,13 +820,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL getCoefficientOrCurve(tempStrings(1), m%LineTypeList(l)%BA, & m%LineTypeList(l)%nBApoints, & m%LineTypeList(l)%dampXs, & - m%LineTypeList(l)%dampYs, ErrStat2, ErrMsg2) + m%LineTypeList(l)%dampYs, ErrStat2, ErrMsg2, p%PriPath) ! process bending stiffness coefficients (which might use lookup tables) CALL getCoefficientOrCurve(tempString3, m%LineTypeList(l)%EI, & m%LineTypeList(l)%nEIpoints, & m%LineTypeList(l)%bstiffXs, & - m%LineTypeList(l)%bstiffYs, ErrStat2, ErrMsg2) + m%LineTypeList(l)%bstiffYs, ErrStat2, ErrMsg2, p%PriPath) ! specify IdNum of line type for error checking m%LineTypeList(l)%IdNum = l @@ -5068,10 +5071,11 @@ SUBROUTINE MD_ReadSyropeWorkingCurves(inputString, owcPath, wcFormula, k1, k2, E CHARACTER(1024) :: NextLine CHARACTER(64) :: OptString CHARACTER(256) :: OptValue - LOGICAL :: FoundOWC = .false. - LOGICAL :: FoundWCFormula = .false. - LOGICAL :: FoundK1 = .false. - LOGICAL :: FoundK2 = .false. + CHARACTER(256) :: Words(2) + LOGICAL :: FoundOWC + LOGICAL :: FoundWCFormula + LOGICAL :: FoundK1 + LOGICAL :: FoundK2 CHARACTER(*), PARAMETER :: RoutineName = 'MD_ReadSyropeWorkingCurves' INTEGER(IntKi) :: i, ios @@ -5085,6 +5089,10 @@ SUBROUTINE MD_ReadSyropeWorkingCurves(inputString, owcPath, wcFormula, k1, k2, E k2 = 0.0_DbKi ErrStat3 = ErrID_None ErrMsg3 = '' + FoundOWC = .false. + FoundWCFormula = .false. + FoundK1 = .false. + FoundK2 = .false. ! read file @@ -5104,13 +5112,19 @@ SUBROUTINE MD_ReadSyropeWorkingCurves(inputString, owcPath, wcFormula, k1, k2, E NextLine = TRIM(FileInfo%Lines(i)) IF (LEN_TRIM(NextLine) == 0) CYCLE - READ(NextLine, *, IOSTAT=ios) OptValue, OptString - IF (ios /= 0) THEN + OptValue = '' + OptString = '' + Words = '' + CALL GetWords(NextLine, Words, 2) + IF (LEN_TRIM(Words(1)) == 0 .OR. LEN_TRIM(Words(2)) == 0) THEN CALL SetErrStat(ErrID_Fatal, & - 'Failed to read options in Syrope working curve file ' // TRIM(inputString) // '.', & + 'Failed to read options in Syrope working curve file ' // TRIM(inputString) // & + '. Expected lines in the form " ". Line was: ' // TRIM(NextLine), & ErrStat3, ErrMsg3, RoutineName) RETURN END IF + OptValue = TRIM(Words(1)) + OptString = TRIM(Words(2)) CALL Conv2UC(OptString) diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index 3061b350f9..0f5fff9c20 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -241,9 +241,10 @@ END SUBROUTINE setupBathymetry ! read in stiffness/damping coefficient or load nonlinear data file if applicable - SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, LineProp_Xs, LineProp_Ys, ErrStat3, ErrMsg3) + SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, LineProp_Xs, LineProp_Ys, & + ErrStat3, ErrMsg3, PriPath) - CHARACTER(40), INTENT(IN ) :: inputString + CHARACTER(*), INTENT(IN ) :: inputString REAL(DbKi), INTENT(INOUT) :: LineProp_c INTEGER(IntKi), INTENT( OUT) :: LineProp_nPoints REAL(DbKi), INTENT( OUT) :: LineProp_Xs (MD_MaxNCoef) ! MD_MaxNCoef set in registry @@ -251,6 +252,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line INTEGER(IntKi), INTENT( OUT) :: ErrStat3 ! Error status of the operation CHARACTER(*), INTENT( OUT) :: ErrMsg3 ! Error message if ErrStat /= ErrID_None + CHARACTER(*), OPTIONAL, INTENT(IN) :: PriPath ! Path to the primary MoorDyn input file INTEGER(IntKi) :: I INTEGER(IntKi) :: UnCoef ! unit number for coefficient input file @@ -259,6 +261,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line INTEGER(IntKi) :: ErrStat4 CHARACTER(120) :: ErrMsg4 CHARACTER(120) :: Line2 + CHARACTER(1024) :: FileName if (SCAN(inputString, "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ") == 0) then ! "eE" are exluded as they're used for scientific notation! @@ -273,9 +276,14 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line LineProp_c = 0.0 ! load lookup table data from file + IF ( PRESENT(PriPath) .AND. PathIsRelative( inputString ) ) THEN + FileName = TRIM(PriPath)//TRIM(inputString) + ELSE + FileName = TRIM(inputString) + END IF CALL GetNewUnit( UnCoef ) - CALL OpenFInpFile( UnCoef, TRIM(inputString), ErrStat4, ErrMsg4 ) ! add error handling? + CALL OpenFInpFile( UnCoef, TRIM(FileName), ErrStat4, ErrMsg4 ) ! add error handling? IF (ErrStat4 == ErrID_Fatal) then ErrStat3 = ErrStat4 ErrMsg3 = ErrMsg4 @@ -303,7 +311,7 @@ SUBROUTINE getCoefficientOrCurve(inputString, LineProp_c, LineProp_npoints, Line if (I < 2) then ErrStat3 = ErrID_Fatal - ErrMsg3 = "Less than the minimum of 2 data lines found in file "//TRIM(inputString)//" (first 3 lines are headers)." + ErrMsg3 = "Less than the minimum of 2 data lines found in file "//TRIM(FileName)//" (first 3 lines are headers)." LineProp_npoints = 0 Close (UnCoef) RETURN @@ -322,7 +330,7 @@ SUBROUTINE SplitByBars(instring, n, outstrings) CHARACTER(*), INTENT(INOUT) :: instring INTEGER(IntKi), INTENT( OUT) :: n - CHARACTER(40), INTENT(INOUT) :: outstrings(6) ! array of output strings. Up to 6 strings can be read + CHARACTER(*), INTENT(INOUT) :: outstrings(6) ! array of output strings. Up to 6 strings can be read INTEGER :: pos1, pos2