Skip to content

Commit ec9b169

Browse files
andrew-plattandrew-platt
authored
Merge pull request #2782 from luwang00/f/SDMultiRotor
SubDyn upgrade to support multiple transition pieces
2 parents 06d2505 + 08633ac commit ec9b169

File tree

13 files changed

+33370
-31686
lines changed

13 files changed

+33370
-31686
lines changed

‎modules/openfast-library/src/FAST_Mapping.f90‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1144,8 +1144,8 @@ subroutine InitMappings_ED(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg)
11441144
case (Module_SD)
11451145

11461146
call MapLoadMesh(Turbine, Mappings, SrcMod=SrcMod, DstMod=DstMod, &
1147-
SrcDL=DatLoc(SD_y_Y1Mesh), & ! SD%y%Y1mesh, &
1148-
SrcDispDL=DatLoc(SD_u_TPMesh), & ! SD%u%TPMesh
1147+
SrcDL=DatLoc(SD_y_Y1Mesh,DstMod%Ins), & ! SD%y%Y1mesh, &
1148+
SrcDispDL=DatLoc(SD_u_TPMesh,DstMod%Ins), & ! SD%u%TPMesh
11491149
DstDL=DatLoc(ED_u_PlatformPtMesh), & ! ED%u%PlatformPtMesh
11501150
DstDispDL=DatLoc(ED_y_PlatformPtMesh), & ! ED%y%PlatformPtMesh
11511151
ErrStat=ErrStat2, ErrMsg=ErrMsg2)
@@ -1808,8 +1808,8 @@ subroutine InitMappings_SD(Mappings, SrcMod, DstMod, Turbine, ErrStat, ErrMsg)
18081808
case (Module_ED)
18091809

18101810
call MapMotionMesh(Turbine, Mappings, SrcMod=SrcMod, DstMod=DstMod, &
1811-
SrcDL=DatLoc(ED_y_PlatformPtMesh), & ! ED%y%PlatformPtMesh
1812-
DstDL=DatLoc(SD_u_TPMesh), & ! SD%u%TPMesh
1811+
SrcDL=DatLoc(ED_y_PlatformPtMesh), & ! ED%y%PlatformPtMesh
1812+
DstDL=DatLoc(SD_u_TPMesh, SrcMod%Ins), & ! SD%u%TPMesh
18131813
ErrStat=ErrStat2, ErrMsg=ErrMsg2); if(Failed()) return
18141814

18151815
case (Module_FEAM)

