diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 43d5c2a3a5..9a580ec981 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -76,6 +76,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! CHARACTER(1024) :: priPath ! The path to the primary MoorDyn input file REAL(DbKi) :: t ! instantaneous time, to be used during IC generation INTEGER(IntKi) :: l ! index + INTEGER(IntKi) :: il ! index + INTEGER(IntKi) :: iil ! index + INTEGER(IntKi) :: Success ! flag for checking whether line is attached to failure point INTEGER(IntKi) :: I ! Current line number of input file INTEGER(IntKi) :: J ! index INTEGER(IntKi) :: K ! index @@ -556,7 +559,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ALLOCATE(m%BodyList( p%nBodies ), STAT = ErrStat2 ); if(AllocateFailed("BodyList" )) return ALLOCATE(m%RodList( p%nRods ), STAT = ErrStat2 ); if(AllocateFailed("RodList" )) return - ALLOCATE(m%PointList( p%nPoints ), STAT = ErrStat2 ); if(AllocateFailed("PointList" )) return + ALLOCATE(m%PointList( p%nPointsExtra), STAT = ErrStat2 ); if(AllocateFailed("PointList" )) return ALLOCATE(m%LineList( p%nLines ), STAT = ErrStat2 ); if(AllocateFailed("LineList" )) return ALLOCATE(m%FailList( p%nFails ), STAT = ErrStat2 ); if(AllocateFailed("FailList" )) return @@ -565,16 +568,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! Allocate associated index arrays (note: some are allocated larger than will be used, for simplicity) ALLOCATE(m%BodyStateIs1(p%nBodies ), m%BodyStateIsN(p%nBodies ), STAT=ErrStat2); if(AllocateFailed("BodyStateIs1/N")) return ALLOCATE(m%RodStateIs1(p%nRods ), m%RodStateIsN(p%nRods ), STAT=ErrStat2); if(AllocateFailed("RodStateIs1/N" )) return - ALLOCATE(m%PointStateIs1(p%nPoints), m%PointStateIsN(p%nPoints), STAT=ErrStat2); if(AllocateFailed("PointStateIs1/N" )) return + ALLOCATE(m%PointStateIs1(p%nPointsExtra), m%PointStateIsN(p%nPointsExtra), STAT=ErrStat2); if(AllocateFailed("PointStateIs1/N" )) return ALLOCATE(m%LineStateIs1(p%nLines) , m%LineStateIsN(p%nLines) , STAT=ErrStat2); if(AllocateFailed("LineStateIs1/N")) return ALLOCATE(m%FreeBodyIs( p%nBodies ), STAT=ErrStat2); if(AllocateFailed("FreeBodyIs")) return ALLOCATE(m%FreeRodIs( p%nRods ), STAT=ErrStat2); if(AllocateFailed("FreeRodIs")) return - ALLOCATE(m%FreePointIs( p%nPoints), STAT=ErrStat2); if(AllocateFailed("FreePointIs")) return + ALLOCATE(m%FreePointIs(p%nPointsExtra), STAT=ErrStat2); if(AllocateFailed("FreePointIs")) return ALLOCATE(m%CpldBodyIs(p%nBodies , p%nTurbines), STAT=ErrStat2); if(AllocateFailed("CpldBodyIs")) return ALLOCATE(m%CpldRodIs( p%nRods , p%nTurbines), STAT=ErrStat2); if(AllocateFailed("CpldRodIs")) return - ALLOCATE(m%CpldPointIs(p%nPoints, p%nTurbines), STAT=ErrStat2); if(AllocateFailed("CpldPointIs")) return + ALLOCATE(m%CpldPointIs(p%nPointsExtra, p%nTurbines), STAT=ErrStat2); if(AllocateFailed("CpldPointIs")) return ! ---------------------- now go through again and process file contents -------------------- @@ -811,10 +814,30 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er READ(tempString4, *) m%BodyList(l)%BodyCdA(1) m%BodyList(l)%BodyCdA(2) = m%BodyList(l)%BodyCdA(1) m%BodyList(l)%BodyCdA(3) = m%BodyList(l)%BodyCdA(1) + m%BodyList(l)%BodyCdA(4) = m%BodyList(l)%BodyCdA(1) + m%BodyList(l)%BodyCdA(5) = m%BodyList(l)%BodyCdA(1) + m%BodyList(l)%BodyCdA(6) = m%BodyList(l)%BodyCdA(1) + else if (N ==2) then + READ(tempStrings(1), *) m%BodyList(l)%BodyCdA(1) + READ(tempStrings(4), *) m%BodyList(l)%BodyCdA(4) + m%BodyList(l)%BodyCdA(2) = m%BodyList(l)%BodyCdA(1) + m%BodyList(l)%BodyCdA(3) = m%BodyList(l)%BodyCdA(1) + m%BodyList(l)%BodyCdA(5) = m%BodyList(l)%BodyCdA(4) + m%BodyList(l)%BodyCdA(6) = m%BodyList(l)%BodyCdA(4) else if (N==3) then ! all three coordinates provided READ(tempStrings(1), *) m%BodyList(l)%BodyCdA(1) READ(tempStrings(2), *) m%BodyList(l)%BodyCdA(2) READ(tempStrings(3), *) m%BodyList(l)%BodyCdA(3) + READ(tempStrings(1), *) m%BodyList(l)%BodyCdA(4) + READ(tempStrings(2), *) m%BodyList(l)%BodyCdA(5) + READ(tempStrings(3), *) m%BodyList(l)%BodyCdA(6) + else if (N==6) then + READ(tempStrings(1), *) m%BodyList(l)%BodyCdA(1) + READ(tempStrings(2), *) m%BodyList(l)%BodyCdA(2) + READ(tempStrings(3), *) m%BodyList(l)%BodyCdA(3) + READ(tempStrings(4), *) m%BodyList(l)%BodyCdA(4) + READ(tempStrings(5), *) m%BodyList(l)%BodyCdA(5) + READ(tempStrings(6), *) m%BodyList(l)%BodyCdA(6) else CALL SetErrStat( ErrID_Fatal, 'Body '//trim(Num2LStr(l))//' CdA entry (col 13) must have 1 or 3 numbers.' , ErrStat, ErrMsg, RoutineName ) end if @@ -877,7 +900,29 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%FreeBodyIs(p%nFreeBodies) = l - ! TODO: add option for body coupling to different turbines in FAST.Farm <<< + else if ((let1 == "TURBINE") .or. (let1 == "T")) then ! turbine-coupled in FAST.Farm case + + if (len_trim(num1) > 0) then + READ(num1, *) J ! convert to int, representing turbine index + + if ((J <= p%nTurbines) .and. (J > 0)) then + + m%BodyList(l)%typeNum = -1 ! set as coupled type + p%nCpldBodies(J)=p%nCpldBodies(J)+1 ! increment counter for the appropriate turbine + m%CpldBodyIs(p%nCpldBodies(J),J) = l + CALL WrScr(' added Body '//TRIM(int2lstr(l))//' for turbine '//trim(int2lstr(J))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Added Body '//TRIM(int2lstr(l))//' for turbine '//trim(int2lstr(J)) + end if + + else + CALL SetErrStat( ErrID_Fatal, "Turbine ID out of bounds for Body "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + return + end if + else + CALL SetErrStat( ErrID_Fatal, "No number provided for Body "//trim(Num2LStr(l))//" Turbine attachment.", ErrStat, ErrMsg, RoutineName ) + return + end if else if (let1 == "FREE") then ! if a free body m%BodyList(l)%typeNum = 0 @@ -1058,7 +1103,28 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%CpldRodIs(p%nCpldRods(1),1) = l m%FreeRodIs(p%nFreeRods) = l - ! TODO: add option for body coupling to different turbines in FAST.Farm <<< + else if ((let1 == "TURBINE") .or. (let1 == "T")) then ! turbine-coupled in FAST.Farm case + + if (len_trim(num1) > 0) then + READ(num1, *) J ! convert to int, representing turbine index + + if ((J <= p%nTurbines) .and. (J > 0)) then + m%RodList(l)%typeNum = -2 ! set as coupled type + p%nCpldRods(J)=p%nCpldRods(J)+1 ! increment counter for the appropriate turbine + m%CpldRodIs(p%nCpldRods(J),J) = l + CALL WrScr(' added Rod '//TRIM(int2lstr(l))//' for turbine '//trim(int2lstr(J))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Added Rod '//TRIM(int2lstr(l))//' for turbine '//trim(int2lstr(J)) + end if + + else + CALL SetErrStat( ErrID_Fatal, "Turbine ID out of bounds for Rod "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + return + end if + else + CALL SetErrStat( ErrID_Fatal, "No number provided for Rod "//trim(Num2LStr(l))//" Turbine attachment.", ErrStat, ErrMsg, RoutineName ) + return + end if else if ((let1 == "ROD") .or. (let1 == "R") .or. (let1 == "FREE")) then m%RodList(l)%typeNum = 0 @@ -1390,7 +1456,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er return end if else - CALL SetErrStat( ErrID_Fatal, "Error: rod point ID out of bounds for line "//trim(Num2LStr(l))//" end A attachment.", ErrStat, ErrMsg, RoutineName ) + CALL SetErrStat( ErrID_Fatal, " Rod ID out of bounds for line "//trim(Num2LStr(l))//" end A attachment.", ErrStat, ErrMsg, RoutineName ) return end if @@ -1431,7 +1497,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er return end if else - CALL SetErrStat( ErrID_Fatal, "Error: rod connection ID out of bounds for line "//trim(Num2LStr(l))//" end B attachment.", ErrStat, ErrMsg, RoutineName ) + CALL SetErrStat( ErrID_Fatal, " Rod ID out of bounds for line "//trim(Num2LStr(l))//" end B attachment.", ErrStat, ErrMsg, RoutineName ) return end if @@ -1526,6 +1592,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! count commas to determine how many line IDs specified for this channel N = count(transfer(Line, 'a', len(Line)) == ",") + 1 ! number of line IDs given + + ! check for correct number of columns in current line (CountWords splits by comma and space, so 2 columns means number of line ID's plus one more for the control channel) + IF ( CountWords( Line ) /= N+1 ) THEN + CALL SetErrStat( ErrID_Fatal, ' Unable to parse controls '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be 2 columns.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF ! parse out entries: CtrlChan, LineIdNums read(Line, *) Itemp, TempIDnums(1:N) ! parse out each line ID @@ -1539,16 +1612,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er write(p%UnLog,'(A)') 'Assigned Line '//TRIM(Int2LStr(TempIDnums(J)))//' to control channel '//TRIM(Int2LStr(Itemp)) end if else - CALL WrScr('Error: Line '//TRIM(Int2LStr(TempIDnums(J)))//' already is assigned to control channel '//TRIM(Int2LStr(m%LineList( TempIDnums(J) )%CtrlChan))//' so cannot also be assigned to channel '//TRIM(Int2LStr(Itemp))) - if (p%writeLog > 0) then - write(p%UnLog,'(A)') 'Error: Line '//TRIM(Int2LStr(TempIDnums(J)))//' already is assigned to control channel '//TRIM(Int2LStr(m%LineList( TempIDnums(J) )%CtrlChan))//' so cannot also be assigned to channel '//TRIM(Int2LStr(Itemp)) - end if + CALL SetErrStat( ErrID_Fatal, ' Line '//TRIM(Int2LStr(TempIDnums(J)))//' already is assigned to control channel '//TRIM(Int2LStr(m%LineList( TempIDnums(J) )%CtrlChan))//' so cannot also be assigned to channel '//TRIM(Int2LStr(Itemp)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return end if else - CALL WrScr('Error: Line ID '//TRIM(Int2LStr(TempIDnums(J)))//' of CtrlChan '//TRIM(Int2LStr(Itemp))//' is out of range') - if (p%writeLog > 0) then - write(p%UnLog,'(A)') 'Error: Line ID '//TRIM(Int2LStr(TempIDnums(J)))//' of CtrlChan '//TRIM(Int2LStr(Itemp))//' is out of range' - end if + CALL SetErrStat( ErrID_Fatal, ' Line ID '//TRIM(Int2LStr(TempIDnums(J)))//' of CtrlChan '//TRIM(Int2LStr(Itemp))//' is out of range', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return end if END DO @@ -1558,12 +1629,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er !------------------------------------------------------------------------------------------- else if (INDEX(Line, "FAILURE") > 0) then ! if failure conditions header - - CALL WrScr(" Warning: Failure capabilities are not yet implemented in MoorDyn.") - if (p%writeLog > 0) then - write(p%UnLog,'(A)') " Warning: Failure capabilities are not yet implemented in MoorDyn." - end if + IF (wordy > 0) then + CALL WrScr(" Reading failure inputs") + endif + + ! TODO: allocate fail list size (we need to do this before though right?) + + ! skip following two lines (label line and unit line) Line = NextLine(i) Line = NextLine(i) @@ -1573,13 +1646,169 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er !read into a line Line = NextLine(i) + + ! count commas to determine how many line IDs specified for this channel + N = count(transfer(Line, 'a', len(Line)) == ",") + 1 ! number of line IDs given + ! check for correct number of columns in current line (CountWords splits by comma and space, so 2 columns means number of line ID's plus 4 more for the other 4 channels) - ! TODO: Failure capabilities still need to be completed - READ(Line,*,IOSTAT=ErrStat2) m%LineList(l)%IdNum, tempString1, m%LineList(l)%UnstrLen, & - m%LineList(l)%N, tempString2, tempString3, LineOutString - - END DO - + IF ( CountWords( Line ) /= N+4 ) 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: FailID Point Lines FailTime FailTen + IF (ErrStat2 == 0) THEN + + READ(Line,*,IOSTAT=ErrStat2) m%FailList(l)%IdNum, TempString1, TempIDnums(1:N), m%FailList(l)%failTime, m%FailList(l)%failTen + + ! check for duplicate failure ID's + ! check for sequential IdNums + IF ( m%FailList(l)%IdNum .NE. l ) THEN + CALL SetErrStat( ErrID_Fatal, 'Failure ID numbers must be sequential starting from 1.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + + CALL Conv2UC(TempString1) ! convert to uppercase so that matching is not case-sensitive + + call DecomposeString(TempString1, let1, num1, let2, num2, let3) ! divided failPoint into letters and numbers + + if (len_trim(num1)<1) then + CALL SetErrStat( ErrID_Fatal, "Error: no point number provided for line failure "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + end if + + READ(num1, *) m%FailList(l)%attachID ! convert to int + + ! if id starts with an "R" or "Rod" + if ((let1 == "R") .OR. (let1 == "ROD")) then + if ((m%FailList(l)%attachID <= p%nRods) .AND. (m%FailList(l)%attachID > 0)) then + if (let2 == "A") then + m%FailList(l)%isRod = 1 + else if (let2 == "B") then + m%FailList(l)%isRod = 2 + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//'. Rod end must be A or B.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//'. Rod number out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + endif + + if ((len_trim(let1)<1) .OR. (let1 == "P") .OR. (let1 == "POINT")) then + if ((m%FailList(l)%attachID <= p%nPoints) .AND. (m%FailList(l)%attachID > 0)) then + m%FailList(l)%isRod = 0 + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//'. Point number out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + endif + + ! get lines + m%FailList(l)%nLinesToDetach = N + + DO il = 1, m%FailList(l)%nLinesToDetach + if (TempIDnums(il) <= p%nLines) then ! ensure line ID is in range + m%FailList(l)%lineIDs(il) = TempIDnums(il) + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse line failure '//trim(Num2LStr(l))//'. Line number '//TRIM(Int2LStr(TempIDnums(il)))//' out of bounds.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + ! check whether line is attached to fail point at fairlead or anchor and assing line tops + if (m%FailList(l)%isRod == 0) then ! point + + Success = 0 + DO iil = 1, m%PointList(m%FailList(l)%attachID)%nAttached ! find index of line + if (m%PointList(m%FailList(l)%attachID)%Attached(iil) == m%FailList(l)%lineIDs(il)) then + m%FailList(l)%lineTops(il) = m%PointList(m%FailList(l)%attachID)%Top(iil) + Success = 1 + exit + endif + ENDDO + + if (Success == 0) then + CALL SetErrStat( ErrID_Fatal, " Line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" not attached to point "//trim(num2lstr(m%FailList(l)%attachID))//" for failure "//trim(num2lstr(m%FailList(l)%IdNum)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + elseif (m%FailList(l)%isRod == 1) then ! Rod end A + + Success = 0 + DO iil = 1, m%RodList(m%FailList(l)%attachID)%nAttachedA ! find index of line + if (m%RodList(m%FailList(l)%attachID)%AttachedA(iil) == m%FailList(l)%lineIDs(il)) then + m%FailList(l)%lineTops(il) = m%RodList(m%FailList(l)%attachID)%TopA(iil) + Success = 1 + exit + endif + ENDDO + + if (Success == 0) then + CALL SetErrStat( ErrID_Fatal, " Line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" not attached to R"//trim(num2lstr(m%FailList(l)%attachID))//"A for failure "//trim(num2lstr(m%FailList(l)%IdNum)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + elseif (m%FailList(l)%isRod == 2) then ! Rod end B + + Success = 0 + DO iil = 1, m%RodList(m%FailList(l)%attachID)%nAttachedB ! find index of line + if (m%RodList(m%FailList(l)%attachID)%AttachedB(iil) == m%FailList(l)%lineIDs(il)) then + m%FailList(l)%lineTops(il) = m%RodList(m%FailList(l)%attachID)%TopB(iil) + Success = 1 + exit + endif + ENDDO + + if (Success == 0) then + CALL SetErrStat( ErrID_Fatal, " Line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" not attached to R"//trim(num2lstr(m%FailList(l)%attachID))//"B for failure "//trim(num2lstr(m%FailList(l)%IdNum)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + + else + CALL SetErrStat( ErrID_Fatal, " isRod out of range for failure "//trim(num2lstr(m%FailList(l)%IdNum)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + return + endif + ENDDO + + ! cant have both time and tension conditions, time is prioritized + if ((m%FailList(l)%failTime > 0) .AND. (m%FailList(l)%failTen > 0)) then + CALL SetErrStat( ErrID_Info, ' MoorDyn failure condition checks time before tension. If time reached before tension, failure '//trim(Num2LStr(m%FailList(l)%IdNum))//' will trigger.', ErrStat, ErrMsg, RoutineName ) + endif + + if ((m%PointList(m%FailList(l)%attachID)%typeNum == 0) .AND. (m%PointList(m%FailList(l)%attachID)%nAttached == m%FailList(l)%nLinesToDetach)) then + + ! if X lines called to fail from a free point with only X lines attached + Call SetErrStat(ErrID_Warn, trim(num2lstr(m%FailList(l)%nLinesToDetach))//" lines called to fail from a free point with only "//trim(num2lstr(m%FailList(l)%nLinesToDetach))//" lines attached. Failure "//trim(num2lstr(l))//" ignored.", ErrStat, ErrMsg, RoutineName ) + m%FailList(l)%failStatus = 2 + + elseif ((m%FailList(l)%failTime == 0) .AND. (m%FailList(l)%failTen == 0)) then + + CALL SetErrStat( ErrID_Warn, ' MoorDyn failure condition must have non-zero time or tension. Failure condition '//trim(Num2LStr(m%FailList(l)%IdNum))//' ignored.', ErrStat, ErrMsg, RoutineName ) + m%FailList(l)%failStatus = 2 + + else + + m%FailList(l)%failStatus = 0; ! initialize as unfailed + + endif + + endif + enddo + !------------------------------------------------------------------------------------------- else if (INDEX(Line, "OUTPUT") > 0) then ! if output header @@ -1788,13 +2017,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! prepare state vector etc. !------------------------------------------------------------------------------------ - ! the number of states is Nx + ! the number of states is Nx and Nxtra includes additional states for potential line failures m%Nx = Nx + m%Nxtra = m%Nx + 6*2*p%nLines - IF (wordy > 0) print *, "allocating state vectors to size ", Nx + IF (wordy > 0) print *, "allocating state vectors to size ", m%Nxtra ! allocate state vector and temporary state vectors based on size just calculated - ALLOCATE ( x%states(m%Nx), m%xTemp%states(m%Nx), m%xdTemp%states(m%Nx), STAT = ErrStat2 ) + ALLOCATE ( x%states(m%Nxtra), m%xTemp%states(m%Nxtra), m%xdTemp%states(m%Nxtra), STAT = ErrStat2 ) IF ( ErrStat2 /= ErrID_None ) THEN ErrMsg = ' Error allocating state vectors.' !CALL CleanUp() @@ -2498,7 +2728,7 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er INTEGER(IntKi) , INTENT(IN ) :: n TYPE(MD_InputType) , INTENT(INOUT) :: u(:) ! INTENT(INOUT) ! had to change this to INOUT REAL(DbKi) , INTENT(IN ) :: t_array(:) - TYPE(MD_ParameterType) , INTENT(IN ) :: p ! INTENT(IN ) + TYPE(MD_ParameterType) , INTENT(INOUT) :: p ! INTENT(IN ) TYPE(MD_ContinuousStateType) , INTENT(INOUT) :: x ! INTENT(INOUT) TYPE(MD_DiscreteStateType) , INTENT(INOUT) :: xd ! INTENT(INOUT) TYPE(MD_ConstraintStateType) , INTENT(INOUT) :: z ! INTENT(INOUT) @@ -2520,6 +2750,8 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er INTEGER(IntKi) :: NdtM ! number of time steps to integrate through with RK2 INTEGER(IntKi) :: I INTEGER(IntKi) :: J + INTEGER(IntKi) :: l, li, lii, il ! index + INTEGER(IntKi) :: tension ! tension at line attachment to failure point nTime = size(u) ! the number of times of input data provided? <<<<<<< not used @@ -2574,22 +2806,15 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er IF (Is_NaN(x%states(J))) THEN ErrStat = ErrID_Fatal ErrMsg = ' NaN state detected.' + IF (wordy > 1) THEN + print *, ". Here is the state vector: " + print *, x%states + END IF + CALL CheckError(ErrStat, ErrMsg) EXIT END IF END DO - IF (ErrStat == ErrID_Fatal) THEN - CALL WrScr("NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn.") - if (p%writeLog > 0) then - write(p%UnLog,'(A)') "NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn." - end if - IF (wordy > 1) THEN - print *, ". Here is the state vector: " - print *, x%states - END IF - EXIT - END IF - END DO ! I time steps @@ -2599,6 +2824,7 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er CALL MD_DestroyInput(u_interp, ErrStat, ErrMsg) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error destroying dxdt or x2.' + CALL CheckError(ErrStat, ErrMsg) END IF ! CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MD_UpdateStates') @@ -2608,23 +2834,181 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er IF (Is_NaN(x%states(J))) THEN ErrStat = ErrID_Fatal ErrMsg = ' NaN state detected.' + IF (wordy > 1) THEN + print *, ". Here is the state vector: " + print *, x%states + END IF + CALL CheckError(ErrStat, ErrMsg) EXIT END IF END DO - - IF (ErrStat == ErrID_Fatal) THEN - CALL WrScr("NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn.") - if (p%writeLog > 0) then - write(p%UnLog,'(A)') "NaN detected at time "//TRIM(Num2LStr(t2))//" in MoorDyn." - end if - IF (wordy > 1) THEN - print *, ". Here is the state vector: " - print *, x%states - END IF - END IF + + ! do we want to check failures here (at the coupling step level? Or at the dtM level?) + ! --------------- check for line failures (detachments!) ---------------- + DO l= 1,p%nFails + + if (m%FailList(l)%failStatus == 0) then + + if ((t >= m%FailList(l)%failTime) .AND. (m%FailList(l)%failTime .NE. 0.0)) then + + ! step 1: check for time-triggered failures + + CALL WrScr("Failure number "//trim(Num2LStr(l))//" triggered by t = "//trim(Num2LStr(t))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "Failure number "//trim(Num2LStr(l))//" triggered by t = "//trim(Num2LStr(t)) + end if + + m%FailList(l)%failStatus = 1; ! set status to failed so it's not checked again + CALL DetachLines(m%FailList(l)%attachID, m%FailList(l)%isRod, m%FailList(l)%lineIDs, m%FailList(l)%lineTops, m%FailList(l)%nLinesToDetach, t) + + elseif (m%FailList(l)%failTen .NE. 0.0) then + + ! step 2: check for tension-triggered failures (this will require specifying max tension things) + + DO il = 1,m%FailList(l)%nLinesToDetach ! for each line in the failure, check the tension at the attachment + + ! check line ID is right + if (m%FailList(l)%lineIDs(il) .NE. m%LineList(m%FailList(l)%lineIDs(il))%IdNum) then + CALL CheckError(ErrID_Fatal, " Line ID's dont match for failure "//trim(num2lstr(m%FailList(l)%IdNum))) + endif + + ! if fairlead else anchor + if (m%FailList(l)%lineTops(il) == 1) then + tension = Line_GetNodeTen(m%LineList(m%FailList(l)%lineIDs(il)), m%LineList(m%FailList(l)%lineIDs(il))%N) + else + tension = Line_GetNodeTen(m%LineList(m%FailList(l)%lineIDs(il)), 0) + endif + + if (tension >= m%FailList(l)%failTen) then + CALL WrScr("Failure number "//trim(Num2LStr(l))//" triggered by line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" tension = "//trim(Num2LStr(tension))//" at time = "//trim(Num2LStr(t))) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "Failure number "//trim(Num2LStr(l))//" triggered by line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" tension = "//trim(Num2LStr(tension))//" at time = "//trim(Num2LStr(t)) + end if + + m%FailList(l)%failStatus = 1; ! set status to failed so it's not checked again + CALL DetachLines(m%FailList(l)%attachID, m%FailList(l)%isRod, m%FailList(l)%lineIDs, m%FailList(l)%lineTops, m%FailList(l)%nLinesToDetach, t) + exit ! dont need to check the other lines becasue failure already detected + + endif + + ENDDO ! il = 1,m%FailList(l)%nLinesToDetach + + endif ! end checking time and tension thresholds non-zero + + if (m%FailList(l)%failStatus == 1) then + + ! if a failure is triggered, remove all lines from that failure from all other failures + DO li = 1, p%nFails + if (m%FailList(li)%IdNum .NE. m%FailList(l)%IdNum) then ! if not this failure + if ((m%FailList(l)%attachID == m%FailList(li)%attachID) .AND. (m%FailList(l)%isRod == m%FailList(li)%isRod)) then ! if failures are for the same point + + DO il = 1, m%FailList(l)%nLinesToDetach ! loop through lines of this failure + DO lii = 1, m%FailList(li)%nLinesToDetach ! loop through lines of the other failure + + if (m%FailList(l)%lineIDs(il) == m%FailList(li)%lineIDs(lii)) then ! if this failure's line IDs are found in any of the other failure's line IDs + + m%FailList(li)%nLinesToDetach = m%FailList(li)%nLinesToDetach - 1 ! reduce the size of nLinesToDetach of the other failure + m%FailList(li)%lineIDs(lii) = m%FailList(li)%lineIDs(lii+1) ! move subsequent line ID's forward one spot in the list to eliminate this line ID + CALL CheckError(ErrID_Warn, " Line "//trim(num2lstr(m%FailList(l)%lineIDs(il)))//" removed from Failure "//trim(num2lstr(li))//" becasue it already failed by Failure "//trim(num2lstr(l))) + + endif + + ENDDO + ENDDO + + if (m%FailList(li)%nLinesToDetach == 0) then + ! invalid failure + m%FailList(li)%failStatus = 2 + CALL CheckError (ErrID_Warn, " Failure "//trim(num2lstr(li))//" is a duplicate of Failure "//trim(num2lstr(l))//" and will be ignored.") + endif + + endif + endif + ENDDO !li = 1, p%nFails + endif ! m%FailList(l)%failStatus == 1 + + endif ! m%FailList(l)%failStatus == 0 + + ENDDO ! l= 0,nFails CONTAINS + SUBROUTINE DetachLines (attachID, isRod, lineIDs, lineTops, nLinesToDetach, time) + INTEGER(IntKi), INTENT(IN) :: attachID + INTEGER(IntKi), INTENT(IN) :: isRod + INTEGER(IntKi), INTENT(IN ) :: lineIDs(:) + INTEGER(IntKi), INTENT( OUT) :: lineTops(:) + INTEGER(IntKi), INTENT(IN) :: nLinesToDetach + REAL(DbKi), INTENT(IN ) :: time + INTEGER(IntKi) :: k ! index + REAL(DbKi) :: dummyPointState(6) = 0.0_DbKi ! dummy state array to hold kinematics of old attachment point (format in terms of part of point state vector: r[J] = X[3 + J]; rd[J] = X[J]; ) + + ! add point to list of free ones and add states for it + p%nPoints = p%nPoints + 1 ! add 1 to the number of points (this is now the number of the new point) + p%nFreePoints=p%nFreePoints+1 + m%FreePointIs(p%nFreePoints) = p%nPoints + m%PointStateIs1(p%nFreePoints) = m%Nx+1 ! assign start index of this point's states + m%PointStateIsN(p%nFreePoints) = m%Nx+6 + m%Nx = m%Nx + 6 ! add 6 state variables for each point + + ! note: for the DetachLines subroutine, p%nPoints == m%FreePointIs(p%nFreePoints) and can be used interchangeably for indexing. p%nPoints is used to make things easier to read + + ! check to make sure we haven't gone beyond the extra size allotted to the state arrays or the points list <<<< really should throw an error here + if (p%nPoints > p%nPointsExtra) then + CALL CheckError(ErrID_Fatal, " DetachLines: p%nPoints > p%nPointsExtra") + endif + if (m%Nx > m%Nxtra) then + CALL CheckError(ErrID_Fatal, " DetachLines: nX > m%Nx") + endif + + ! create new massless point for detached end(s) of line(s) + m%PointList(p%nPoints)%IdNum = p%nPoints + m%PointList(p%nPoints)%r = 0.0_DbKi ! will be set by Point_SetState + m%PointList(p%nPoints)%rd = 0.0_DbKi ! will be set by Point_SetState + m%PointList(p%nPoints)%pointM = 0.0_DbKi + m%PointList(p%nPoints)%pointV = 0.0_DbKi + m%PointList(p%nPoints)%pointCa = 0.0_DbKi + m%PointList(p%nPoints)%pointCda = 0.0_DbKi + m%PointList(p%nPoints)%typeNum = 0_IntKi ! free point + ! not used + m%PointList(p%nPoints)%pointFX = 0.0_DbKi + m%PointList(p%nPoints)%pointFY = 0.0_DbKi + m%PointList(p%nPoints)%pointFZ = 0.0_DbKi + CALL Point_Initialize(m%PointList(p%nPoints), x%states(m%PointStateIs1(p%nFreePoints) : m%pointStateIsN(p%nFreePoints)), m) + + ! detach lines from old Rod or Point, and get kinematics of the old attachment point + + DO k=1,nLinesToDetach + + if (isRod==1) then ! end A + CALL Rod_RemoveLine(m%RodList(attachID), lineIDs(k), lineTops(k), 0, dummyPointState(4:6), dummyPointState(1:3)) + elseif (isRod==2) then ! end B + CALL Rod_RemoveLine(m%RodList(attachID), lineIDs(k), lineTops(k), 1, dummyPointState(4:6), dummyPointState(1:3)) + elseif (isRod==0) then + CALL Point_RemoveLine(m%PointList(attachID), lineIDs(k), lineTops(k), dummyPointState(4:6), dummyPointState(1:3)) + else + CALL CheckError(ErrID_Fatal, " DetachLines: Failure doesn't have a valid isRod value of 0, 1, or 2.") + endif + + ENDDO + + ! attach lines to new point + DO k=1,nLinesToDetach ! for each relevant line + CALL Point_AddLine(m%PointList(p%nPoints), lineIDs(k), lineTops(k)) + ENDDO + + ! update point kinematics to match old line attachment point kinematics and set positions of attached line ends + CALL Point_SetState(m%PointList(p%nPoints),dummyPointState, time, m) + + ! now make the state vector up to date! + DO k=1,6 + x%states(m%PointStateIs1(p%nFreePoints)+(k-1)) = dummyPointState(k) + ENDDO + + IF (wordy > 0) print *, "Set up new Point ", p%nPoints, " of type ", m%PointList(p%nPoints)%typeNum + + END SUBROUTINE DetachLines + SUBROUTINE CheckError(ErrId, Msg) ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev @@ -2638,14 +3022,16 @@ SUBROUTINE CheckError(ErrId, Msg) ErrMsg = TRIM(ErrMsg)//' MD_UpdateStates:'//TRIM(Msg) ! add current error message ErrStat = MAX(ErrStat, ErrID) - CALL WrScr( ErrMsg ) ! do this always or only if warning level? - if (p%writeLog > 0) then - write(p%UnLog,'(A)') ErrMsg - end if + ! if (ErrStat <= ErrID_Warn) then + ! CALL WrScr( ErrMsg ) ! do this always or only if warning level? + ! endif IF( ErrStat > ErrID_Warn ) THEN + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ErrMsg + end if ! CALL MD_DestroyInput( u_interp, ErrStat, ErrMsg ) - RETURN + RETURN END IF END IF @@ -2884,8 +3270,8 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er INTEGER(IntKi) :: J ! index INTEGER(IntKi) :: K ! index INTEGER(IntKi) :: iTurb ! index -! INTEGER(IntKi) :: Istart ! start index of line/connect in state vector -! INTEGER(IntKi) :: Iend ! end index of line/connect in state vector +! INTEGER(IntKi) :: Istart ! start index of line/point in state vector +! INTEGER(IntKi) :: Iend ! end index of line/point in state vector ! REAL(DbKi) :: temp(3) ! temporary for passing kinematics diff --git a/modules/moordyn/src/MoorDyn_Body.f90 b/modules/moordyn/src/MoorDyn_Body.f90 index acf6f92098..743c3c5791 100644 --- a/modules/moordyn/src/MoorDyn_Body.f90 +++ b/modules/moordyn/src/MoorDyn_Body.f90 @@ -286,7 +286,7 @@ SUBROUTINE Body_SetState(Body, X, t, m) CALL Body_SetDependentKin(Body, t, m) else - print *, "Error: Body::setState called for a non-free Body type in MoorDyn" ! <<< + Call WrScr("Error: Body::setState called for a non-free Body type in MoorDyn") ! <<< end if END SUBROUTINE Body_SetState @@ -407,7 +407,7 @@ SUBROUTINE Body_GetStateDeriv(Body, Xd, m, p) ! check for NaNs (should check all state derivatives, not just first 6) DO J = 1, 6 IF (Is_NaN(Xd(J))) THEN - CALL WrScr("NaN detected at time "//trim(Num2LStr(Body%time))//" in Body "//trim(Int2LStr(Body%IdNum))//"in MoorDyn,") + CALL WrScr("NaN detected at time "//trim(Num2LStr(Body%time))//" in Body "//trim(Int2LStr(Body%IdNum))//" in MoorDyn,") IF (wordy > 0) print *, "state derivatives:" IF (wordy > 0) print *, Xd EXIT @@ -436,6 +436,9 @@ SUBROUTINE Body_DoRHS(Body, m, p) Real(DbKi) :: vi(6) ! relative water velocity (last 3 terms are rotatonal and will be set to zero Real(DbKi) :: F6_i(6) ! net force and moments from an attached object Real(DbKi) :: M6_i(6,6) ! mass and inertia from an attached object + Real(DbKi) :: cda(6) ! body drag coefficients + Real(DbKi) :: cda_t(3,3) = 0.0 ! matrix with translational drag coefficients as diagonals + Real(DbKi) :: cda_r(3,3) = 0.0 ! matrix with rotational drag coefficients as diagonals ! Initialize variables U = 0.0_DbKi ! Set to zero for now @@ -465,8 +468,18 @@ SUBROUTINE Body_DoRHS(Body, m, p) vi(1:3) = U - Body%v6(1:3) ! relative flow velocity over body ref point vi(4:6) = - Body%v6(4:6) ! for rotation, this is just the negative of the body's rotation for now (not allowing flow rotation) - Body%F6net = Body%F6net + 0.5*p%rhoW * vi * abs(vi) * Body%bodyCdA - ! <<< NOTE, for body this should be fixed to account for orientation!! <<< what about drag in rotational DOFs??? <<<<<<<<<<<<<< + cda_t(1,1) = Body%bodyCdA(1) + cda_t(2,2) = Body%bodyCdA(2) + cda_t(3,3) = Body%bodyCdA(3) + cda_r(1,1) = Body%bodyCdA(4) + cda_r(2,2) = Body%bodyCdA(5) + cda_r(3,3) = Body%bodyCdA(6) + + cda(1:3) = MATMUL( MATMUL( MATMUL(Body%OrMat,cda_t) , transpose(Body%OrMat) ) , vi(1:3) * norm2(vi(1:3)) ); + cda(4:6) = MATMUL( MATMUL( MATMUL(Body%OrMat,cda_r) , transpose(Body%OrMat) ) , vi(4:6) * norm2(vi(4:6)) ); + Body%F6net = Body%F6net + 0.5*p%rhoW*cda + + @@ -477,7 +490,7 @@ SUBROUTINE Body_DoRHS(Body, m, p) CALL Point_GetNetForceAndMass( m%PointList(Body%attachedC(l)), Body%r6(1:3), F6_i, M6_i, m, p) if (ABS(F6_i(5)) > 1.0E12) then - print *, "Warning: extreme pitch moment from body-attached Point ", l + Call WrScr( "Warning: extreme pitch moment from body-attached Point "//trim(num2lstr(l))) end if ! sum quantitites @@ -492,7 +505,7 @@ SUBROUTINE Body_DoRHS(Body, m, p) CALL Rod_GetNetForceAndMass(m%RodList(Body%attachedR(l)), Body%r6(1:3), F6_i, M6_i, m, p) if (ABS(F6_i(5)) > 1.0E12) then - print *, "Warning: extreme pitch moment from body-attached Rod ", l + Call WrScr("Warning: extreme pitch moment from body-attached Rod "//trim(num2lstr(l))) end if ! sum quantitites @@ -546,7 +559,7 @@ SUBROUTINE Body_GetCoupledForce(Body, Fnet_out, m, p) Fnet_out = Body%F6net else - print *, "ERROR, Body_GetCoupledForce called for wrong (non-coupled) body type in MoorDyn!" + Call WrScr( "ERROR, Body_GetCoupledForce called for wrong (non-coupled) body type in MoorDyn!") end if END SUBROUTINE Body_GetCoupledForce diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index 5a5f319a09..82ec10f977 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -118,6 +118,7 @@ MODULE MoorDyn_IO PUBLIC :: MDIO_CloseOutput PUBLIC :: MDIO_ProcessOutList PUBLIC :: MDIO_WriteOutputs + PUBLIC :: Line_GetNodeTen CONTAINS @@ -1365,11 +1366,11 @@ SUBROUTINE MDIO_WriteOutputs( Time, p, m, y, ErrStat, ErrMsg ) CASE (FZ) y%WriteOutput(I) = m%LineList(p%OutParam(I)%ObjID)%Fnet(3,p%OutParam(I)%NodeID) ! node force in z CASE (Ten) - y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), p%OutParam(I)%NodeID, p) ! this is actually the segment tension ( 1 < NodeID < N ) Should deal with properly! + y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), p%OutParam(I)%NodeID) ! this is actually the segment tension ( 1 < NodeID < N ) Should deal with properly! CASE (TenA) - y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), 0, p) + y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), 0) CASE (TenB) - y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), m%LineList(p%OutParam(I)%ObjID)%N, p) + y%WriteOutput(I) = Line_GetNodeTen(m%LineList(p%OutParam(I)%ObjID), m%LineList(p%OutParam(I)%ObjID)%N) CASE DEFAULT y%WriteOutput(I) = 0.0_ReKi ErrStat = ErrID_Warn @@ -1886,11 +1887,10 @@ END SUBROUTINE MDIO_WriteOutputs ! get tension at any node including fairlead or anchor (accounting for weight in these latter cases) !-------------------------------------------------------------- - FUNCTION Line_GetNodeTen(Line, i, p) result(NodeTen) + FUNCTION Line_GetNodeTen(Line, i) result(NodeTen) TYPE(MD_Line), INTENT(IN ) :: Line ! label for the current line, for convenience INTEGER(IntKi), INTENT(IN ) :: i ! node index to get tension at - TYPE(MD_ParameterType), INTENT(IN ) :: p ! Parameters REAL(DbKi) :: NodeTen ! returned calculation of tension at node INTEGER(IntKi) :: J diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 190cc4d7eb..a6be2b217e 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -357,12 +357,8 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) Line%r(2,J) = Line%r(2,0) + (Line%r(2,N) - Line%r(2,0))*REAL(J, DbKi)/REAL(N, DbKi) Line%r(3,J) = Line%r(3,0) + (Line%r(3,N) - Line%r(3,0))*REAL(J, DbKi)/REAL(N, DbKi) - ! print*, Line%r(:,J) ENDDO - ! print*,"FYI line end A and B node coords are" - ! print*, Line%r(:,0) - ! print*, Line%r(:,N) ENDIF ENDIF @@ -1464,7 +1460,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p) !, FairFtot, FairMtot, AnchFtot, ! check for NaNs DO J = 1, 6*(N-1) IF (Is_NaN(Xd(J))) THEN - print *, "NaN detected at time ", Line%time, " in Line ", Line%IdNum, " in MoorDyn." + Call WrScr( "NaN detected at time "//trim(num2lstr(Line%time))//" in Line "//trim(num2lstr(Line%IdNum))//" in MoorDyn.") IF (wordy > 1) THEN print *, "state derivatives:" print *, Xd diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index c54dea1f2e..cfc82ed6f4 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -1500,6 +1500,10 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ! Close the inputs file CLOSE ( UnElev ) + + IF (WaveTimeIn(1) .NE. 0.0) THEN + CALL SetErrStat( ErrID_Warn, ' MoorDyn WaveElev time series should start at t = 0 seconds. First two lines are read as headers.',ErrStat, ErrMsg, RoutineName); return + ENDIF call WrScr( "Read "//trim(num2lstr(ntIn))//" time steps from input file." ) diff --git a/modules/moordyn/src/MoorDyn_Point.f90 b/modules/moordyn/src/MoorDyn_Point.f90 index 9328f4805d..8f5d87f07f 100644 --- a/modules/moordyn/src/MoorDyn_Point.f90 +++ b/modules/moordyn/src/MoorDyn_Point.f90 @@ -97,24 +97,18 @@ SUBROUTINE Point_SetKinematics(Point, r_in, rd_in, a_in, t, m) Point%time = t - ! if (Point%typeNum==0) THEN ! anchor ( <<< to be changed/expanded) ... in MoorDyn F also used for coupled points - - ! set position and velocity - Point%r = r_in - Point%rd = rd_in - Point%a = a_in - - ! pass latest kinematics to any attached lines - DO l=1,Point%nAttached - CALL Line_SetEndKinematics(m%LineList(Point%attached(l)), Point%r, Point%rd, t, Point%Top(l)) - END DO - - ! else - ! - ! PRINT*,"Error: setKinematics called for wrong Point type. Point ", Point%IdNum, " type ", Point%typeNum - - ! END IF - + + ! set position and velocity + Point%r = r_in + Point%rd = rd_in + Point%a = a_in + + ! pass latest kinematics to any attached lines + DO l=1,Point%nAttached + CALL Line_SetEndKinematics(m%LineList(Point%attached(l)), Point%r, Point%rd, t, Point%Top(l)) + END DO + + END SUBROUTINE Point_SetKinematics !-------------------------------------------------------------- @@ -379,6 +373,7 @@ SUBROUTINE Point_RemoveLine(Point, lineID, TopOfLine, rEnd, rdEnd) REAL(DbKi), INTENT(INOUT) :: rdEnd(3) Integer(IntKi) :: l,m,J + Integer(IntKi) :: found = 0 DO l = 1,Point%nAttached ! look through attached lines @@ -386,7 +381,7 @@ SUBROUTINE Point_RemoveLine(Point, lineID, TopOfLine, rEnd, rdEnd) TopOfLine = Point%Top(l); ! record which end of the line was attached - DO m = l,Point%nAttached-1 + DO m = l,Point%nAttached Point%Attached(m) = Point%Attached(m+1) ! move subsequent line links forward one spot in the list to eliminate this line link Point%Top( m) = Point%Top(m+1) @@ -404,13 +399,15 @@ SUBROUTINE Point_RemoveLine(Point, lineID, TopOfLine, rEnd, rdEnd) EXIT END DO - IF (l == Point%nAttached) THEN ! detect if line not found - print *, "Error: failed to find line to remove during removeLineFromPoint call to point ", Point%IdNum, ". Line ", lineID - END IF + found = 1 END IF END DO + + IF (found == 0) THEN ! detect if line not found TODO: fix this, its wrong. If pointNnattached is oprginally 2, then it will be 1 after one run of the loop and l will also be 1 + CALL WrScr("Error: failed to find line to remove during RemoveLine call to Point "//trim(num2lstr(Point%IdNum))//". Line "//trim(num2lstr(lineID))) + END IF END SUBROUTINE Point_RemoveLine diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 283fca0bed..f2b509f9cf 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -265,9 +265,16 @@ typedef ^ ^ DbKi EndMomentB {3} typedef ^ ^ IntKi LineUnOut - - - "unit number of line output file" typedef ^ ^ DbKi LineWrOutput {:} - - "one row of output data for this line" -# this is the Fail type, which holds data for possible line failure descriptors TO BE FILLED IN LATER -typedef ^ MD_Fail IntKi IdNum - - - "integer identifier of this failure" - +# this is the Fail type, which holds data for possible line failure descriptors +typedef ^ MD_Fail IntKi IdNum - - - "integer identifier of this failure" "-" +typedef ^ ^ IntKi attachID - - - "ID of connection or Rod the lines are attached to" "-" +typedef ^ ^ IntKi isRod - - - "1 Rod end A, 2 Rod end B, 0 if point" "-" +typedef ^ ^ IntKi lineIDs {30} - - "array of one or more lines to detach (starting from 1...)" "-" +typedef ^ ^ IntKi lineTops {30} - - "an array that will be FILLED IN to return which end of each line was disconnected ... 1 = top/fairlead(end B), 0 = bottom/anchor(end A)" "-" +typedef ^ ^ IntKi nLinesToDetach - - - "how many lines to dettach" "-" +typedef ^ ^ DbKi failTime - - - "time of failure" "s" +typedef ^ ^ DbKi failTen - - - "tension threshold of failure" "N" +typedef ^ ^ IntKi failStatus - - - "0 not failed yet, 1 failed, 2 invalid" "-" # this is the MDOutParmType - a less literal alternative of the NWTC OutParmType for MoorDyn (to avoid huge lists of possible output channel permutations) typedef ^ MD_OutParmType CHARACTER(10) Name - - - "name of output channel" @@ -333,6 +340,7 @@ typedef ^ ^ IntKi RodStateIsN {:} typedef ^ ^ IntKi BodyStateIs1 {:} - - "starting index of each body's states in state vector" "" typedef ^ ^ IntKi BodyStateIsN {:} - - "ending index of each body's states in state vector" "" typedef ^ ^ IntKi Nx - - - "number of states and size of state vector" "" +typedef ^ ^ IntKi Nxtra - - - "number of states and size of state vector including points for potential line failures" "" typedef ^ ^ IntKi WaveTi - - - "current interpolation index for wave time series data" "" typedef ^ ^ MD_ContinuousStateType xTemp - - - "contains temporary state vector used in integration (put here so it's only allocated once)" typedef ^ ^ MD_ContinuousStateType xdTemp - - - "contains temporary state derivative vector used in integration (put here so it's only allocated once)" diff --git a/modules/moordyn/src/MoorDyn_Rod.f90 b/modules/moordyn/src/MoorDyn_Rod.f90 index f5f718198a..50bd217a0c 100644 --- a/modules/moordyn/src/MoorDyn_Rod.f90 +++ b/modules/moordyn/src/MoorDyn_Rod.f90 @@ -259,7 +259,7 @@ SUBROUTINE Rod_SetKinematics(Rod, r6_in, v6_in, a6_in, t, m) ! handled, along with passing kinematics to dependent lines, by separate call to setState else - print *, "Error: Rod_SetKinematics called for a free Rod in MoorDyn." ! <<< + Call WrScr("Error: Rod_SetKinematics called for a free Rod in MoorDyn. Rod number"//trim(num2lstr(Rod%IdNum))) ! <<< end if @@ -324,7 +324,7 @@ SUBROUTINE Rod_SetState(Rod, X, t, m) CALL Rod_SetDependentKin(Rod, t, m, .FALSE.) else - print *, "Error: Rod::setState called for a non-free rod type in MoorDyn" ! <<< + Call WrScr("Error: Rod::setState called for a non-free rod type in MoorDyn") ! <<< end if ! update Rod direction unit vector (simply equal to last three entries of r6) @@ -1018,7 +1018,7 @@ SUBROUTINE Rod_GetCoupledForce(Rod, Fnet_out, m, p) Rod%F6net(4:6) = 0.0_DbKi Fnet_out = Rod%F6net else - print *, "ERROR, Rod_GetCoupledForce called for wrong (non-coupled) rod type!" + Call WrScr("ERROR, Rod_GetCoupledForce called for wrong (non-coupled) rod type!") end if END SUBROUTINE Rod_GetCoupledForce @@ -1113,6 +1113,7 @@ SUBROUTINE Rod_RemoveLine(Rod, lineID, TopOfLine, endB, rEnd, rdEnd) REAL(DbKi), INTENT(INOUT) :: rdEnd(3) Integer(IntKi) :: l,m,J + Integer(IntKi) :: foundA, foundB = 0 if (endB==1) then ! attaching to end B @@ -1122,7 +1123,7 @@ SUBROUTINE Rod_RemoveLine(Rod, lineID, TopOfLine, endB, rEnd, rdEnd) TopOfLine = Rod%TopB(l); ! record which end of the line was attached - DO m = l,Rod%nAttachedB-1 + DO m = l,Rod%nAttachedB Rod%AttachedB(m) = Rod%AttachedB(m+1) ! move subsequent line links forward one spot in the list to eliminate this line link Rod%TopB( m) = Rod%TopB(m+1) @@ -1135,17 +1136,19 @@ SUBROUTINE Rod_RemoveLine(Rod, lineID, TopOfLine, endB, rEnd, rdEnd) rdEnd(J) = Rod%rd(J,Rod%N) END DO - call WrScr( "Detached line "//trim(num2lstr(lineID))//" from Rod "//trim(num2lstr(Rod%IdNum))//" end B") + CALL WrScr( "Detached line "//trim(num2lstr(lineID))//" from Rod "//trim(num2lstr(Rod%IdNum))//" end B") EXIT END DO - - IF (l == Rod%nAttachedB) THEN ! detect if line not found - print *, "Error: failed to find line to remove during RemoveLine call to Rod ", Rod%IdNum, ". Line ", lineID - END IF + + foundB = 1 + END IF END DO - + IF (foundB == 0) THEN ! detect if line not found + CALL WrScr("Error: failed to find line to remove during RemoveLine call to Rod "//trim(num2lstr(Rod%IdNum))//" end B. Line "//trim(num2lstr(lineID))) + END IF + else ! attaching to end A DO l = 1,Rod%nAttachedA ! look through attached lines @@ -1154,7 +1157,7 @@ SUBROUTINE Rod_RemoveLine(Rod, lineID, TopOfLine, endB, rEnd, rdEnd) TopOfLine = Rod%TopA(l); ! record which end of the line was attached - DO m = l,Rod%nAttachedA-1 + DO m = l,Rod%nAttachedA Rod%AttachedA(m) = Rod%AttachedA(m+1) ! move subsequent line links forward one spot in the list to eliminate this line link Rod%TopA( m) = Rod%TopA(m+1) @@ -1167,16 +1170,19 @@ SUBROUTINE Rod_RemoveLine(Rod, lineID, TopOfLine, endB, rEnd, rdEnd) rdEnd(J) = Rod%rd(J,0) END DO - call WrScr( "Detached line "//trim(num2lstr(lineID))//" from Rod "//trim(num2lstr(Rod%IdNum))//" end A") + CALL WrScr( "Detached line "//trim(num2lstr(lineID))//" from Rod "//trim(num2lstr(Rod%IdNum))//" end A") EXIT END DO - - IF (l == Rod%nAttachedA) THEN ! detect if line not found - print *, "Error: failed to find line to remove during RemoveLine call to Rod ", Rod%IdNum, ". Line ", lineID - END IF + + foundA = 1 + END IF END DO + + IF (foundA == 0) THEN ! detect if line not found + CALL WrScr("Error: failed to find line to remove during RemoveLine call to Rod "//trim(num2lstr(Rod%IdNum))//" end A. Line "//trim(num2lstr(lineID))) + END IF end if diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index e862234ef1..a9361da552 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -290,6 +290,14 @@ MODULE MoorDyn_Types ! ========= MD_Fail ======= TYPE, PUBLIC :: MD_Fail INTEGER(IntKi) :: IdNum = 0_IntKi !< integer identifier of this failure [-] + INTEGER(IntKi) :: attachID = 0_IntKi !< ID of connection or Rod the lines are attached to [-] + INTEGER(IntKi) :: isRod = 0_IntKi !< 1 Rod end A, 2 Rod end B, 0 if point [-] + INTEGER(IntKi) , DIMENSION(1:30) :: lineIDs = 0_IntKi !< array of one or more lines to detach (starting from 1...) [-] + INTEGER(IntKi) , DIMENSION(1:30) :: lineTops = 0_IntKi !< an array that will be FILLED IN to return which end of each line was disconnected ... 1 = top/fairlead(end B), 0 = bottom/anchor(end A) [-] + INTEGER(IntKi) :: nLinesToDetach = 0_IntKi !< how many lines to dettach [-] + REAL(DbKi) :: failTime = 0.0_R8Ki !< time of failure [s] + REAL(DbKi) :: failTen = 0.0_R8Ki !< tension threshold of failure [N] + INTEGER(IntKi) :: failStatus = 0_IntKi !< 0 not failed yet, 1 failed, 2 invalid [-] END TYPE MD_Fail ! ======================= ! ========= MD_OutParmType ======= @@ -368,6 +376,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: BodyStateIs1 !< starting index of each body's states in state vector [] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: BodyStateIsN !< ending index of each body's states in state vector [] INTEGER(IntKi) :: Nx = 0_IntKi !< number of states and size of state vector [] + INTEGER(IntKi) :: Nxtra = 0_IntKi !< number of states and size of state vector including points for potential line failures [] INTEGER(IntKi) :: WaveTi = 0_IntKi !< current interpolation index for wave time series data [] TYPE(MD_ContinuousStateType) :: xTemp !< contains temporary state vector used in integration (put here so it's only allocated once) [-] TYPE(MD_ContinuousStateType) :: xdTemp !< contains temporary state derivative vector used in integration (put here so it's only allocated once) [-] @@ -2173,6 +2182,14 @@ subroutine MD_CopyFail(SrcFailData, DstFailData, CtrlCode, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = '' DstFailData%IdNum = SrcFailData%IdNum + DstFailData%attachID = SrcFailData%attachID + DstFailData%isRod = SrcFailData%isRod + DstFailData%lineIDs = SrcFailData%lineIDs + DstFailData%lineTops = SrcFailData%lineTops + DstFailData%nLinesToDetach = SrcFailData%nLinesToDetach + DstFailData%failTime = SrcFailData%failTime + DstFailData%failTen = SrcFailData%failTen + DstFailData%failStatus = SrcFailData%failStatus end subroutine subroutine MD_DestroyFail(FailData, ErrStat, ErrMsg) @@ -2190,6 +2207,14 @@ subroutine MD_PackFail(RF, Indata) character(*), parameter :: RoutineName = 'MD_PackFail' if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%IdNum) + call RegPack(RF, InData%attachID) + call RegPack(RF, InData%isRod) + call RegPack(RF, InData%lineIDs) + call RegPack(RF, InData%lineTops) + call RegPack(RF, InData%nLinesToDetach) + call RegPack(RF, InData%failTime) + call RegPack(RF, InData%failTen) + call RegPack(RF, InData%failStatus) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -2199,6 +2224,14 @@ subroutine MD_UnPackFail(RF, OutData) character(*), parameter :: RoutineName = 'MD_UnPackFail' if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%IdNum); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%attachID); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%isRod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%lineIDs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%lineTops); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%nLinesToDetach); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%failTime); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%failTen); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%failStatus); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine MD_CopyOutParmType(SrcOutParmTypeData, DstOutParmTypeData, CtrlCode, ErrStat, ErrMsg) @@ -3016,6 +3049,7 @@ subroutine MD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) DstMiscData%BodyStateIsN = SrcMiscData%BodyStateIsN end if DstMiscData%Nx = SrcMiscData%Nx + DstMiscData%Nxtra = SrcMiscData%Nxtra DstMiscData%WaveTi = SrcMiscData%WaveTi call MD_CopyContState(SrcMiscData%xTemp, DstMiscData%xTemp, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -3313,6 +3347,7 @@ subroutine MD_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%BodyStateIs1) call RegPackAlloc(RF, InData%BodyStateIsN) call RegPack(RF, InData%Nx) + call RegPack(RF, InData%Nxtra) call RegPack(RF, InData%WaveTi) call MD_PackContState(RF, InData%xTemp) call MD_PackContState(RF, InData%xdTemp) @@ -3443,6 +3478,7 @@ subroutine MD_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%BodyStateIs1); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BodyStateIsN); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Nx); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Nxtra); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WaveTi); if (RegCheckErr(RF, RoutineName)) return call MD_UnpackContState(RF, OutData%xTemp) ! xTemp call MD_UnpackContState(RF, OutData%xdTemp) ! xdTemp diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 07b182092d..769810a590 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -468,6 +468,7 @@ seast_regression("seastate_wavemod5" "seastate") # pla # MoorDyn regression tests md_regression("md_5MW_OC4Semi" "moordyn") +md_regression("md_lineFail" "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_Node_Check_N20" "moordyn") diff --git a/reg_tests/r-test b/reg_tests/r-test index 3210b20016..c2c1f07c99 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 3210b20016c5319ff1b0c8fa7b41cdf4af5da871 +Subproject commit c2c1f07c99eaeb4572c9a51579514302f3212a5b