Skip to content

Line Failures added to MoorDyn #2214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrew-platt merged 11 commits into from
May 17, 2024
Merged

Line Failures added to MoorDyn #2214

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
c26f800
Adding bodies and rods to FF
RyanDavies19 Apr 2, 2024
06815b0
Merge remote-tracking branch 'openfast/dev' into dev
RyanDavies19 Apr 22, 2024
462d6d5
Added rotational drag and fixed orientation for drag coefficients
RyanDavies19 Apr 22, 2024
41ecd9f
Line loss/failures added to MoorDyn
RyanDavies19 May 10, 2024
b1f509a
updated registry
RyanDavies19 May 10, 2024
5776123
Updated registry and types, added warning for wave elevation profile
RyanDavies19 May 13, 2024
327b61e
Merge remote-tracking branch 'openfast/dev' into failures
RyanDavies19 May 13, 2024
e5fe808
Cleaning up print statements
RyanDavies19 May 14, 2024
6888ad0
Merge remote-tracking branch 'openfast/dev' into failures
RyanDavies19 May 16, 2024
250ff38
Updated Types
RyanDavies19 May 16, 2024
dddeaf6
MD: add line failure regression test
andrew-platt May 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
508 changes: 447 additions & 61 deletions modules/moordyn/src/MoorDyn.f90
View file
Open in desktop

Large diffs are not rendered by default.

27 changes: 20 additions & 7 deletions modules/moordyn/src/MoorDyn_Body.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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





Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions modules/moordyn/src/MoorDyn_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,6 +118,7 @@ MODULE MoorDyn_IO
PUBLIC :: MDIO_CloseOutput
PUBLIC :: MDIO_ProcessOutList
PUBLIC :: MDIO_WriteOutputs
PUBLIC :: Line_GetNodeTen


CONTAINS
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
6 changes: 1 addition & 5 deletions modules/moordyn/src/MoorDyn_Line.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions modules/moordyn/src/MoorDyn_Misc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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." )

Expand Down
41 changes: 19 additions & 22 deletions modules/moordyn/src/MoorDyn_Point.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
!--------------------------------------------------------------
Expand Down Expand Up @@ -379,14 +373,15 @@ SUBROUTINE Point_RemoveLine(Point, lineID, TopOfLine, rEnd, rdEnd)
REAL(DbKi), INTENT(INOUT) :: rdEnd(3)

Integer(IntKi) :: l,m,J
Integer(IntKi) :: found = 0
Copy link
Contributor

@bjonkman bjonkman Aug 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In Fortran, initializing a variable in the type declaration statement automatically gives it a "SAVE" attribute, which means it retains its value each time it is called (i.e., it's set to 0 only the first time the routine is called). Since this is called from an UpdateStates routine, I'd be surprised if this is how it was intended to be used.... and if that was the intention, it wouldn't work properly in a restart situation.

I'd recommend either

  1. rewriting like this:
Integer(IntKi) :: found

found = 0

or

  1. putting this variable in a MiscVar type, if you do intend for it to have the SAVE attribute.

Copy link
Contributor

@bjonkman bjonkman Aug 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, using a LOGICAL here might make it easier to read, unless you have plans for checking something other than found vs not-found .

Copy link
Contributor Author

@RyanDavies19 RyanDavies19 Aug 19, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @bjonkman for the review. @andrew-platt if you have a PR or something merging soon can you add the found = 0 and the logical fix there? Otherwise I can open a PR. It's intent is a boolean flag that prints an error message (I probably should've made it throw a warning, because if it fails the simulation continues with the line attached)

Copy link
Collaborator

@andrew-platt andrew-platt Aug 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll include this into a minor PR that I'll post tomorrow.

Copy link
Contributor Author

@RyanDavies19 RyanDavies19 Aug 19, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks!

Copy link
Collaborator

@andrew-platt andrew-platt Aug 29, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See #2400


DO l = 1,Point%nAttached ! look through attached lines

IF (Point%Attached(l) == lineID) THEN ! if this is the line's entry in the attachment list

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)
Expand All @@ -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

Expand Down
14 changes: 11 additions & 3 deletions modules/moordyn/src/MoorDyn_Registry.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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)"
Expand Down
38 changes: 22 additions & 16 deletions modules/moordyn/src/MoorDyn_Rod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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

Expand Down
Loading