diff --git a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 index 1ea74acce7..1dc5f5e45e 100644 --- a/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 +++ b/glue-codes/fast-farm/src/FAST_Farm_Subs.f90 @@ -813,6 +813,9 @@ SUBROUTINE Farm_InitMD( farm, ErrStat, ErrMsg ) ALLOCATE( farm%MD%Input( 2 ), farm%MD%InputTimes( 2 ), STAT = ErrStat2 ) if (Failed0("MD%Input and MD%InputTimes.")) return; + ! Assign the SS pointer of the first SS instance (turbine 1) to MD. Because MD in FF mode will only pull frequency info, instance of SS doesn't matter (error will be thrown by MD if user asks for SS grid). + MD_InitInp%WaveField => farm%FWrap(1)%m%Turbine%SeaSt%p%WaveField ! this is the same wave field as Init%OutData_SeaSt%WaveField in FAST_subs.f90 (as set by line 278 in SeaSt_Init). Cant use Init%OutData_SeaSt%WaveField because Init is a local variable to FAST_InitializeAll + ! initialize MoorDyn CALL MD_Init( MD_InitInp, farm%MD%Input(1), farm%MD%p, farm%MD%x, farm%MD%xd, farm%MD%z, & farm%MD%OtherSt, farm%MD%y, farm%MD%m, farm%p%DT_mooring, MD_InitOut, ErrStat2, ErrMsg2 ) diff --git a/modules/moordyn/CMakeLists.txt b/modules/moordyn/CMakeLists.txt index 7cb0cf82e5..5ec97ef21c 100644 --- a/modules/moordyn/CMakeLists.txt +++ b/modules/moordyn/CMakeLists.txt @@ -28,7 +28,7 @@ add_library(moordynlib STATIC src/MoorDyn_Rod.f90 src/MoorDyn_Types.f90 ) -target_link_libraries(moordynlib nwtclibs) +target_link_libraries(moordynlib seastlib nwtclibs) # Driver add_executable(moordyn_driver @@ -40,7 +40,7 @@ target_link_libraries(moordyn_driver moordynlib versioninfolib) add_library(moordyn_c_binding SHARED src/MoorDyn_C_Binding.f90 ) -target_link_libraries(moordyn_c_binding moordynlib versioninfolib) +target_link_libraries(moordyn_c_binding moordynlib seastlib versioninfolib) if(APPLE OR UNIX) target_compile_definitions(moordyn_c_binding PRIVATE IMPLICIT_DLLEXPORT) endif() diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index e4850c9020..12f381927e 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -204,6 +204,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er DepthValue = "" ! Start off as empty string, to only be filled if MD setting is specified (otherwise InitInp%WtrDepth is used) ! DepthValue and InitInp%WtrDepth are processed later by setupBathymetry. WaterKinValue = "" + + ! Read in the SeaState wave field pointer for wave kinematics (regardless if kinematics are enabled in MD or not) + p%WaveField => InitInp%WaveField m%PtfmInit = InitInp%PtfmInit(:,1) ! is this copying necssary in case this is an individual instance in FAST.Farm? @@ -693,6 +696,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er 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 @@ -715,6 +719,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (N > m%LineTypeList(l)%ElasticMod) then CALL SetErrStat( ErrID_Fatal, 'A line type BA entry cannot have more (bar-separated) values than its EA entry.', ErrStat, ErrMsg, RoutineName ) CALL CleanUp() + RETURN else if (N==2) then ! visco-elastic case when two BA values provided read(tempStrings(2), *) m%LineTypeList(l)%BA_D else if (m%LineTypeList(l)%ElasticMod > 1) then ! case where there is no dynamic damping for viscoelastic model (will it work)? @@ -925,7 +930,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF ( ErrStat2 /= 0 ) THEN CALL WrScr(' Unable to parse Body '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file.') ! Specific screen output because errors likely CALL WrScr(' Ensure row has all 13 columns needed in MDv2 input file (13th Dec 2021).') - CALL SetErrStat( ErrID_Fatal, 'Failed to read bodies.' , ErrStat, ErrMsg, RoutineName ) + CALL SetErrStat( ErrID_Fatal, 'Failed to read bodies.' , ErrStat, ErrMsg, RoutineName ) if (p%writeLog > 0) then write(p%UnLog,'(A)') ' Unable to parse Body '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file.' write(p%UnLog,'(A)') ' Ensure row has all 13 columns needed in MDv2 input file (13th Dec 2021).' @@ -983,11 +988,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Turbine ID out of bounds for Body "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "No number provided for Body "//trim(Num2LStr(l))//" Turbine attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else if (let1 == "FREE") then ! if a free body @@ -1005,6 +1012,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Unidentified Body type string for Body "//trim(Num2LStr(l))//": "//trim(tempString1), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1075,6 +1083,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%RodList(l)%N, LineOutString END IF + ! check p%numRodTypes is greater than 0, if not users are missing Rod Types section of input file + IF (p%nRodTypes == 0) THEN + CALL SetErrStat( ErrID_Fatal, 'No rod types defined in input file. Please define rod types before defining rods.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + END IF + ! find Rod properties index DO J = 1,p%nRodTypes IF (trim(tempString1) == trim(m%RodTypeList(J)%name)) THEN @@ -1083,6 +1098,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END IF IF (J == p%nRodTypes) THEN ! call an error if there is no match CALL SetErrStat( ErrID_Fatal, 'Unable to find matching rod type name for Rod '//trim(Num2LStr(l))//": "//trim(tempString1), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() RETURN END IF END DO @@ -1136,17 +1152,20 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Unidentified Type/BodyID for Rod "//trim(Num2LStr(l))//": "//trim(tempString2), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "Body ID out of bounds for Rod "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "No number provided for Rod "//trim(Num2LStr(l))//" Body attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else if ((let1 == "VESSEL") .or. (let1 == "VES") .or. (let1 == "COUPLED") .or. (let1 == "CPLD")) then ! if a rigidly coupled rod, add to list and add @@ -1185,11 +1204,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Turbine ID out of bounds for Rod "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "No number provided for Rod "//trim(Num2LStr(l))//" Turbine attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else if ((let1 == "ROD") .or. (let1 == "R") .or. (let1 == "FREE")) then @@ -1206,6 +1227,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Unidentified Type/BodyID for Rod "//trim(Num2LStr(l))//": "//trim(tempString2), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1321,11 +1343,11 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF ( ErrStat2 /= 0 ) THEN CALL WrScr(' Unable to parse Point '//trim(Num2LStr(l))//' row in input file.') ! Specific screen output because errors likely CALL WrScr(' Ensure row has all 9 columns, including CdA and Ca.') ! to be caused by non-updated input file formats. - CALL SetErrStat( ErrID_Fatal, 'Failed to read points.' , ErrStat, ErrMsg, RoutineName ) ! would be nice to specify which line <<<<<<<<< - if (p%writeLog > 0) then - write(p%UnLog,'(A)') ' Unable to parse Point '//trim(Num2LStr(l))//' row in input file.' - write(p%UnLog,'(A)') ' Ensure row has all 9 columns, including CdA and Ca.' - end if + CALL SetErrStat( ErrID_Fatal, 'Failed to read points.' , ErrStat, ErrMsg, RoutineName ) ! would be nice to specify which line <<<<<<<<< + if (p%writeLog > 0) then + write(p%UnLog,'(A)') ' Unable to parse Point '//trim(Num2LStr(l))//' row in input file.' + write(p%UnLog,'(A)') ' Ensure row has all 9 columns, including CdA and Ca.' + end if CALL CleanUp() RETURN END IF @@ -1354,11 +1376,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Body ID out of bounds for Point "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "No number provided for Point "//trim(Num2LStr(l))//" Body attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else if ((let1 == "VESSEL") .or. (let1 == "VES") .or. (let1 == "COUPLED") .or. (let1 == "CPLD")) then ! if a fairlead, add to list and add @@ -1397,15 +1421,18 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er else CALL SetErrStat( ErrID_Fatal, "Turbine ID out of bounds for Point "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "No number provided for Point "//trim(Num2LStr(l))//" Turbine attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if else CALL SetErrStat( ErrID_Fatal, "Unidentified Type/BodyID for Point "//trim(Num2LStr(l))//": "//trim(tempString1), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1482,6 +1509,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END IF IF (J == p%nLineTypes) THEN ! call an error if there is no match CALL SetErrStat( ErrID_Fatal, 'Unable to find matching line type name for Line '//trim(Num2LStr(l)), ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() RETURN END IF END DO @@ -1510,6 +1538,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (len_trim(num1)<1) then CALL SetErrStat( ErrID_Fatal, "Error: no number provided for line "//trim(Num2LStr(l))//" end A attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1525,10 +1554,12 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL Rod_AddLine(m%RodList(J), l, 0, 1) ! add line l (end A, denoted by 0) to rod J (end B, denoted by 1) else CALL SetErrStat( ErrID_Fatal, "Error: rod end (A or B) must be specified for line "//trim(Num2LStr(l))//" end A attachment. Instead seeing "//let2, ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else CALL SetErrStat( ErrID_Fatal, " Rod ID out of bounds for line "//trim(Num2LStr(l))//" end A attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1539,6 +1570,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL Point_AddLine(m%PointList(J), l, 0) ! add line l (end A, denoted by 0) to point J else CALL SetErrStat( ErrID_Fatal, "Error: point out of bounds for line "//trim(Num2LStr(l))//" end A attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1551,6 +1583,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (len_trim(num1)<1) then CALL SetErrStat( ErrID_Fatal, "Error: no number provided for line "//trim(Num2LStr(l))//" end B attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1566,10 +1599,12 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL Rod_AddLine(m%RodList(J), l, 1, 1) ! add line l (end B, denoted by 1) to rod J (end B, denoted by 1) else CALL SetErrStat( ErrID_Fatal, "Error: rod end (A or B) must be specified for line "//trim(Num2LStr(l))//" end B attachment. Instead seeing "//let2, ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return end if else CALL SetErrStat( ErrID_Fatal, " Rod ID out of bounds for line "//trim(Num2LStr(l))//" end B attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1580,6 +1615,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL Point_AddLine(m%PointList(J), l, 1) ! add line l (end B, denoted by 1) to point J else CALL SetErrStat( ErrID_Fatal, "Error: point out of bounds for line "//trim(Num2LStr(l))//" end B attachment.", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return end if @@ -1793,11 +1829,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er END IF ELSE CALL SetErrStat( ErrID_Fatal, "Body ID out of bounds for External Load "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return END IF ELSE CALL SetErrStat( ErrID_Fatal, "No number provided for External Load "//trim(Num2LStr(l))//" BODY attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return END IF ELSEIF (let1 == "POINT" .OR. let1 == "P") THEN IF (len_trim(num1) > 0) THEN @@ -1808,11 +1846,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%PointList(J)%Bquad = m%PointList(J)%Bquad + m%ExtLdList(l)%Bquad ELSE CALL SetErrStat( ErrID_Fatal, "Point ID out of bounds for External Load "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return END IF ELSE CALL SetErrStat( ErrID_Fatal, "No number provided for External Load "//trim(Num2LStr(l))//" POINT attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return END IF ELSEIF (let1 == "ROD" .OR. let1 == "R") THEN IF (len_trim(num1) > 0) THEN @@ -1823,11 +1863,13 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er m%RodList(J)%Bquad = m%RodList(J)%Bquad(1:2) + m%ExtLdList(l)%Bquad(1:2) ! rods only have axial and transverse ELSE CALL SetErrStat( ErrID_Fatal, "Rod ID out of bounds for External Load "//trim(Num2LStr(l))//".", ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() return END IF ELSE CALL SetErrStat( ErrID_Fatal, "No number provided for External Load "//trim(Num2LStr(l))//" ROD attachment.", ErrStat, ErrMsg, RoutineName ) - return + CALL CleanUp() + return END IF END IF @@ -2693,8 +2735,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er end if ! boost drag coefficient of each line type <<<<<<<< does this actually do anything or do lines hold these coefficients??? + m%IC_gen = .True. ! turn on IC_gen flag DO I = 1, p%nLines - m%LineList(I)%IC_gen = .True. ! turn on IC_gen flag for Line VIV model m%LineList(I)%Cdn = m%LineList(I)%Cdn * InputFileDat%CdScaleIC m%LineList(I)%Cdt = m%LineList(I)%Cdt * InputFileDat%CdScaleIC END DO @@ -2852,9 +2894,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL MD_DestroyInput( u_array(1), ErrStat2, ErrMsg2 ) - ! UNboost drag coefficient of each line type <<< + ! Unboost drag coefficient of each line type <<< + m%IC_gen = .False. ! turn off IC_gen flag DO I = 1, p%nLines - m%LineList(I)%IC_gen = .False. ! turn off IC_gen flag for Line VIV model m%LineList(I)%Cdn = m%LineList(I)%Cdn / InputFileDat%CdScaleIC m%LineList(I)%Cdt = m%LineList(I)%Cdt / InputFileDat%CdScaleIC END DO @@ -2949,7 +2991,7 @@ SUBROUTINE CheckError(ErrID,Msg) CHARACTER(*), INTENT(INOUT) :: Msg ! The error message (ErrMsg) INTEGER(IntKi) :: ErrStat3 ! The error identifier (ErrStat) - CHARACTER(1024) :: ErrMsg3 ! The error message (ErrMsg) + CHARACTER(ErrMsgLen) :: ErrMsg3 ! The error message (ErrMsg) ! Set error status/message; IF ( ErrID /= ErrID_None ) THEN @@ -3745,7 +3787,7 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er ! give Lines latest state variable values for internal nodes DO l = 1,p%nLines - CALL Line_SetState(m%LineList(l), x%states(m%LineStateIs1(l):m%LineStateIsN(l)), t) + CALL Line_SetState(m%LineList(l), x%states(m%LineStateIs1(l):m%LineStateIsN(l)), t, m) END DO ! calculate dynamics of free objects (will also calculate forces (doRHS()) from any child/dependent objects)... diff --git a/modules/moordyn/src/MoorDyn_Body.f90 b/modules/moordyn/src/MoorDyn_Body.f90 index 263a3156f3..3fbf2d6bac 100644 --- a/modules/moordyn/src/MoorDyn_Body.f90 +++ b/modules/moordyn/src/MoorDyn_Body.f90 @@ -82,7 +82,7 @@ SUBROUTINE Body_Setup( Body, tempArray, p, ErrStat, ErrMsg) Body%BquadL = 0.0_DbKi !also set number of attached rods and points to zero initially - Body%nAttachedC = 0 + Body%nAttachedP = 0 Body%nAttachedR = 0 ! set up body initial mass matrix (excluding any rods or attachements) @@ -327,7 +327,7 @@ SUBROUTINE Body_SetDependentKin(Body, t, m) Body%OrMat = TRANSPOSE( EulerConstruct( Body%r6(4:6) ) ) ! full Euler angle approach <<<< need to check order ! set kinematics of any dependent points - do l = 1,Body%nAttachedC + do l = 1,Body%nAttachedP CALL transformKinematics(Body%rPointRel(:,l), Body%r6, Body%OrMat, Body%v6, rPoint, rdPoint) !<<< should double check this function @@ -522,7 +522,7 @@ SUBROUTINE Body_DoRHS(Body, m, p) ! Get contributions from any dependent points - do l = 1,Body%nAttachedC + do l = 1,Body%nAttachedP ! get net force and mass from Point on body ref point (global orientation) CALL Point_GetNetForceAndMass( m%PointList(Body%attachedC(l)), Body%r6(1:3), Body%v6(4:6), F6_i, M6_i, m, p) @@ -622,10 +622,10 @@ SUBROUTINE Body_AddPoint(Body, pointID, coords) IF (wordy > 0) Print*, "P", pointID, "->B", Body%IdNum - IF(Body%nAttachedC < 30) THEN ! this is currently just a maximum imposed by a fixed array size. could be improved. - Body%nAttachedC = Body%nAttachedC + 1 ! increment the number pointed - Body%AttachedC(Body%nAttachedC) = pointID - Body%rPointRel(:,Body%nAttachedC) = coords ! store relative position of point on body + IF(Body%nAttachedP < 30) THEN ! this is currently just a maximum imposed by a fixed array size. could be improved. + Body%nAttachedP = Body%nAttachedP + 1 ! increment the number pointed + Body%AttachedC(Body%nAttachedP) = pointID + Body%rPointRel(:,Body%nAttachedP) = coords ! store relative position of point on body ELSE call WrScr("too many Points attached to Body "//trim(num2lstr(Body%IdNum))//" in MoorDyn!") END IF diff --git a/modules/moordyn/src/MoorDyn_C_Binding.f90 b/modules/moordyn/src/MoorDyn_C_Binding.f90 index baf3742772..445348389d 100644 --- a/modules/moordyn/src/MoorDyn_C_Binding.f90 +++ b/modules/moordyn/src/MoorDyn_C_Binding.f90 @@ -144,7 +144,7 @@ SUBROUTINE MD_C_Init( & REAL(C_FLOAT) , INTENT(IN ) :: G_C REAL(C_FLOAT) , INTENT(IN ) :: RHO_C REAL(C_FLOAT) , INTENT(IN ) :: DEPTH_C - REAL(C_FLOAT) , INTENT(IN ) :: PtfmInit_C(6) + REAL(C_FLOAT) , INTENT(IN ) :: PtfmInit_C(6) ! TODO: make this more flexible, can we not have 6 DOF only coupling? INTEGER(C_INT) , INTENT(IN ) :: InterpOrder_C INTEGER(C_INT) , INTENT( OUT) :: NumChannels_C CHARACTER(KIND=C_CHAR) , INTENT( OUT) :: OutputChannelNames_C(100000) @@ -729,4 +729,4 @@ SUBROUTINE Set_OutputLoadArray() tmpForces(4:6,1) = MD_LoadMesh%Moment(1:3,1) END SUBROUTINE Set_OutputLoadArray -END MODULE +END MODULE \ No newline at end of file diff --git a/modules/moordyn/src/MoorDyn_Driver.f90 b/modules/moordyn/src/MoorDyn_Driver.f90 index 1e7e5950f7..6e2d4326c5 100644 --- a/modules/moordyn/src/MoorDyn_Driver.f90 +++ b/modules/moordyn/src/MoorDyn_Driver.f90 @@ -22,13 +22,14 @@ PROGRAM MoorDyn_Driver USE MoorDyn_Types USE MoorDyn + USE SeaState_Types + USE SeaState USE NWTC_Library USE VersionInfo IMPLICIT NONE TYPE MD_Drvr_InitInput - ! LOGICAL :: Echo REAL(DbKi) :: Gravity REAL(DbKi) :: rhoW REAL(DbKi) :: WtrDepth @@ -37,6 +38,8 @@ PROGRAM MoorDyn_Driver CHARACTER(1024) :: OutRootName REAL(DbKi) :: TMax REAL(DbKi) :: dtC + + CHARACTER(1024) :: SeaStateInputFile INTEGER :: FarmSize REAL(DbKi) :: FarmPositions(8,40) @@ -48,10 +51,10 @@ PROGRAM MoorDyn_Driver INTEGER(IntKi) :: ErrStat ! Status of error message - CHARACTER(1024) :: ErrMsg ! Error message if ErrStat /= ErrID_None + CHARACTER(ErrMsgLen) :: ErrMsg ! Error message if ErrStat /= ErrID_None INTEGER(IntKi) :: ErrStat2 ! Status of error message - CHARACTER(1024) :: ErrMsg2 ! Error message if ErrStat /= ErrID_None + CHARACTER(ErrMsgLen) :: ErrMsg2 ! Error message if ErrStat /= ErrID_None CHARACTER(1024) :: drvrFilename ! Filename and path for the driver input file. This is passed in as a command line argument when running the Driver exe. TYPE(MD_Drvr_InitInput) :: drvrInitInp ! Initialization data for the driver program @@ -72,6 +75,19 @@ PROGRAM MoorDyn_Driver TYPE (MD_OutputType) :: MD_y ! Output file identifier + ! SeaState types + TYPE(SeaSt_InitInputType) :: InitInData_SeaSt ! Input data for initialization + TYPE(SeaSt_InitOutputType) :: InitOutData_SeaSt ! Output data from initialization + type(SeaSt_ContinuousStateType) :: x_SeaSt ! Continuous states + type(SeaSt_DiscreteStateType) :: xd_SeaSt ! Discrete states + type(SeaSt_ConstraintStateType) :: z_SeaSt ! Constraint states + type(SeaSt_OtherStateType) :: OtherState_SeaSt ! Other states + type(SeaSt_MiscVarType) :: m_SeaSt ! Misc/optimization variables + type(SeaSt_ParameterType) :: p_SeaSt ! Parameters + type(SeaSt_InputType) :: u_SeaSt(1) ! System inputs + type(SeaSt_OutputType) :: y_SeaSt ! System outputs + LOGICAL :: SeaState_Initialized = .FALSE. + ! Motion file parsing type(FileInfoType) :: FileInfo_PrescribeMtn !< The derived type for holding the prescribed forces input file for parsing -- we may pass this in the future integer(IntKi) :: CurLine !< current entry in FileInfo_In%Lines array @@ -153,7 +169,6 @@ PROGRAM MoorDyn_Driver ! do any initializing and allocating needed in prep for calling MD_Init ! set the input file name and other environment terms - !MD_InitInp%NStepWave = 1 ! an arbitrary number > 0 (to set the size of the wave data, which currently contains all zero values) MD_InitInp%Tmax = drvrInitInp%TMax MD_InitInp%g = drvrInitInp%Gravity MD_InitInp%rhoW = drvrInitInp%rhoW @@ -161,9 +176,6 @@ PROGRAM MoorDyn_Driver MD_InitInp%FileName = drvrInitInp%MDInputFile MD_InitInp%RootName = drvrInitInp%OutRootName MD_InitInp%UsePrimaryInputFile = .TRUE. - !MD_InitInp%PassedPrimaryInputData = - ! MD_InitInp%Echo = drvrInitInp%Echo - !MD_InitInp%OutList = <<<< never used? MD_InitInp%Linearize = .FALSE. TMax = drvrInitInp%TMax @@ -211,29 +223,37 @@ PROGRAM MoorDyn_Driver ! -------------------------------- ----------------------------------- + + IF (LEN_TRIM(drvrInitInp%SeaStateInputFile) > 0 ) THEN ! If SeaState input file path in driver input file is not empty. Error checks for Null pointer in MD_Init -> setupWaterKin + ! Initialize the SeaState module + InitInData_SeaSt%hasIce = .FALSE. + InitInData_SeaSt%Gravity = MD_InitInp%g + InitInData_SeaSt%defWtrDens = MD_InitInp%rhoW + InitInData_SeaSt%defWtrDpth = MD_InitInp%WtrDepth + InitInData_SeaSt%defMSL2SWL = 0.0_DbKi ! MoorDyn does not allow for a sea level offset + InitInData_SeaSt%UseInputFile = .TRUE. + InitInData_SeaSt%InputFile = drvrInitInp%SeaStateInputFile + InitInData_SeaSt%OutRootName = trim(MD_InitInp%RootName)//'.SEA' + InitInData_SeaSt%TMax = MD_InitInp%TMax + InitInData_SeaSt%Linearize = MD_InitInp%Linearize + + CALL SeaSt_Init( InitInData_SeaSt, u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, dtC, InitOutData_SeaSt, ErrStat2, ErrMsg2 ); call AbortIfFailed() + SeaState_Initialized = .TRUE. + + IF ( dtC /= drvrInitInp%dtC) THEN + ErrMsg2 = 'The SeaState Module attempted to change the coupling timestep, but this is not allowed. The SeaState Module must use the Driver coupling timestep.' + ErrStat2 = ErrID_Fatal + CALL AbortIfFailed() + ENDIF - ! fill in the hydrodynamics data - !ALLOCATE( MD_InitInp%WaveVel (2,200,3)) - !ALLOCATE( MD_InitInp%WaveAcc (2,200,3)) - !ALLOCATE( MD_InitInp%WavePDyn(2,200) ) - !ALLOCATE( MD_InitInp%WaveElev(2,200) ) - !ALLOCATE( MD_InitInp%WaveTime(2) ) - !MD_InitInp%WaveVel = 0.0_ReKi - !MD_InitInp%WaveAcc = 0.0_ReKi - !MD_InitInp%WavePDyn = 0.0_ReKi - !MD_InitInp%WaveElev = 0.0_ReKi - !MD_InitInp%WaveTime = 0.0_ReKi - !DO I = 1,SIZE(MD_InitInp%WaveTime) - ! MD_InitInp%WaveTime(I) = 600.0*I - !END DO - - ! open driver output file >>> not yet used <<< - !CALL GetNewUnit( Un ) - !OPEN(Unit=Un,FILE='MD.out',STATUS='UNKNOWN') + ! pass the pointer + MD_InitInp%WaveField => InitOutData_SeaSt%WaveField + + END IF ! call the initialization routine CALL MD_Init( MD_InitInp, MD_u(1), MD_p, MD_x , MD_xd, MD_xc, MD_xo, MD_y, MD_m, dtC, MD_InitOut, ErrStat2, ErrMsg2 ); call AbortIfFailed() - + CALL MD_DestroyInitInput ( MD_InitInp , ErrStat2, ErrMsg2 ); call AbortIfFailed() CALL MD_DestroyInitOutput ( MD_InitOut , ErrStat2, ErrMsg2 ); call AbortIfFailed() @@ -673,6 +693,9 @@ PROGRAM MoorDyn_Driver CALL RunTimes( ProgStrtTime, ProgStrtCPU, SimStrtTime, SimStrtCPU, t ) ! Destroy all objects + IF (SeaState_Initialized) THEN + CALL SeaSt_End( u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, ErrStat2, ErrMsg2); call AbortIfFailed() + ENDIF CALL MD_End( MD_u(1), MD_p, MD_x, MD_xd, MD_xc , MD_xo, MD_y, MD_m, ErrStat2, ErrMsg2 ); call AbortIfFailed() do j = 2,MD_interp_order+1 @@ -683,7 +706,7 @@ PROGRAM MoorDyn_Driver CALL WrScr1( "Errors: " ) CALL WrScr( trim(GetErrStr(ErrStat))//': '//trim(ErrMsg) ) endif - + !close (un) call CleanUp() CALL NormStop() @@ -696,6 +719,16 @@ SUBROUTINE AbortIfFailed() call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver') if (ErrStat >= AbortErrLev) then + if (SeaState_Initialized) then + call SeaSt_End( u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, ErrStat2, ErrMsg2) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' ) + end if + + CALL SeaSt_DestroyInitOutput( InitOutData_SeaSt, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' ) + CALL SeaSt_DestroyInitInput( InitInData_SeaSt, ErrStat2, ErrMsg2 ) + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' ) + call CleanUp() Call ProgAbort(trim(ErrMsg)) elseif ( ErrStat2 /= ErrID_None ) THEN @@ -724,8 +757,9 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp) TYPE(MD_Drvr_InitInput), INTENT( OUT ) :: InitInp ! Local variables INTEGER :: J ! generic integer for counting + INTEGER :: i ! generic integer for counting - ! CHARACTER(1024) :: EchoFile ! Name of MoorDyn echo file + CHARACTER(1024) :: tmpString ! temporary string CHARACTER(1024) :: FileName ! Name of MoorDyn input file CHARACTER(1024) :: FilePath ! Name of path to MoorDyn input file @@ -739,23 +773,20 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp) call AbortIfFailed() CALL WrScr( 'Opening MoorDyn Driver input file: '//FileName ) - - ! Read until "echo" - CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 1', ErrStat2, ErrMsg2); call AbortIfFailed() - CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 2', ErrStat2, ErrMsg2); call AbortIfFailed() - ! CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo Input', ErrStat2, ErrMsg2); call AbortIfFailed() - ! ! If we echo, we rewind - ! IF ( InitInp%Echo ) THEN - ! EchoFile = TRIM(FileName)//'.echo' - ! CALL GetNewUnit( UnEcho ) - ! CALL OpenEcho ( UnEcho, EchoFile, ErrStat2, ErrMsg2 ); call AbortIfFailed() - ! REWIND(UnIn) - ! CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 1', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - ! CALL ReadCom( UnIn, FileName, 'MoorDyn Driver input file header line 2', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - ! CALL ReadVar ( UnIn, FileName, InitInp%Echo, 'Echo', 'Echo the input file data', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - ! END IF + + ! Read through the header lines until hitting --- + DO I=1,30 ! max of 30 header lines + CALL ReadVar( UnIn, FileName, tmpString, '', 'MoorDyn Driver input file header line', ErrStat2, ErrMsg2); call AbortIfFailed() + IF (INDEX(tmpString, '---') > 0) EXIT ! exit the loop if we hit the end of the header + ENDDO + ! make sure the user didn't give more than 30 lines of header text + IF (I == 30) THEN + ErrStat2 = ErrID_Fatal + ErrMsg2 = ' The MoorDyn Driver input file can have a maximum of 30 header lines.' + CALL AbortIfFailed() + END IF !---------------------- ENVIRONMENTAL CONDITIONS ------------------------------------------------- - CALL ReadCom( UnIn, FileName, 'Environmental conditions header', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() + ! The Environmental conditions header is read at the end of the above loop. CALL ReadVar( UnIn, FileName, InitInp%Gravity, 'Gravity', 'Gravity', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() CALL ReadVar( UnIn, FileName, InitInp%rhoW , 'rhoW', 'water density', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() CALL ReadVar( UnIn, FileName, InitInp%WtrDepth, 'WtrDepth', 'water depth', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() @@ -768,7 +799,17 @@ SUBROUTINE ReadDriverInputFile( inputFile, InitInp) CALL ReadVar( UnIn, FileName, InitInp%InputsMod , 'InputsMode', 'Mode for the inputs - zero/steady/time-series', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() CALL ReadVar( UnIn, FileName, InitInp%InputsFile , 'InputsFile', 'Filename for the MoorDyn inputs', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() CALL ReadVar( UnIn, FileName, InitInp%FarmSize , 'NumTurbines', 'number of turbines in FAST.Farm', ErrStat2, ErrMsg2, UnEcho); call AbortIfFailed() - CALL ReadCom( UnIn, FileName, 'Initial positions header', ErrStat2, ErrMsg2); call AbortIfFailed() + CALL ReadVar( UnIn, FileName, tmpString , 'SeaStateFile', 'Filename for the SeaState inputs', ErrStat2, ErrMsg2, UnEcho) + ! Check if SeaState path is given. If not provided, then initialize it as an empty string. This keeps things backwards compatible + IF (INDEX(tmpString, '---') > 0) THEN + InitInp%SeaStateInputFile = '' + CALL WrScr('No SeaState input file specified in the MoorDyn driver. SeaState will not be initialized.') + ELSE + InitInp%SeaStateInputFile = tmpString + CALL ReadCom( UnIn, FileName, 'Initial positions header', ErrStat2, ErrMsg2); call AbortIfFailed() ! skip the inital positions header if SeaState path exists (need to read an extra line) + END IF + !---------------------- Initial Positions -------------------------------------------------------- + ! The Initial Positions conditions header is read by the above SeaState path logic CALL ReadCom( UnIn, FileName, 'Initial positions table header line 1', ErrStat2, ErrMsg2); call AbortIfFailed() CALL ReadCom( UnIn, FileName, 'Initial positions table header line 2', ErrStat2, ErrMsg2); call AbortIfFailed() do J=1,MAX(1,InitInp%FarmSize) diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index 773cf03566..12764b381a 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -1305,11 +1305,6 @@ SUBROUTINE MDIO_CloseOutput ( p, m, ErrStat, ErrMsg ) IF (ALLOCATED(m%MDWrOutput)) THEN DEALLOCATE(m%MDWrOutput) ENDIF - DO I=1,p%NLines - IF (ALLOCATED(m%LineList(I)%LineWrOutput)) THEN - DEALLOCATE(m%LineList(I)%LineWrOutput) ! this may be unnecessary and handled by Line destructor - ENDIF - END DO END SUBROUTINE MDIO_CloseOutput !----------------------------------------------------------------------------------------============ diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index f40d482a75..345a29702a 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -1033,11 +1033,12 @@ END SUBROUTINE Line_Initialize !-------------------------------------------------------------- - SUBROUTINE Line_SetState(Line, X, t) + SUBROUTINE Line_SetState(Line, X, t, m) TYPE(MD_Line), INTENT(INOUT) :: Line ! the current Line object Real(DbKi), INTENT(IN ) :: X(:) ! state vector section for this line Real(DbKi), INTENT(IN ) :: t ! instantaneous time + TYPE(MD_MiscVarType), INTENT(INOUT) :: m ! passing along all mooring objects INTEGER(IntKi) :: i ! index of segments or nodes along line INTEGER(IntKi) :: J ! index @@ -1063,7 +1064,7 @@ SUBROUTINE Line_SetState(Line, X, t) end if ! if using the viv mdodel, also set the lift force phase - if (Line%Cl > 0 .AND. (.NOT. Line%IC_gen) .AND. t > 0) then ! not needed in IC_gen, and t=0 should be skipped to avoid setting these all to zero. Initialize as distribution on 0-2pi + if (Line%Cl > 0 .AND. (.NOT. m%IC_gen) .AND. t > 0) then ! not needed in IC_gen, and t=0 should be skipped to avoid setting these all to zero. Initialize as distribution on 0-2pi do I=0, Line%N if (Line%ElasticMod > 1) then ! if both additional states are included then N-1 entries after internal node states and visco segment states Line%phi(I) = X( 7*Line%N-6 + I+1) - (2 * Pi * floor(X( 7*Line%N-6 + I+1) / (2*Pi))) ! Map integrated phase to 0-2Pi range. Is this necessary? sin (a-b) is the same if b is 100 pi or 2pi @@ -1167,6 +1168,12 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Real(DbKi) :: phi_dot ! frequency of lift force (rad/s) Real(DbKi) :: f_hat ! non-dimensional frequency + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + CHARACTER(120) :: RoutineName = 'Line_GetStateDeriv' + + ErrStat = ErrID_None + ErrMsg = "" N = Line%N ! for convenience d = Line%d @@ -1215,7 +1222,8 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! apply wave kinematics (if there are any) DO i=0,N - CALL getWaterKin(p, Line%r(1,i), Line%r(2,i), Line%r(3,i), Line%time, m%WaveTi, Line%U(:,i), Line%Ud(:,i), Line%zeta(i), Line%PDyn(i)) + CALL getWaterKin(p, m, Line%r(1,i), Line%r(2,i), Line%r(3,i), Line%time, Line%U(:,i), Line%Ud(:,i), Line%zeta(i), Line%PDyn(i), ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) END DO ! --------- calculate line partial submergence (Line::calcSubSeg from MD-C) --------- @@ -1554,7 +1562,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Vortex Induced Vibration (VIV) cross-flow lift force Line%Lf(:,I) = 0.0_DbKi ! Zero lift force - IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen + IF ((Line%Cl > 0.0) .AND. (.NOT. m%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen ! Note: This logic is slightly different than MD-C, but equivalent. MD-C runs the VIV model for only the internal ! nodes. That means in MD-F the state vector has N+1 extra states when using the VIV model while the MD-C state @@ -1746,7 +1754,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai endif ! ! for checking rdd_old - ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen) then + ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. m%IC_gen) then ! print*, "rdd_old at t = ", Line%time ! DO I = 0, 4 ! print*, "I =", I, "rdd_old =", Line%rdd_old(:,I) diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index bf26a7ab1b..2a3bd59bf6 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -21,6 +21,8 @@ MODULE MoorDyn_Misc USE MoorDyn_Types + USE SeaSt_WaveField + USE Current USE NWTC_Library USE NWTC_FFTPACK @@ -75,7 +77,7 @@ SUBROUTINE UnitVector( r1, r2, u, Length ) ! note: order of parameters u = r2 - r1 length = TwoNorm(u) - if ( .NOT. EqualRealNos(length, 0.0_DbKi ) ) THEN + IF ( .NOT. EqualRealNos(length, 0.0_DbKi ) ) THEN u = u / Length END IF @@ -98,11 +100,11 @@ SUBROUTINE ScaleVector( u_in, newlength, u_out ) length_squared = length_squared + u_in(J)*u_in(J) END DO - if (length_squared > 0) then + IF (length_squared > 0) THEN scaler = newlength/sqrt(length_squared) - else ! if original vector is zero, return zero + ELSE ! if original vector is zero, return zero scaler = 0.0_DbKi - end if + END IF DO J=1,3 u_out(J) = u_in(J) * scaler @@ -115,63 +117,63 @@ END SUBROUTINE ScaleVector ! convenience function to calculate curvature based on adjacent segments' direction vectors and their combined length function GetCurvature(length, q1, q2) - real(DbKi), intent(in ) :: length - real(DbKi), intent(in ) :: q1(3) - real(DbKi), intent(in ) :: q2(3) - real(DbKi) :: GetCurvature + Real(DbKi), intent(in ) :: length + Real(DbKi), intent(in ) :: q1(3) + Real(DbKi), intent(in ) :: q2(3) + Real(DbKi) :: GetCurvature - real(DbKi) :: q1_dot_q2 + Real(DbKi) :: q1_dot_q2 ! note "length" here is combined from both segments q1_dot_q2 = dot_product( q1, q2 ) - if (q1_dot_q2 > 1.0) then ! this is just a small numerical error, so set q1_dot_q2 to 1 + IF (q1_dot_q2 > 1.0) THEN ! this is just a small numerical error, so set q1_dot_q2 to 1 GetCurvature = 0.0_DbKi ! this occurs when there's no curvature, so return zero curvature - !else if (q1_dot_q2 < 0) ! this is a bend of more than 90 degrees, too much, call an error! + !ELSE IF (q1_dot_q2 < 0) ! this is a bend of more than 90 degrees, too much, call an error! - else ! normal case + ELSE ! normal case GetCurvature = 4.0/length * sqrt(0.5*(1.0 - q1_dot_q2)) ! this is the normal curvature calculation - end if + END IF - return - end function GetCurvature + RETURN + END function GetCurvature ! calculate orientation angles of a direction vector !----------------------------------------------------------------------- subroutine GetOrientationAngles(vec, phi, sinPhi, cosPhi, tanPhi, beta, sinBeta, cosBeta, k_hat) - real(DbKi), intent(in ) :: vec(3) !p1(3),p2(3) - real(DbKi), intent( out) :: phi, sinPhi, cosPhi, tanPhi, beta, sinBeta, cosBeta, k_hat(3) + Real(DbKi), intent(in ) :: vec(3) !p1(3),p2(3) + Real(DbKi), intent( out) :: phi, sinPhi, cosPhi, tanPhi, beta, sinBeta, cosBeta, k_hat(3) - real(DbKi) :: vecLen, vecLen2D + Real(DbKi) :: vecLen, vecLen2D vecLen = SQRT(Dot_Product(vec,vec)) vecLen2D = SQRT(vec(1)**2+vec(2)**2) - if ( vecLen < 0.000001 ) then - if (wordy > 0) then + IF ( vecLen < 0.000001 ) THEN + IF (wordy > 0) THEN print *, "ERROR in GetOrientationAngles in MoorDyn. Supplied vector is near zero" print *, vec - endif + ENDIF k_hat = NaN ! 1.0/0.0 - else + ELSE k_hat = vec / vecLen phi = atan2(vecLen2D, vec(3)) ! incline angle - end if - if ( phi < 0.000001) then + END IF + IF ( phi < 0.000001) THEN beta = 0.0_ReKi - else + ELSE beta = atan2(vec(2), vec(1)) ! heading of incline - endif + ENDIF sinPhi = sin(phi) cosPhi = cos(phi) tanPhi = tan(phi) sinBeta = sin(beta) cosBeta = cos(beta) - end subroutine GetOrientationAngles + END subroutine GetOrientationAngles !----------------------------------------------------------------------- @@ -541,11 +543,11 @@ SUBROUTINE LUsolve(n, A, LU, b, y, x) sum = sum + LU(k,p)*LU(p,j) END DO - if ( EqualRealNos( LU(k,k), 0.0_DbKi) ) then - LU(k,j) = 0.0_DbKi ! avoid divide by zero <<< numerator likely zero too. check if this is safe <<< - else + IF ( EqualRealNos( LU(k,k), 0.0_DbKi) ) THEN + LU(k,j) = 0.0_DbKi ! avoid divide by zero <<< numerator likely zero too. check IF this is safe <<< + ELSE LU(k,j) = (A(k,j)-sum)/LU(k,k) - end if + END IF END DO !j END DO !K @@ -558,11 +560,11 @@ SUBROUTINE LUsolve(n, A, LU, b, y, x) sum = sum + LU(i,k)*y(k); END DO - if ( EqualRealNos( LU(i,i), 0.0_DbKi) ) then - y(i) = 0.0_DbKi ! avoid divide by zero <<< numerator likely zero too. check if this is safe <<< - else + IF ( EqualRealNos( LU(i,i), 0.0_DbKi) ) THEN + y(i) = 0.0_DbKi ! avoid divide by zero <<< numerator likely zero too. check IF this is safe <<< + ELSE y(i) = (b(i)-sum)/LU(i,i) - end if + END IF END DO DO j=1,n ! this is actually for looping through i in reverse @@ -600,15 +602,15 @@ SUBROUTINE getInterpNums(xlist, xin, istart, i, fout) nx = SIZE(xlist) - if (xin <= xlist(1)) THEN ! below lowest data point + IF (xin <= xlist(1)) THEN ! below lowest data point i = 1_IntKi fout = 0.0_DbKi - else if (xlist(nx) <= xin) THEN ! above highest data point + ELSE IF (xlist(nx) <= xin) THEN ! above highest data point i = nx fout = 0.0_DbKi - else ! within the data range + ELSE ! within the data range IF (xlist(min(istart,nx)) < xin) i1 = istart ! if istart is below the actual value, start with it instead of starting at 1 to save time, but make sure it doesn't overstep the array @@ -638,15 +640,15 @@ SUBROUTINE getInterpNumsSiKi(xlist, xin, istart, i, fout) nx = SIZE(xlist) - if (xin <= xlist(1)) THEN ! below lowest data point + IF (xin <= xlist(1)) THEN ! below lowest data point i = 1_IntKi fout = 0.0_SiKi - else if (xlist(nx) <= xin) THEN ! above highest data point + ELSE IF (xlist(nx) <= xin) THEN ! above highest data point i = nx fout = 0.0_SiKi - else ! within the data range + ELSE ! within the data range IF (xlist(min(istart,nx)) < xin) i1 = istart ! if istart is below the actual value, start with it instead of starting at 1 to save time, but make sure it doesn't overstep the array @@ -672,29 +674,29 @@ SUBROUTINE calculate4Dinterpolation(f, ix0, iy0, iz0, it0, fx, fy, fz, ft, c) REAL(SiKi) :: c00, c01, c10, c11, c0, c1 ! handle end case conditions - if (fx == 0) then + IF (fx == 0) THEN ix1 = ix0 - else + ELSE ix1 = min(ix0+1,size(f,4)) ! don't overstep bounds - end if + END IF - if (fy == 0) then + IF (fy == 0) THEN iy1 = iy0 - else + ELSE iy1 = min(iy0+1,size(f,3)) ! don't overstep bounds - end if + END IF - if (fz == 0) then + IF (fz == 0) THEN iz1 = iz0 - else + ELSE iz1 = min(iz0+1,size(f,2)) ! don't overstep bounds - end if + END IF - if (ft == 0) then + IF (ft == 0) THEN it1 = it0 - else + ELSE it1 = min(it0+1,size(f,1)) ! don't overstep bounds - end if + END IF c000 = f(it0,iz0,iy0,ix0)*(1.0-ft) + f(it1,iz0,iy0,ix0)*ft c001 = f(it0,iz1,iy0,ix0)*(1.0-ft) + f(it1,iz1,iy0,ix0)*ft @@ -732,23 +734,23 @@ SUBROUTINE calculate3Dinterpolation(f, ix0, iy0, iz0, fx, fy, fz, c) ! note that "z" could also be "t" - dimension names are arbitrary ! handle end case conditions - if (fx == 0) then + IF (fx == 0) THEN ix1 = ix0 - else + ELSE ix1 = min(ix0+1,size(f,3)) ! don't overstep bounds - end if + END IF - if (fy == 0) then + IF (fy == 0) THEN iy1 = iy0 - else + ELSE iy1 = min(iy0+1,size(f,2)) ! don't overstep bounds - end if + END IF - if (fz == 0) then + IF (fz == 0) THEN iz1 = iz0 - else + ELSE iz1 = min(iz0+1,size(f,1)) ! don't overstep bounds - end if + END IF c000 = f(iz0,iy0,ix0) c001 = f(iz1,iy0,ix0) @@ -880,16 +882,16 @@ SUBROUTINE getDepthFromBathymetry(BathymetryGrid, BathGrid_Xs, BathGrid_Ys, Line ! get local slope dx = BathGrid_Xs(ix1) - BathGrid_Xs(ix0) dy = BathGrid_Ys(iy1) - BathGrid_Ys(iy0) - if ( dx > 0.0 ) then + IF ( dx > 0.0 ) THEN dc_dx = (c1y-c0y)/dx - else + ELSE dc_dx = 0.0_DbKi ! maybe this should raise an error - end if - if ( dy > 0.0 ) then + END IF + IF ( dy > 0.0 ) THEN dc_dy = (cx1-cx0)/dy - else + ELSE dc_dy = 0.0_DbKi ! maybe this should raise an error - end if + END IF tempVector(1) = dc_dx tempVector(2) = dc_dy @@ -902,377 +904,135 @@ END SUBROUTINE getDepthFromBathymetry ! :::::::::::::::::::::::::: wave and current subroutines ::::::::::::::::::::::::::::::: - ! master function to get wave/water kinematics at a given point -- called by each object from grid-based data - SUBROUTINE getWaterKin(p, x, y, z, t, tindex, U, Ud, zeta, PDyn) + ! master function to get wave/water kinematics at a given point -- called by each object + SUBROUTINE getWaterKin(p, m, x, y, z, t, U, Ud, zeta, PDyn, ErrStat, ErrMsg) ! This whole approach assuems that px, py, and pz are in increasing order. - ! Wheeler stretching is now built in. + ! Wheeler stretching is built in for WaterKin 1 and 2. - TYPE(MD_ParameterType),INTENT (IN ) :: p ! MoorDyn parameters (contains the wave info for now) - Real(DbKi), INTENT (IN ) :: x - Real(DbKi), INTENT (IN ) :: y - Real(DbKi), INTENT (IN ) :: z - Real(DbKi), INTENT (IN ) :: t - INTEGER(IntKi), INTENT (INOUT) :: tindex ! pass time index to try starting from, returns identified time index - Real(DbKi), INTENT (INOUT) :: U(3) - Real(DbKi), INTENT (INOUT) :: Ud(3) - Real(DbKi), INTENT (INOUT) :: zeta - Real(DbKi), INTENT (INOUT) :: PDyn - - + TYPE(MD_ParameterType), INTENT (IN ) :: p ! MoorDyn parameters + TYPE(MD_MiscVarType), INTENT (INOUT) :: m ! MoorDyn misc data + Real(DbKi), INTENT (IN ) :: x ! node position + Real(DbKi), INTENT (IN ) :: y ! node position + Real(DbKi), INTENT (IN ) :: z ! node position + Real(DbKi), INTENT (IN ) :: t ! time + Real(DbKi), INTENT (INOUT) :: U(3) ! fluid speed + Real(DbKi), INTENT (INOUT) :: Ud(3) ! fluid acceleration + Real(DbKi), INTENT (INOUT) :: zeta ! water height above node + Real(DbKi), INTENT (INOUT) :: PDyn ! dynamic fluid pressure + INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message IF ErrStat /= ErrID_None + + ! outputs of WaveField_GetNodeWaveKin not in above list (single precision required by wave grid) + Real(ReKi) :: xyz_sp(3) ! single precision + Real(SiKi) :: WaveElev1 ! Unused by MD. per SeaSt_WaveField: zeta = WaveElev = WaveElev1 + WaveElev2 + Real(SiKi) :: WaveElev2 ! Unused by MD. per SeaSt_WaveField: zeta = WaveElev = WaveElev1 + WaveElev2 + Real(SiKi) :: zeta_sp ! single precison + Real(SiKi) :: U_sp(3) ! single precision + Real(SiKi) :: Ud_sp(3) ! single precision + Real(SiKi) :: FAMCF(3) ! Unused by MD + Real(SiKi) :: PDyn_sp ! single precision + INTEGER(IntKi) :: nodeInWater ! Unused by MD + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + + ! interpolation variables for Waterkin = 1, 2 INTEGER(IntKi) :: ix, iy, iz, it ! indices for interpolation INTEGER(IntKi) :: iz0, iz1 ! special indices for currrent interpolation -! INTEGER(IntKi) :: N ! number of rod elements for convenience Real(SiKi) :: fx, fy, fz, ft ! interpolation fractions Real(DbKi) :: zp ! zprime coordinate used for Wheeler stretching + + ! local + CHARACTER(120) :: RoutineName = 'getWaterKin' - - ! if wave kinematics enabled, get interpolated values from grid - if (p%WaveKin > 0) then - - ! find time interpolation indices and coefficients - !CALL getInterpNums(p%tWave, t, tindex, it, ft) - it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 - ft = (t - (it-1)*p%dtWave)/p%dtWave - tindex = it - - ! find x-y interpolation indices and coefficients - CALL getInterpNumsSiKi(p%pxWave , REAL(x,SiKi), 1, ix, fx) - CALL getInterpNumsSiKi(p%pyWave , REAL(y,SiKi), 1, iy, fy) + ErrStat = ErrID_None + ErrMsg = "" + + IF (p%WaterKin == 3 .AND. (.NOT. m%IC_gen)) THEN ! disable wavekin 3 during IC_gen, otherwise will never find stead state (becasue of waves) - ! interpolate wave elevation - CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) + ! SeaState throws warning when queried location is out of bounds from the SeaState grid, so no need to handle here + + ! Pack all MD inputs to WaveGrid input data types (double to single) + ! (only pos needed becasue time is double in wave field, all other are outputs that will be set by WaveField_GetNodeWaveKin) + xyz_sp = REAL((/ x, y, z /),SiKi) + + ! for now we will force the node to be in the water (forceNodeInWater = True). Rods handle partial submergence seperately so they need to get information from SeaState + CALL WaveField_GetNodeWaveKin(p%WaveField, m%WaveField_m, t, xyz_sp, .TRUE., nodeInWater, WaveElev1, WaveElev2, zeta_sp, PDyn_sp, U_sp, Ud_sp, FAMCF, ErrStat2, ErrMsg2 ) ! outputs: nodeInWater, WaveElev1, WaveElev2, FAMCF all unused + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + + ! Unpack all WaveGrid outputs to MD output data types (single to double) + U = REAL(U_sp,DbKi) + Ud = REAL(Ud_sp,DbKi) + zeta = REAL(zeta_sp,DbKi) + PDyn = REAL(PDyn_sp,DbKi) + + ELSEIF (p%WaterKin == 1 .OR. p%WaterKin == 2) THEN ! old or hybrid approach. SeaState contributions handeled in setupWaterKin, just proceed using old method + + ! If wave kinematics enabled, get interpolated values from grid + IF (p%WaveKin > 0) THEN + + ! find time interpolation indices and coefficients + !CALL getInterpNums(p%tWave, t, tindex, it, ft) + it = floor(t/ p%dtWave) + 1 ! add 1 because Fortran indexing starts at 1 + ft = (t - (it-1)*p%dtWave)/p%dtWave + m%WaveTi = it + + ! find x-y interpolation indices and coefficients + CALL getInterpNumsSiKi(p%pxWave , REAL(x,SiKi), 1, ix, fx) ! wave grid + CALL getInterpNumsSiKi(p%pyWave , REAL(y,SiKi), 1, iy, fy) ! wave grid + + ! interpolate wave elevation + CALL calculate3Dinterpolation(p%zeta, ix, iy, it, fx, fy, ft, zeta) + + ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching + zp = ( z - zeta ) * p%WtrDpth/( p%WtrDpth + zeta ) + + CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) ! wave grid + + ! interpolate everything else + CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) + CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1) ) + CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2) ) + CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3) ) + CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1) ) + CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2) ) + CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3) ) - ! compute modified z coordinate to be used for interpolating velocities and accelerations with Wheeler stretching - zp = ( z - zeta ) * p%WtrDpth/( p%WtrDpth + zeta ) + ELSE ! set things to zero if wave kinematics not enabled + U = 0.0_DbKi + Ud = 0.0_DbKi + zeta = 0.0_DbKi + PDyn = 0.0_DbKi + + ENDIF + + ! IF current kinematics enabled, add interpolated current values from profile + IF (p%Current > 0) THEN - CALL getInterpNumsSiKi(p%pzWave , REAL(zp,SiKi), 1, iz, fz) - - ! interpolate everything else - CALL calculate4Dinterpolation(p%PDyn , ix, iy, iz, it, fx, fy, fz, ft, PDyn) - CALL calculate4Dinterpolation(p%uxWave, ix, iy, iz, it, fx, fy, fz, ft, U(1) ) - CALL calculate4Dinterpolation(p%uyWave, ix, iy, iz, it, fx, fy, fz, ft, U(2) ) - CALL calculate4Dinterpolation(p%uzWave, ix, iy, iz, it, fx, fy, fz, ft, U(3) ) - CALL calculate4Dinterpolation(p%axWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(1) ) - CALL calculate4Dinterpolation(p%ayWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(2) ) - CALL calculate4Dinterpolation(p%azWave, ix, iy, iz, it, fx, fy, fz, ft, Ud(3) ) - else + CALL getInterpNumsSiKi(p%pzCurrent, REAL(z,SiKi), 1, iz0, fz) + + IF (fz == 0) THEN ! handle end case conditions + iz1 = iz0 + ELSE + iz1 = min(iz0+1,size(p%pzCurrent)) ! don't overstep bounds + END IF + + ! Add the current velocities to the wave velocities (if any) + U(1) = U(1) + (1.0-fz)*p%uxCurrent(iz0) + fz*p%uxCurrent(iz1) + U(2) = U(2) + (1.0-fz)*p%uyCurrent(iz0) + fz*p%uyCurrent(iz1) + END IF + + ELSEIF (p%WaterKin > 3) THEN + CALL SetErrStat(ErrID_Fatal, "Invalid p%WaterKin value found in getWaterKin", ErrStat, ErrMsg, RoutineName) + + ELSE ! set things to zero if Water Kinematics not enabled U = 0.0_DbKi Ud = 0.0_DbKi zeta = 0.0_DbKi PDyn = 0.0_DbKi - end if - - - ! if current kinematics enabled, add interpolated current values from profile - if (p%Current > 0) then - - CALL getInterpNumsSiKi(p%pzCurrent, REAL(z,SiKi), 1, iz0, fz) - - IF (fz == 0) THEN ! handle end case conditions - iz1 = iz0 - ELSE - iz1 = min(iz0+1,size(p%pzCurrent)) ! don't overstep bounds - END IF - - U(1) = U(1) + (1.0-fz)*p%uxCurrent(iz0) + fz*p%uxCurrent(iz1) - U(2) = U(2) + (1.0-fz)*p%uyCurrent(iz0) + fz*p%uyCurrent(iz1) - end if - - END SUBROUTINE getWaterKin - - - !! ! unused routine with old code for taking wave kinematic grid inputs from HydroDyn - !! SUBROUTINE CopyWaterKinFromHydroDyn(p, InitInp) - !! - !! TYPE(MD_InitInputType), INTENT(IN ) :: InitInp ! INTENT(INOUT) : Input data for initialization routine - !! TYPE(MD_ParameterType), INTENT( OUT) :: p ! INTENT( OUT) : Parameters - !! - !! INTEGER(IntKi) :: I, J, K, Itemp - !! - !! - !! ! ----------------------------- Arrays for wave kinematics ----------------------------- - !! - !! - !!! :::::::::::::: BELOW WILL BE USED EVENTUALLY WHEN WAVE INFO IS AN INPUT :::::::::::::::::: - !!! ! The rAll array contains all nodes or reference points in the system - !!! ! (x,y,z global coordinates for each) in the order of bodies, rods, points, internal line nodes. - !!! - !!! ! count the number of nodes to use for passing wave kinematics - !!! J=0 - !!! ! Body reference point coordinates - !!! J = J + p%nBodies - !!! ! Rod node coordinates (including ends) - !!! DO l = 1, p%nRods - !!! J = J + (m%RodList(l)%N + 1) - !!! END DO - !!! ! Point reference point coordinates - !!! J = J + p%nConnects - !!! ! Line internal node coordinates - !!! DO l = 1, p%nLines - !!! J = J + (m%LineList(l)%N - 1) - !!! END DO - !!! - !!! ! allocate all relevant arrays - !!! ! allocate state vector and temporary state vectors based on size just calculated - !!! ALLOCATE ( y%rAll(3,J), u%U(3,J), u%Ud(3,J), u%zeta(J), u%PDyn(J), STAT = ErrStat ) - !!! IF ( ErrStat /= ErrID_None ) THEN - !!! ErrMsg = ' Error allocating wave kinematics vectors.' - !!! RETURN - !!! END IF - !!! - !!! - !!! ! go through the nodes and fill in the data (this should maybe be turned into a global function) - !!! J=0 - !!! ! Body reference point coordinates - !!! DO I = 1, p%nBodies - !!! J = J + 1 - !!! y%rAll(:,J) = m%BodyList(I)%r6(1:3) - !!! END DO - !!! ! Rod node coordinates - !!! DO I = 1, p%nRods - !!! DO K = 0,m%RodList(I)%N - !!! J = J + 1 - !!! y%rAll(:,J) = m%RodList(I)%r(:,K) - !!! END DO - !!! END DO - !!! ! Point reference point coordinates - !!! DO I = 1, p%nConnects - !!! J = J + 1 - !!! y%rAll(:,J) = m%ConnectList(I)%r - !!! END DO - !!! ! Line internal node coordinates - !!! DO I = 1, p%nLines - !!! DO K = 1,m%LineList(I)%N-1 - !!! J = J + 1 - !!! y%rAll(:,J) = m%LineList(I)%r(:,K) - !!! END DO - !!! END DO - !! ! :::::::::::::::: the above might be used eventually. For now, let's store wave info grids within this module ::::::::::::::::: - !! - !! - !! ! ----- copy wave grid data over from HydroDyn (as was done in USFLOWT branch) ----- - !! - !! ! get grid and time info (currently this is hard-coded to match what's in HydroDyn_Input - !! ! DO I=1,p%nzWave - !! ! p%pz(I) = 1.0 - 2.0**(p%nzWave-I) ! -127, -63, -31, -15, -7, -3, -1, 0 - !! ! END DO - !! ! DO J = 1,p%nyWave - !! ! p%py(J) = WaveGrid_y0 + WaveGrid_dy*(J-1) - !! ! END DO - !! ! DO K = 1,p%nxWave - !! ! p%px(K) = WaveGrid_x0 + WaveGrid_dx*(K-1) - !! ! END DO - !! ! - !! ! p%tWave = InitInp%WaveTime - !! - !! DO I=1,p%nzWave - !! DO J = 1,p%nyWave - !! DO K = 1,p%nxWave - !! Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node on 3D grid - !! - !! p%uxWave (:,I,J,K) = InitInp%WaveVel( :,Itemp,1) ! note: indices are t, z, y, x - !! p%uyWave (:,I,J,K) = InitInp%WaveVel( :,Itemp,2) - !! p%uzWave (:,I,J,K) = InitInp%WaveVel( :,Itemp,3) - !! p%axWave (:,I,J,K) = InitInp%WaveAcc( :,Itemp,1) - !! p%ayWave (:,I,J,K) = InitInp%WaveAcc( :,Itemp,2) - !! p%azWave (:,I,J,K) = InitInp%WaveAcc( :,Itemp,3) - !! p%PDyn( :,I,J,K) = InitInp%WavePDyn(:,Itemp) - !! END DO - !! END DO - !! END DO - !! - !! DO J = 1,p%nyWave - !! DO K = 1,p%nxWave - !! Itemp = (J-1)*p%nxWave + K ! index of actual node on surface 2D grid - !! p%zeta(:,J,K) = InitInp%WaveElev(:,Itemp) - !! END DO - !! END DO - !! - !! END SUBROUTINE CopyWaterKinFromHydroDyn - - - ! ----- write wave grid spacing to output file ----- - SUBROUTINE WriteWaveGrid(p, ErrStat, ErrMsg) - - TYPE(MD_ParameterType), INTENT(INOUT) :: p ! Parameters - - INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None - - INTEGER(IntKi) :: ErrStat2 -! CHARACTER(120) :: ErrMsg2 - - CHARACTER(120) :: Frmt - INTEGER(IntKi) :: UnOut ! for outputing wave kinematics data - INTEGER(IntKi) :: I - - - CALL GetNewUnit( UnOut) - - CALL OpenFOutFile ( UnOut, "waves.txt", ErrStat, ErrMsg ) - IF ( ErrStat > ErrID_None ) THEN - ErrMsg = ' Error opening wave grid file: '//TRIM(ErrMsg) - ErrStat = ErrID_Fatal - RETURN - END IF - - WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( 'MoorDyn v2 wave/current kinematics grid file' ) - WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( '---------------------------------------------' ) - WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( 'The following 6 lines (4-9) specify the input type then the inputs for x, then, y, then z coordinates.' ) - - WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - X input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) - Frmt = '('//TRIM(Int2LStr(5))//'(A1,e10.4))' - WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pxWave(I))), I=1,p%nxWave ) - - WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - Y input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) - Frmt = '('//TRIM(Int2LStr(5))//'(A1,e10.4))' - WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pyWave(I))), I=1,p%nyWave ) - - WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - Z input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) - Frmt = '('//TRIM(Int2LStr(8))//'(A1,e10.4))' - WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pzWave(I))), I=1,p%nzWave ) - - CLOSE(UnOut, IOSTAT = ErrStat2 ) - IF ( ErrStat2 /= 0 ) THEN - ErrStat = ErrID_Severe - ErrMsg = 'Error closing wave grid file' - END IF - - END SUBROUTINE WriteWaveGrid - - - ! ----- write wave kinematics grid data to output file ----- - SUBROUTINE WriteWaveData(p, ErrStat, ErrMsg) - - TYPE(MD_ParameterType), INTENT(INOUT) :: p ! Parameters - - INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None - - INTEGER(IntKi) :: ErrStat2 -! CHARACTER(120) :: ErrMsg2 - - INTEGER(IntKi) :: UnOut ! for outputing wave kinematics data - INTEGER(IntKi) :: I,J,K, l, Itemp - - CALL GetNewUnit( UnOut) - - CALL OpenFOutFile ( UnOut, "wave data.txt", ErrStat, ErrMsg ) - IF ( ErrStat > ErrID_None ) THEN - ErrMsg = ' Error opening wave grid file: '//TRIM(ErrMsg) - ErrStat = ErrID_Fatal - RETURN - END IF - - ! write channel labels - - - ! time - WRITE(UnOut,"(A10)", IOSTAT=ErrStat2, advance="no") "Time" - - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ze0", Num2Lstr(J+10*K) - END DO - END DO - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ux", Num2Lstr(I+10*J+100*K) - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " uy", Num2Lstr(I+10*J+100*K) - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " uz", Num2Lstr(I+10*J+100*K) - END DO - END DO - END DO - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ax", Num2Lstr(I+10*J+100*K) - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ay", Num2Lstr(I+10*J+100*K) - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " az", Num2Lstr(I+10*J+100*K) - END DO - END DO - END DO - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " pd", Num2Lstr(I+10*J+100*K) - END DO - END DO - END DO - - ! end the line - WRITE(UnOut, "(A1)", IOSTAT=ErrStat2) " " - - - - DO l=1, p%ntWave ! loop through all time steps - - ! time - WRITE(UnOut,"(F10.4)", IOSTAT=ErrStat2, advance="no") p%dtWave*(l-1) - !WRITE(UnOut,"(F10.4)", IOSTAT=ErrStat2, advance="no") InitInp%WaveTime(l) - - ! wave elevation (all slices for now, to check) - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - Itemp = (J-1)*p%nxWave + K ! index of actual node - - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%zeta(l,J,K) - END DO - END DO - - ! wave velocities - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node - - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uxWave(l,I,J,K) - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uyWave(l,I,J,K) - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uzWave(l,I,J,K) - END DO - END DO - END DO - - ! wave accelerations - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node - - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%axWave(l,I,J,K) - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%ayWave(l,I,J,K) - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%azWave(l,I,J,K) - END DO - END DO - END DO - - ! dynamic pressure - DO I=1,p%nzWave !z - DO J = 1,p%nyWave !y - DO K = 1,p%nxWave !x - Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node - - WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%PDyn(l,I,J,K) - END DO - END DO - END DO - - ! end the line - WRITE(UnOut, "(A1)", IOSTAT=ErrStat2) " " - - - END DO - - - CLOSE(UnOut, IOSTAT = ErrStat ) - IF ( ErrStat /= 0 ) THEN - ErrMsg = 'Error closing wave grid file' - END IF - - END SUBROUTINE WriteWaveData + ENDIF + END SUBROUTINE getWaterKin ! ----- process WaterKin input value, potentially reading wave inputs and generating wave field ----- SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) @@ -1287,10 +1047,12 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) INTEGER(IntKi) :: ntIn ! number of time series inputs from file INTEGER(IntKi) :: UnIn ! unit number for coefficient input file INTEGER(IntKi) :: UnEcho - REAL(SiKi) :: pzCurrentTemp(100) ! current depth increments read in from input file (positive-down at this stage) + REAL(SiKi), ALLOCATABLE :: pzCurrentTemp(:) ! current depth increments read in from input file (positive-down at this stage) REAL(SiKi) :: uxCurrentTemp(100) REAL(SiKi) :: uyCurrentTemp(100) + TYPE(Current_InitOutputType) :: Current_InitOutput ! Current outputs from SS subroutine Current_Init (for current mod = 2) + CHARACTER(120) :: tmpString CHARACTER(120) :: WaveKinFile INTEGER(IntKi) :: UnElev ! unit number for coefficient input file REAL(SiKi), ALLOCATABLE :: WaveTimeIn(:) ! temporarily holds wave input time series @@ -1307,7 +1069,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) INTEGER(IntKi) :: NStepWave ! INTEGER(IntKi) :: NStepWave2 ! REAL(SiKi) :: WaveTMax ! max wave elevation time series duration after optimizing lenght for FFT - REAL(SiKi) :: WaveDOmega + REAL(SiKi) :: WaveDOmega ! frequency step REAL(SiKi) :: SinWaveDir ! SIN( WaveDirArr(I) ) -- Each wave frequency has a unique wave direction. REAL(SiKi) :: CosWaveDir ! COS( WaveDirArr(I) ) -- Each wave frequency has a unique wave direction. @@ -1327,384 +1089,690 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) COMPLEX(SiKi), ALLOCATABLE :: WaveAccCHy(:) ! Discrete Fourier transform of the instantaneous horizontal acceleration in y-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2) COMPLEX(SiKi), ALLOCATABLE :: WaveAccCV( :) ! Discrete Fourier transform of the instantaneous vertical acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2) COMPLEX(SiKi), ALLOCATABLE :: WaveDynPC( :) ! Discrete Fourier transform of the instantaneous dynamic pressure of incident waves before applying stretching at the zi-coordinates for points (N/m^2) - COMPLEX(SiKi), ALLOCATABLE :: WaveVelCHx(:) ! Discrete Fourier transform of the instantaneous horizontal velocity of incident waves before applying stretching at the zi-coordinates for points (m/s) - COMPLEX(SiKi), ALLOCATABLE :: WaveVelCHy(:) ! Discrete Fourier transform of the instantaneous horizontal velocity in x-direction of incident waves before applying stretching at the zi-coordinates for points (m/s) - COMPLEX(SiKi), ALLOCATABLE :: WaveVelCV( :) ! Discrete Fourier transform of the instantaneous vertical velocity in y-direction of incident waves before applying stretching at the zi-coordinates for points (m/s) + COMPLEX(SiKi), ALLOCATABLE :: WaveVelCHx(:) ! Discrete Fourier transform of the instantaneous horizontal velocity in x-direction of incident waves before applying stretching at the zi-coordinates for points (m/s) + COMPLEX(SiKi), ALLOCATABLE :: WaveVelCHy(:) ! Discrete Fourier transform of the instantaneous horizontal velocity in y-direction of incident waves before applying stretching at the zi-coordinates for points (m/s) + COMPLEX(SiKi), ALLOCATABLE :: WaveVelCV( :) ! Discrete Fourier transform of the instantaneous vertical velocity of incident waves before applying stretching at the zi-coordinates for points (m/s) ! COMPLEX(SiKi) :: WGNC ! Discrete Fourier transform of the realization of a White Gaussian Noise (WGN) time series process with unit variance for the current frequency component (-) - INTEGER(IntKi) :: ErrStatTmp INTEGER(IntKi) :: ErrStat2 - CHARACTER(120) :: ErrMsg2 + CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(120) :: RoutineName = 'SetupWaveKin' - - ErrStatTmp = ErrID_None ! TODO: get rid of redundancy <<< ErrStat2 = ErrID_None ErrMsg2 = "" IF (LEN_TRIM(WaterKinString) == 0) THEN ! If the input is empty (not provided), there are no water kinematics to be included - p%WaveKin = 0 - p%Current = 0 - return + p%WaveKin = 0 + p%Current = 0 + p%WaterKin = 0 + RETURN ELSE IF (SCAN(WaterKinString, "abcdfghijklmnopqrstuvwxyzABCDFGHIJKLMNOPQRSTUVWXYZ") == 0) THEN ! If the input has no letters, let's assume it's a number - call WrScr( "ERROR WaveKin option does not currently support numeric entries. It must be a filename." ) - p%WaveKin = 0 - p%Current = 0 - return + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR WaveKin option does not currently support numeric entries. It must be a filename." + ENDIF + p%WaveKin = 0 + p%Current = 0 + p%WaterKin = 0 + CALL SetErrStat( ErrID_Fatal, "ERROR WaveKin option does not currently support numeric entries. It must be a filename.", ErrStat, ErrMsg, RoutineName); RETURN END IF + tmpString = WaterKinString + CALL Conv2UC(tmpString) - ! otherwise interpret the input as a file name to load the bathymetry lookup data from - call WrScr( " The waterKin input contains letters so will load a water kinematics input file" ) - - - ! -------- load water kinematics input file ------------- - - IF ( PathIsRelative( WaterKinString ) ) THEN ! properly handle relative path <<< - !CALL GetPath( TRIM(InitInp%InputFile), TmpPath ) - FileName = TRIM(p%PriPath)//TRIM(WaterKinString) - ELSE - FileName = trim(WaterKinString) - END IF - - - - UnEcho=-1 - CALL GetNewUnit( UnIn ) - CALL OpenFInpFile( UnIn, FileName, ErrStat2, ErrMsg2); if(Failed()) return - - - CALL ReadCom( UnIn, FileName, 'MoorDyn water kinematics input file header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadCom( UnIn, FileName, 'MoorDyn water kinematics input file header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - ! ----- waves ----- - CALL ReadCom( UnIn, FileName, 'waves header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadVar( UnIn, FileName, p%WaveKin , 'WaveKinMod' , 'WaveKinMod' , ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadVar( UnIn, FileName, WaveKinFile, 'WaveKinFile', 'WaveKinFile' , ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadVar( UnIn, FileName, p%dtWave , 'dtWave', 'time step for waves', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadVar( UnIn, FileName, WaveDir , 'WaveDir' , 'wave direction', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - ! X grid points - READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type - READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed - CALL gridAxisCoords(coordtype, entries2, p%pxWave, p%nxWave, ErrStat2, ErrMsg2) - Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') - ! Y grid points - READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type - READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed - CALL gridAxisCoords(coordtype, entries2, p%pyWave, p%nyWave, ErrStat2, ErrMsg2) - Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') - ! Z grid points - READ(UnIn,*, IOSTAT=ErrStat2) coordtype ! get the entry type - READ(UnIn,'(A)', IOSTAT=ErrStat2) entries2 ! get entries as string to be processed - CALL gridAxisCoords(coordtype, entries2, p%pzWave, p%nzWave, ErrStat2, ErrMsg2) - Call SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, 'MD_getWaterKin') - ! ----- current ----- - CALL ReadCom( UnIn, FileName, 'current header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadVar( UnIn, FileName, p%Current, 'CurrentMod', 'CurrentMod', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadCom( UnIn, FileName, 'current profile header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - CALL ReadCom( UnIn, FileName, 'current profile header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return - ! current profile table... (read through to end of file or ---) - DO I=1,100 - READ(UnIn, *, IOSTAT=ErrStat2) pzCurrentTemp(i), uxCurrentTemp(i), uyCurrentTemp(i) ! read into a line - if (ErrStat2 /= 0) then - p%nzCurrent = i-1 ! save number of valid current depth points in profile - EXIT ! break out of the loop if it couldn't read the line (i.e. if at end of file) - end if - if (i == 100) then - call WrScr("WARNING: MD can handle a maximum of 100 current profile points") - exit - end if - END DO - + ! Error check for SeaState rhoW and MoorDyn rhoW, SeaState gravity and MoorDyn gravity - CLOSE(UnIn) + IF (tmpString == "SEASTATE") THEN ! if its the SeaState Keyword, then we are using the SeaState method (WaterKin = 3) - - ! ------------------- start with wave kinematics ----------------------- - - ! WaveKin options: 0 - none or set externally during the sim (Waves object not needed unless there's current) [default] - ! 1 - set externally for each node in each object (Waves object not needed unless there's current) (TBD) - ! 2 - set from inputted wave elevation FFT, grid approach* (TBD) - ! 3 - set from inputted wave elevation time series, grid approach* [supported] - ! 4 - set from inputted wave elevation FFT, node approach (TBD) - ! 5 - set from inputted wave elevation time series, node approach (TBD) - ! 6 - set from inputted velocity, acceleration, and wave elevation grid data (TBD)** - - ! Current options: 0 - no currents or set externally (as part of WaveKin =0 or 1 approach) [default] - ! 1 - read in steady current profile, grid approach (current_profile.txt)** [supported] - ! 2 - read in dynamic current profile, grid approach (current_profile_dynamic.txt)** (TBD) - ! 3 - read in steady current profile, node approach (current_profile.txt) (TBD) - ! 4 - read in dynamic current profile, node approach (current_profile_dynamic.txt) (TBD) - - ! * the first call to any of these will attempt to load water_grid.txt to define the grid to put things on - ! ** if a grid has already been set, these will interpolate onto it, otherwise they'll make a new grid based on their provided coordinates + ! Error check to make sure we are not running in FAST.Farm + IF (p%nTurbines > 1) THEN ! if running in FAST.Farm, cannot use the SeaState coupling method because of multiple SeaState instances and only a single MoorDyn instance + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " MoorDyn fully coupled with SeaState is not compatible with FAST.Farm. Please use the WaterKin 1 or 2 methods." + ENDIF + CALL SetErrStat(ErrID_Fatal, "MoorDyn fully coupled with SeaState is not compatible with FAST.Farm. Please use the WaterKin 1 or 2 methods.", ErrStat, ErrMsg, RoutineName); RETURN + END IF - ! NOTE: lots of partial code is available from MD-C for supporting various wave kinematics input options - - ! WaveKin and Current compatibility check could go here in future - - - ! --------------------- set from inputted wave elevation time series, grid approach ------------------- - if (p%WaveKin == 3) then + ! Error check to make sure the wave field pointer is not null + IF (.NOT. ASSOCIATED(p%WaveField)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR WaveField pointer is null. SeaState method requires SeaState to be enabled. Please check input files." + ENDIF + CALL SetErrStat(ErrID_Fatal, "WaveField pointer is null. SeaState method requires SeaState to be enabled. Please check input files.", ErrStat, ErrMsg, RoutineName); RETURN + END IF - call WrScr( 'Setting up WaveKin 3 option: read wave elevation time series from file' ) + ! Warning check to make sure SeaState and MoorDyn have the same water depth + IF (p%WaveField%WtrDpth /= p%WtrDpth) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING SeaState water depth does not match MoorDyn water depth. Using SeaState values for water kinematics." + ENDIF + CALL SetErrStat(ErrID_Warn, "SeaState water depth does not match MoorDyn water depth. Using SeaState values for water kinematics.", ErrStat, ErrMsg, RoutineName) + END IF - IF ( LEN_TRIM( WaveKinFile ) == 0 ) THEN - CALL SetErrStat( ErrID_Fatal,'WaveKinFile must not be an empty string.',ErrStat, ErrMsg, RoutineName); return - RETURN + ! Error check to make sure SeaState and MoorDyn have the same water density and gravity. This also checks that the pointer is valid + IF (p%WaveField%RhoXg /= REAL((p%rhoW * p%g), SiKi)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR SeaState (water density * gravity) ["//trim(num2lstr(p%WaveField%RhoXg))//"] does not match MoorDyn water (density * gravity) ["//trim(num2lstr(REAL((p%rhoW * p%g), SiKi)))//"]. The SeaState pointer may be corrupted." + ENDIF + CALL SetErrStat(ErrID_Fatal, " ERROR SeaState (water density * gravity) ["//trim(num2lstr(p%WaveField%RhoXg))//"] does not match MoorDyn water (density * gravity) ["//trim(num2lstr(REAL((p%rhoW * p%g), SiKi)))//"]. The SeaState pointer may be corrupted.", ErrStat, ErrMsg, RoutineName) END IF - IF ( PathIsRelative( WaveKinFile ) ) THEN ! properly handle relative path <<< + ! Check for if SeaState grid does not match water depth + IF (p%WaveField%GridParams%Z_Depth /= p%WtrDpth) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " INFO SeaState grid depth does not match MoorDyn water depth." + ENDIF + END IF + + ! turn on flags, for getWaterKin to know to just pull data from the WaveField pointer. Nothing more is needed for setup + p%WaveKin = 3 + p%Current = 3 + p%WaterKin = 3 + + ELSE ! we must be using the old or hybrid methods. Start setting up the user grids by reading the WaveKin File + + ! otherwise interpret the input as a file name to load the bathymetry lookup data from + CALL WrScr( " The waterKin input contains a filename so will load a water kinematics input file" ) + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " The waterKin input contains a filename so will load a water kinematics input file" + ENDIF + + + ! -------- load water kinematics input file ------------- + + IF ( PathIsRelative( WaterKinString ) ) THEN ! properly handle relative path <<< !CALL GetPath( TRIM(InitInp%InputFile), TmpPath ) - WaveKinFile = TRIM(p%PriPath)//TRIM(WaveKinFile) + FileName = TRIM(p%PriPath)//TRIM(WaterKinString) + ELSE + FileName = trim(WaterKinString) END IF - ! note: following is adapted from MoorDyn_Driver - CALL GetNewUnit( UnElev ) - - CALL OpenFInpFile ( UnElev, WaveKinFile, ErrStat2, ErrMsg2 ); if(Failed()) return - - call WrScr( 'Reading wave elevation data from '//trim(WaveKinFile) ) - ! Read through length of file to find its length - i = 0 ! start line counter - numHdrLn = 0 ! start header-line counter - dataBegin = .FALSE. ! started reading the data section - DO - READ(UnElev,'(A)',IOSTAT=ErrStat2) Line !read into a line - IF (ErrStat2 /= 0) EXIT ! break out of the loop if it couldn't read the line (i.e. if at end of file) - i = i+1 - READ(Line,*,IOSTAT=ErrStatTmp) tmpReal - IF (ErrStatTmp/=0) THEN ! Not a number - IF (dataBegin) THEN - CALL SetErrStat( ErrID_Fatal,' Non-data line detected in WaveKinFile past the header lines.',ErrStat, ErrMsg, RoutineName); return + UnEcho=-1 + CALL GetNewUnit( UnIn ) + CALL OpenFInpFile( UnIn, FileName, ErrStat2, ErrMsg2); IF(Failed()) RETURN + + + CALL ReadCom( UnIn, FileName, 'MoorDyn water kinematics input file header', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + CALL ReadCom( UnIn, FileName, 'MoorDyn water kinematics input file header', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + ! ----- waves ----- + CALL ReadCom( UnIn, FileName, 'waves header', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + CALL ReadVar( UnIn, FileName, p%WaveKin , 'WaveKinMod' , 'WaveKinMod' , ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + ! log the method being used + IF (p%writeLog > 0) THEN + IF (p%WaveKin == 2) THEN + WRITE(p%UnLog, '(A)' ) " WaveKinMod = 2. Reading in the user provided wave grid and using SeaState for frequency calculation." + ELSE IF (p%WaveKin == 1) THEN + WRITE(p%UnLog, '(A)' ) " WaveKinMod = 1. Reading in the user provided wave grid and frequency information." + ELSE IF (p%WaveKin == 0) THEN + WRITE(p%UnLog, '(A)' ) " WaveKinMod = 0. No wave kinematics enabled." + ELSE + WRITE(p%UnLog, '(A)' ) " Invalid value for WaveKinMod" + Call SetErrStat(ErrID_Fatal, "Invalid value for WaveKinMod", ErrStat, ErrMsg, RoutineName); RETURN + ENDIF + ENDIF + CALL ReadVar( UnIn, FileName, WaveKinFile, 'WaveKinFile', 'WaveKinFile' , ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + CALL ReadVar( UnIn, FileName, p%dtWave , 'dtWave', 'time step for waves', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + CALL ReadVar( UnIn, FileName, WaveDir , 'WaveDir' , 'wave direction', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + ! X grid points + CALL ReadVar( UnIn, FileName, coordtype , 'coordtype' , '', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN ! get the entry type + READ(UnIn, '(A)', IOSTAT=ErrStat2) entries2; IF(ErrStat2 /= ErrID_None) ErrMsg2 = "There was an error reading in the grid description"; IF(Failed()) RETURN ! get entries as string to be processed + CALL gridAxisCoords(coordtype, entries2, p%pxWave, p%nxWave); IF(Failed()) RETURN + ! Y grid points + CALL ReadVar( UnIn, FileName, coordtype , 'coordtype' , '', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN ! get the entry type + READ(UnIn, '(A)', IOSTAT=ErrStat2) entries2; IF(ErrStat2 /= ErrID_None) ErrMsg2 = "There was an error reading in the grid description"; IF(Failed()) RETURN ! get entries as string to be processed + CALL gridAxisCoords(coordtype, entries2, p%pyWave, p%nyWave); IF(Failed()) RETURN + ! Z grid points + CALL ReadVar( UnIn, FileName, coordtype , 'coordtype' , '', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN ! get the entry type + READ(UnIn, '(A)', IOSTAT=ErrStat2) entries2; IF(ErrStat2 /= ErrID_None) ErrMsg2 = "There was an error reading in the grid description"; IF(Failed()) RETURN ! get entries as string to be processed + CALL gridAxisCoords(coordtype, entries2, p%pzWave, p%nzWave); IF(Failed()) RETURN + ! ----- current ----- + CALL ReadCom( UnIn, FileName, 'current header', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + CALL ReadVar( UnIn, FileName, p%Current, 'CurrentMod', 'CurrentMod', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN + ! read in the current profile + IF (p%Current == 2) THEN + + ! log the method being used + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " CurrentMod = 2. Reading in the user provided current grid and using SeaState for profile calculation." + ENDIF + + ! Z grid points + CALL ReadVar( UnIn, FileName, coordtype , 'coordtype' , '', ErrStat2, ErrMsg2, UnEcho); IF(Failed()) RETURN ! get the entry type + READ(UnIn, '(A)', IOSTAT=ErrStat2) entries2; IF(ErrStat2 /= ErrID_None) ErrMsg2 = "There was an error reading in the grid description"; IF(Failed()) RETURN ! get entries as string to be processed + CALL gridAxisCoords(coordtype, entries2, pzCurrentTemp, p%nzCurrent); IF(Failed()) RETURN ! max size of 100 because gridAxisCoords has a 100 element temporary array used for processing entries2 + uxCurrentTemp = 0.0 ! set these to zero to avoid unitialized values. This will be set later in this routine by SeaState + uyCurrentTemp = 0.0 ! set these to zero to avoid unitialized values. This will be set later in this routine by SeaState + + ELSE IF (p%Current == 1) THEN + + ! read two header lines. Old file format has 2 header lines, new file format has four (two info lines for CurrentMod = 2 and then the two headers). + CALL ReadCom( UnIn, FileName, 'current profile header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return + CALL ReadCom( UnIn, FileName, 'current profile header', ErrStat2, ErrMsg2, UnEcho); if(Failed()) return + + CALL AllocAry(pzCurrentTemp, 100, 'pzCurrentTemp', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN ! allocate pzCurrentTemp if reading in user provided current table + + DO I=1,4 ! ignore lines until table is found, or exit after 3 lines. Old file format has 2 header lines, new file format has four (two info lines for CurrentMod = 2 and then the two headers). + READ(UnIn, *, IOSTAT=ErrStat2) pzCurrentTemp(1), uxCurrentTemp(1), uyCurrentTemp(1) ! try to read into a line to first elements in the array + IF (ErrStat2 == 0) THEN + EXIT ! break out of the loop if successfully reads first table line END IF - numHdrLn = numHdrLn + 1 - ELSE - dataBegin = .TRUE. - END IF - END DO + ENDDO + + IF (I == 4) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR: Could not read the current profile table from the input file. Check the file format." + ENDIF + CALL SetErrStat(ErrID_Fatal, "Could not read the current profile table from the input file. Check the file format.", ErrStat, ErrMsg, RoutineName); RETURN + ENDIF + + ! log the first line read + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " CurrentMod = 1. Reading in the user provided current profile." + IF (p%writeLog > 1) THEN + WRITE(p%UnLog, '(A)' ) " The first line read from the current table is:" + WRITE(p%UnLog, '(A)' ) " "//trim(num2lstr(pzCurrentTemp(1)))//" "//trim(num2lstr(uxCurrentTemp(1)))//" "//trim(num2lstr(uyCurrentTemp(1))) + ENDIF + ENDIF + + ! current profile table... (read until no more rows in table) + DO I=2,100 ! start at second row because first is already read + READ(UnIn, *, IOSTAT=ErrStat2) pzCurrentTemp(i), uxCurrentTemp(i), uyCurrentTemp(i) ! read into a line + IF (ErrStat2 /= ErrID_None) THEN + p%nzCurrent = i-1 ! save number of valid current depth points in profile + EXIT ! break out of the loop if it couldn't read the line (i.e. if at end of file) + END IF + IF (i == 100) THEN + p%nzCurrent = 100 + CALL WrScr("WARNING: MD can handle a maximum of 100 current profile points") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING: MD can handle a maximum of 100 current profile points" + ENDIF + EXIT + END IF + END DO - ! rewind to start of input file to re-read things now that we know how long it is - REWIND(UnElev) + ! Check that all z values are below the water line + DO I=1,p%nzCurrent + IF (pzCurrentTemp(I) > 0) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING: MoorDyn current profile z values are above the water line." + ENDIF + CALL SetErrStat(ErrID_Warn, "Mooring current profile z values are above the water line.", ErrStat, ErrMsg, RoutineName) + END IF + END DO - ntIn = i-numHdrLn ! save number of lines of file - + ! Check that order is valid (deepest to shallowest depths), flip array if necessary + IF (pzCurrentTemp(1) > pzCurrentTemp(p%nzCurrent)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " INFO: Current profile is decreasing depths. MoorDyn needs increasing depths. Flipping arrays" + ENDIF + DO I=1,p%nzCurrent/2 + tmpReal = pzCurrentTemp(I) + pzCurrentTemp(I) = pzCurrentTemp(p%nzCurrent-I+1) + pzCurrentTemp(p%nzCurrent-I+1) = tmpReal + tmpReal = uxCurrentTemp(I) + uxCurrentTemp(I) = uxCurrentTemp(p%nzCurrent-I+1) + uxCurrentTemp(p%nzCurrent-I+1) = tmpReal + tmpReal = uyCurrentTemp(I) + uyCurrentTemp(I) = uyCurrentTemp(p%nzCurrent-I+1) + uyCurrentTemp(p%nzCurrent-I+1) = tmpReal + END DO + END IF - ! allocate space for input wave elevation array (including time column) - CALL AllocAry(WaveTimeIn, ntIn, 'WaveTimeIn', ErrStat2, ErrMsg2 ); if(Failed()) return - CALL AllocAry(WaveElevIn, ntIn, 'WaveElevIn', ErrStat2, ErrMsg2 ); if(Failed()) return + ! Check for valid array (strictly increasing z values) + DO I=1,p%nzCurrent-1 + IF (pzCurrentTemp(I) > pzCurrentTemp(I+1)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR: Current profile z values are not strictly increasing. Check the file format." + ENDIF + CALL SetErrStat(ErrID_Fatal, "Current profile z values are not strictly increasing. Check the file format.", ErrStat, ErrMsg, RoutineName); RETURN + END IF + END DO - ! read the data in from the file - DO i = 1, numHdrLn - READ(UnElev,'(A)',IOSTAT=ErrStat2) Line ! skip header lines - END DO + ELSE IF (p%Current == 0) THEN + ! log the method being used + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " CurrentMod = 0. No currents will be simulated." + ENDIF + ELSE + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " Invalid value for CurrentMod" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Invalid value for CurrentMod", ErrStat, ErrMsg, RoutineName); RETURN + ENDIF ! if p%Current != 1 or 2, then no current - DO i = 1, ntIn - READ (UnElev, *, IOSTAT=ErrStat2) WaveTimeIn(i), WaveElevIn(i) - - IF ( ErrStat2 /= 0 ) THEN - CALL SetErrStat( ErrID_Fatal,'Error reading WaveElev input file.',ErrStat, ErrMsg, RoutineName); return - END IF - END DO - ! Close the inputs file - CLOSE ( UnElev ) + CLOSE(UnIn) - IF (WaveTimeIn(1) .NE. 0.0) THEN - CALL SetErrStat( ErrID_Fatal, ' MoorDyn WaveElev time series should start at t = 0 seconds.',ErrStat, ErrMsg, RoutineName); return + IF (p%Current > 2 .OR. p%WaveKin > 2) THEN + CALL SetErrStat( ErrID_Fatal,'WaveKinMod and CurrentMod must be less than 3 if using user defined MoorDyn grid',ErrStat, ErrMsg, RoutineName); RETURN + RETURN + END IF + + ! set p%WaterKin, note that p%WaterKin is not used in this routine, but is necessary for getWaterKin + IF (p%Current == p%WaveKin) THEN ! if they are the same, set WaveKin as that value + p%WaterKin = p%WaveKin + ELSEIF (p%Current == 0) THEN ! if one is zero, use the other for water kin + p%WaterKin = p%WaveKin + ELSEIF (p%WaveKin == 0) THEN ! if one is zero, use the other for water kin + p%WaterKin = p%Current + ELSE + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR WaveKinMod and CurrentMod must be equal or one must be zero" + ENDIF + CALL SetErrStat( ErrID_Fatal,'WaveKinMod and CurrentMod must be equal or one must be zero',ErrStat, ErrMsg, RoutineName); RETURN ! TODO: can we find a way to enable wave mod = 2 and current mod = 1? ENDIF + + ! ------------------- start with wave kinematics ----------------------- - call WrScr( "Read "//trim(num2lstr(ntIn))//" time steps from input file." ) + IF (p%WaveKin > 0) THEN + + ! Check that all wave grid z values are below the water line, otherwise COSHNumOvrCOSHDen calcs will nan + DO I=1,p%nzWave + IF (p%pzWave(I) > 0) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR: MoorDyn wave grid z values are above the water line. The wave grid cannot contain values greater than 0. This may be caused by interpolation over integer z-grid spacing." + ENDIF + CALL SetErrStat(ErrID_Fatal, "Mooring wave grid z values are above the water line. The wave grid cannot contain values greater than 0.", ErrStat, ErrMsg, RoutineName); RETURN + END IF + END DO + + IF (p%WaveKin == 2) THEN ! must be using the hybrid approach + + ! Error check to make sure the wave field pointer is not null + IF (.NOT. ASSOCIATED(p%WaveField)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR WaveField pointer is null. Hybrid method requires SeaState to be enabled. Please check input files." + ENDIF + CALL SetErrStat(ErrID_Fatal, "WaveField pointer is null. Hybrid method requires SeaState to be enabled. Please check input files.", ErrStat, ErrMsg, RoutineName); RETURN + END IF + + ! Warning check to make sure SeaState and MoorDyn have the same water depth + IF (p%WaveField%WtrDpth /= p%WtrDpth) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING SeaState water depth does not match MoorDyn water depth. Using MoorDyn values for interpolating SeaState data to MoorDyn grid." + ENDIF + CALL SetErrStat(ErrID_Warn, "SeaState water depth does not match MoorDyn water depth. Using MoorDyn values for interpolating SeaState data to MoorDyn grid.", ErrStat, ErrMsg, RoutineName) + END IF + + ! Error check to make sure SeaState and MoorDyn have the same water density and gravity. This also checks the pointer is valid. + IF (p%WaveField%RhoXg /= REAL((p%rhoW * p%g), SiKi)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR SeaState (water density * gravity) ["//trim(num2lstr(p%WaveField%RhoXg))//"] does not match MoorDyn water (density * gravity) ["//trim(num2lstr(REAL((p%rhoW * p%g), SiKi)))//"]. The SeaState pointer may be corrupted." + ENDIF + CALL SetErrStat(ErrID_Fatal, " ERROR SeaState (water density * gravity) ["//trim(num2lstr(p%WaveField%RhoXg))//"] does not match MoorDyn water (density * gravity) ["//trim(num2lstr(REAL((p%rhoW * p%g), SiKi)))//"]. The SeaState pointer may be corrupted.", ErrStat, ErrMsg, RoutineName) + END IF - ! if (WaveTimeIn(ntIn) < TMax) then <<<< need to handle if time series is too short? - - ! specify stepping details - p%ntWave = CEILING(Tmax/p%dtWave) ! number of wave time steps + ! Warning check to make sure SeaState and MoorDyn have the same wave dir. For now, no wave spreading. This can be updated + IF (p%WaveField%WaveDir /= WaveDir) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING SeaState WaveDir does not match MoorDyn WaveDir. Using MoorDyn values for interpolating SeaState data to MoorDyn grid." + ENDIF + CALL SetErrStat(ErrID_Warn, "SeaState WaveDir does not match MoorDyn WaveDir.Using MoorDyn values for interpolating SeaState data to MoorDyn grid.", ErrStat, ErrMsg, RoutineName) + END IF - - ! allocate space for processed reference wave elevation time series - ALLOCATE ( WaveElev0( 0:p%ntWave ), STAT=ErrStatTmp ) ! this has an extra entry of zero in case it needs to be padded to be even - IF (ErrStatTmp /= 0) CALL SetErrStat(ErrID_Fatal,'Cannot allocate array WaveElev0.',ErrStat,ErrMsg,RoutineName) - WaveElev0 = 0.0_SiKi - - ! go through and interpolate (should replace with standard function) - DO i = 1, p%ntWave - t = p%dtWave*(i-1) - - ! interpolation routine - DO iIn = 1,ntIn-1 - IF (WaveTimeIn(iIn+1) > t) THEN ! find the right two points to interpolate between (remember that the first column of PtfmMotIn is time) - frac = (t - WaveTimeIn(iIn) )/( WaveTimeIn(iIn+1) - WaveTimeIn(iIn) ) ! interpolation fraction (0-1) between two interpolation points - WaveElev0(i-1) = WaveElevIn(iIn) + frac*(WaveElevIn(iIn+1) - WaveElevIn(iIn)) ! get interpolated wave elevation - EXIT ! break out of the loop for this time step once we've done its interpolation + ! Info check for if MoorDyn dtWave is non-zero. Users may set this accidentially, or could be left over in an input file + IF (p%dtWave > 0) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " MoorDyn dtWave is ignored when using WaveKinMod = 2 becasue wave frequency information is supplied by SeaState" + ENDIF END IF - END DO - END DO - - ! note: following is adapted from UserWaves.v90 UserWaveElevations_Init - - - - ! Set new value for NStepWave so that the FFT algorithms are efficient. We will use the values passed in rather than what is read from the file - - IF ( MOD(p%ntWave,2) == 1 ) p%ntWave = p%ntWave + 1 ! Set NStepWave to an even integer - NStepWave2 = MAX( p%ntWave/2, 1 ) ! Make sure that NStepWave is an even product of small factors (PSF) that is - NStepWave = 2*PSF ( NStepWave2, 9 ) ! greater or equal to WaveTMax/WaveDT to ensure that the FFT is efficient. - NStepWave2 = NStepWave/2 ! Update the value of NStepWave2 based on the value needed for NStepWave. - WaveTMax = NStepWave*p%dtWave ! Update the value of WaveTMax based on the value needed for NStepWave. - WaveDOmega = TwoPi/TMax ! Compute the frequency step for incident wave calculations. - p%ntWave = NStepWave - - - - ! Allocate array to hold the wave elevations for calculation of FFT. - ALLOCATE ( TmpFFTWaveElev( 0:NStepWave-1 ), STAT=ErrStatTmp ) - IF (ErrStatTmp /= 0) CALL SetErrStat(ErrID_Fatal,'Cannot allocate array TmpFFTWaveElev.',ErrStat,ErrMsg,RoutineName) + ! check SS only is being run w/ no wave spreading, because MoorDyn is not compatiable with those(for now) + IF (p%WaveField%WaveMultiDir) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR MoorDyn WaveKinMod2 does not support wave spreading. Please use WaveDirMod = 0 in SeaState." + ENDIF + CALL SetErrStat(ErrID_Fatal, "MoorDyn WaveKinMod2 does not support wave spreading. Please use WaveDirMod = 0 in SeaState.", ErrStat, ErrMsg, RoutineName); RETURN + END IF - ! Allocate frequency array for the wave elevation information in frequency space - ALLOCATE ( WaveElevC0(2, 0:NStepWave2 ), STAT=ErrStatTmp ) - IF (ErrStatTmp /= 0) CALL SetErrStat(ErrID_Fatal,'Cannot allocate array WaveElevC0.',ErrStat,ErrMsg,RoutineName) - + ! check SS only is being run w/ first order waves, because MoorDyn is not compatiable with those(for now) + IF (ALLOCATED(p%WaveField%WaveElev2)) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR MoorDyn WaveKinMod2 does not support second order waves. Please use disable in SeaState." + ENDIF + CALL SetErrStat(ErrID_Fatal, "MoorDyn WaveKinMod2 does not support second order waves. Please use disable in SeaState.", ErrStat, ErrMsg, RoutineName); RETURN + END IF - ! Now check if all the allocations worked properly - IF ( ErrStat >= AbortErrLev ) THEN - CALL CleanUp() - RETURN - END IF - - ! Set the values - TmpFFTWaveElev = 0.0_DbKi - WaveElevC0(:,:) = 0.0_DbKi + ! set dtWave to WaveField dtWave (for interpolation use in getWaterKin) + p%dtWave = p%WaveField%WaveTime(2)-p%WaveField%WaveTime(1) ! from Waves.f90 (line 1040), DO I = 0,WaveField%NStepWave; WaveField%WaveTime(I) = I*REAL(InitInp%WaveDT,SiKi) + + ! Interpolations need the following from SeaState: WaveElevC0, NStepWave2, NStepWave, WaveDOmega. In p%WaveKin = 1, these values are extracted from the wave elevation time series + WaveElevC0 = p%WaveField%WaveElevC0 + WaveDOmega = p%WaveField%WaveDOmega + NStepWave2 = p%WaveField%NStepWave2 + NStepWave = p%WaveField%NStepWave + p%ntWave = NStepWave ! set ntWave to NStepWave + + ELSEIF (p%WaveKin == 1) THEN ! must be a filepath therefore read wave elevations from timeseries + + ! NOTE: there is a decent ammount of code duplication (intentional for now) with what is in SeaState that eventually + ! we will want to synchronize with a standard library at some point + + ! --------------------- set from inputted wave elevation time series, grid approach ------------------- + CALL WrScr( ' WaveKinMod = 1. Reading wave elevation time series from file' ) + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WaveKinMod = 1. Reading wave elevation time series from file" + ENDIF + + IF ( LEN_TRIM( WaveKinFile ) == 0 ) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR WaveKinFile must not be an empty string." + ENDIF + CALL SetErrStat( ErrID_Fatal,'WaveKinFile must not be an empty string.',ErrStat, ErrMsg, RoutineName); RETURN + RETURN + END IF + + IF ( PathIsRelative( WaveKinFile ) ) THEN ! properly handle relative path <<< + !CALL GetPath( TRIM(InitInp%InputFile), TmpPath ) + WaveKinFile = TRIM(p%PriPath)//TRIM(WaveKinFile) + END IF + + ! note: following is adapted from MoorDyn_Driver + + CALL GetNewUnit( UnElev ) + + CALL OpenFInpFile ( UnElev, WaveKinFile, ErrStat2, ErrMsg2 ); IF(Failed()) RETURN + + CALL WrScr( ' Reading wave elevation data from '//trim(WaveKinFile) ) + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) ' Reading wave elevation data from '//trim(WaveKinFile) + ENDIF + + ! Read through length of file to find its length + i = 0 ! start line counter + numHdrLn = 0 ! start header-line counter + dataBegin = .FALSE. ! started reading the data section + DO + READ(UnElev,'(A)',IOSTAT=ErrStat2) Line !read into a line + IF (ErrStat2 /= 0) EXIT ! break out of the loop if it couldn't read the line (i.e. if at end of file) + i = i+1 + READ(Line,*,IOSTAT=ErrStat2) tmpReal + IF (ErrStat2/=0) THEN ! Not a number + IF (dataBegin) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR Non-data line detected in WaveKinFile past the header lines." + ENDIF + CALL SetErrStat( ErrID_Fatal,' Non-data line detected in WaveKinFile past the header lines.',ErrStat, ErrMsg, RoutineName); RETURN + END IF + numHdrLn = numHdrLn + 1 + ELSE + dataBegin = .TRUE. + END IF + END DO + + ! rewind to start of input file to re-read things now that we know how long it is + REWIND(UnElev) + + ntIn = i-numHdrLn ! save number of lines of file + + ! allocate space for input wave elevation array (including time column) + CALL AllocAry(WaveTimeIn, ntIn, 'WaveTimeIn', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN + CALL AllocAry(WaveElevIn, ntIn, 'WaveElevIn', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN - ! Copy values over - DO I=0, MIN(SIZE(WaveElev0), NStepWave)-1 - TmpFFTWaveElev(I) = WaveElev0(I) - ENDDO + ! read the data in from the file + DO i = 1, numHdrLn + READ(UnElev,'(A)',IOSTAT=ErrStat2) Line ! skip header lines + END DO + + DO i = 1, ntIn + READ (UnElev, *, IOSTAT=ErrStat2) WaveTimeIn(i), WaveElevIn(i) ! read wave elevation time series + + IF ( ErrStat2 /= 0 ) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR reading WaveElev input file." + ENDIF + CALL SetErrStat( ErrID_Fatal,'Error reading WaveElev input file.',ErrStat, ErrMsg, RoutineName); RETURN + END IF + END DO + + ! Close the inputs file + CLOSE ( UnElev ) + + IF (WaveTimeIn(1) .NE. 0.0) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)') "ERROR MoorDyn WaveElev time series should start at t = 0 seconds." + ENDIF + CALL SetErrStat( ErrID_Fatal, ' MoorDyn WaveElev time series should start at t = 0 seconds.',ErrStat, ErrMsg, RoutineName); RETURN + ENDIF + + CALL WrScr( " Read "//trim(num2lstr(ntIn))//" time steps from input file." ) + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " Read "//trim(num2lstr(ntIn))//" time steps from input file." + ENDIF - ! Initialize the FFT - CALL InitFFT ( NStepWave, FFT_Data, .FALSE., ErrStatTmp ) - CALL SetErrStat(ErrStatTmp,'Error occured while initializing the FFT.',ErrStat,ErrMsg,RoutineName); if(Failed()) return + ! IF (WaveTimeIn(ntIn) < TMax) THEN <<<< need to handle if time series is too short? + + ! specify stepping details + p%ntWave = CEILING(Tmax/p%dtWave) ! number of wave time steps - ! Apply the forward FFT to get the real and imaginary parts of the frequency information. - CALL ApplyFFT_f ( TmpFFTWaveElev(:), FFT_Data, ErrStatTmp ) ! Note that the TmpFFTWaveElev now contains the real and imaginary bits. - CALL SetErrStat(ErrStatTmp,'Error occured while applying the forwards FFT to TmpFFTWaveElev array.',ErrStat,ErrMsg,RoutineName); if(Failed()) return + + ! allocate space for processed reference wave elevation time series + ALLOCATE ( WaveElev0( 0:p%ntWave ), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate array WaveElev0.'; IF(Failed0()) RETURN ! this has an extra entry of zero in case it needs to be padded to be even + WaveElev0 = 0.0_SiKi + + ! go through and interpolate (should replace with standard function) + DO i = 1, p%ntWave + t = p%dtWave*(i-1) + + ! interpolation routine to dtWave spacing + DO iIn = 1,ntIn-1 + IF (WaveTimeIn(iIn+1) > t) THEN ! find the right two points to interpolate between (remember that the first column of PtfmMotIn is time) + frac = (t - WaveTimeIn(iIn) )/( WaveTimeIn(iIn+1) - WaveTimeIn(iIn) ) ! interpolation fraction (0-1) between two interpolation points + WaveElev0(i-1) = WaveElevIn(iIn) + frac*(WaveElevIn(iIn+1) - WaveElevIn(iIn)) ! get interpolated wave elevation + EXIT ! break out of the loop for this time step once we've DOne its interpolation + END IF + END DO + END DO + + ! note: following is adapted from UserWaves.v90 UserWaveElevations_Init + + + + ! Set new value for NStepWave so that the FFT algorithms are efficient. We will use the values passed in rather than what is read from the file + + IF ( MOD(p%ntWave,2) == 1 ) p%ntWave = p%ntWave + 1 ! Set NStepWave to an even integer + NStepWave2 = MAX( p%ntWave/2, 1 ) ! Make sure that NStepWave is an even product of small factors (PSF) that is + NStepWave = 2*PSF ( NStepWave2, 9 ) ! greater or equal to WaveTMax/WaveDT to ensure that the FFT is efficient. + NStepWave2 = NStepWave/2 ! Update the value of NStepWave2 based on the value needed for NStepWave. + WaveTMax = NStepWave*p%dtWave ! Update the value of WaveTMax based on the value needed for NStepWave. + WaveDOmega = TwoPi/TMax ! Compute the frequency step for incident wave calculations. + p%ntWave = NStepWave + + + - ! Copy the resulting TmpFFTWaveElev(:) data over to the WaveElevC0 array - DO I=1,NStepWave2-1 - WaveElevC0 (1,I) = TmpFFTWaveElev(2*I-1) - WaveElevC0 (2,I) = TmpFFTWaveElev(2*I) - ENDDO - WaveElevC0(:,NStepWave2) = 0.0_SiKi + ! Allocate array to hold the wave elevations for calculation of FFT. + ALLOCATE ( TmpFFTWaveElev( 0:NStepWave-1), STAT=ErrStat2 ); ErrMsg2 = 'Cannot allocate array TmpFFTWaveElev.'; IF (Failed0()) RETURN - CALL ExitFFT(FFT_Data, ErrStatTmp) - CALL SetErrStat(ErrStatTmp,'Error occured while cleaning up after the FFTs.', ErrStat,ErrMsg,RoutineName); if(Failed()) return + ! Allocate frequency array for the wave elevation information in frequency space + ALLOCATE ( WaveElevC0(2, 0:NStepWave2), STAT=ErrStat2 ); ErrMsg2 = 'Cannot allocate array WaveElevC0.'; IF (Failed0()) RETURN + + ! Now check if all the allocations worked properly + IF ( ErrStat >= AbortErrLev ) THEN + CALL CleanUp() + RETURN + END IF + + ! Set the values + TmpFFTWaveElev = 0.0_DbKi + WaveElevC0(:,:) = 0.0_DbKi - IF (ALLOCATED( WaveElev0 )) DEALLOCATE( WaveElev0 , STAT=ErrStatTmp) - IF (ALLOCATED( TmpFFTWaveElev )) DEALLOCATE( TmpFFTWaveElev, STAT=ErrStatTmp) + ! Copy values over + DO I=0, MIN(SIZE(WaveElev0), NStepWave)-1 + TmpFFTWaveElev(I) = WaveElev0(I) + ENDDO - - ! note: following is a very streamlined adaptation from from Waves.v90 VariousWaves_Init - - ! allocate all the wave kinematics FFT arrays - ALLOCATE( WaveNmbr (0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveNmbr. ',ErrStat,ErrMsg,RoutineName) - ALLOCATE( tmpComplex(0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate tmpComplex.',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveElevC (0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveElevC .',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveDynPC (0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveDynPC .',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveVelCHx(0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveVelCHx.',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveVelCHy(0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveVelCHy.',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveVelCV (0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveVelCV .',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveAccCHx(0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveAccCHx.',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveAccCHy(0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveAccCHy.',ErrStat,ErrMsg,RoutineName) - ALLOCATE( WaveAccCV (0:NStepWave2), STAT=ErrStatTmp); CALL SetErrStat(ErrStatTmp,'Cannot allocate WaveAccCV .',ErrStat,ErrMsg,RoutineName) - - ! allocate time series grid data arrays (now that we know the number of time steps coming from the IFFTs) - CALL allocateKinematicsArrays() - - - ! Set the CosWaveDir and SinWaveDir values - CosWaveDir=COS(D2R*WaveDir) - SinWaveDir=SIN(D2R*WaveDir) - - ! get wave number array once - DO I = 0, NStepWave2 - WaveNmbr(i) = WaveNumber ( REAL(I*WaveDOmega, R8Ki), p%g, p%WtrDpth ) - tmpComplex(I) = CMPLX(WaveElevC0(1,I), WaveElevC0(2,I)) - END DO - - ! set up FFTer for doing IFFTs - CALL InitFFT ( NStepWave, FFT_Data, .TRUE., ErrStatTmp ) - CALL SetErrStat(ErrStatTmp,'Error occured while initializing the FFT.', ErrStat, ErrMsg, routineName); if(Failed()) return - - ! Loop through all points where the incident wave kinematics will be computed - do ix = 1,p%nxWave - do iy = 1,p%nyWave - do iz = 1,p%nzWave - - ! Compute the discrete Fourier transform of the incident wave kinematics - do i = 0, NStepWave2 ! Loop through the positive frequency components (including zero) of the discrete Fourier transforms - - Omega = i*WaveDOmega - ImagOmega = ImagNmbr*Omega - - WaveElevC (i) = tmpComplex(i) * EXP( -ImagNmbr*WaveNmbr(i)*( p%pxWave(ix)*CosWaveDir + p%pyWave(iy)*SinWaveDir )) - WaveDynPC (i) = p%rhoW*p%g* WaveElevC(i) * COSHNumOvrCOSHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) - WaveVelCHx(i) = Omega*WaveElevC(i) * COSHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) *CosWaveDir - WaveVelCHy(i) = Omega*WaveElevC(i) * COSHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) *SinWaveDir - WaveVelCV (i) = ImagOmega*WaveElevC(i) * SINHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) - WaveAccCHx(i) = ImagOmega*WaveVelCHx(i) - WaveAccCHy(i) = ImagOmega*WaveVelCHy(i) - WaveAccCV (i) = ImagOmega*WaveVelCV (i) - end do ! I, frequencies + ! Initialize the FFT + CALL InitFFT ( NStepWave, FFT_Data, .FALSE., ErrStat2 ); ErrMsg2 = 'Error occured while initializing the FFT.'; IF(Failed()) RETURN + + ! Apply the forward FFT on the wave elevation timeseries to get the real and imaginary parts of the frequency information. + CALL ApplyFFT_f ( TmpFFTWaveElev(:), FFT_Data, ErrStat2 ); ErrMsg2 = 'Error occured while applying the forwards FFT to TmpFFTWaveElev array.'; IF(Failed()) RETURN ! Note that the TmpFFTWaveElev now contains the real and imaginary bits. + + ! Copy the resulting TmpFFTWaveElev(:) data over to the WaveElevC0 array + DO I=1,NStepWave2-1 + WaveElevC0 (1,I) = TmpFFTWaveElev(2*I-1) ! all the odd indicies in array (real part ?) + WaveElevC0 (2,I) = TmpFFTWaveElev(2*I) ! all the even indicies in array (imaginary part ?) + ENDDO + WaveElevC0(:,NStepWave2) = 0.0_SiKi + + CALL ExitFFT(FFT_Data, ErrStat2); ErrMsg2 = 'Error occured while cleaning up after the FFTs.'; IF(Failed()) RETURN + + IF (ALLOCATED( WaveElev0 )) DEALLOCATE( WaveElev0 , STAT=ErrStat2); ErrMsg2 = 'Cannot deallocate WaveElev0.' ; IF (Failed0()) RETURN + IF (ALLOCATED( TmpFFTWaveElev )) DEALLOCATE( TmpFFTWaveElev, STAT=ErrStat2); ErrMsg2 = 'Cannot deallocate TmpFFTWaveElev.'; IF (Failed0()) RETURN + + ENDIF ! End getting frequency data either from time series or WaveField. The below needs to happen for both old and hybrid approach + + ! NOTE: there is a decent ammount of code duplication (intentional for now) with what is in SeaState that eventually + ! we will want to synchronize with a standard library at some point + + ! allocate all the wave kinematics FFT arrays + ALLOCATE( WaveNmbr (0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveNmbr.' ; IF (Failed0()) RETURN + ALLOCATE( tmpComplex(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate tmpComplex.'; IF (Failed0()) RETURN + ALLOCATE( WaveElevC (0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveElevC.' ; IF (Failed0()) RETURN + ALLOCATE( WaveDynPC (0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveDynPC.' ; IF (Failed0()) RETURN + ALLOCATE( WaveVelCHx(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveVelCHx.'; IF (Failed0()) RETURN + ALLOCATE( WaveVelCHy(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveVelCHy.'; IF (Failed0()) RETURN + ALLOCATE( WaveVelCV (0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveVelCV.' ; IF (Failed0()) RETURN + ALLOCATE( WaveAccCHx(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveAccCHx.'; IF (Failed0()) RETURN + ALLOCATE( WaveAccCHy(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveAccCHy.'; IF (Failed0()) RETURN + ALLOCATE( WaveAccCV (0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate WaveAccCV .'; IF (Failed0()) RETURN + + ! allocate time series grid data arrays (now that we know the number of time steps coming from the IFFTs) + CALL allocateKinematicsArrays() + + ! Set the CosWaveDir and SinWaveDir values + CosWaveDir=COS(D2R*WaveDir) + SinWaveDir=SIN(D2R*WaveDir) + + ! get wave number array once + DO I = 0, NStepWave2 + WaveNmbr(i) = WaveNumber ( REAL(I*WaveDOmega, R8Ki), p%g, p%WtrDpth ) + tmpComplex(I) = CMPLX(WaveElevC0(1,I), WaveElevC0(2,I)) ! 1 are real, 2 are imaginary + END DO + + ! set up FFTer for doing IFFTs + CALL InitFFT ( NStepWave, FFT_Data, .TRUE., ErrStat2 ); ErrMsg2 = 'Error occured while initializing the FFT.'; IF(Failed()) RETURN + + ! Loop through all points where the incident wave kinematics will be computed + DO ix = 1,p%nxWave + DO iy = 1,p%nyWave + DO iz = 1,p%nzWave - ! now IFFT all the wave kinematics except surface elevation and save it into the grid of data - CALL ApplyFFT_cx( p%PDyn (:,iz,iy,ix), WaveDynPC , FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveDynP.', ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%uxWave(:,iz,iy,ix), WaveVelCHx, FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveVelHx.',ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%uyWave(:,iz,iy,ix), WaveVelCHy, FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveVelHy.',ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%uzWave(:,iz,iy,ix), WaveVelCV , FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveVelV.', ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%axWave(:,iz,iy,ix), WaveAccCHx, FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveAccHx.',ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%ayWave(:,iz,iy,ix), WaveAccCHy, FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveAccHy.',ErrStat,ErrMsg,RoutineName) - CALL ApplyFFT_cx( p%azWave(:,iz,iy,ix), WaveAccCV , FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveAccV.', ErrStat,ErrMsg,RoutineName) - - end do ! iz - - ! IFFT wave elevation here because it's only at the surface - CALL ApplyFFT_cx( p%zeta(:,iy,ix) , WaveElevC , FFT_Data, ErrStatTmp ); CALL SetErrStat(ErrStatTmp,'Error IFFTing WaveElev.', ErrStat,ErrMsg,RoutineName) - end do ! iy - end do ! ix - - ! could also reproduce the wave elevation at 0,0,0 on a separate channel for verification... + ! Compute the discrete Fourier transform of the incident wave kinematics + DO i = 0, NStepWave2 ! Loop through the positive frequency components (including zero) of the discrete Fourier transforms + + Omega = i*WaveDOmega + ImagOmega = ImagNmbr*Omega + + WaveElevC (i) = tmpComplex(i) * EXP( -ImagNmbr*WaveNmbr(i)*( p%pxWave(ix)*CosWaveDir + p%pyWave(iy)*SinWaveDir )) ! Discrete Fourier transform of the instantaneous elevation of incident waves at the ref point (meters) + WaveDynPC (i) = p%rhoW*p%g* WaveElevC(i) * COSHNumOvrCOSHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) ! Discrete Fourier transform of the instantaneous dynamic pressure of incident waves before applying stretching at the zi-coordinates for points (N/m^2) + WaveVelCHx(i) = Omega*WaveElevC(i) * COSHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) *CosWaveDir ! Discrete Fourier transform of the instantaneous horizontal velocity in x-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2) + WaveVelCHy(i) = Omega*WaveElevC(i) * COSHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) *SinWaveDir ! Discrete Fourier transform of the instantaneous horizontal velocity in y-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2) + WaveVelCV (i) = ImagOmega*WaveElevC(i) * SINHNumOvrSINHDen( WaveNmbr(i), p%WtrDpth, REAL(p%pzWave(iz), R8Ki) ) ! Discrete Fourier transform of the instantaneous vertical velocity of incident waves before applying stretching at the zi-coordinates for points (m/s) + WaveAccCHx(i) = ImagOmega*WaveVelCHx(i) ! Discrete Fourier transform of the instantaneous horizontal acceleration in x-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2) + WaveAccCHy(i) = ImagOmega*WaveVelCHy(i) ! Discrete Fourier transform of the instantaneous horizontal acceleration in y-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2) + WaveAccCV (i) = ImagOmega*WaveVelCV (i) ! Discrete Fourier transform of the instantaneous vertical acceleration of incident waves before applying stretching at the zi-coordinates for points (m/s^2) + END DO ! I, frequencies + + ! now IFFT all the wave kinematics except surface elevation and save it into the grid of data + CALL ApplyFFT_cx( p%PDyn (:,iz,iy,ix), WaveDynPC , FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveDynPC.'; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%uxWave(:,iz,iy,ix), WaveVelCHx, FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveVelHx.'; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%uyWave(:,iz,iy,ix), WaveVelCHy, FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveVelHy.'; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%uzWave(:,iz,iy,ix), WaveVelCV , FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveVelV.' ; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%axWave(:,iz,iy,ix), WaveAccCHx, FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveAccHx.'; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%ayWave(:,iz,iy,ix), WaveAccCHy, FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveAccHy.'; IF (Failed()) RETURN + CALL ApplyFFT_cx( p%azWave(:,iz,iy,ix), WaveAccCV , FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveAccV.' ; IF (Failed()) RETURN + + END DO ! iz + + ! IFFT wave elevation here because it's only at the surface + CALL ApplyFFT_cx( p%zeta(:,iy,ix) , WaveElevC , FFT_Data, ErrStat2 ); ErrMsg2 = 'Error IFFTing WaveElev.'; IF (Failed()) RETURN + END DO ! iy + END DO ! ix + + ! could also reproduce the wave elevation at 0,0,0 on a separate channel for verIFication... + + CALL ExitFFT(FFT_Data, ErrStat2); ErrMsg2 = 'Error occured while cleaning up after the IFFTs.'; IF(Failed()) RETURN - CALL ExitFFT(FFT_Data, ErrStatTmp) - CALL SetErrStat(ErrStatTmp,'Error occured while cleaning up after the IFFTs.', ErrStat,ErrMsg,RoutineName); if(Failed()) return - - end if ! p%WaveKin == 3 + END IF ! p%WaveKin > 0 - ! --------------------------------- now do currents -------------------------------- - if (p%Current == 1) then - - ! allocate current profile arrays to correct size - CALL AllocAry( p%pzCurrent, p%nzCurrent, 'pzCurrent', ErrStat2, ErrMsg2 ); if(Failed()) return - CALL AllocAry( p%uxCurrent, p%nzCurrent, 'uxCurrent', ErrStat2, ErrMsg2 ); if(Failed()) return - CALL AllocAry( p%uyCurrent, p%nzCurrent, 'uyCurrent', ErrStat2, ErrMsg2 ); if(Failed()) return + ! --------------------------------- now do currents -------------------------------- + + IF (p%Current > 0) THEN - ! copy over data, flipping sign of depth values (to be positive-up) and reversing order - do i = 1,p%nzCurrent - p%pzCurrent(i) = -pzCurrentTemp(p%nzCurrent + 1 - i) ! flip sign so depth is positive-up - p%uxCurrent(i) = uxCurrentTemp(p%nzCurrent + 1 - i) - p%uyCurrent(i) = uyCurrentTemp(p%nzCurrent + 1 - i) - end do - - end if ! p%Current == 1 + IF (p%Current == 2) THEN + + ! Overwrite the SeaState current grid to use the MoorDyn grid + p%WaveField%Current_InitInput%WaveKinGridzi = pzCurrentTemp + p%WaveField%Current_InitInput%NGridPts = p%nzCurrent + + ! Calculate the current profile in the MD grid with SS inputs + CALL Current_Init(p%WaveField%Current_InitInput, Current_InitOutput, ErrStat2, ErrMsg2); IF(Failed()) RETURN + + ! Check output current arrays are the right size + IF (SIZE(Current_InitOutput%CurrVxi) /= p%nzCurrent) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR Current_Init output size does not match the MoorDyn grid size. "//trim(num2lstr(SIZE(Current_InitOutput%CurrVxi)))//" vs "//trim(num2lstr(REAL(p%nzCurrent, SiKi))) + ENDIF + CALL SetErrStat( ErrID_Fatal,'Current_Init output size does not match the MoorDyn grid size.',ErrStat, ErrMsg, RoutineName); RETURN + ENDIF + + ! extract the currents from outputs + DO I = 1, p%nzCurrent + uxCurrentTemp(I) = Current_InitOutput%CurrVxi(I) + uyCurrentTemp(I) = Current_InitOutput%CurrVyi(I) + END DO + + ENDIF + + ! allocate current profile arrays to correct size + CALL AllocAry( p%pzCurrent, p%nzCurrent, 'pzCurrent', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN + CALL AllocAry( p%uxCurrent, p%nzCurrent, 'uxCurrent', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN + CALL AllocAry( p%uyCurrent, p%nzCurrent, 'uyCurrent', ErrStat2, ErrMsg2 ); IF(Failed()) RETURN + + ! copy over data + DO i = 1,p%nzCurrent + p%pzCurrent(i) = pzCurrentTemp(i) + p%uxCurrent(i) = uxCurrentTemp(i) + p%uyCurrent(i) = uyCurrentTemp(i) + END DO + + ENDIF ! p%Current >0 + ENDIF ! (tmpString == "SEASTATE") + + ! some status messages: + IF (p%WaterKin == 0) THEN + CALL WrScr(" No water kinematics set up") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " No water kinematics set up" + ENDIF + ELSEIF (p%WaterKin == 1) THEN + CALL WrScr(" Water kinematics will be simulated using the old method (user provided grid and water kinematic data)") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " Water kinematics will be simulated using the old method (user provided grid and water kinematic data)" + ENDIF + ELSEIF (p%WaterKin == 2) THEN + CALL WrScr(" Water kinematics will be simulated using the hybrid method (user provided grid with SeaState water kinematics data)") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " Water kinematics will be simulated using the hybrid method (user provided grid with SeaState water kinematics data)" + ENDIF + ELSEIF (p%WaterKin == 3) THEN + CALL WrScr(" Water kinematics will be simulated using the SeaState method (SeaState provided grid and water kinematics data). These will be disabled during IC solve.") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " Water kinematics will be simulated using the SeaState method (SeaState provided grid and water kinematics data). These will be disabled during IC solve." + ENDIF + ELSE + CALL SetErrStat( ErrID_Fatal,"Invalid value for WaterKin",ErrStat, ErrMsg, RoutineName); RETURN + ENDIF ! ------------------------------ clean up and finished --------------------------- CALL cleanup() @@ -1714,76 +1782,91 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) ! get grid axis coordinates, initialize/record in array, and return size - SUBROUTINE gridAxisCoords(coordtype, entries, coordarray, n, ErrStat, ErrMsg) + SUBROUTINE gridAxisCoords(coordtype, entries, coordarray, n) INTEGER(IntKi), INTENT(IN ) :: coordtype CHARACTER(*), INTENT(INOUT) :: entries REAL(SiKi), ALLOCATABLE, INTENT(INOUT) :: coordarray(:) INTEGER(IntKi), INTENT( OUT) :: n - - - INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None - + REAL(ReKi) :: tempArray (100) REAL(ReKi) :: dx INTEGER(IntKi) :: nEntries, I IF (len(trim(entries)) == len(entries)) THEN - call WrScr("Warning: Only 120 characters read from wave grid coordinates") + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING in gridAxisCoords: Only "//trim(num2lstr(len(entries)))//" characters read from wave grid coordinates" + ENDIF + CALL SetErrStat(ErrID_Warn, "Only "//trim(num2lstr(len(entries)))//" characters read from wave grid coordinates", ErrStat, ErrMsg, RoutineName) END IF IF (entries(len(entries):len(entries)) == ',') THEN - ErrStat = ErrID_Fatal - ErrMsg = 'Last character of wave grid coordinate list cannot be comma' + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in gridAxisCoords: Last character of wave grid coordinate list cannot be comma" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Last character of wave grid coordinate list cannot be comma", ErrStat, ErrMsg, RoutineName); RETURN ELSE ! get array of coordinate entries - CALL stringToArray(entries, nEntries, tempArray) + CALL stringToArray(entries, nEntries, tempArray); IF(Failed()) RETURN ! set number of coordinates - if ( coordtype==0) then ! 0: not used - make one grid point at zero + IF ( coordtype==0) THEN ! 0: not used - make one grid point at zero n = 1; - else if (coordtype==1) then ! 1: list values in ascending order + ELSE IF (coordtype==1) THEN ! 1: list values in ascending order n = nEntries - else if (coordtype==2) then ! 2: uniform specified by -xlim, xlim, num + ELSE IF (coordtype==2) THEN ! 2: uniform specified by -xlim, xlim, num n = int(tempArray(3)) - else - call WrScr("Error: invalid coordinate type specified to gridAxisCoords") - end if + + ! Check that tmpArray range / n is not an integer value, warn in log file if so + IF (MOD((tempArray(2)-tempArray(1)),REAL(n,ReKi)) == 0.0_ReKi) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " WARNING in gridAxisCoords: grid range / num grid points is an integer. This may cause interpolation issues. For integer grid spacing, add 1 to num grid points." + ENDIF + ENDIF + + ELSE + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in gridAxisCoords: invalid coordinate type specified to gridAxisCoords. Check WaterKin input file format for missing lines" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Invalid coordinate type. Check WaterKin input file format for missing lines", ErrStat, ErrMsg, RoutineName); RETURN + END IF ! allocate coordinate array CALL AllocAry(coordarray, n, 'x,y, or z grid points' , ErrStat, ErrMsg) !ALLOCATE ( coordarray(n), STAT=ErrStat) ! fill in coordinates - if ( coordtype==0) then + IF ( coordtype==0) THEN coordarray(1) = 0.0_ReKi - else if (coordtype==1) then + ELSE IF (coordtype==1) THEN coordarray(1:n) = tempArray(1:n) - else if (coordtype==2) then + ELSE IF (coordtype==2) THEN coordarray(1) = tempArray(1) coordarray(n) = tempArray(2) dx = (coordarray(n)-coordarray(1))/REAL(n-1) - do i=2,n + DO i=2,n coordarray(i) = coordarray(i-1) + dx - end do + END DO - else - call WrScr("Error: invalid coordinate type specified to gridAxisCoords") - end if + ELSE + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in gridAxisCoords: invalid coordinate type specified to gridAxisCoords. Check WaterKin input file format for missing lines" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Invalid coordinate type specified to gridAxisCoords. Check WaterKin input file format for missing lines", ErrStat, ErrMsg, RoutineName); RETURN + END IF ! print *, "Set water grid coordinates to :" ! DO i=1,n ! print *, " ", coordarray(i) - ! end do + ! END DO END IF END SUBROUTINE gridAxisCoords - ! Extract an array of numbers out of a string with comma-separated numbers (this could go in a more general location) + ! Extract an array of numbers out of a string with comma-separated numbers (cannot use NWTC function ReadArr for now becasue len not known) SUBROUTINE stringToArray(instring, n, outarray) CHARACTER(*), INTENT(INOUT) :: instring @@ -1802,14 +1885,29 @@ SUBROUTINE stringToArray(instring, n, outarray) pos2 = INDEX(instring(pos1:), ",") ! find index of next comma IF (pos2 == 0) THEN ! if there isn't another comma, read the last entry and call it done (this could be the only entry if no commas) n = n + 1 - READ(instring(pos1:), *) outarray(n) + READ(instring(pos1:), *, IOSTAT=ErrStat2) outarray(n) + IF (ErrStat2 /= ErrID_None) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in stringToArray: invalid value in string" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Invalid value in string", ErrStat, ErrMsg, RoutineName); RETURN + ENDIF EXIT END IF n = n + 1 - if (n > 100) then - call WrScr( "ERROR - stringToArray cannot do more than 100 entries") - end if - READ(instring(pos1:pos1+pos2-2), *) outarray(n) + IF (n > 100) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in stringToArray: cannot do more than 100 entries" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Cannot do more than 100 entries", ErrStat, ErrMsg, RoutineName); RETURN + END IF + READ(instring(pos1:pos1+pos2-2), *, IOSTAT=ErrStat2) outarray(n) + IF (ErrStat2 /= ErrID_None) THEN + IF (p%writeLog > 0) THEN + WRITE(p%UnLog, '(A)' ) " ERROR in stringToArray: invalid value in string" + ENDIF + CALL SetErrStat(ErrID_Fatal, "Invalid value in string", ErrStat, ErrMsg, RoutineName); RETURN + ENDIF pos1 = pos2+pos1 END DO @@ -1821,45 +1919,57 @@ END SUBROUTINE stringToArray SUBROUTINE allocateKinematicsArrays() ! error check print *, "Error in Waves::makeGrid, a time or space array is size zero." << endl; - ALLOCATE ( p%uxWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%uyWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%uzWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%axWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%ayWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%azWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%PDyn ( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStatTmp) - ALLOCATE ( p%zeta ( p%ntWave,p%nyWave,p%nxWave), STAT = ErrStatTmp ) ! 2D grid over x and y only + ALLOCATE ( p%uxWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate uxWave'; IF (Failed0()) RETURN + ALLOCATE ( p%uyWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate uyWave'; IF (Failed0()) RETURN + ALLOCATE ( p%uzWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate uzWave'; IF (Failed0()) RETURN + ALLOCATE ( p%axWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate axWave'; IF (Failed0()) RETURN + ALLOCATE ( p%ayWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate ayWave'; IF (Failed0()) RETURN + ALLOCATE ( p%azWave( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate azWave'; IF (Failed0()) RETURN + ALLOCATE ( p%PDyn ( p%ntWave,p%nzWave,p%nyWave,p%nxWave), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate PDyn' ; IF (Failed0()) RETURN + ALLOCATE ( p%zeta ( p%ntWave,p%nyWave,p%nxWave) , STAT=ErrStat2); ErrMsg2 = 'Cannot allocate zeta' ; IF (Failed0()) RETURN ! 2D grid over x and y only END SUBROUTINE allocateKinematicsArrays ! compact way to set the right error status and check if an abort is needed (and do cleanup if so) LOGICAL FUNCTION Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SetupWaterKin') + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'SetupWaterKin') Failed = ErrStat >= AbortErrLev - if (Failed) call CleanUp() + IF (Failed) CALL CleanUp() + CALL SetErrStat(ErrStat2, 'Cleanup was not succcessful', ErrStat, ErrMsg, 'SetupWaterKin') ! check that deallocation was successful END FUNCTION Failed + + ! check for failed where /= 0 is fatal + LOGICAL FUNCTION Failed0() + IF (ErrStat2 /= ErrID_None) THEN + ErrStat2 = ErrID_Fatal + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + ENDIF + Failed0 = ErrStat >= AbortErrLev + IF (Failed0) CALL CleanUp() + CALL SetErrStat(ErrStat2, 'Cleanup was not succcessful', ErrStat, ErrMsg, 'SetupWaterKin') ! check that deallocation was successful + END FUNCTION Failed0 SUBROUTINE CleanUp - !IF (ALLOCATED( WaveElev )) DEALLOCATE( WaveElev, STAT=ErrStatTmp) - !IF (ALLOCATED( WaveTime )) DEALLOCATE( WaveTime, STAT=ErrStatTmp) - IF (ALLOCATED( TmpFFTWaveElev )) DEALLOCATE( TmpFFTWaveElev, STAT=ErrStatTmp) - IF (ALLOCATED( WaveElevC0 )) DEALLOCATE( WaveElevC0, STAT=ErrStatTmp) + !IF (ALLOCATED( WaveElev )) DEALLOCATE( WaveElev, STAT=ErrStat2) + !IF (ALLOCATED( WaveTime )) DEALLOCATE( WaveTime, STAT=ErrStat2) + IF (ALLOCATED( TmpFFTWaveElev )) DEALLOCATE( TmpFFTWaveElev, STAT=ErrStat2) + IF (ALLOCATED( WaveElevC0 )) DEALLOCATE( WaveElevC0, STAT=ErrStat2) ! >>> missing some things <<< - IF (ALLOCATED( WaveNmbr )) DEALLOCATE( WaveNmbr , STAT=ErrStatTmp) - IF (ALLOCATED( tmpComplex )) DEALLOCATE( tmpComplex , STAT=ErrStatTmp) - IF (ALLOCATED( WaveElevC )) DEALLOCATE( WaveElevC , STAT=ErrStatTmp) - IF (ALLOCATED( WaveDynPC )) DEALLOCATE( WaveDynPC , STAT=ErrStatTmp) - IF (ALLOCATED( WaveVelCHx )) DEALLOCATE( WaveVelCHx , STAT=ErrStatTmp) - IF (ALLOCATED( WaveVelCHy )) DEALLOCATE( WaveVelCHy , STAT=ErrStatTmp) - IF (ALLOCATED( WaveVelCV )) DEALLOCATE( WaveVelCV , STAT=ErrStatTmp) - IF (ALLOCATED( WaveAccCHx )) DEALLOCATE( WaveAccCHx , STAT=ErrStatTmp) - IF (ALLOCATED( WaveAccCHy )) DEALLOCATE( WaveAccCHy , STAT=ErrStatTmp) - IF (ALLOCATED( WaveAccCV )) DEALLOCATE( WaveAccCV , STAT=ErrStatTmp) + IF (ALLOCATED( WaveNmbr )) DEALLOCATE( WaveNmbr , STAT=ErrStat2) + IF (ALLOCATED( tmpComplex )) DEALLOCATE( tmpComplex , STAT=ErrStat2) + IF (ALLOCATED( WaveElevC )) DEALLOCATE( WaveElevC , STAT=ErrStat2) + IF (ALLOCATED( WaveDynPC )) DEALLOCATE( WaveDynPC , STAT=ErrStat2) + IF (ALLOCATED( WaveVelCHx )) DEALLOCATE( WaveVelCHx , STAT=ErrStat2) + IF (ALLOCATED( WaveVelCHy )) DEALLOCATE( WaveVelCHy , STAT=ErrStat2) + IF (ALLOCATED( WaveVelCV )) DEALLOCATE( WaveVelCV , STAT=ErrStat2) + IF (ALLOCATED( WaveAccCHx )) DEALLOCATE( WaveAccCHx , STAT=ErrStat2) + IF (ALLOCATED( WaveAccCHy )) DEALLOCATE( WaveAccCHy , STAT=ErrStat2) + IF (ALLOCATED( WaveAccCV )) DEALLOCATE( WaveAccCV , STAT=ErrStat2) END SUBROUTINE CleanUp @@ -2133,9 +2243,188 @@ FUNCTION SINHNumOvrSINHDen ( k, h, z ) END FUNCTION SINHNumOvrSINHDen END SUBROUTINE setupWaterKin + + ! Unused water kinematics file writing subroutines: + + ! ----- write wave grid spacing to output file ----- + SUBROUTINE WriteWaveGrid(p, ErrStat, ErrMsg) + TYPE(MD_ParameterType), INTENT(INOUT) :: p ! Parameters + + INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None + + INTEGER(IntKi) :: ErrStat2 + ! CHARACTER(ErrMsgLen) :: ErrMsg2 - - + CHARACTER(120) :: Frmt + INTEGER(IntKi) :: UnOut ! for outputing wave kinematics data + INTEGER(IntKi) :: I + + + CALL GetNewUnit( UnOut) + + CALL OpenFOutFile ( UnOut, "waves.txt", ErrStat, ErrMsg ) + IF ( ErrStat > ErrID_None ) THEN + ErrMsg = ' Error opening wave grid file: '//TRIM(ErrMsg) + ErrStat = ErrID_Fatal + RETURN + END IF + + WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( 'MoorDyn v2 wave/current kinematics grid file' ) + WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( '---------------------------------------------' ) + WRITE(UnOut, *, IOSTAT=ErrStat2) TRIM( 'The following 6 lines (4-9) specify the input type then the inputs for x, then, y, then z coordinates.' ) + + WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - X input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) + Frmt = '('//TRIM(Int2LStr(5))//'(A1,e10.4))' + WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pxWave(I))), I=1,p%nxWave ) + + WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - Y input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) + Frmt = '('//TRIM(Int2LStr(5))//'(A1,e10.4))' + WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pyWave(I))), I=1,p%nyWave ) + + WRITE(UnOut,*, IOSTAT=ErrStat2) TRIM( '1 - Z input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)' ) + Frmt = '('//TRIM(Int2LStr(8))//'(A1,e10.4))' + WRITE(UnOut,*, IOSTAT=ErrStat2) ( " ", TRIM(Num2LStr(p%pzWave(I))), I=1,p%nzWave ) + + CLOSE(UnOut, IOSTAT = ErrStat2 ) + IF ( ErrStat2 /= 0 ) THEN + ErrStat = ErrID_Severe + ErrMsg = 'Error closing wave grid file' + END IF + + END SUBROUTINE WriteWaveGrid + + + ! ----- write wave kinematics grid data to output file ----- + SUBROUTINE WriteWaveData(p, ErrStat, ErrMsg) + + TYPE(MD_ParameterType), INTENT(INOUT) :: p ! Parameters + + INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None + + INTEGER(IntKi) :: ErrStat2 + ! CHARACTER(ErrMsgLen) :: ErrMsg2 + + INTEGER(IntKi) :: UnOut ! for outputing wave kinematics data + INTEGER(IntKi) :: I,J,K, l, Itemp + + CALL GetNewUnit( UnOut) + + CALL OpenFOutFile ( UnOut, "wave data.txt", ErrStat, ErrMsg ) + IF ( ErrStat > ErrID_None ) THEN + ErrMsg = ' Error opening wave grid file: '//TRIM(ErrMsg) + ErrStat = ErrID_Fatal + RETURN + END IF + + ! write channel labels + + + ! time + WRITE(UnOut,"(A10)", IOSTAT=ErrStat2, advance="no") "Time" + + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ze0", Num2Lstr(J+10*K) + END DO + END DO + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ux", Num2Lstr(I+10*J+100*K) + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " uy", Num2Lstr(I+10*J+100*K) + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " uz", Num2Lstr(I+10*J+100*K) + END DO + END DO + END DO + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ax", Num2Lstr(I+10*J+100*K) + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " ay", Num2Lstr(I+10*J+100*K) + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " az", Num2Lstr(I+10*J+100*K) + END DO + END DO + END DO + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + WRITE(UnOut,"(A3,A8)", IOSTAT=ErrStat2, advance="no") " pd", Num2Lstr(I+10*J+100*K) + END DO + END DO + END DO + + ! end the line + WRITE(UnOut, "(A1)", IOSTAT=ErrStat2) " " + + + + DO l=1, p%ntWave ! loop through all time steps + + ! time + WRITE(UnOut,"(F10.4)", IOSTAT=ErrStat2, advance="no") p%dtWave*(l-1) + !WRITE(UnOut,"(F10.4)", IOSTAT=ErrStat2, advance="no") InitInp%WaveTime(l) + + ! wave elevation (all slices for now, to check) + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + Itemp = (J-1)*p%nxWave + K ! index of actual node + + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%zeta(l,J,K) + END DO + END DO + + ! wave velocities + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node + + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uxWave(l,I,J,K) + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uyWave(l,I,J,K) + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%uzWave(l,I,J,K) + END DO + END DO + END DO + + ! wave accelerations + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node + + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%axWave(l,I,J,K) + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%ayWave(l,I,J,K) + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%azWave(l,I,J,K) + END DO + END DO + END DO + + ! dynamic pressure + DO I=1,p%nzWave !z + DO J = 1,p%nyWave !y + DO K = 1,p%nxWave !x + Itemp = (I-1)*p%nxWave*p%nyWave + (J-1)*p%nxWave + K ! index of actual node + + WRITE(UnOut,"(A1,e10.3)", IOSTAT=ErrStat2, advance="no") " ", p%PDyn(l,I,J,K) + END DO + END DO + END DO + + ! end the line + WRITE(UnOut, "(A1)", IOSTAT=ErrStat2) " " + + + END DO + + + CLOSE(UnOut, IOSTAT = ErrStat ) + IF ( ErrStat /= 0 ) THEN + ErrMsg = 'Error closing wave grid file' + END IF + + END SUBROUTINE WriteWaveData END MODULE MoorDyn_Misc diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 60a4ad7c8e..48b39a8c41 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -11,7 +11,10 @@ # column 10 Units # Keyword ModuleName/ModName Derived data type Field type Variable name variable dimension Initial value for mix language, not used Description Units +# ...... Include files (definitions from NWTC Library) ............................................................................ +# make sure that the file name does not have any trailing white spaces! include Registry_NWTC_Library.txt +usefrom SeaSt_WaveField.txt ## ====== some data read from the input file, but not needed after init ====== typedef MoorDyn/MD MD_InputFileType DbKi DTIC - 0.5 - "convergence check time step for IC generation" "[s]" @@ -20,26 +23,22 @@ typedef ^ ^ ReKi CdScaleIC - 1 typedef ^ ^ ReKi threshIC - 0.01 - "convergence tolerance for ICs (0.01 means 1%)" "[]" ## ============================== Define initialization input types here: ============================================================================================================================= -typedef MoorDyn/MD InitInputType ReKi g - -999.9 - "gravity constant" "[m/s^2]" -typedef ^ ^ ReKi rhoW - -999.9 - "sea density" "[kg/m^3]" -typedef ^ ^ ReKi WtrDepth - -999.9 - "depth of water" "[m]" -typedef ^ ^ ReKi PtfmInit {:}{:} - - "initial position of platform(s) shape: 6, nTurbines" - -typedef ^ ^ IntKi FarmSize - 0 - "Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >0" - -typedef ^ ^ ReKi TurbineRefPos {:}{:} - - "reference position of turbines in farm, shape: 3, nTurbines" - -typedef ^ ^ ReKi Tmax - - - "simulation duration" "[s]" -typedef ^ ^ CHARACTER(1024) FileName - "" - "MoorDyn input file" -typedef ^ ^ CHARACTER(1024) RootName - - - "RootName for writing output files" - -typedef ^ ^ LOGICAL UsePrimaryInputFile - .TRUE. - "Read input file instead of passed data" - -typedef ^ ^ FileInfoType PassedPrimaryInputData - - - "Primary input file as FileInfoType (set by driver/glue code) -- String array with metadata" - -typedef ^ ^ LOGICAL Echo - "" - "echo parameter - do we want to echo the header line describing the input file?" -typedef ^ ^ CHARACTER(ChanLen) OutList {:} "" - "string containing list of output channels requested in input file" -typedef ^ ^ Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." - -typedef ^ ^ Logical VisMeshes - .FALSE. - "Glue code requesting visualization meshes" - - -#typedef ^ ^ DbKi UGrid {:}{:}{:} - - "water velocities time series at each grid point" - -#typedef ^ ^ DbKi UdGrid {:}{:}{:} - - "water accelerations time series at each grid point" - -#typedef ^ ^ DbKi zetaGrid {:}{:} - - "water surface elevations time series at each grid point" - -#typedef ^ ^ DbKi PDynGrid {:}{:} - - "water dynamic pressure time series at each grid point" - +typedef MoorDyn/MD InitInputType ReKi g - -999.9 - "gravity constant" "[m/s^2]" +typedef ^ ^ ReKi rhoW - -999.9 - "sea density" "[kg/m^3]" +typedef ^ ^ ReKi WtrDepth - -999.9 - "depth of water" "[m]" +typedef ^ ^ ReKi PtfmInit {:}{:} - - "initial position of platform(s) shape: 6, nTurbines" - +typedef ^ ^ IntKi FarmSize - 0 - "Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >1" - +typedef ^ ^ ReKi TurbineRefPos {:}{:} - - "reference position of turbines in farm, shape: 3, nTurbines" - +typedef ^ ^ ReKi Tmax - - - "simulation duration" "[s]" +typedef ^ ^ CHARACTER(1024) FileName - "" - "MoorDyn input file" +typedef ^ ^ CHARACTER(1024) RootName - - - "RootName for writing output files" - +typedef ^ ^ LOGICAL UsePrimaryInputFile - .TRUE. - "Read input file instead of passed data" - +typedef ^ ^ FileInfoType PassedPrimaryInputData - - - "Primary input file as FileInfoType (set by driver/glue code) -- String array with metadata" - +typedef ^ ^ LOGICAL Echo - "" - "echo parameter - do we want to echo the header line describing the input file?" +typedef ^ ^ CHARACTER(ChanLen) OutList {:} "" - "string containing list of output channels requested in input file" +typedef ^ ^ Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." - +typedef ^ ^ Logical VisMeshes - .FALSE. - "Glue code requesting visualization meshes" - +typedef ^ ^ SeaSt_WaveFieldType *WaveField - - - "Pointer to SeaState wave field" - # nvm # Farm-level simulation inputs - these are passed by FAST.Farm - the arrays are populated from the individual turbine-level MoorDyn instances # nvm typedef ^ ^ MeshType FarmCoupledKinematics {:} - - "array of input kinematics meshes from each of the turbine-level MoorDyn instances" "[m, m/s]" @@ -95,7 +94,7 @@ typedef ^ MD_Body IntKi IdNum - typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=free, 1=fixed, -1=coupled, 2=coupledpinned" typedef ^ ^ IntKi AttachedC {30} - - "list of IdNums of points attached to this body" typedef ^ ^ IntKi AttachedR {30} - - "list of IdNums of rods attached to this body" -typedef ^ ^ IntKi nAttachedC - - - "number of attached points" +typedef ^ ^ IntKi nAttachedP - - - "number of attached points" typedef ^ ^ IntKi nAttachedR - - - "number of attached rods" typedef ^ ^ DbKi rPointRel {3}{30} - - "relative position of point on body" typedef ^ ^ DbKi r6RodRel {6}{30} - - "relative position and orientation of rod on body" @@ -294,7 +293,6 @@ typedef ^ ^ IntKi LineUnOut - typedef ^ ^ DbKi LineWrOutput {:} - - "one row of output data for this line" # NOTE: Below are vars used in the line VIV model typedef ^ ^ DbKi phi {:} - - "phase of lift force" -typedef ^ ^ LOGICAL IC_gen - .FALSE. - "boolean to indicate dynamic relaxation occuring" "-" typedef ^ ^ DbKi n_m - 500 - "Num timesteps for rolling RMS of crossflow velocity phase" "-" typedef ^ ^ DbKi phi_yd - - - "The crossflow motion phase" "-" typedef ^ ^ DbKi t_old - - - "old t" "s" @@ -399,6 +397,8 @@ typedef ^ ^ DbKi BathymetryGrid {:}{:} typedef ^ ^ DbKi BathGrid_Xs {:} - - "array of x-coordinates in the bathymetry grid" typedef ^ ^ DbKi BathGrid_Ys {:} - - "array of y-coordinates in the bathymetry grid" typedef ^ ^ IntKi BathGrid_npoints {:} - - "number of grid points to describe the bathymetry grid" +typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ LOGICAL IC_gen - .FALSE. - "boolean to indicate dynamic relaxation occuring" "-" ## ============================== Parameters ============================================================================================================================================ @@ -439,9 +439,11 @@ typedef ^ ^ CHARACTER(1024) PriPath - typedef ^ ^ IntKi writeLog - -1 - "Switch for level of log file output" #NOTE: there may be an issue with start/restart with the UnLog stored in parameters. We'll ignore this for now -- ADP typedef ^ ^ IntKi UnLog - -1 - "Unit number of log file" -typedef ^ ^ IntKi WaveKin - - - "Flag for whether or how to consider water kinematics" -typedef ^ ^ IntKi Current - - - "Flag for whether or how to consider water kinematics" -typedef ^ ^ IntKi nTurbines - - - "Number of turbines if MoorDyn is performing an array-level simulation with FAST.Farm, otherwise 0" +typedef ^ ^ IntKi WaveKin - 0 - "Flag for whether or how to consider wave kinematics. (0: no waves, 1: old method, 2: hybrid method, 3: SeaState method)" +typedef ^ ^ IntKi Current - 0 - "Flag for whether or how to consider currents. (0: no currents, 1: old method, 2: hybrid method, 3: SeaState method)" +typedef ^ ^ IntKi WaterKin - 0 - "Flag for whether or how to consider water kinematics. (0: no kinematics, 1: old method, 2: hybrid method, 3: SeaState method)" +typedef ^ ^ SeaSt_WaveFieldType *WaveField - - - "Pointer to SeaState wave field" - +typedef ^ ^ IntKi nTurbines - - - "Number of turbines if MoorDyn is performing an array-level simulation with FAST.Farm, otherwise 1" typedef ^ ^ ReKi TurbineRefPos {:}{:} - - "reference position of turbines in farm, shape: 3, nTurbines" - typedef ^ ^ DbKi mu_kT - - - "transverse kinetic friction coefficient" "(-)" typedef ^ ^ DbKi mu_kA - - - "axial kinetic friction coefficient" "(-)" diff --git a/modules/moordyn/src/MoorDyn_Rod.f90 b/modules/moordyn/src/MoorDyn_Rod.f90 index c01c5b6ae5..fda235bedf 100644 --- a/modules/moordyn/src/MoorDyn_Rod.f90 +++ b/modules/moordyn/src/MoorDyn_Rod.f90 @@ -88,7 +88,7 @@ SUBROUTINE Rod_Setup(Rod, RodProp, endCoords, p, ErrStat, ErrMsg) ! allocate segment scalar quantities if (Rod%N == 0) then ! special case of zero-length Rod - ALLOCATE(Rod%l(1), Rod%V(N), STAT=ErrStat2); if(AllocateFailed("Rod: l and V")) return + ALLOCATE(Rod%l(1), Rod%V(1), STAT=ErrStat2); if(AllocateFailed("Rod: l and V")) return else ! normal case ALLOCATE(Rod%l(N), Rod%V(N), STAT=ErrStat2); if(AllocateFailed("Rod: l and V")) return end if @@ -148,7 +148,7 @@ SUBROUTINE Rod_Setup(Rod, RodProp, endCoords, p, ErrStat, ErrMsg) Rod%OrMat = CalcOrientation(phi, beta, 0.0_DbKi) ! get rotation matrix to put things in global rather than rod-axis orientations - IF (wordy > 0) print *, "Set up Rod ",Rod%IdNum, ", type ", Rod%typeNum + IF (wordy > 0) CALL WrScr("Set up Rod "//trim(num2lstr(Rod%IdNum))//", type "//trim(num2lstr(Rod%typeNum))) ! need to add cleanup sub <<< @@ -184,7 +184,7 @@ SUBROUTINE Rod_Initialize(Rod, states, m) ! REAL(DbKi) :: rRef(3) ! reference position of mesh node ! REAL(DbKi) :: OrMat(3,3) ! DCM for body orientation based on r6_in - IF (wordy > 0) print *, "initializing Rod ", Rod%idNum + IF (wordy > 0) CALL WrScr("initializing Rod "//trim(num2lstr(Rod%idNum))) ! the r6 and v6 vectors should have already been set ! r and rd of ends have already been set by setup function or by parent object <<<<< right? <<<<< @@ -587,6 +587,9 @@ SUBROUTINE Rod_DoRHS(Rod, m, p) Real(DbKi) :: depth ! local interpolated depth from bathymetry grid [m] Real(DbKi) :: nvec(3) ! local seabed surface normal vector (positive out) + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + N = Rod%N @@ -619,7 +622,9 @@ SUBROUTINE Rod_DoRHS(Rod, m, p) ! apply wave kinematics (if there are any) DO i=0,N - CALL getWaterKin(p, Rod%r(1,i), Rod%r(2,i), Rod%r(3,i), Rod%time, m%WaveTi, Rod%U(:,i), Rod%Ud(:,i), Rod%zeta(i), Rod%PDyn(i)) + CALL getWaterKin(p, m, Rod%r(1,i), Rod%r(2,i), Rod%r(3,i), Rod%time, Rod%U(:,i), Rod%Ud(:,i), Rod%zeta(i), Rod%PDyn(i), ErrStat2, ErrMsg2) + ! TODO: Handle error messages. Roads broadly needs error flags to be supported + !F(i) = 1.0 ! set VOF value to one for now (everything submerged - eventually this should be element-based!!!) <<<< ! <<<< currently F is not being used and instead a VOF variable is used within the node loop END DO @@ -634,7 +639,7 @@ SUBROUTINE Rod_DoRHS(Rod, m, p) else if ((Rod%r(3,N) < zeta) .and. (Rod%r(3,0) > zeta)) then ! check if it's crossing the water plane but upside down Rod%h0 = -(zeta - Rod%r(3,0))/Rod%q(3) ! negative distance along rod centerline from end A to the waterplane else - Rod%h0 = 0.0_DbKi ! fully unsubmerged case (ever applicable?) + Rod%h0 = 0.0_DbKi ! fully unsubmerged case (ever applicable?). (Rod%r(3,0) > zeta) .and. (Rod%r(3,N) > zeta) end if @@ -879,6 +884,22 @@ SUBROUTINE Rod_DoRHS(Rod, m, p) Rod%Fnet(:,I) = Rod%W(:,I) + Rod%Bo(:,I) + Rod%Dp(:,I) + Rod%Dq(:,I) & + Rod%Ap(:,I) + Rod%Aq(:,I) + Rod%Pd(:,I) + Rod%B(:,I) + Rod%Bp(:,I) + Rod%Bq(:,I) + DO J=1,3 + IF (Is_NaN(Rod%Fnet(J,I)) .AND. wordy>1) THEN + CALL WrScr("NaN detected in Rod%Fnet at node "//trim(num2lstr(I))) + CALL WrScr("Rod%IDNum: "//trim(Int2LStr(Rod%IdNum))) + CALL WrScr("Rod%W("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%W(J,I)))) + CALL WrScr("Rod%Bo("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Bo(J,I)))) + CALL WrScr("Rod%Dp("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Dp(J,I)))) + CALL WrScr("Rod%Dq("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Dq(J,I)))) + CALL WrScr("Rod%Ap("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Ap(J,I)))) + CALL WrScr("Rod%Aq("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Aq(J,I)))) + CALL WrScr("Rod%Pd("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Pd(J,I)))) + CALL WrScr("Rod%B("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%B(J,I)))) + CALL WrScr("Rod%Bp("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Bp(J,I)))) + CALL WrScr("Rod%Bq("//trim(num2lstr(J))//","//trim(num2lstr(I))//"): "//trim(num2lstr(Rod%Bq(J,I)))) + END IF + ENDDO END DO ! I - done looping through nodes @@ -985,8 +1006,25 @@ SUBROUTINE Rod_DoRHS(Rod, m, p) Rod%F6net(1:3) = Rod%F6net(1:3) + Fcentripetal Rod%F6net(4:6) = Rod%F6net(4:6) + Mcentripetal + Rod%Mext + DO J=1,3 + IF ((Is_NaN(Rod%F6net(J)) .OR. Is_NaN(Rod%F6net(J+3))) .AND. wordy>1) THEN ! convoluted logic to check 1:6 with a 1:3 loop + CALL WrScr("NaN detected in Rod%F6net("//trim(num2lstr(J))//") after adding centripetal force/moment") + CALL WrScr("Rod%IDNum: "//trim(num2lstr(Rod%IdNum))) + CALL WrScr("Fcentripetal("//trim(num2lstr(J))//"): "//trim(num2lstr(Fcentripetal(J)))) + CALL WrScr("Mcentripetal("//trim(num2lstr(J))//"): "//trim(num2lstr(Mcentripetal(J)))) + CALL WrScr("Rod%Mext("//trim(num2lstr(J))//"): "//trim(num2lstr(Rod%Mext(J)))) + END IF + ENDDO + ! add in user defined external forces to end A Rod%F6net(1:3) = Rod%F6net(1:3) + Rod%FextU + DO J=1,3 + IF (Is_NaN(Rod%F6net(J)) .AND. wordy>1) THEN + CALL WrScr("NaN detected in Rod%F6net("//trim(num2lstr(J))//") after adding Rod%FextU") + CALL WrScr("Rod%IDNum: "//trim(num2lstr(Rod%IdNum))) + CALL WrScr("Rod%FextU("//trim(num2lstr(J))//"): "//trim(num2lstr(Rod%FextU(J)))) + END IF + ENDDO ! Note: F6net saves the Rod's net forces and moments (excluding inertial ones) for use in later output ! (this is what the rod will apply to whatever it's attached to, so should be zero moments if pinned). @@ -1119,7 +1157,7 @@ SUBROUTINE Rod_AddLine(Rod, lineID, TopOfLine, endB) if (endB==1) then ! attaching to end B - IF (wordy > 0) Print*, "L", lineID, "->R", Rod%IdNum , "b" + IF (wordy > 0) CALL WrScr("L"//trim(num2lstr(lineID))//"->R"//trim(num2lstr(Rod%IdNum))//"b") IF (Rod%nAttachedB <10) THEN ! this is currently just a maximum imposed by a fixed array size. could be improved. Rod%nAttachedB = Rod%nAttachedB + 1 ! add the line to the number connected @@ -1131,7 +1169,7 @@ SUBROUTINE Rod_AddLine(Rod, lineID, TopOfLine, endB) else ! attaching to end A - IF (wordy > 0) Print*, "L", lineID, "->R", Rod%IdNum , "a" + IF (wordy > 0) CALL WrScr("L"//trim(num2lstr(lineID))//"->R"//trim(num2lstr(Rod%IdNum))//"a") IF (Rod%nAttachedA <10) THEN ! this is currently just a maximum imposed by a fixed array size. could be improved. Rod%nAttachedA = Rod%nAttachedA + 1 ! add the line to the number connected diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index cd85475cf6..8278b06d80 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -31,6 +31,7 @@ !! unpack routines associated with each defined data type. This code is automatically generated by the FAST Registry. MODULE MoorDyn_Types !--------------------------------------------------------------------------------------------------------------------------------- +USE SeaSt_WaveField_Types USE NWTC_Library IMPLICIT NONE ! ========= MD_InputFileType ======= @@ -47,7 +48,7 @@ MODULE MoorDyn_Types REAL(ReKi) :: rhoW = -999.9 !< sea density [[kg/m^3]] REAL(ReKi) :: WtrDepth = -999.9 !< depth of water [[m]] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: PtfmInit !< initial position of platform(s) shape: 6, nTurbines [-] - INTEGER(IntKi) :: FarmSize = 0 !< Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >0 [-] + INTEGER(IntKi) :: FarmSize = 0 !< Indicates normal FAST module mode if 0, FAST.Farm coupled mode and =nTurbines if >1 [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: TurbineRefPos !< reference position of turbines in farm, shape: 3, nTurbines [-] REAL(ReKi) :: Tmax = 0.0_ReKi !< simulation duration [[s]] CHARACTER(1024) :: FileName !< MoorDyn input file [-] @@ -58,6 +59,7 @@ MODULE MoorDyn_Types CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: OutList !< string containing list of output channels requested in input file [-] LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-] LOGICAL :: VisMeshes = .FALSE. !< Glue code requesting visualization meshes [-] + TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to SeaState wave field [-] END TYPE MD_InitInputType ! ======================= ! ========= MD_LineProp ======= @@ -112,7 +114,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: typeNum = 0_IntKi !< integer identifying the type. 0=free, 1=fixed, -1=coupled, 2=coupledpinned [-] INTEGER(IntKi) , DIMENSION(1:30) :: AttachedC = 0_IntKi !< list of IdNums of points attached to this body [-] INTEGER(IntKi) , DIMENSION(1:30) :: AttachedR = 0_IntKi !< list of IdNums of rods attached to this body [-] - INTEGER(IntKi) :: nAttachedC = 0_IntKi !< number of attached points [-] + INTEGER(IntKi) :: nAttachedP = 0_IntKi !< number of attached points [-] INTEGER(IntKi) :: nAttachedR = 0_IntKi !< number of attached rods [-] REAL(DbKi) , DIMENSION(1:3,1:30) :: rPointRel = 0.0_R8Ki !< relative position of point on body [-] REAL(DbKi) , DIMENSION(1:6,1:30) :: r6RodRel = 0.0_R8Ki !< relative position and orientation of rod on body [-] @@ -313,7 +315,6 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: LineUnOut = 0_IntKi !< unit number of line output file [-] REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: LineWrOutput !< one row of output data for this line [-] REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: phi !< phase of lift force [-] - LOGICAL :: IC_gen = .FALSE. !< boolean to indicate dynamic relaxation occuring [-] REAL(DbKi) :: n_m = 500 !< Num timesteps for rolling RMS of crossflow velocity phase [-] REAL(DbKi) :: phi_yd = 0.0_R8Ki !< The crossflow motion phase [-] REAL(DbKi) :: t_old = 0.0_R8Ki !< old t [s] @@ -434,6 +435,8 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_Xs !< array of x-coordinates in the bathymetry grid [-] REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_Ys !< array of y-coordinates in the bathymetry grid [-] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_npoints !< number of grid points to describe the bathymetry grid [-] + TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + LOGICAL :: IC_gen = .FALSE. !< boolean to indicate dynamic relaxation occuring [-] END TYPE MD_MiscVarType ! ======================= ! ========= MD_ParameterType ======= @@ -474,9 +477,11 @@ MODULE MoorDyn_Types CHARACTER(1024) :: PriPath !< The path to the primary MoorDyn input file, used if looking for additional input files [-] INTEGER(IntKi) :: writeLog = -1 !< Switch for level of log file output [-] INTEGER(IntKi) :: UnLog = -1 !< Unit number of log file [-] - INTEGER(IntKi) :: WaveKin = 0_IntKi !< Flag for whether or how to consider water kinematics [-] - INTEGER(IntKi) :: Current = 0_IntKi !< Flag for whether or how to consider water kinematics [-] - INTEGER(IntKi) :: nTurbines = 0_IntKi !< Number of turbines if MoorDyn is performing an array-level simulation with FAST.Farm, otherwise 0 [-] + INTEGER(IntKi) :: WaveKin = 0 !< Flag for whether or how to consider wave kinematics. (0: no waves, 1: old method, 2: hybrid method, 3: SeaState method) [-] + INTEGER(IntKi) :: Current = 0 !< Flag for whether or how to consider currents. (0: no currents, 1: old method, 2: hybrid method, 3: SeaState method) [-] + INTEGER(IntKi) :: WaterKin = 0 !< Flag for whether or how to consider water kinematics. (0: no kinematics, 1: old method, 2: hybrid method, 3: SeaState method) [-] + TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to SeaState wave field [-] + INTEGER(IntKi) :: nTurbines = 0_IntKi !< Number of turbines if MoorDyn is performing an array-level simulation with FAST.Farm, otherwise 1 [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: TurbineRefPos !< reference position of turbines in farm, shape: 3, nTurbines [-] REAL(DbKi) :: mu_kT = 0.0_R8Ki !< transverse kinetic friction coefficient [(-)] REAL(DbKi) :: mu_kA = 0.0_R8Ki !< axial kinetic friction coefficient [(-)] @@ -645,6 +650,7 @@ subroutine MD_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrSta end if DstInitInputData%Linearize = SrcInitInputData%Linearize DstInitInputData%VisMeshes = SrcInitInputData%VisMeshes + DstInitInputData%WaveField => SrcInitInputData%WaveField end subroutine subroutine MD_DestroyInitInput(InitInputData, ErrStat, ErrMsg) @@ -667,12 +673,14 @@ subroutine MD_DestroyInitInput(InitInputData, ErrStat, ErrMsg) if (allocated(InitInputData%OutList)) then deallocate(InitInputData%OutList) end if + nullify(InitInputData%WaveField) end subroutine subroutine MD_PackInitInput(RF, Indata) type(RegFile), intent(inout) :: RF type(MD_InitInputType), intent(in) :: InData character(*), parameter :: RoutineName = 'MD_PackInitInput' + logical :: PtrInIndex if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%g) call RegPack(RF, InData%rhoW) @@ -689,6 +697,13 @@ subroutine MD_PackInitInput(RF, Indata) call RegPackAlloc(RF, InData%OutList) call RegPack(RF, InData%Linearize) call RegPack(RF, InData%VisMeshes) + call RegPack(RF, associated(InData%WaveField)) + if (associated(InData%WaveField)) then + call RegPackPointer(RF, c_loc(InData%WaveField), PtrInIndex) + if (.not. PtrInIndex) then + call SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, InData%WaveField) + end if + end if if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -699,6 +714,8 @@ subroutine MD_UnPackInitInput(RF, OutData) integer(B4Ki) :: LB(2), UB(2) integer(IntKi) :: stat logical :: IsAllocAssoc + integer(B8Ki) :: PtrIdx + type(c_ptr) :: Ptr if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%g); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%rhoW); if (RegCheckErr(RF, RoutineName)) return @@ -715,6 +732,24 @@ subroutine MD_UnPackInitInput(RF, OutData) call RegUnpackAlloc(RF, OutData%OutList); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Linearize); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%VisMeshes); if (RegCheckErr(RF, RoutineName)) return + if (associated(OutData%WaveField)) deallocate(OutData%WaveField) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackPointer(RF, Ptr, PtrIdx); if (RegCheckErr(RF, RoutineName)) return + if (c_associated(Ptr)) then + call c_f_pointer(Ptr, OutData%WaveField) + else + allocate(OutData%WaveField,stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%WaveField.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + RF%Pointers(PtrIdx) = c_loc(OutData%WaveField) + call SeaSt_WaveField_UnpackSeaSt_WaveFieldType(RF, OutData%WaveField) ! WaveField + end if + else + OutData%WaveField => null() + end if end subroutine subroutine MD_CopyLineProp(SrcLinePropData, DstLinePropData, CtrlCode, ErrStat, ErrMsg) @@ -914,7 +949,7 @@ subroutine MD_CopyBody(SrcBodyData, DstBodyData, CtrlCode, ErrStat, ErrMsg) DstBodyData%typeNum = SrcBodyData%typeNum DstBodyData%AttachedC = SrcBodyData%AttachedC DstBodyData%AttachedR = SrcBodyData%AttachedR - DstBodyData%nAttachedC = SrcBodyData%nAttachedC + DstBodyData%nAttachedP = SrcBodyData%nAttachedP DstBodyData%nAttachedR = SrcBodyData%nAttachedR DstBodyData%rPointRel = SrcBodyData%rPointRel DstBodyData%r6RodRel = SrcBodyData%r6RodRel @@ -962,7 +997,7 @@ subroutine MD_PackBody(RF, Indata) call RegPack(RF, InData%typeNum) call RegPack(RF, InData%AttachedC) call RegPack(RF, InData%AttachedR) - call RegPack(RF, InData%nAttachedC) + call RegPack(RF, InData%nAttachedP) call RegPack(RF, InData%nAttachedR) call RegPack(RF, InData%rPointRel) call RegPack(RF, InData%r6RodRel) @@ -1002,7 +1037,7 @@ subroutine MD_UnPackBody(RF, OutData) call RegUnpack(RF, OutData%typeNum); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%AttachedC); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%AttachedR); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%nAttachedC); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%nAttachedP); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nAttachedR); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%rPointRel); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%r6RodRel); if (RegCheckErr(RF, RoutineName)) return @@ -2113,7 +2148,6 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) end if DstLineData%phi = SrcLineData%phi end if - DstLineData%IC_gen = SrcLineData%IC_gen DstLineData%n_m = SrcLineData%n_m DstLineData%phi_yd = SrcLineData%phi_yd DstLineData%t_old = SrcLineData%t_old @@ -2342,7 +2376,6 @@ subroutine MD_PackLine(RF, Indata) call RegPack(RF, InData%LineUnOut) call RegPackAlloc(RF, InData%LineWrOutput) call RegPackAlloc(RF, InData%phi) - call RegPack(RF, InData%IC_gen) call RegPack(RF, InData%n_m) call RegPack(RF, InData%phi_yd) call RegPack(RF, InData%t_old) @@ -2431,7 +2464,6 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpack(RF, OutData%LineUnOut); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%LineWrOutput); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%phi); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%IC_gen); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%n_m); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%phi_yd); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%t_old); if (RegCheckErr(RF, RoutineName)) return @@ -3457,6 +3489,10 @@ subroutine MD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%BathGrid_npoints = SrcMiscData%BathGrid_npoints end if + call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + DstMiscData%IC_gen = SrcMiscData%IC_gen end subroutine subroutine MD_DestroyMisc(MiscData, ErrStat, ErrMsg) @@ -3607,6 +3643,8 @@ subroutine MD_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%BathGrid_npoints)) then deallocate(MiscData%BathGrid_npoints) end if + call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine MD_PackMisc(RF, Indata) @@ -3717,6 +3755,8 @@ subroutine MD_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%BathGrid_Xs) call RegPackAlloc(RF, InData%BathGrid_Ys) call RegPackAlloc(RF, InData%BathGrid_npoints) + call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call RegPack(RF, InData%IC_gen) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3862,6 +3902,8 @@ subroutine MD_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%BathGrid_Xs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BathGrid_Ys); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BathGrid_npoints); if (RegCheckErr(RF, RoutineName)) return + call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call RegUnpack(RF, OutData%IC_gen); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) @@ -3963,6 +4005,8 @@ subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%UnLog = SrcParamData%UnLog DstParamData%WaveKin = SrcParamData%WaveKin DstParamData%Current = SrcParamData%Current + DstParamData%WaterKin = SrcParamData%WaterKin + DstParamData%WaveField => SrcParamData%WaveField DstParamData%nTurbines = SrcParamData%nTurbines if (allocated(SrcParamData%TurbineRefPos)) then LB(1:2) = lbound(SrcParamData%TurbineRefPos) @@ -4257,6 +4301,7 @@ subroutine MD_DestroyParam(ParamData, ErrStat, ErrMsg) end do deallocate(ParamData%OutParam) end if + nullify(ParamData%WaveField) if (allocated(ParamData%TurbineRefPos)) then deallocate(ParamData%TurbineRefPos) end if @@ -4331,6 +4376,7 @@ subroutine MD_PackParam(RF, Indata) character(*), parameter :: RoutineName = 'MD_PackParam' integer(B4Ki) :: i1, i2, i3, i4 integer(B4Ki) :: LB(4), UB(4) + logical :: PtrInIndex if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%nLineTypes) call RegPack(RF, InData%nRodTypes) @@ -4378,6 +4424,14 @@ subroutine MD_PackParam(RF, Indata) call RegPack(RF, InData%UnLog) call RegPack(RF, InData%WaveKin) call RegPack(RF, InData%Current) + call RegPack(RF, InData%WaterKin) + call RegPack(RF, associated(InData%WaveField)) + if (associated(InData%WaveField)) then + call RegPackPointer(RF, c_loc(InData%WaveField), PtrInIndex) + if (.not. PtrInIndex) then + call SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, InData%WaveField) + end if + end if call RegPack(RF, InData%nTurbines) call RegPackAlloc(RF, InData%TurbineRefPos) call RegPack(RF, InData%mu_kT) @@ -4436,6 +4490,8 @@ subroutine MD_UnPackParam(RF, OutData) integer(B4Ki) :: LB(4), UB(4) integer(IntKi) :: stat logical :: IsAllocAssoc + integer(B8Ki) :: PtrIdx + type(c_ptr) :: Ptr if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%nLineTypes); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nRodTypes); if (RegCheckErr(RF, RoutineName)) return @@ -4487,6 +4543,25 @@ subroutine MD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%UnLog); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WaveKin); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Current); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%WaterKin); if (RegCheckErr(RF, RoutineName)) return + if (associated(OutData%WaveField)) deallocate(OutData%WaveField) + call RegUnpack(RF, IsAllocAssoc); if (RegCheckErr(RF, RoutineName)) return + if (IsAllocAssoc) then + call RegUnpackPointer(RF, Ptr, PtrIdx); if (RegCheckErr(RF, RoutineName)) return + if (c_associated(Ptr)) then + call c_f_pointer(Ptr, OutData%WaveField) + else + allocate(OutData%WaveField,stat=stat) + if (stat /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating OutData%WaveField.', RF%ErrStat, RF%ErrMsg, RoutineName) + return + end if + RF%Pointers(PtrIdx) = c_loc(OutData%WaveField) + call SeaSt_WaveField_UnpackSeaSt_WaveFieldType(RF, OutData%WaveField) ! WaveField + end if + else + OutData%WaveField => null() + end if call RegUnpack(RF, OutData%nTurbines); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%TurbineRefPos); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%mu_kT); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90 index 9d9ec54e55..d5cdfc82a5 100644 --- a/modules/openfast-library/src/FAST_Subs.f90 +++ b/modules/openfast-library/src/FAST_Subs.f90 @@ -1197,6 +1197,9 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, SED, BD, S Init%InData_MD%Linearize = p_FAST%Linearize if (p_FAST%WrVTK /= VTK_None) Init%InData_MD%VisMeshes=.true. + ! Assign the seastate pointer here + Init%InData_MD%WaveField => Init%OutData_SeaSt%WaveField + CALL MD_Init( Init%InData_MD, MD%Input(1), MD%p, MD%x(STATE_CURR), MD%xd(STATE_CURR), MD%z(STATE_CURR), & MD%OtherSt(STATE_CURR), MD%y, MD%m, p_FAST%dt_module( MODULE_MD ), Init%OutData_MD, ErrStat2, ErrMsg2 ) CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) diff --git a/modules/seastate/src/SeaSt_WaveField.txt b/modules/seastate/src/SeaSt_WaveField.txt index f5730bc868..4dc9ba0732 100644 --- a/modules/seastate/src/SeaSt_WaveField.txt +++ b/modules/seastate/src/SeaSt_WaveField.txt @@ -1,3 +1,5 @@ +# ...... Include files ..... +usefrom Current.txt #--------------------------------------------------------------------------------------------------------------------------------------------------------- # Data structures for representing wave fields. # @@ -72,3 +74,4 @@ typedef ^ ^ INTEGER WaveMod typedef ^ ^ INTEGER NStepWave - - - "Total number of frequency components = total number of time steps in the incident wave" - typedef ^ ^ INTEGER NStepWave2 - - - "NStepWave / 2" - +typedef ^ ^ Current_InitInputType Current_InitInput - - - "InitInputs in the Current Module. For coupling with MD." - diff --git a/modules/seastate/src/SeaSt_WaveField_Types.f90 b/modules/seastate/src/SeaSt_WaveField_Types.f90 index 4654a04040..6db5590929 100644 --- a/modules/seastate/src/SeaSt_WaveField_Types.f90 +++ b/modules/seastate/src/SeaSt_WaveField_Types.f90 @@ -31,6 +31,7 @@ !! unpack routines associated with each defined data type. This code is automatically generated by the FAST Registry. MODULE SeaSt_WaveField_Types !--------------------------------------------------------------------------------------------------------------------------------- +USE Current_Types USE NWTC_Library IMPLICIT NONE INTEGER(IntKi), PUBLIC, PARAMETER :: WaveDirMod_None = 0 ! WaveDirMod = 0 [Directional spreading function is NONE] [-] @@ -103,6 +104,7 @@ MODULE SeaSt_WaveField_Types INTEGER(IntKi) :: WaveMod = 0_IntKi !< Incident wave kinematics model: See valid values in SeaSt_WaveField module parameters. [-] INTEGER(IntKi) :: NStepWave = 0_IntKi !< Total number of frequency components = total number of time steps in the incident wave [-] INTEGER(IntKi) :: NStepWave2 = 0_IntKi !< NStepWave / 2 [-] + TYPE(Current_InitInputType) :: Current_InitInput !< InitInputs in the Current Module. For coupling with MD. [-] END TYPE SeaSt_WaveFieldType ! ======================= CONTAINS @@ -420,6 +422,9 @@ subroutine SeaSt_WaveField_CopySeaSt_WaveFieldType(SrcSeaSt_WaveFieldTypeData, D DstSeaSt_WaveFieldTypeData%WaveMod = SrcSeaSt_WaveFieldTypeData%WaveMod DstSeaSt_WaveFieldTypeData%NStepWave = SrcSeaSt_WaveFieldTypeData%NStepWave DstSeaSt_WaveFieldTypeData%NStepWave2 = SrcSeaSt_WaveFieldTypeData%NStepWave2 + call Current_CopyInitInput(SrcSeaSt_WaveFieldTypeData%Current_InitInput, DstSeaSt_WaveFieldTypeData%Current_InitInput, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return end subroutine subroutine SeaSt_WaveField_DestroySeaSt_WaveFieldType(SeaSt_WaveFieldTypeData, ErrStat, ErrMsg) @@ -478,6 +483,8 @@ subroutine SeaSt_WaveField_DestroySeaSt_WaveFieldType(SeaSt_WaveFieldTypeData, E if (allocated(SeaSt_WaveFieldTypeData%WaveDirArr)) then deallocate(SeaSt_WaveFieldTypeData%WaveDirArr) end if + call Current_DestroyInitInput(SeaSt_WaveFieldTypeData%Current_InitInput, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine subroutine SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, Indata) @@ -522,6 +529,7 @@ subroutine SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, Indata) call RegPack(RF, InData%WaveMod) call RegPack(RF, InData%NStepWave) call RegPack(RF, InData%NStepWave2) + call Current_PackInitInput(RF, InData%Current_InitInput) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -570,6 +578,7 @@ subroutine SeaSt_WaveField_UnPackSeaSt_WaveFieldType(RF, OutData) call RegUnpack(RF, OutData%WaveMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NStepWave); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NStepWave2); if (RegCheckErr(RF, RoutineName)) return + call Current_UnpackInitInput(RF, OutData%Current_InitInput) ! Current_InitInput end subroutine END MODULE SeaSt_WaveField_Types !ENDOFREGISTRYGENERATEDFILE diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90 index 10dca7ab7d..af5d0b15bf 100644 --- a/modules/seastate/src/SeaState.f90 +++ b/modules/seastate/src/SeaState.f90 @@ -153,6 +153,7 @@ SUBROUTINE SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, Init ! Initialize Current module CALL Current_Init(InputFileData%Current, Current_InitOut, ErrStat2, ErrMsg2 ); if(Failed()) return; + p%WaveField%Current_InitInput = InputFileData%Current ! Save the current input data for later use by MD ! Move initialization output data from Current module into the initialization input data for the Waves module IF (ALLOCATED(Current_InitOut%CurrVxi)) CALL Move_Alloc( Current_InitOut%CurrVxi, InputFileData%Waves%CurrVxi ) @@ -275,7 +276,7 @@ SUBROUTINE SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, Init DEALLOCATE ( InitOut%WriteOutputHdr ) END IF - InitOut%WaveField => p%WaveField + InitOut%WaveField => p%WaveField ! Tell HydroDyn if state-space wave excitation is not allowed: InitOut%InvalidWithSSExctn = InputFileData%WaveMod == WaveMod_ExtFull .or. & ! 'Externally generated full wave-kinematics time series cannot be used with state-space wave excitations. Set WaveMod 0, 1, 1P#, 2, 3, 4, or 5.' diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index ba0421447a..4c02606695 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -3452,14 +3452,19 @@ def read_WaterKin(self,WaterKin_file): self.fst_vt['WaterKin']['Z_Grid'] = read_array(f,None,split_val='-',array_type=float) f.readline() self.fst_vt['WaterKin']['CurrentMod'] = int_read(f.readline().split()[0]) - f.readline() - f.readline() - data_line = readline_filterComments(f).split() - while data_line[0] and data_line[0][:3] != '---': # OpenFAST searches for ---, so we'll do the same - self.fst_vt['WaterKin']['z-depth'].append(float(data_line[0])) - self.fst_vt['WaterKin']['x-current'].append(float(data_line[1])) - self.fst_vt['WaterKin']['y-current'].append(float(data_line[2])) - data_line = readline_filterComments(f).split() + # depending on CurrentMod, the rest of the WaterKin input file changes + if self.fst_vt['WaterKin']['CurrentMod'] == 2: # user provided depths + self.fst_vt['WaterKin']['current_Z_type'] = int_read(f.readline().split()[0]) + self.fst_vt['WaterKin']['current_Z_Grid'] = read_array(f,None,split_val='-',array_type=float) + elif self.fst_vt['WaterKin']['CurrentMod'] == 1: # user provided depths and current speeds + f.readline() + f.readline() + data_line = readline_filterComments(f).split() + while data_line[0] and data_line[0][:3] != '---': # OpenFAST searches for ---, so we'll do the same + self.fst_vt['WaterKin']['z-depth'].append(float(data_line[0])) + self.fst_vt['WaterKin']['x-current'].append(float(data_line[1])) + self.fst_vt['WaterKin']['y-current'].append(float(data_line[2])) + data_line = readline_filterComments(f).split() f.close() def read_NonLinearEA(self,Stiffness_file): # read and return the nonlinear line stiffness lookup table for a given line diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index c1cd955de7..30e594181c 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -2471,7 +2471,7 @@ def write_MAP(self): os.fsync(f) f.close() - def write_MoorDyn(self): # TODO: see TODO's in FAST_reader/read_MoorDyn.py and make corresponding changes here + def write_MoorDyn(self): self.fst_vt['Fst']['MooringFile'] = self.FAST_namingOut + '_MoorDyn.dat' moordyn_file = os.path.join(self.FAST_runDirectory, self.fst_vt['Fst']['MooringFile']) @@ -2657,10 +2657,10 @@ def write_WaterKin(self,WaterKin_file): f.write('MoorDyn v2 Waves and Currents input file set up\n') f.write('This file was written by FAST_writer.py, comments from seed file have been dropped.\n') f.write('--------------------------- WAVES -------------------------------------\n') - f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveKinMod'], 'WaveKinMod', '- type of wave input {0 no waves; 3 set up grid of wave data based on time series}\n')) - f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveKinFile'], 'WaveKinFile', '- file containing wave elevation time series at 0,0,0\n')) - f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['dtWave'], 'dtWave', '- time step to use in setting up wave kinematics grid (s)\n')) - f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveDir'], 'WaveDir', '- wave heading (deg)\n')) + f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveKinMod'], 'WaveKinMod', '- flag for waves (0 no waves; 1 compute waves using provided elevation timeseries; 2 use wave elevation data from seastate)\n')) + f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveKinFile'], 'WaveKinFile', '- file containing wave elevation time series at 0,0,0. This is ignored if WaveKinMod is not 1\n')) + f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['dtWave'], 'dtWave', '- time step to use in setting up wave kinematics grid (s). This is ignored if WaveKinMod is not 1\n')) + f.write('{:<22} {:<11} {:}'.format(self.fst_vt['WaterKin']['WaveDir'], 'WaveDir', '- wave heading (deg). This is ignored if WaveKinMod is not 1\n')) f.write('{:<22} {:}'.format(self.fst_vt['WaterKin']['X_Type'], '- X wave input type (0: not used; 1: list values in ascending order; 2: uniform specified by -xlim, xlim, num)\n')) f.write('{:<22} {:}'.format(', '.join(['{:.3f}'.format(i) for i in self.fst_vt['WaterKin']['X_Grid']]), '- X wave grid point data separated by commas\n')) f.write('{:<22} {:}'.format(self.fst_vt['WaterKin']['Y_Type'], '- Y wave input type (0: not used; 1: list values in ascending order; 2: uniform specified by -Ylim, Ylim, num)\n')) @@ -2668,12 +2668,17 @@ def write_WaterKin(self,WaterKin_file): f.write('{:<22} {:}'.format(self.fst_vt['WaterKin']['Z_Type'], '- Z wave input type (0: not used; 1: list values in ascending order; 2: uniform specified by -Zlim, Zlim, num)\n')) f.write('{:<22} {:}'.format(', '.join(['{:.3f}'.format(i) for i in self.fst_vt['WaterKin']['Z_Grid']]), '- Z wave grid point data separated by commas\n')) f.write('--------------------------- CURRENT -------------------------------------\n') - f.write('{:} CurrentMod - type of current input (0 no current; 1 steady current profile described below) \n'.format(self.fst_vt['WaterKin']['CurrentMod'])) - f.write('z-depth x-current y-current\n') - f.write('(m) (m/s) (m/s)\n') - if self.fst_vt['WaterKin']['z-depth']: - for i in range(len(self.fst_vt['WaterKin']['z-depth'])): - f.write('{:.3f} {:.3f} {:.3f}'.format(self.fst_vt['WaterKin']['z-depth'][i], self.fst_vt['WaterKin']['x-current'][i], self.fst_vt['WaterKin']['y-current'][i])) + f.write('{:} CurrentMod - flag for currents (0 no current; 1 use depth and current described below; 2 use current speed data from SeaState and depth spacing defined below) \n'.format(self.fst_vt['WaterKin']['CurrentMod'])) + # depending on CurrentMod, the rest of the WaterKin input file changes + if self.fst_vt['WaterKin']['CurrentMod'] == 2: # user provided depths + f.write('{:<22} {:}'.format(self.fst_vt['WaterKin']['current_Z_Type'], '- Z current input type (0: not used; 1: list values in ascending order; 2: uniform specified by -Zlim, Zlim, num). This is ignored unless CurrentMod = 2\n')) + f.write('{:<22} {:}'.format(', '.join(['{:.3f}'.format(i) for i in self.fst_vt['WaterKin']['current_Z_Grid']]), '- Z current grid point data separated by commas. This is ignored unless CurrentMod = 2\n')) + elif self.fst_vt['WaterKin']['CurrentMod'] == 1: # user provided depths and current speeds + f.write('z-depth x-current y-current - This table is ignored unless CurrentMod = 1\n') + f.write('(m) (m/s) (m/s)\n') + if self.fst_vt['WaterKin']['z-depth']: + for i in range(len(self.fst_vt['WaterKin']['z-depth'])): + f.write('{:.3f} {:.3f} {:.3f}'.format(self.fst_vt['WaterKin']['z-depth'][i], self.fst_vt['WaterKin']['x-current'][i], self.fst_vt['WaterKin']['y-current'][i])) f.write('--------------------- need this line ------------------\n') f.flush() diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index b5ca00e1fe..99c76c0bc0 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -524,6 +524,9 @@ md_regression("md_horizontal" "moordyn") md_regression("md_no_line" "moordyn") md_regression("md_vertical" "moordyn") md_regression("md_BdyExtLdDmpg" "moordyn") +md_regression("md_VIV" "moordyn") +md_regression("md_waterkin2" "moordyn") +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") diff --git a/reg_tests/r-test b/reg_tests/r-test index 1b388a27ae..a18c4d72d6 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 1b388a27aec6c744840df3eb347679a5fefef8bc +Subproject commit a18c4d72d643ac8a8aeebba6e02465671c14c5b0