‎modules/openfast-library/src/FAST_Subs.f90‎

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -924,7 +924,11 @@ SUBROUTINE FAST_InitializeAll( t_initial, m_Glue, p_FAST, y_FAST, m_FAST, ED, SE
924924
Init%InData_SD%g = p_FAST%Gravity
925925
Init%InData_SD%SDInputFile = p_FAST%SubFile
926926
Init%InData_SD%RootName = p_FAST%OutFileRoot
927-
Init%InData_SD%TP_RefPoint = ED%y(iED)%PlatformPtMesh%Position(:,1) ! "Interface point" where loads will be transferred to
927+
! TODO: Below is only a temporary patch to keep things working as before. Need to be updated to support multirotor capability.
928+
Init%InData_SD%nTP = 1
929+
Allocate(Init%InData_SD%TP_RefPoint(3,Init%InData_SD%nTP), stat=ErrStat2)
930+
CALL SetErrStat(ErrStat2,'Allocating TP reference points',ErrStat,ErrMsg,RoutineName)
931+
Init%InData_SD%TP_RefPoint(:,iED) = ED%y(iED)%PlatformPtMesh%Position(:,1) ! "Interface points" where loads and motion will be transferred to and from
928932
Init%InData_SD%SubRotateZ = 0.0 ! Used by driver to rotate structure around z
929933

930934
CALL SD_Init( Init%InData_SD, SD%Input(1), SD%p, SD%x(STATE_CURR), SD%xd(STATE_CURR), SD%z(STATE_CURR), &
@@ -5902,8 +5906,9 @@ SUBROUTINE WrVTK_AllMeshes(p_FAST, y_FAST, ED, SED, BD, AD, IfW, ExtInfw, HD, SD
59025906
!call MeshWrVTK(p_FAST%TurbinePos, SD%Input(1)%TPMesh, trim(p_FAST%VTK_OutFileRoot)//'.SD_TPMesh_motion', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth )
59035907
call MeshWrVTK(p_FAST%TurbinePos, SD%Input(1)%LMesh, trim(p_FAST%VTK_OutFileRoot)//'.SD_LMesh_y2Mesh', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, SD%y%y2Mesh )
59045908
call MeshWrVTK(p_FAST%TurbinePos, SD%Input(1)%LMesh, trim(p_FAST%VTK_OutFileRoot)//'.SD_LMesh_y3Mesh', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, SD%y%y3Mesh )
5905-
5906-
call MeshWrVTK(p_FAST%TurbinePos, SD%y%y1Mesh, trim(p_FAST%VTK_OutFileRoot)//'.SD_y1Mesh_TPMesh', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, SD%Input(1)%TPMesh )
5909+
DO j = 1,size(SD%y%y1Mesh)
5910+
call MeshWrVTK(p_FAST%TurbinePos, SD%y%y1Mesh(j), trim(p_FAST%VTK_OutFileRoot)//'.SD_y1Mesh_TPMesh'//trim(num2lstr(j)), y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, SD%Input(1)%TPMesh(j) )
5911+
END DO
59075912
!call MeshWrVTK(p_FAST%TurbinePos, SD%y%y3Mesh, trim(p_FAST%VTK_OutFileRoot)//'.SD_y3Mesh_motion', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth )
59085913
ELSE IF ( p_FAST%CompSub == Module_ExtPtfm .and. allocated(ExtPtfm%Input)) THEN
59095914
call MeshWrVTK(p_FAST%TurbinePos, ExtPtfm%y%PtfmMesh, trim(p_FAST%VTK_OutFileRoot)//'.ExtPtfm', y_FAST%VTK_count, p_FAST%VTK_fields, ErrStat2, ErrMsg2, p_FAST%VTK_tWidth, ExtPtfm%Input(1)%PtfmMesh )
@@ -6460,7 +6465,9 @@ SUBROUTINE WriteInputMeshesToFile(u_ED, u_AD, u_SD, u_HD, u_MAP, u_BD, FileName,
64606465
CALL MeshWrBin( unOut, u_ED(J_local)%TowerPtLoads, ErrStat, ErrMsg )
64616466
CALL MeshWrBin( unOut, u_ED(J_local)%PlatformPtMesh, ErrStat, ErrMsg )
64626467
end do
6463-
CALL MeshWrBin( unOut, u_SD%TPMesh, ErrStat, ErrMsg )
6468+
do J_local = 1,size(u_SD%TPMesh)
6469+
CALL MeshWrBin( unOut, u_SD%TPMesh(J_local), ErrStat, ErrMsg )
6470+
end do
64646471
CALL MeshWrBin( unOut, u_SD%LMesh, ErrStat, ErrMsg )
64656472
CALL MeshWrBin( unOut, u_HD%Morison%Mesh, ErrStat, ErrMsg )
64666473
CALL MeshWrBin( unOut, u_HD%WAMITMesh, ErrStat, ErrMsg )
@@ -6550,7 +6557,9 @@ SUBROUTINE WriteMotionMeshesToFile(time, y_ED, u_SD, y_SD, u_HD, u_MAP, y_BD, u_
65506557
CALL MeshWrBin( unOut, y_ED(J_local)%TowerLn2Mesh, ErrStat, ErrMsg )
65516558
CALL MeshWrBin( unOut, y_ED(J_local)%PlatformPtMesh, ErrStat, ErrMsg )
65526559
end do
6553-
CALL MeshWrBin( unOut, u_SD%TPMesh, ErrStat, ErrMsg )
6560+
do J_local = 1,size(u_SD%TPMesh)
6561+
CALL MeshWrBin( unOut, u_SD%TPMesh(J_local), ErrStat, ErrMsg )
6562+
end do
65546563
CALL MeshWrBin( unOut, y_SD%y2Mesh, ErrStat, ErrMsg )
65556564
CALL MeshWrBin( unOut, y_SD%y3Mesh, ErrStat, ErrMsg )
65566565
CALL MeshWrBin( unOut, u_HD%Morison%Mesh, ErrStat, ErrMsg )

‎modules/subdyn/src/FEM.f90‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1586,7 +1586,6 @@ SUBROUTINE PseudoInverse(A, Ainv, ErrStat, ErrMsg)
15861586
end do
15871587
! Compute Ainv = 1.0*V^t * U^t + 0.0*Ainv V*(inv(S))*U'
15881588
!call DGEMM( 'T', 'T', N, M, K, 1.0, V, K, U, M, 0.0, Ainv, N)
1589-
print*,'8'
15901589
call LAPACK_GEMM( 'T', 'T', 1.0_FEKi, Vt, U, 0.0_FEKi, Ainv, ErrStat, ErrMsg)
15911590
! --- Compute rank
15921591
!tol=maxval(shape(A))*epsilon(maxval(S))

‎modules/subdyn/src/SD_FEM.f90‎

Lines changed: 101 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ MODULE SD_FEM
2525

2626
INTEGER(IntKi), PARAMETER :: MaxMemJnt = 20 ! Maximum number of members at one joint
2727
INTEGER(IntKi), PARAMETER :: MaxOutChs = 2000 ! Max number of Output Channels to be read in
28-
INTEGER(IntKi), PARAMETER :: nDOFL_TP = 6 !TODO rename me ! 6 degrees of freedom (length of u subarray [UTP])
28+
INTEGER(IntKi), PARAMETER :: nDOFL_TP = 6 !TODO rename me ! 6 degrees of freedom (length of u subarray [UTP]), current used for output only. Need to update this to reflect multiple transition pieces.
2929

3030
! values of these parameters are ordered by their place in SubDyn input file:
3131
INTEGER(IntKi), PARAMETER :: JointsCol = 9 ! Number of columns in Joints (JointID, JointXss, JointYss, JointZss, JointType, JointDirX JointDirY JointDirZ JointStiff)
@@ -78,8 +78,8 @@ MODULE SD_FEM
7878
! Types of Guyan Damping
7979
INTEGER(IntKi), PARAMETER :: idGuyanDamp_None = 0
8080
INTEGER(IntKi), PARAMETER :: idGuyanDamp_Rayleigh = 1
81-
INTEGER(IntKi), PARAMETER :: idGuyanDamp_66 = 2
82-
INTEGER(IntKi) :: idGuyanDamp_Valid(3) = (/idGuyanDamp_None, idGuyanDamp_Rayleigh, idGuyanDamp_66 /)
81+
INTEGER(IntKi), PARAMETER :: idGuyanDamp_Matrix = 2
82+
INTEGER(IntKi) :: idGuyanDamp_Valid(3) = (/idGuyanDamp_None, idGuyanDamp_Rayleigh, idGuyanDamp_Matrix /)
8383

8484
INTEGER(IntKi), PARAMETER :: SDMaxInpCols = MAX(JointsCol,InterfCol,MembersCol,PropSetsBCCol,PropSetsBRCol,PropSetsXCol,PropSetsSCol,COSMsCol,CMassCol)
8585

@@ -91,7 +91,7 @@ MODULE SD_FEM
9191
! Implementation Flags
9292
LOGICAL, PARAMETER :: DEV_VERSION = .false.
9393
LOGICAL, PARAMETER :: BC_Before_CB = .true.
94-
LOGICAL, PARAMETER :: ANALYTICAL_LIN = .true.
94+
LOGICAL, PARAMETER :: ANALYTICAL_LIN = .false. ! Analytical linearization is no longer valid in general
9595
LOGICAL, PARAMETER :: GUYAN_RIGID_FLOATING = .true.
9696

9797

@@ -243,17 +243,18 @@ END FUNCTION NodeHasRigidElem
243243
!! Typically called to get:
244244
!! - the transformation from the interface points to the TP point
245245
!! - the transformation from the bottom nodes to SubDyn origin (0,0,)
246-
SUBROUTINE RigidTrnsf(Init, p, RefPoint, DOF, nDOF, T_ref, ErrStat, ErrMsg)
246+
SUBROUTINE RigidTrnsf(Init, p, RefPoint, DOF, nDOF, nTP, T_ref, ErrStat, ErrMsg)
247247
TYPE(SD_InitType), INTENT(IN ) :: Init ! Input data for initialization routine
248248
TYPE(SD_ParameterType), INTENT(IN ) :: p
249-
REAL(ReKi), INTENT(IN ) :: RefPoint(3) ! Coordinate of the reference point
249+
REAL(ReKi), INTENT(IN ) :: RefPoint(3,*) ! Coordinate of the reference point
250250
INTEGER(IntKi), INTENT(IN ) :: nDOF ! Number of DOFS
251251
INTEGER(IntKi), INTENT(IN ) :: DOF(nDOF) ! DOF indices that are used to create the transformation matrix
252-
REAL(ReKi), INTENT( OUT) :: T_ref(nDOF,6) ! matrix that relates the subset of DOFs to the reference point
252+
INTEGER(IntKi), INTENT(IN ) :: nTP ! Number of transition pieces
253+
REAL(ReKi), INTENT( OUT) :: T_ref(nDOF,6*nTP) ! matrix that relates the subset of DOFs to the reference point
253254
INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation
254255
CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None
255256
! local variables
256-
INTEGER :: I, iDOF, iiDOF, iNode, nDOFPerNode
257+
INTEGER :: I, iDOF, iiDOF, iNode, nDOFPerNode, iTP
257258
REAL(ReKi) :: dx, dy, dz
258259
REAL(ReKi), dimension(6) :: Line
259260
ErrStat = ErrID_None
@@ -264,6 +265,11 @@ SUBROUTINE RigidTrnsf(Init, p, RefPoint, DOF, nDOF, T_ref, ErrStat, ErrMsg)
264265
iNode = p%DOFred2Nodes(iDOF,1) ! First column is node
265266
nDOFPerNode = p%DOFred2Nodes(iDOF,2) ! Second column is number of DOF per node
266267
iiDOF = p%DOFred2Nodes(iDOF,3) ! Third column is dof index for this joint (1-6 for cantilever)
268+
if (nTP==1) then
269+
iTP = 1
270+
else
271+
iTP = p%TPIdx(FINDLOCI(p%Nodes_I(:,1),iNode))
272+
endif
267273

268274
if ((iiDOF<1) .or. (iiDOF>6)) then
269275
ErrMsg = 'RigidTrnsf, node DOF number is not valid. DOF:'//trim(Num2LStr(iDOF))//' Node:'//trim(Num2LStr(iNode))//' iiDOF:'//trim(Num2LStr(iiDOF)); ErrStat = ErrID_Fatal
@@ -273,13 +279,17 @@ SUBROUTINE RigidTrnsf(Init, p, RefPoint, DOF, nDOF, T_ref, ErrStat, ErrMsg)
273279
ErrMsg = 'RigidTrnsf, node doesnt have 6 DOFs. DOF:'//trim(Num2LStr(iDOF))//' Node:'//trim(Num2LStr(iNode))//' nDOF:'//trim(Num2LStr(nDOFPerNode)); ErrStat = ErrID_Fatal
274280
return
275281
endif
282+
if (iTP<1) then
283+
ErrMsg = 'RigidTrnsf, interface node not attached to any transition piece. Node:'//trim(Num2LStr(iNode)); ErrStat = ErrID_Fatal
284+
return
285+
endif
276286

277-
dx = Init%Nodes(iNode, 2) - RefPoint(1)
278-
dy = Init%Nodes(iNode, 3) - RefPoint(2)
279-
dz = Init%Nodes(iNode, 4) - RefPoint(3)
287+
dx = Init%Nodes(iNode, 2) - RefPoint(1,iTP)
288+
dy = Init%Nodes(iNode, 3) - RefPoint(2,iTP)
289+
dz = Init%Nodes(iNode, 4) - RefPoint(3,iTP)
280290

281291
CALL RigidTransformationLine(dx,dy,dz,iiDOF,Line) !returns Line
282-
T_ref(I, 1:6) = Line
292+
T_ref(I, (6*iTP-5):(6*iTP)) = Line
283293
ENDDO
284294
END SUBROUTINE RigidTrnsf
285295

@@ -870,7 +880,7 @@ END SUBROUTINE SetNewPropBR
870880
END SUBROUTINE SD_Discrt
871881

872882

873-
!> Store relative vector between nodes and TP point, to later compute Guyan rigid body motion
883+
!> Store relative vector between nodes and the first TP point, to later compute Guyan rigid body motion
874884
subroutine StoreNodesRelPos(Init, p, ErrStat, ErrMsg)
875885
type(SD_InitType), intent(in ) :: Init
876886
type(SD_ParameterType), intent(inout) :: p
@@ -886,9 +896,9 @@ subroutine StoreNodesRelPos(Init, p, ErrStat, ErrMsg)
886896
call AllocAry(p%DP0, 3, size(Init%Nodes,1), 'DP0', ErrStat2, ErrMsg2); if(Failed()) return
887897

888898
do i = 1, size(Init%Nodes,1)
889-
p%DP0(1, i) = Init%Nodes(i, 2) - Init%TP_RefPoint(1)
890-
p%DP0(2, i) = Init%Nodes(i, 3) - Init%TP_RefPoint(2)
891-
p%DP0(3, i) = Init%Nodes(i, 4) - Init%TP_RefPoint(3)
899+
p%DP0(1, i) = Init%Nodes(i, 2) - Init%TP_RefPoint(1,1)
900+
p%DP0(2, i) = Init%Nodes(i, 3) - Init%TP_RefPoint(2,1)
901+
p%DP0(3, i) = Init%Nodes(i, 4) - Init%TP_RefPoint(3,1)
892902
enddo
893903

894904
contains
@@ -1347,13 +1357,22 @@ SUBROUTINE CheckBCs(p, ErrStat, ErrMsg)
13471357
END SUBROUTINE CheckBCs
13481358

13491359
!> Check interface inputs, and remap 0s and 1s
1350-
SUBROUTINE CheckIntf(p, ErrStat, ErrMsg)
1360+
SUBROUTINE CheckIntf(p, TPIdxInput, RB_RefJoint, ErrStat, ErrMsg)
13511361
TYPE(SD_ParameterType),INTENT(INOUT) :: p
1362+
INTEGER(IntKi), INTENT(INOUT) :: TPIdxInput(:)
1363+
INTEGER(IntKi), INTENT( OUT) :: RB_RefJoint
13521364
INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation
13531365
CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None
1354-
INTEGER(IntKi) :: I, J, iNode
1366+
INTEGER(IntKi) :: I, J, iNode, iFound, TPIDOffset
1367+
LOGICAL, ALLOCATABLE :: TPExist(:)
1368+
LOGICAL :: TP1NodeFound
13551369
ErrMsg = ""
13561370
ErrStat = ErrID_None
1371+
1372+
call AllocAry(TPExist,p%nNodes_I,'TPExist',ErrStat,ErrMsg); if (ErrStat/=ErrID_None) return;
1373+
TPExist = .FALSE.
1374+
p%nTP = -1;
1375+
13571376
DO I = 1, p%nNodes_I
13581377
iNode = p%Nodes_I(I,1) ! Node index
13591378
DO J = 1, 6 ! ItfTDXss ItfTDYss ItfTDZss ItfRDXss ItfRDYss ItfRDZss
@@ -1369,6 +1388,70 @@ SUBROUTINE CheckIntf(p, ErrStat, ErrMsg)
13691388
endif
13701389
ENDDO
13711390
ENDDO
1391+
1392+
TPIDOffset = 0
1393+
DO I = 1, p%nNodes_I
1394+
IF ( TPIdxInput(I) == 0 ) THEN
1395+
TPIDOffset = 1
1396+
EXIT
1397+
END IF
1398+
ENDDO
1399+
TPIdxInput = TPIdxInput + TPIDOffset
1400+
1401+
DO I = 1, p%nNodes_I
1402+
iNode = p%Nodes_I(I,1) ! Node index
1403+
if ( (TPIdxInput(I)<=0) .or. (TPIdxInput(I)>p%nNodes_I) ) then
1404+
ErrStat = ErrID_Fatal
1405+
ErrMsg = 'Transition-piece index must be sequential starting from 0 (floating with multiple transition pieces) or 1 (all other cases), check interface node:'//trim(Num2LStr(iNode))
1406+
return
1407+
else
1408+
if (TPIdxInput(I)>p%nTP) p%nTP = TPIdxInput(I)
1409+
TPExist(TPIdxInput(I)) = .TRUE.
1410+
p%TPIdx(I) = TPIdxInput(I)
1411+
endif
1412+
ENDDO
1413+
1414+
DO I = 1,p%nTP
1415+
if (.not.TPExist(I)) then
1416+
ErrStat = ErrID_Fatal
1417+
ErrMsg = 'Transition-piece index must be sequential starting from 0 (floating with multiple transition pieces) or 1 (all other cases), missing reference to transition piece '//trim(Num2LStr(I-TPIDOffset))
1418+
return
1419+
end if
1420+
ENDDO
1421+
1422+
p%nDOFRB = 0
1423+
if (p%Floating .and. (p%nTP>1)) then
1424+
p%TP1IsRBRefPt = .true.
1425+
p%nDOFRB = 6
1426+
if (TPIDOffset == 0) then
1427+
ErrStat = ErrID_Fatal
1428+
ErrMsg = 'For a floating structure with more than one transition pieces, must have one and only one interface joint assigned to the dummy transition piece (TPID=0) used to represent rigid-body motion. '
1429+
return
1430+
end if
1431+
else
1432+
p%TP1IsRBRefPt = .false.
1433+
end if
1434+
1435+
TP1NodeFound = .false.
1436+
if (p%TP1IsRBRefPt) then
1437+
do I = 1, p%nNodes_I
1438+
if (p%TPIdx(I) == 1) then
1439+
if (.not.TP1NodeFound) then
1440+
TP1NodeFound = .true.
1441+
RB_RefJoint = p%Nodes_I(I,1)
1442+
else
1443+
ErrStat = ErrID_Fatal
1444+
ErrMsg = ' Only one joint can be assigned to the dummy transition piece (TPID=0) serving as the rigid-body reference point for the floating structure. '
1445+
return
1446+
end if
1447+
end if
1448+
end do
1449+
else
1450+
RB_RefJoint = -1
1451+
endif
1452+
1453+
deallocate(TPExist)
1454+
13721455
END SUBROUTINE CheckIntf
13731456

13741457

0 commit comments

Comments
 (0)