diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 67191d4609..3d6892479c 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,77 @@ 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 >= 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)) + + 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 +1768,64 @@ 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 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 + + ! parse out entries: SyropeLineID Tmax Tmean + 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) >= 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 + 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 +2806,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 +5052,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 diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 50e93350bf..a491e67596 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), STAT = ErrStat ) + 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,7 +325,6 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) REAL(DbKi), ALLOCATABLE :: LNodesZ(:) INTEGER(IntKi) :: N - N = Line%N ! for convenience ! try to calculate initial line profile using catenary routine (from FAST v.7) @@ -412,10 +435,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 +446,27 @@ 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 ) + Line%dl_1(I) = Line%dl_1(I) * Line%l(I) + 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 +1216,8 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(120) :: RoutineName = 'Line_GetStateDeriv' + Real(DbKi) :: dl_s ! static stretch (both slow and fast components) + ErrStat = ErrID_None ErrMsg = "" @@ -1366,7 +1412,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 +1466,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 +2169,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 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 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 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