From b3c13d23dcb5bb6131848932436bf2412eba0755 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Mon, 10 Feb 2025 17:13:16 -0700 Subject: [PATCH 01/10] MD: RK4 in MoorDyn. Tested by hardcoding default to RK4 and running reg tests --- modules/moordyn/src/MoorDyn.f90 | 158 +++++++++++++++++++++-- modules/moordyn/src/MoorDyn_Registry.txt | 2 + modules/moordyn/src/MoorDyn_Types.f90 | 12 ++ 3 files changed, 161 insertions(+), 11 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 80d21a7810..338845696d 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -117,6 +117,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CHARACTER(20) :: LineOutString ! String to temporarially hold characters specifying line output options CHARACTER(20) :: OptString ! String to temporarially hold name of option variable CHARACTER(40) :: OptValue ! String to temporarially hold value of options variable input + CHARACTER(40) :: tSchemeString = 'RK2'! String to temporarially hold value of tScheme variable input (initial string sets RK2 as default if not provided by user) CHARACTER(40) :: DepthValue ! Temporarily stores the optional WtrDpth setting for MD, which could be a number or a filename CHARACTER(40) :: WaterKinValue ! Temporarily stores the optional WaterKin setting for MD, which is typically a filename INTEGER(IntKi) :: nOpts ! number of options lines in input file @@ -443,6 +444,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er end if else if ( OptString == 'DTM') THEN read (OptValue,*) p%dtM0 + else if ( OptString == 'TSCHEME') THEN + read (OptValue,*) tSchemeString else if ( OptString == 'G') then read (OptValue,*) p%g else if (( OptString == 'RHOW') .or. ( OptString == 'RHO')) then @@ -491,6 +494,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (p%writeLog > 1) then write(p%UnLog, '(A)' ) " - Options List:" write(p%UnLog, '(A17,f12.4)') " dtm : ", p%dtM0 + write(p%UnLog, '(A17,A)' ) " tScheme : ", tSchemeString write(p%UnLog, '(A17,f12.4)') " g : ", p%g write(p%UnLog, '(A17,f12.4)') " rhoW : ", p%rhoW write(p%UnLog, '(A17,A)' ) " Depth : ", DepthValue ! water depth input read in as a string to be processed by setupBathymetry @@ -544,6 +548,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! set up wave and current kinematics CALL setupWaterKin(WaterKinValue, p, InitInp%Tmax, ErrStat2, ErrMsg2); if(Failed()) return + ! set up time integration method + IF (tSchemeString == 'RK2') THEN + p%tScheme = 0 + ELSEIF (tSchemeString == 'RK4') THEN + p%tScheme = 1 + ELSE + CALL SetErrStat( ErrID_Fatal, 'Unrecognized tScheme option: '//tSchemeString//' Only RK2 and RK4 supported in MD-F', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + ENDIF ! ----------------------------- misc checks to be sorted ----------------------------- @@ -2238,15 +2252,22 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! allocate state vector and temporary state vectors based on size just calculated ALLOCATE ( x%states(m%Nxtra), m%xTemp%states(m%Nxtra), m%xdTemp%states(m%Nxtra), STAT = ErrStat2 ) - IF ( ErrStat2 /= ErrID_None ) THEN - ErrMsg = ' Error allocating state vectors.' - !CALL CleanUp() - RETURN - END IF + x%states = 0.0_DbKi m%xTemp%states = 0.0_DbKi m%xdTemp%states = 0.0_DbKi + ! Allocate kSum if using RK4 method. Not needed for RK2 becasue slopes not summed + IF (p%tScheme == 1_Intki) THEN + ALLOCATE (m%kSum%states(m%Nxtra), STAT = ErrStat2) + m%kSum%states = 0.0_DbKi + END IF + + IF ( ErrStat2 /= ErrID_None ) THEN + ErrMsg = ' Error allocating state vectors.' + !CALL CleanUp() + RETURN + END IF ! ================================ initialize system ================================ @@ -2677,7 +2698,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er t = 0.0_DbKi ! start time at zero ! because TimeStep wants an array... - call MD_CopyInput( u, u_array(1), MESH_NEWCOPY, ErrStat2, ErrMsg2 ) ! make a size=1 array of inputs (since MD_RK2 expects an array to InterpExtrap) + call MD_CopyInput( u, u_array(1), MESH_NEWCOPY, ErrStat2, ErrMsg2 ) ! make a size=1 array of inputs (since MD_RK2 and MD_RK4 expects an array to InterpExtrap) call MD_CopyInput( u, u_interp, MESH_NEWCOPY, ErrStat2, ErrMsg2 ) ! also make an inputs object to interpExtrap to t_array(1) = t ! fill in the times "array" for u_array @@ -2687,7 +2708,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er !loop through line integration time steps DO J = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) - CALL MD_RK2(t, dtM, u_interp, u_array, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) + Call MD_Step(t, dtM, u_interp, u_array, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) + IF ( ErrStat2 /= ErrID_None ) THEN + CALL CheckError(ErrStat2, ErrMsg2) + END IF ! check for NaNs - is this a good place/way to do it? DO K = 1, m%Nx @@ -3012,8 +3036,10 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er !loop through line integration time steps DO I = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) - CALL MD_RK2(t2, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) - + Call MD_Step(t2, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) + IF ( ErrStat2 /= ErrID_None ) THEN + CALL CheckError(ErrStat2, ErrMsg2) + END IF ! check for NaNs - is this a good place/way to do it? DO J = 1, m%Nx @@ -3855,7 +3881,36 @@ END SUBROUTINE MD_End !----------------------------------------------------------------------------------------================== - ! RK2 integrater (part of what was in TimeStep) + ! Advances the state one time step using the integration scheme specified in the input file + SUBROUTINE MD_Step ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg ) + + REAL(DbKi) , INTENT(INOUT) :: t ! intial time (s) for this integration step + REAL(DbKi) , INTENT(IN ) :: dtM ! single time step size (s) for this integration step + TYPE( MD_InputType ) , INTENT(INOUT) :: u_interp ! interpolated instantaneous input values to be calculated for each mooring time step + TYPE( MD_InputType ) , INTENT(INOUT) :: u(:) ! INTENT(IN ) + REAL(DbKi) , INTENT(IN ) :: t_array(:) ! times corresponding to elements of u(:)? + TYPE( MD_ParameterType ) , INTENT(IN ) :: p ! INTENT(IN ) + TYPE( MD_ContinuousStateType ) , INTENT(INOUT) :: x + TYPE( MD_DiscreteStateType ) , INTENT(IN ) :: xd ! INTENT(IN ) + TYPE( MD_ConstraintStateType ) , INTENT(IN ) :: z ! INTENT(IN ) + TYPE( MD_OtherStateType ) , INTENT(IN ) :: other ! INTENT(INOUT) + TYPE(MD_MiscVarType) , INTENT(INOUT) :: m ! INTENT(INOUT) + INTEGER(IntKi) , INTENT( OUT) :: ErrStat + CHARACTER(*) , INTENT( OUT) :: ErrMsg + + IF (p%tScheme == 0_Intki) THEN + CALL MD_RK2(t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg) + ELSEIF (p%tScheme == 1_Intki) THEN + CALL MD_RK4(t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg) + ELSE + ErrStat = ErrID_Fatal + ErrMsg = ' Unrecognized tScheme option in MD_Step' + RETURN + ENDIF + + END SUBROUTINE MD_Step + + ! RK2 integrator (part of what was in TimeStep) !-------------------------------------------------------------- SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg ) @@ -3887,8 +3942,9 @@ SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + ! k0 = m%xdTemp DO J = 1, m%Nx - m%xTemp%states(J) = x%states(J) + 0.5*dtM*m%xdTemp%states(J) !x1 = x0 + dt*f0/2.0; + m%xTemp%states(J) = x%states(J) + 0.5_DbKi*dtM*m%xdTemp%states(J) !x1 = x0 + dt*k0/2.0; END DO ! step 2 @@ -3896,6 +3952,7 @@ SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t+0.5*dtM) CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) !called with updated states x2 and time = t + dt/2.0 + ! k1 = m%xdTemp DO J = 1, m%Nx x%states(J) = x%states(J) + dtM*m%xdTemp%states(J) END DO @@ -3907,6 +3964,85 @@ SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat END SUBROUTINE MD_RK2 !-------------------------------------------------------------- + ! RK4 integrator + !-------------------------------------------------------------- + SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg ) + + REAL(DbKi) , INTENT(INOUT) :: t ! intial time (s) for this integration step + REAL(DbKi) , INTENT(IN ) :: dtM ! single time step size (s) for this integration step + TYPE( MD_InputType ) , INTENT(INOUT) :: u_interp ! interpolated instantaneous input values to be calculated for each mooring time step + TYPE( MD_InputType ) , INTENT(INOUT) :: u(:) ! INTENT(IN ) + REAL(DbKi) , INTENT(IN ) :: t_array(:) ! times corresponding to elements of u(:)? + TYPE( MD_ParameterType ) , INTENT(IN ) :: p ! INTENT(IN ) + TYPE( MD_ContinuousStateType ) , INTENT(INOUT) :: x + TYPE( MD_DiscreteStateType ) , INTENT(IN ) :: xd ! INTENT(IN ) + TYPE( MD_ConstraintStateType ) , INTENT(IN ) :: z ! INTENT(IN ) + TYPE( MD_OtherStateType ) , INTENT(IN ) :: other ! INTENT(INOUT) + TYPE(MD_MiscVarType) , INTENT(INOUT) :: m ! INTENT(INOUT) + INTEGER(IntKi) , INTENT( OUT) :: ErrStat + CHARACTER(*) , INTENT( OUT) :: ErrMsg + + + INTEGER(IntKi) :: I ! counter + INTEGER(IntKi) :: J ! counter + + ! ------------------------------------------------------------------------------- + ! RK4 integrator written here, calling CalcContStateDeriv + !-------------------------------------------------------------------------------- + + ! k0 + CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) + + CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + ! k0 = m%xdTemp + + ! k1 + DO J = 1, m%Nx + m%kSum%states(J) = m%xdTemp%states(J) ! k0 - In loop to avoid unecessary extra loop + m%xTemp%states(J) = x%states(J) + 0.5_DbKi*dtM*m%xdTemp%states(J) !x1 = x0 + dt*k0/2.0 = m%xTemp; + END DO + + CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) + + CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + ! k1 = m%xdTemp + + ! k2 + DO J = 1, m%Nx + m%kSum%states(J) = m%kSum%states(J) + 2.0_DbKi * m%xdTemp%states(J) ! 2 * k1 - In loop to avoid unecessary extra loop + m%xTemp%states(J) = x%states(J) + 0.5_DbKi*dtM*m%xdTemp%states(J) !x2 = x0 + dt*k1/2.0 = m%xTemp; + END DO + + CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) TODO: is this needed, it is already called for k1 + + CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + ! k2 = m%xdTemp + + ! k3 + DO J = 1, m%Nx + m%kSum%states(J) = m%kSum%states(J) + 2.0_DbKi * m%xdTemp%states(J) ! 2 * k2 - In loop to avoid unecessary extra loop + m%xTemp%states(J) = x%states(J) + dtM*m%xdTemp%states(J) !x3 = x0 + dt*k2 = m%xTemp; + END DO + + CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) + + CALL MD_CalcContStateDeriv( (t + dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + ! k3 = m%xdTemp + + ! Apply + DO J = 1, m%Nx + m%kSum%states(J) = m%kSum%states(J) + m%xdTemp%states(J) ! k3 - In loop to avoid unecessary extra loop + + ! x(t+dtM) = x(t) + dtM * kSum / 6 = x(t) + dtM * (k0 + 2*k1 + 2*k2 + k3) / 6 + x%states(J) = x%states(J) + dtM*(m%kSum%states(J)) / 6.0_DbKi + END DO + + t = t + dtM ! update time + + !TODO error check? <<<< + + END SUBROUTINE MD_RK4 + !-------------------------------------------------------------- !----------------------------------------------------------------------------------------================ ! this would do a full (coupling) time step and is no longer used diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index a6c9beb996..a1bab2859a 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -373,6 +373,7 @@ typedef ^ ^ IntKi Nxtra - typedef ^ ^ IntKi WaveTi - - - "current interpolation index for wave time series data" "" typedef ^ ^ MD_ContinuousStateType xTemp - - - "contains temporary state vector used in integration (put here so it's only allocated once)" typedef ^ ^ MD_ContinuousStateType xdTemp - - - "contains temporary state derivative vector used in integration (put here so it's only allocated once)" +typedef ^ ^ MD_ContinuousStateType kSum - - - "Sum of RK4 slope estimates: k0 + 2*k1 + 2*k2 + k3" typedef ^ ^ DbKi zeros6 {6} - - "array of zeros for convenience" typedef ^ ^ DbKi MDWrOutput {:} - - "Data from time step to be written to a MoorDyn output file" typedef ^ ^ DbKi LastOutTime - - - "Time of last writing to MD output files" @@ -403,6 +404,7 @@ typedef ^ ^ IntKi nCpldPoints {:} typedef ^ ^ IntKi NConns - 0 - "number of Connect type Points - not to be confused with NPoints" "" typedef ^ ^ IntKi NAnchs - 0 - "number of Anchor type Points" "" typedef ^ ^ DbKi Tmax - - - "simulation duration" "[s]" +typedef ^ ^ IntKi tScheme - 0 - "Time integration scheme (0 = RK2, 1 = RK4). Default is RK2" - typedef ^ ^ DbKi g - 9.81 - "gravitational constant (positive)" "[m/s^2]" typedef ^ ^ DbKi rhoW - 1025 - "density of seawater" "[kg/m^3]" typedef ^ ^ DbKi WtrDpth - - - "water depth" "[m]" diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index 601ff63433..cba62a494a 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -410,6 +410,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: WaveTi = 0_IntKi !< current interpolation index for wave time series data [] TYPE(MD_ContinuousStateType) :: xTemp !< contains temporary state vector used in integration (put here so it's only allocated once) [-] TYPE(MD_ContinuousStateType) :: xdTemp !< contains temporary state derivative vector used in integration (put here so it's only allocated once) [-] + TYPE(MD_ContinuousStateType) :: kSum !< Sum of RK4 slope estimates: k0 + 2*k1 + 2*k2 + k3 [-] REAL(DbKi) , DIMENSION(1:6) :: zeros6 = 0.0_R8Ki !< array of zeros for convenience [-] REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: MDWrOutput !< Data from time step to be written to a MoorDyn output file [-] REAL(DbKi) :: LastOutTime = 0.0_R8Ki !< Time of last writing to MD output files [-] @@ -441,6 +442,7 @@ MODULE MoorDyn_Types INTEGER(IntKi) :: NConns = 0 !< number of Connect type Points - not to be confused with NPoints [] INTEGER(IntKi) :: NAnchs = 0 !< number of Anchor type Points [] REAL(DbKi) :: Tmax = 0.0_R8Ki !< simulation duration [[s]] + INTEGER(IntKi) :: tScheme = 0 !< Time integration scheme (0 = RK2, 1 = RK4). Default is RK2 [-] REAL(DbKi) :: g = 9.81 !< gravitational constant (positive) [[m/s^2]] REAL(DbKi) :: rhoW = 1025 !< density of seawater [[kg/m^3]] REAL(DbKi) :: WtrDpth = 0.0_R8Ki !< water depth [[m]] @@ -3258,6 +3260,9 @@ subroutine MD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) call MD_CopyContState(SrcMiscData%xdTemp, DstMiscData%xdTemp, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return + call MD_CopyContState(SrcMiscData%kSum, DstMiscData%kSum, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return DstMiscData%zeros6 = SrcMiscData%zeros6 if (allocated(SrcMiscData%MDWrOutput)) then LB(1:1) = lbound(SrcMiscData%MDWrOutput) @@ -3454,6 +3459,8 @@ subroutine MD_DestroyMisc(MiscData, ErrStat, ErrMsg) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call MD_DestroyContState(MiscData%xdTemp, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call MD_DestroyContState(MiscData%kSum, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (allocated(MiscData%MDWrOutput)) then deallocate(MiscData%MDWrOutput) end if @@ -3570,6 +3577,7 @@ subroutine MD_PackMisc(RF, Indata) call RegPack(RF, InData%WaveTi) call MD_PackContState(RF, InData%xTemp) call MD_PackContState(RF, InData%xdTemp) + call MD_PackContState(RF, InData%kSum) call RegPack(RF, InData%zeros6) call RegPackAlloc(RF, InData%MDWrOutput) call RegPack(RF, InData%LastOutTime) @@ -3714,6 +3722,7 @@ subroutine MD_UnPackMisc(RF, OutData) call RegUnpack(RF, OutData%WaveTi); if (RegCheckErr(RF, RoutineName)) return call MD_UnpackContState(RF, OutData%xTemp) ! xTemp call MD_UnpackContState(RF, OutData%xdTemp) ! xdTemp + call MD_UnpackContState(RF, OutData%kSum) ! kSum call RegUnpack(RF, OutData%zeros6); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%MDWrOutput); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%LastOutTime); if (RegCheckErr(RF, RoutineName)) return @@ -3789,6 +3798,7 @@ subroutine MD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%NConns = SrcParamData%NConns DstParamData%NAnchs = SrcParamData%NAnchs DstParamData%Tmax = SrcParamData%Tmax + DstParamData%tScheme = SrcParamData%tScheme DstParamData%g = SrcParamData%g DstParamData%rhoW = SrcParamData%rhoW DstParamData%WtrDpth = SrcParamData%WtrDpth @@ -4209,6 +4219,7 @@ subroutine MD_PackParam(RF, Indata) call RegPack(RF, InData%NConns) call RegPack(RF, InData%NAnchs) call RegPack(RF, InData%Tmax) + call RegPack(RF, InData%tScheme) call RegPack(RF, InData%g) call RegPack(RF, InData%rhoW) call RegPack(RF, InData%WtrDpth) @@ -4312,6 +4323,7 @@ subroutine MD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%NConns); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NAnchs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Tmax); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%tScheme); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%g); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%rhoW); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%WtrDpth); if (RegCheckErr(RF, RoutineName)) return From 1f5936ac4b3d751ce02bd1591b41634f43fa9bba Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Thu, 13 Feb 2025 13:13:12 -0700 Subject: [PATCH 02/10] MD: cleaned up spacing, commented out timestep function, added dynamic line damping error check --- modules/moordyn/src/MoorDyn.f90 | 166 ++++++++++---------- modules/moordyn/src/MoorDyn_Line.f90 | 13 ++ modules/moordyn/src/MoorDyn_Registry.txt | 186 +++++++++++------------ 3 files changed, 189 insertions(+), 176 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 338845696d..d681172055 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -4044,117 +4044,117 @@ SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat END SUBROUTINE MD_RK4 !-------------------------------------------------------------- - !----------------------------------------------------------------------------------------================ - ! this would do a full (coupling) time step and is no longer used - SUBROUTINE TimeStep ( t, dtStep, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg ) + ! !----------------------------------------------------------------------------------------================ + ! ! this would do a full (coupling) time step and is no longer used + ! SUBROUTINE TimeStep ( t, dtStep, u, t_array, p, x, xd, z, other, m, ErrStat, ErrMsg ) - REAL(DbKi) , INTENT(INOUT) :: t - REAL(DbKi) , INTENT(IN ) :: dtStep ! how long to advance the time for - TYPE( MD_InputType ) , INTENT(INOUT) :: u(:) ! INTENT(IN ) - REAL(DbKi) , INTENT(IN ) :: t_array(:) ! times corresponding to elements of u(:)? - TYPE( MD_ParameterType ) , INTENT(IN ) :: p ! INTENT(IN ) - TYPE( MD_ContinuousStateType ) , INTENT(INOUT) :: x - TYPE( MD_DiscreteStateType ) , INTENT(IN ) :: xd ! INTENT(IN ) - TYPE( MD_ConstraintStateType ) , INTENT(IN ) :: z ! INTENT(IN ) - TYPE( MD_OtherStateType ) , INTENT(IN ) :: other ! INTENT(INOUT) - TYPE(MD_MiscVarType) , INTENT(INOUT) :: m ! INTENT(INOUT) - INTEGER(IntKi) , INTENT( OUT) :: ErrStat - CHARACTER(*) , INTENT( OUT) :: ErrMsg - - - TYPE(MD_ContinuousStateType) :: dxdt ! time derivatives of continuous states (initialized in CalcContStateDeriv) - TYPE(MD_ContinuousStateType) :: x2 ! temporary copy of continuous states used in RK2 calculations - INTEGER(IntKi) :: NdtM ! the number of time steps to make with the mooring model - Real(DbKi) :: dtM ! the actual time step size to use - INTEGER(IntKi) :: Nx ! size of states vector - INTEGER(IntKi) :: I ! counter - INTEGER(IntKi) :: J ! counter - TYPE(MD_InputType) :: u_interp ! interpolated instantaneous input values to be calculated for each mooring time step - - ! Real(DbKi) :: tDbKi ! double version because that's what MD_Input_ExtrapInterp needs. + ! REAL(DbKi) , INTENT(INOUT) :: t + ! REAL(DbKi) , INTENT(IN ) :: dtStep ! how long to advance the time for + ! TYPE( MD_InputType ) , INTENT(INOUT) :: u(:) ! INTENT(IN ) + ! REAL(DbKi) , INTENT(IN ) :: t_array(:) ! times corresponding to elements of u(:)? + ! TYPE( MD_ParameterType ) , INTENT(IN ) :: p ! INTENT(IN ) + ! TYPE( MD_ContinuousStateType ) , INTENT(INOUT) :: x + ! TYPE( MD_DiscreteStateType ) , INTENT(IN ) :: xd ! INTENT(IN ) + ! TYPE( MD_ConstraintStateType ) , INTENT(IN ) :: z ! INTENT(IN ) + ! TYPE( MD_OtherStateType ) , INTENT(IN ) :: other ! INTENT(INOUT) + ! TYPE(MD_MiscVarType) , INTENT(INOUT) :: m ! INTENT(INOUT) + ! INTEGER(IntKi) , INTENT( OUT) :: ErrStat + ! CHARACTER(*) , INTENT( OUT) :: ErrMsg + + + ! TYPE(MD_ContinuousStateType) :: dxdt ! time derivatives of continuous states (initialized in CalcContStateDeriv) + ! TYPE(MD_ContinuousStateType) :: x2 ! temporary copy of continuous states used in RK2 calculations + ! INTEGER(IntKi) :: NdtM ! the number of time steps to make with the mooring model + ! Real(DbKi) :: dtM ! the actual time step size to use + ! INTEGER(IntKi) :: Nx ! size of states vector + ! INTEGER(IntKi) :: I ! counter + ! INTEGER(IntKi) :: J ! counter + ! TYPE(MD_InputType) :: u_interp ! interpolated instantaneous input values to be calculated for each mooring time step + + ! ! Real(DbKi) :: tDbKi ! double version because that's what MD_Input_ExtrapInterp needs. - ! allocate space for x2 - CALL MD_CopyContState( x, x2, 0, ErrStat, ErrMsg) + ! ! allocate space for x2 + ! CALL MD_CopyContState( x, x2, 0, ErrStat, ErrMsg) - ! create space for arrays/meshes in u_interp ... is it efficient to do this every time step??? - CALL MD_CopyInput(u(1), u_interp, MESH_NEWCOPY, ErrStat, ErrMsg) + ! ! create space for arrays/meshes in u_interp ... is it efficient to do this every time step??? + ! CALL MD_CopyInput(u(1), u_interp, MESH_NEWCOPY, ErrStat, ErrMsg) - Nx = size(x%states) ! <<<< should this be the m%Nx parameter instead? + ! Nx = size(x%states) ! <<<< should this be the m%Nx parameter instead? - ! round dt to integer number of time steps - NdtM = ceiling(dtStep/p%dtM0) ! get number of mooring time steps to do based on desired time step size - dtM = dtStep/REAL(NdtM,DbKi) ! adjust desired time step to satisfy dt with an integer number of time steps + ! ! round dt to integer number of time steps + ! NdtM = ceiling(dtStep/p%dtM0) ! get number of mooring time steps to do based on desired time step size + ! dtM = dtStep/REAL(NdtM,DbKi) ! adjust desired time step to satisfy dt with an integer number of time steps - !loop through line integration time steps - DO I = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) + ! !loop through line integration time steps + ! DO I = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) - ! tDbKi = t ! get DbKi version of current time (why does ExtrapInterp except different time type than UpdateStates?) + ! ! tDbKi = t ! get DbKi version of current time (why does ExtrapInterp except different time type than UpdateStates?) - ! ------------------------------------------------------------------------------- - ! RK2 integrator written here, now calling CalcContStateDeriv - !-------------------------------------------------------------------------------- + ! ! ------------------------------------------------------------------------------- + ! ! RK2 integrator written here, now calling CalcContStateDeriv + ! !-------------------------------------------------------------------------------- - ! step 1 + ! ! step 1 - CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) + ! CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) - CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, dxdt, ErrStat, ErrMsg ) - DO J = 1, Nx - x2%states(J) = x%states(J) + 0.5*dtM*dxdt%states(J) !x1 = x0 + dt*f0/2.0; - END DO + ! CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, dxdt, ErrStat, ErrMsg ) + ! DO J = 1, Nx + ! x2%states(J) = x%states(J) + 0.5*dtM*dxdt%states(J) !x1 = x0 + dt*f0/2.0; + ! END DO - ! step 2 + ! ! step 2 - CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t+0.5*dtM) + ! CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t+0.5*dtM) - CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, x2, xd, z, other, m, dxdt, ErrStat, ErrMsg ) !called with updated states x2 and time = t + dt/2.0 - DO J = 1, Nx - x%states(J) = x%states(J) + dtM*dxdt%states(J) - END DO - - t = t + dtM ! update time - - !---------------------------------------------------------------------------------- - - ! >>> below should no longer be necessary thanks to using ExtrapInterp of u(:) within the mooring time stepping loop.. <<< - ! ! update Fairlead positions by integrating velocity and last position (do this AFTER the processing of the time step rather than before) - ! DO J = 1, p%nCpldPoints - ! DO K = 1, 3 - ! m%PointList(m%CpldPointIs(J))%r(K) = m%PointList(m%CpldPointIs(J))%r(K) + m%PointList(m%CpldPointIs(J))%rd(K)*dtM - ! END DO - ! END DO + ! CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, x2, xd, z, other, m, dxdt, ErrStat, ErrMsg ) !called with updated states x2 and time = t + dt/2.0 + ! DO J = 1, Nx + ! x%states(J) = x%states(J) + dtM*dxdt%states(J) + ! END DO + + ! t = t + dtM ! update time + + ! !---------------------------------------------------------------------------------- + + ! ! >>> below should no longer be necessary thanks to using ExtrapInterp of u(:) within the mooring time stepping loop.. <<< + ! ! ! update Fairlead positions by integrating velocity and last position (do this AFTER the processing of the time step rather than before) + ! ! DO J = 1, p%nCpldPoints + ! ! DO K = 1, 3 + ! ! m%PointList(m%CpldPointIs(J))%r(K) = m%PointList(m%CpldPointIs(J))%r(K) + m%PointList(m%CpldPointIs(J))%rd(K)*dtM + ! ! END DO + ! ! END DO - END DO ! I time steps + ! END DO ! I time steps - ! destroy dxdt and x2, and u_interp - CALL MD_DestroyContState( dxdt, ErrStat, ErrMsg) - CALL MD_DestroyContState( x2, ErrStat, ErrMsg) - CALL MD_DestroyInput(u_interp, ErrStat, ErrMsg) - IF ( ErrStat /= ErrID_None ) THEN - ErrMsg = ' Error destroying dxdt or x2.' - END IF - ! CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MD_UpdateStates') + ! ! destroy dxdt and x2, and u_interp + ! CALL MD_DestroyContState( dxdt, ErrStat, ErrMsg) + ! CALL MD_DestroyContState( x2, ErrStat, ErrMsg) + ! CALL MD_DestroyInput(u_interp, ErrStat, ErrMsg) + ! IF ( ErrStat /= ErrID_None ) THEN + ! ErrMsg = ' Error destroying dxdt or x2.' + ! END IF + ! ! CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MD_UpdateStates') - ! check for NaNs - is this a good place/way to do it? - DO J = 1, Nx - IF (Is_NaN(x%states(J))) THEN - ErrStat = ErrID_Fatal - ErrMsg = ' NaN state detected.' - END IF - END DO + ! ! check for NaNs - is this a good place/way to do it? + ! DO J = 1, Nx + ! IF (Is_NaN(x%states(J))) THEN + ! ErrStat = ErrID_Fatal + ! ErrMsg = ' NaN state detected.' + ! END IF + ! END DO - END SUBROUTINE TimeStep - !-------------------------------------------------------------- + ! END SUBROUTINE TimeStep + ! !-------------------------------------------------------------- diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index fe1131f32e..b978f3e756 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -75,6 +75,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%alphaMBL = LineProp%alphaMBL Line%vbeta = LineProp%vbeta Line%BA_D = LineProp%BA_D + print*, "Line%BA_D", Line%BA_D Line%EI = LineProp%EI !<<< for bending stiffness Line%Can = LineProp%Can @@ -119,6 +120,10 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! - we assume desired damping coefficient is zeta = -LineProp%BA ! - highest axial vibration mode of a segment is wn = sqrt(k/m) = 2N/UnstrLen*sqrt(EA/w) Line%BA = -LineProp%BA * Line%UnstrLen / Line%N * SQRT(LineProp%EA * LineProp%w) + if (p%writeLog > 0) then + write(p%UnLog,'(A)') "Based on -zeta of "//trim(Num2LStr(LineProp%BA))//", BA set to "//trim(Num2LStr(Line%BA)) + end if + IF (wordy > 1) print *, 'Based on zeta, BA set to ', Line%BA IF (wordy > 1) print *, 'Negative BA input detected, treating as -zeta. For zeta = ', -LineProp%BA, ', setting BA to ', Line%BA @@ -129,6 +134,14 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) IF (wordy > 1) print *, 'BA set as input to ', Line%BA, '. Corresponding zeta is ', temp END IF + IF (Line%BA_D < 0) THEN + ErrMsg = ' Line dynamic damping cannot be a ratio' + ErrStat = ErrID_Fatal + !CALL CleanUp() + RETURN + ENDIF + + !temp = 2*Line%N / Line%UnstrLen * sqrt( LineProp%EA / LineProp%w) / TwoPi !print *, 'Segment natural frequency is ', temp, ' Hz' diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index a1bab2859a..289dc4ac13 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -45,7 +45,7 @@ typedef ^ ^ Logical VisMeshes - .FA # nvm typedef ^ ^ MeshType FarmCoupledKinematics {:} - - "array of input kinematics meshes from each of the turbine-level MoorDyn instances" "[m, m/s]" # nvm typedef ^ ^ IntKi FarmNCpldBodies {:} - - "" "" # nvm typedef ^ ^ IntKi FarmNCpldRods {:} - - "" "" -# nvm typedef ^ ^ IntKi FarmNCpldPoints {:} - - "number of Fairlead Points" "" +# nvm typedef ^ ^ IntKi FarmNCpldPoints {:} - - "number of Fairlead Points" "" # ====================================== Internal data types ======================================================================== # line properties from line dictionary input @@ -55,8 +55,8 @@ typedef ^ ^ DbKi d - typedef ^ ^ DbKi w - - - "per-length weight in air" "[kg/m]" typedef ^ ^ DbKi EA - - - "axial stiffness" "[N]" typedef ^ ^ DbKi EA_D - - - "axial stiffness" "[N]" -typedef ^ ^ DbKi alphaMBL - - - "dynamic stiffness constant: Krd alpha term x MBL" "[N]" -typedef ^ ^ DbKi vbeta - - - "dynamic stiffness Lm slope: Krd beta term (to be multiplied by mean load)" "[N]" +typedef ^ ^ DbKi alphaMBL - - - "dynamic stiffness constant: Krd alpha term x MBL" "[N]" +typedef ^ ^ DbKi vbeta - - - "dynamic stiffness Lm slope: Krd beta term (to be multiplied by mean load)" "[N]" typedef ^ ^ DbKi BA - - - "internal damping coefficient times area" "[N-s]" typedef ^ ^ DbKi BA_D - - - "internal damping coefficient times area" "[N-s]" typedef ^ ^ DbKi EI - - - "bending stiffness" "[N-m]" @@ -69,8 +69,8 @@ typedef ^ ^ IntKi nEApoints - typedef ^ ^ DbKi stiffXs {30} - - "x array for stress-strain lookup table (up to nCoef)" typedef ^ ^ DbKi stiffYs {30} - - "y array for stress-strain lookup table" typedef ^ ^ IntKi nBApoints - - - "number of values in stress-strainrate lookup table (0 means using constant c)" -typedef ^ ^ DbKi dampXs {30} - - "x array for stress-strainrate lookup table (up to nCoef)" -typedef ^ ^ DbKi dampYs {30} - - "y array for stress-strainrate lookup table" +typedef ^ ^ DbKi dampXs {30} - - "x array for stress-strainrate lookup table (up to nCoef)" +typedef ^ ^ DbKi dampYs {30} - - "y array for stress-strainrate lookup table" typedef ^ ^ IntKi nEIpoints - - - "number of values in bending stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi bstiffXs {30} - - "x array for stress-strain lookup table (up to nCoef)" typedef ^ ^ DbKi bstiffYs {30} - - "y array for stress-strain lookup table" @@ -94,8 +94,8 @@ typedef ^ ^ IntKi AttachedC {30} typedef ^ ^ IntKi AttachedR {30} - - "list of IdNums of rods attached to this body" typedef ^ ^ IntKi nAttachedC - - - "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" +typedef ^ ^ DbKi rPointRel {3}{30} - - "relative position of point on body" +typedef ^ ^ DbKi r6RodRel {6}{30} - - "relative position and orientation of rod on body" typedef ^ ^ DbKi bodyM - - - "body mass (seperate from attached objects)" "[kg]" typedef ^ ^ DbKi bodyV - - - "body volume (for buoyancy calculation)" "[m^3]" typedef ^ ^ DbKi bodyI {3} - - "body 3x3 inertia matrix diagonals" "[kg-m^2]" @@ -122,19 +122,19 @@ typedef ^ ^ DbKi BlinL {3} typedef ^ ^ DbKi BquadL {3} - - "user-defined quadratic translational damping on the body in the local body-fixed frame" [N/(m/s)^2] # this is the Point type, which holds data for each point object -typedef ^ MD_Point IntKi IdNum - - - "integer identifier of this point" +typedef ^ MD_Point IntKi IdNum - - - "integer identifier of this point" typedef ^ ^ CHARACTER(10) type - - - "type of point: fix, vessel, point" typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 1=fixed, -1=coupled, 0=free" typedef ^ ^ IntKi Attached {10} - - "list of IdNums of lines attached to this point node" typedef ^ ^ IntKi Top {10} - - "list of ints specifying whether each line is attached at 1 = top/fairlead(end B), 0 = bottom/anchor(end A)" typedef ^ ^ IntKi nAttached - - - "number of attached lines" -typedef ^ ^ DbKi pointM - - - "point mass" "[kg]" -typedef ^ ^ DbKi pointV - - - "point volume" "[m^3]" -typedef ^ ^ DbKi pointFX - - - "" -typedef ^ ^ DbKi pointFY - - - "" -typedef ^ ^ DbKi pointFZ - - - "" -typedef ^ ^ DbKi pointCa - - - "added mass coefficient of point" "-" -typedef ^ ^ DbKi pointCdA - - - "product of drag force and frontal area of point" "[m^2]" +typedef ^ ^ DbKi pointM - - - "point mass" "[kg]" +typedef ^ ^ DbKi pointV - - - "point volume" "[m^3]" +typedef ^ ^ DbKi pointFX - - - "" +typedef ^ ^ DbKi pointFY - - - "" +typedef ^ ^ DbKi pointFZ - - - "" +typedef ^ ^ DbKi pointCa - - - "added mass coefficient of point" "-" +typedef ^ ^ DbKi pointCdA - - - "product of drag force and frontal area of point" "[m^2]" typedef ^ ^ DbKi time - - - "current time" "[s]" typedef ^ ^ DbKi r {3} - - "position" typedef ^ ^ DbKi rd {3} - - "velocity" @@ -142,7 +142,7 @@ typedef ^ ^ DbKi a {3} typedef ^ ^ DbKi U {3} - - "water velocity at node" "[m/s]" typedef ^ ^ DbKi Ud {3} - - "water acceleration at node" "[m/s^2]" typedef ^ ^ DbKi zeta - - - "water surface elevation above node" "[m]" -typedef ^ ^ DbKi PDyn {:} - - "water dynamic pressure at node" "[Pa]" +typedef ^ ^ DbKi PDyn {:} - - "water dynamic pressure at node" "[Pa]" typedef ^ ^ DbKi Fnet {3} - - "total force on node (excluding inertial loads)" typedef ^ ^ DbKi M {3}{3} - - "node mass matrix, from attached lines" typedef ^ ^ DbKi Fext {3} - - "vector of user-defined external force on the point always in the global frame" [N] @@ -150,9 +150,9 @@ typedef ^ ^ DbKi Blin {3} typedef ^ ^ DbKi Bquad {3} - - "user-defined quadratic translational damping on the point always in the global frame" [N/(m/s)^2] # this is the Rod type, which holds data for each Rod object -typedef ^ MD_Rod IntKi IdNum - - - "integer identifier of this Line" -typedef ^ ^ CHARACTER(10) type - - - "type of Rod. should match one of RodProp names" -typedef ^ ^ IntKi PropsIdNum - - - "the IdNum of the associated rod properties" - +typedef ^ MD_Rod IntKi IdNum - - - "integer identifier of this Line" +typedef ^ ^ CHARACTER(10) type - - - "type of Rod. should match one of RodProp names" +typedef ^ ^ IntKi PropsIdNum - - - "the IdNum of the associated rod properties" - typedef ^ ^ IntKi typeNum - - - "integer identifying the type. 0=free, 1=pinned, 2=fixed, -1=coupledpinned, -2=coupled" typedef ^ ^ IntKi AttachedA {10} - - "list of IdNums of lines attached to end A" typedef ^ ^ IntKi AttachedB {10} - - "list of IdNums of lines attached to end B" @@ -168,64 +168,64 @@ typedef ^ ^ DbKi UnstrLen - typedef ^ ^ DbKi mass - - - "mass of the rod" "[kg]" typedef ^ ^ DbKi rho - - - "density" "[kg/m3]" typedef ^ ^ DbKi d - - - "volume-equivalent diameter" "[m]" -typedef ^ ^ DbKi Can - - - "" "[-]" -typedef ^ ^ DbKi Cat - - - "" "[-]" -typedef ^ ^ DbKi Cdn - - - "" "[-]" -typedef ^ ^ DbKi Cdt - - - "" "[-]" -typedef ^ ^ DbKi CdEnd - - - "drag coefficient for rod end" "[-]" -typedef ^ ^ DbKi CaEnd - - - "added mass coefficient for rod end" "[-]" -typedef ^ ^ DbKi time - - - "current time" "[s]" -typedef ^ ^ DbKi roll - - - "roll relative to vertical" "[rad]" -typedef ^ ^ DbKi pitch - - - "pitch relative to vertical" "[rad]" -typedef ^ ^ DbKi h0 - - - "submerged length of rod axis, distance along rod centerline from end A to the waterplane (0 <= h0 <= L)" "m" -typedef ^ ^ DbKi r {:}{:} - - "node positions" - -typedef ^ ^ DbKi rd {:}{:} - - "node velocities" - -typedef ^ ^ DbKi q {3} - - "tangent vector for rod as a whole" - -typedef ^ ^ DbKi l {:} - - "segment unstretched length" "[m]" -typedef ^ ^ DbKi V {:} - - "segment volume" "[m^3]" -typedef ^ ^ DbKi U {:}{:} - - "water velocity at node" "[m/s]" -typedef ^ ^ DbKi Ud {:}{:} - - "water acceleration at node" "[m/s^2]" -typedef ^ ^ DbKi zeta {:} - - "water surface elevation above node" "[m]" -typedef ^ ^ DbKi PDyn {:} - - "water dynamic pressure at node" "[Pa]" -typedef ^ ^ DbKi W {:}{:} - - "weight vectors" "[N]" -typedef ^ ^ DbKi Bo {:}{:} - - "buoyancy force vectors" "[N]" -typedef ^ ^ DbKi Pd {:}{:} - - "dynamic pressure force vectors" "[N]" -typedef ^ ^ DbKi Dp {:}{:} - - "node drag (transverse)" "[N]" -typedef ^ ^ DbKi Dq {:}{:} - - "node drag (axial)" "[N]" -typedef ^ ^ DbKi Ap {:}{:} - - "node added mass forcing (transverse)" "[N]" -typedef ^ ^ DbKi Aq {:}{:} - - "node added mass forcing (axial)" "[N]" -typedef ^ ^ DbKi B {:}{:} - - "node bottom contact force" "[N]" -typedef ^ ^ DbKi Bp {:}{:} - - "transverse damping force" "[N]" -typedef ^ ^ DbKi Bq {:}{:} - - "axial damping force" "[N]" -typedef ^ ^ DbKi Fnet {:}{:} - - "total force on node" "[N]" -typedef ^ ^ DbKi M {:}{:}{:} - - "node mass matrix" "[kg]" -typedef ^ ^ DbKi FextA {3} - - "external forces from attached lines on/about end A " - -typedef ^ ^ DbKi FextB {3} - - "external forces from attached lines on/about end A " - -typedef ^ ^ DbKi Mext {3} - - "external moment vector holding sum of any externally applied moments i.e. bending lines" - -typedef ^ ^ DbKi r6 {6} - - "6 DOF position vector" - -typedef ^ ^ DbKi v6 {6} - - "6 DOF velocity vector" - -typedef ^ ^ DbKi a6 {6} - - "6 DOF acceleration vector (only used for coupled Rods)" - -typedef ^ ^ DbKi F6net {6} - - "total force and moment about end A (excluding inertial loads) that Rod may exert on whatever it's attached to" -typedef ^ ^ DbKi M6net {6}{6} - - "total mass matrix about end A of Rod and any attached Points" -typedef ^ ^ DbKi Imat {3}{3} - - "inertia about CG in global frame" -typedef ^ ^ DbKi OrMat {3}{3} - - "DCM for body orientation" -typedef ^ ^ IntKi RodUnOut - - - "unit number of rod output file" -typedef ^ ^ DbKi RodWrOutput {:} - - "one row of output data for this rod" -typedef ^ ^ DbKi FextU {3} - - "vector of user-defined external force on the rod end A always in the local body-fixed frame" "[N]" -typedef ^ ^ DbKi Blin {2} - - "linear damping, transverse damping for rod element always in the local body-fixed frame" "[N/(m/s)]" -typedef ^ ^ DbKi Bquad {2} - - "quadratic damping, transverse damping for rod element always in the local body-fixed frame" "[N/(m/s)^2]" +typedef ^ ^ DbKi Can - - - "" "[-]" +typedef ^ ^ DbKi Cat - - - "" "[-]" +typedef ^ ^ DbKi Cdn - - - "" "[-]" +typedef ^ ^ DbKi Cdt - - - "" "[-]" +typedef ^ ^ DbKi CdEnd - - - "drag coefficient for rod end" "[-]" +typedef ^ ^ DbKi CaEnd - - - "added mass coefficient for rod end" "[-]" +typedef ^ ^ DbKi time - - - "current time" "[s]" +typedef ^ ^ DbKi roll - - - "roll relative to vertical" "[rad]" +typedef ^ ^ DbKi pitch - - - "pitch relative to vertical" "[rad]" +typedef ^ ^ DbKi h0 - - - "submerged length of rod axis, distance along rod centerline from end A to the waterplane (0 <= h0 <= L)" "m" +typedef ^ ^ DbKi r {:}{:} - - "node positions" - +typedef ^ ^ DbKi rd {:}{:} - - "node velocities" - +typedef ^ ^ DbKi q {3} - - "tangent vector for rod as a whole" - +typedef ^ ^ DbKi l {:} - - "segment unstretched length" "[m]" +typedef ^ ^ DbKi V {:} - - "segment volume" "[m^3]" +typedef ^ ^ DbKi U {:}{:} - - "water velocity at node" "[m/s]" +typedef ^ ^ DbKi Ud {:}{:} - - "water acceleration at node" "[m/s^2]" +typedef ^ ^ DbKi zeta {:} - - "water surface elevation above node" "[m]" +typedef ^ ^ DbKi PDyn {:} - - "water dynamic pressure at node" "[Pa]" +typedef ^ ^ DbKi W {:}{:} - - "weight vectors" "[N]" +typedef ^ ^ DbKi Bo {:}{:} - - "buoyancy force vectors" "[N]" +typedef ^ ^ DbKi Pd {:}{:} - - "dynamic pressure force vectors" "[N]" +typedef ^ ^ DbKi Dp {:}{:} - - "node drag (transverse)" "[N]" +typedef ^ ^ DbKi Dq {:}{:} - - "node drag (axial)" "[N]" +typedef ^ ^ DbKi Ap {:}{:} - - "node added mass forcing (transverse)" "[N]" +typedef ^ ^ DbKi Aq {:}{:} - - "node added mass forcing (axial)" "[N]" +typedef ^ ^ DbKi B {:}{:} - - "node bottom contact force" "[N]" +typedef ^ ^ DbKi Bp {:}{:} - - "transverse damping force" "[N]" +typedef ^ ^ DbKi Bq {:}{:} - - "axial damping force" "[N]" +typedef ^ ^ DbKi Fnet {:}{:} - - "total force on node" "[N]" +typedef ^ ^ DbKi M {:}{:}{:} - - "node mass matrix" "[kg]" +typedef ^ ^ DbKi FextA {3} - - "external forces from attached lines on/about end A " - +typedef ^ ^ DbKi FextB {3} - - "external forces from attached lines on/about end A " - +typedef ^ ^ DbKi Mext {3} - - "external moment vector holding sum of any externally applied moments i.e. bending lines" - +typedef ^ ^ DbKi r6 {6} - - "6 DOF position vector" - +typedef ^ ^ DbKi v6 {6} - - "6 DOF velocity vector" - +typedef ^ ^ DbKi a6 {6} - - "6 DOF acceleration vector (only used for coupled Rods)" - +typedef ^ ^ DbKi F6net {6} - - "total force and moment about end A (excluding inertial loads) that Rod may exert on whatever it's attached to" +typedef ^ ^ DbKi M6net {6}{6} - - "total mass matrix about end A of Rod and any attached Points" +typedef ^ ^ DbKi Imat {3}{3} - - "inertia about CG in global frame" +typedef ^ ^ DbKi OrMat {3}{3} - - "DCM for body orientation" +typedef ^ ^ IntKi RodUnOut - - - "unit number of rod output file" +typedef ^ ^ DbKi RodWrOutput {:} - - "one row of output data for this rod" +typedef ^ ^ DbKi FextU {3} - - "vector of user-defined external force on the rod end A always in the local body-fixed frame" "[N]" +typedef ^ ^ DbKi Blin {2} - - "linear damping, transverse damping for rod element always in the local body-fixed frame" "[N/(m/s)]" +typedef ^ ^ DbKi Bquad {2} - - "quadratic damping, transverse damping for rod element always in the local body-fixed frame" "[N/(m/s)^2]" # this is the Line type, which holds data for each line object typedef ^ MD_Line IntKi IdNum - - - "integer identifier of this Line" -#typedef ^ ^ CHARACTER(10) type - - - "type of line. should match one of LineProp names" +#typedef ^ ^ CHARACTER(10) type - - - "type of line. should match one of LineProp names" typedef ^ ^ IntKi PropsIdNum - - - "the IdNum of the associated line properties" - typedef ^ ^ IntKi ElasticMod - - - "Which elasticity model to use: {0 basic, 1 viscoelastic, 2 future SYCOM} " - typedef ^ ^ IntKi OutFlagList {20} - - "array specifying what line quantities should be output (1 vs 0)" - typedef ^ ^ IntKi CtrlChan - - - "index of control channel that will drive line active tensioning (0 for none)" - -typedef ^ ^ IntKi FairPoint - - - "IdNum of Point at fairlead" -typedef ^ ^ IntKi AnchPoint - - - "IdNum of Point at anchor" +typedef ^ ^ IntKi FairPoint - - - "IdNum of Point at fairlead" +typedef ^ ^ IntKi AnchPoint - - - "IdNum of Point at anchor" typedef ^ ^ IntKi N - - - "The number of elements in the line" - typedef ^ ^ IntKi endTypeA - - - "type of connection at end A: 0=pinned to Point, 1=cantilevered to Rod." - typedef ^ ^ IntKi endTypeB - - - "type of connection at end B: 0=pinned to Point, 1=cantilevered to Rod." - @@ -234,15 +234,15 @@ typedef ^ ^ DbKi rho - typedef ^ ^ DbKi d - - - "volume-equivalent diameter" "[m]" typedef ^ ^ DbKi EA - - - "stiffness" "[N]" typedef ^ ^ DbKi EA_D - - - "constant dynamic stiffness when using viscoelastic model" "[N]" -typedef ^ ^ DbKi alphaMBL - - - "load dependent dynamic stiffness constant: Krd alpha term x MBL" "[N]" -typedef ^ ^ DbKi vbeta - - - "load dependent dynamic stiffness Lm slope: Krd beta term (to be multiplied by mean load)" "[N]" +typedef ^ ^ DbKi alphaMBL - - - "load dependent dynamic stiffness constant: Krd alpha term x MBL" "[N]" +typedef ^ ^ DbKi vbeta - - - "load dependent dynamic stiffness Lm slope: Krd beta term (to be multiplied by mean load)" "[N]" typedef ^ ^ DbKi BA - - - "internal damping coefficient times area for this line only" "[N-s]" typedef ^ ^ DbKi BA_D - - - "dynamic internal damping coefficient times area when using viscoelastic model" "[N-s]" typedef ^ ^ DbKi EI - - - "bending stiffness" "[N-m]" -typedef ^ ^ DbKi Can - - - "" "[-]" -typedef ^ ^ DbKi Cat - - - "" "[-]" -typedef ^ ^ DbKi Cdn - - - "" "[-]" -typedef ^ ^ DbKi Cdt - - - "" "[-]" +typedef ^ ^ DbKi Can - - - "" "[-]" +typedef ^ ^ DbKi Cat - - - "" "[-]" +typedef ^ ^ DbKi Cdn - - - "" "[-]" +typedef ^ ^ DbKi Cdt - - - "" "[-]" typedef ^ ^ IntKi nEApoints - - - "number of values in stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi stiffXs {30} - - "x array for stress-strain lookup table (up to nCoef)" typedef ^ ^ DbKi stiffYs {30} - - "y array for stress-strain lookup table" @@ -350,20 +350,20 @@ typedef ^ ^ MD_RodProp RodTypeList {:} typedef ^ ^ MD_Body GroundBody - - - "the single ground body which is the parent of all stationary points" - typedef ^ ^ MD_Body BodyList {:} - - "array of body objects" - typedef ^ ^ MD_Rod RodList {:} - - "array of rod objects" - -typedef ^ ^ MD_Point PointList {:} - - "array of point objects" - +typedef ^ ^ MD_Point PointList {:} - - "array of point objects" - typedef ^ ^ MD_Line LineList {:} - - "array of line objects" - typedef ^ ^ MD_ExtLd ExtLdList {:} - - "array of external load objects" - typedef ^ ^ MD_Fail FailList {:} - - "array of line objects" - -typedef ^ ^ IntKi FreePointIs {:} - - "array of free point indices in PointList vector" "" -typedef ^ ^ IntKi CpldPointIs {:}{:} - - "array of coupled/fairlead point indices in PointList vector" "" +typedef ^ ^ IntKi FreePointIs {:} - - "array of free point indices in PointList vector" "" +typedef ^ ^ IntKi CpldPointIs {:}{:} - - "array of coupled/fairlead point indices in PointList vector" "" typedef ^ ^ IntKi FreeRodIs {:} - - "array of free rod indices in RodList vector" "" -typedef ^ ^ IntKi CpldRodIs {:}{:} - - "array of coupled/fairlead rod indices in RodList vector" "" +typedef ^ ^ IntKi CpldRodIs {:}{:} - - "array of coupled/fairlead rod indices in RodList vector" "" typedef ^ ^ IntKi FreeBodyIs {:} - - "array of free body indices in BodyList vector" "" -typedef ^ ^ IntKi CpldBodyIs {:}{:} - - "array of coupled body indices in BodyList vector" "" +typedef ^ ^ IntKi CpldBodyIs {:}{:} - - "array of coupled body indices in BodyList vector" "" typedef ^ ^ IntKi LineStateIs1 {:} - - "starting index of each line's states in state vector" "" typedef ^ ^ IntKi LineStateIsN {:} - - "ending index of each line's states in state vector" "" -typedef ^ ^ IntKi PointStateIs1 {:} - - "starting index of each point's states in state vector" "" -typedef ^ ^ IntKi PointStateIsN {:} - - "ending index of each point's states in state vector" "" +typedef ^ ^ IntKi PointStateIs1 {:} - - "starting index of each point's states in state vector" "" +typedef ^ ^ IntKi PointStateIsN {:} - - "ending index of each point's states in state vector" "" typedef ^ ^ IntKi RodStateIs1 {:} - - "starting index of each rod's states in state vector" "" typedef ^ ^ IntKi RodStateIsN {:} - - "ending index of each rod's states in state vector" "" typedef ^ ^ IntKi BodyStateIs1 {:} - - "starting index of each body's states in state vector" "" @@ -376,19 +376,19 @@ typedef ^ ^ MD_ContinuousStateType xdTemp - typedef ^ ^ MD_ContinuousStateType kSum - - - "Sum of RK4 slope estimates: k0 + 2*k1 + 2*k2 + k3" typedef ^ ^ DbKi zeros6 {6} - - "array of zeros for convenience" typedef ^ ^ DbKi MDWrOutput {:} - - "Data from time step to be written to a MoorDyn output file" -typedef ^ ^ DbKi LastOutTime - - - "Time of last writing to MD output files" -typedef ^ ^ ReKi PtfmInit {6} - - "initial position of platform for an individual (non-farm) MD instance" - -typedef ^ ^ DbKi BathymetryGrid {:}{:} - - "matrix describing the bathymetry in a grid of x's and y's" -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 ^ ^ DbKi LastOutTime - - - "Time of last writing to MD output files" +typedef ^ ^ ReKi PtfmInit {6} - - "initial position of platform for an individual (non-farm) MD instance" - +typedef ^ ^ DbKi BathymetryGrid {:}{:} - - "matrix describing the bathymetry in a grid of x's and y's" +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" ## ============================== Parameters ============================================================================================================================================ typedef ^ ParameterType IntKi nLineTypes - 0 - "number of line types" "" typedef ^ ^ IntKi nRodTypes - 0 - "number of rod types" "" -typedef ^ ^ IntKi nPoints - 0 - "number of Point objects" "" -typedef ^ ^ IntKi nPointsExtra - 0 - "number of Point objects including space for extra ones that could arise from line failures" "" +typedef ^ ^ IntKi nPoints - 0 - "number of Point objects" "" +typedef ^ ^ IntKi nPointsExtra - 0 - "number of Point objects including space for extra ones that could arise from line failures" "" typedef ^ ^ IntKi nBodies - 0 - "number of Body objects" "" typedef ^ ^ IntKi nRods - 0 - "number of Rod objects" "" typedef ^ ^ IntKi nLines - 0 - "number of Line objects" "" @@ -397,10 +397,10 @@ typedef ^ ^ IntKi nCtrlChans - typedef ^ ^ IntKi nFails - 0 - "number of failure conditions" "" typedef ^ ^ IntKi nFreeBodies - 0 - "" "" typedef ^ ^ IntKi nFreeRods - 0 - "" "" -typedef ^ ^ IntKi nFreePoints - 0 - "" "" +typedef ^ ^ IntKi nFreePoints - 0 - "" "" typedef ^ ^ IntKi nCpldBodies {:} - - "number of coupled bodies (for FAST.Farm, size>1 with an entry for each turbine)" "" typedef ^ ^ IntKi nCpldRods {:} - - "number of coupled rods (for FAST.Farm, size>1 with an entry for each turbine)" "" -typedef ^ ^ IntKi nCpldPoints {:} - - "number of coupled points (for FAST.Farm, size>1 with an entry for each turbine)" "" +typedef ^ ^ IntKi nCpldPoints {:} - - "number of coupled points (for FAST.Farm, size>1 with an entry for each turbine)" "" typedef ^ ^ IntKi NConns - 0 - "number of Connect type Points - not to be confused with NPoints" "" typedef ^ ^ IntKi NAnchs - 0 - "number of Anchor type Points" "" typedef ^ ^ DbKi Tmax - - - "simulation duration" "[s]" @@ -477,10 +477,10 @@ typedef ^ ^ ReKi DeltaLdot {:} ## ============================== Outputs ============================================================================================================================================ -typedef ^ OutputType MeshType CoupledLoads {:} - - "array of point meshes for mooring reaction forces (and moments) at coupling points" "[N]" +typedef ^ OutputType MeshType CoupledLoads {:} - - "array of point meshes for mooring reaction forces (and moments) at coupling points" "[N]" typedef ^ ^ ReKi WriteOutput {:} - - "output vector returned to glue code" "" # should CoupledLoads be an array? -#typedef ^ ^ DbKi rAll {:}{:} - - "Mesh of all point positions: bodies, rods, points, line internal nodes" - +#typedef ^ ^ DbKi rAll {:}{:} - - "Mesh of all point positions: bodies, rods, points, line internal nodes" - typedef ^ ^ MeshType VisLinesMesh {:} - - "Line2 mesh for visualizing mooring lines" - typedef ^ ^ MeshType VisRodsMesh {:} - - "Line2 mesh for visualizing mooring rods" - typedef ^ ^ MeshType VisBodiesMesh {:} - - "Point mesh for visualizing mooring bodies" - From 96fd04a00015699308ddf8d7de8a9ab673c2a814 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Thu, 13 Feb 2025 13:15:44 -0700 Subject: [PATCH 03/10] MD: removed forgotten print statement --- modules/moordyn/src/MoorDyn_Line.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index b978f3e756..20c90417b8 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -75,7 +75,6 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%alphaMBL = LineProp%alphaMBL Line%vbeta = LineProp%vbeta Line%BA_D = LineProp%BA_D - print*, "Line%BA_D", Line%BA_D Line%EI = LineProp%EI !<<< for bending stiffness Line%Can = LineProp%Can From 3201d74b16d93c318deca63777681bf46ecf88b6 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Mon, 24 Feb 2025 11:29:57 -0700 Subject: [PATCH 04/10] MD: VIV attempt 1, broken --- modules/moordyn/src/MoorDyn.f90 | 125 ++++++----- modules/moordyn/src/MoorDyn_IO.f90 | 50 +++-- modules/moordyn/src/MoorDyn_Line.f90 | 251 +++++++++++++++++++++-- modules/moordyn/src/MoorDyn_Registry.txt | 21 +- modules/moordyn/src/MoorDyn_Types.f90 | 132 +++++++++++- 5 files changed, 497 insertions(+), 82 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index d681172055..4fc13f5b30 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -36,7 +36,7 @@ MODULE MoorDyn TYPE(ProgDesc), PARAMETER :: MD_ProgDesc = ProgDesc( 'MoorDyn', 'v2.2.2', '2024-01-16' ) - INTEGER(IntKi), PARAMETER :: wordy = 0 ! verbosity level. >1 = more console output + INTEGER(IntKi), PARAMETER :: wordy = 1 ! verbosity level. >1 = more console output PUBLIC :: MD_Init PUBLIC :: MD_UpdateStates @@ -474,7 +474,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er read (OptValue,*) p%mu_kA else if ( OptString == 'MC') then read (OptValue,*) p%mc - else if ( OptString == 'CV') then + else if (( OptString == 'CV') .or. (OptString == 'FRICDAMP')) then read (OptValue,*) p%cv else if ( OptString == 'INERTIALF') then read (OptValue,*) p%inertialF @@ -549,6 +549,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er CALL setupWaterKin(WaterKinValue, p, InitInp%Tmax, ErrStat2, ErrMsg2); if(Failed()) return ! set up time integration method + CALL Conv2UC(tSchemeString) IF (tSchemeString == 'RK2') THEN p%tScheme = 0 ELSEIF (tSchemeString == 'RK4') THEN @@ -646,16 +647,34 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er Line = NextLine(i) ! check for correct number of columns in current line - IF ( CountWords( Line ) /= 10 ) THEN - CALL SetErrStat( ErrID_Fatal, ' Unable to parse Line type '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be 10 columns.', ErrStat, ErrMsg, RoutineName ) + IF ( CountWords( Line ) < 10 ) THEN + CALL SetErrStat( ErrID_Fatal, ' Unable to parse Line type '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must be at least 10 columns.', ErrStat, ErrMsg, RoutineName ) CALL CleanUp() RETURN END IF - - ! parse out entries: Name Diam MassDenInAir EA cIntDamp EI Cd Ca CdAx CaAx + + ! parse out entries: Name Diam MassDenInAir EA cIntDamp EI Cd Ca CdAx CaAx Cl dF cF + if (CountWords( Line ) == 10) then READ(Line,*,IOSTAT=ErrStat2) m%LineTypeList(l)%name, m%LineTypeList(l)%d, & m%LineTypeList(l)%w, tempString1, tempString2, tempString3, & m%LineTypeList(l)%Cdn, m%LineTypeList(l)%Can, m%LineTypeList(l)%Cdt, m%LineTypeList(l)%Cat + elseif (CountWords( Line ) == 11) then + READ(Line,*,IOSTAT=ErrStat2) m%LineTypeList(l)%name, m%LineTypeList(l)%d, & + m%LineTypeList(l)%w, tempString1, tempString2, tempString3, & + m%LineTypeList(l)%Cdn, m%LineTypeList(l)%Can, m%LineTypeList(l)%Cdt, m%LineTypeList(l)%Cat, m%LineTypeList(l)%Cl + m%LineTypeList(l)%dF = 0.08 ! set to default Thorsen synchronization range if not provided + m%LineTypeList(l)%cF = 0.18 ! set to default Thorsen synchronization centering if not provided + elseif (CountWords( Line ) == 13) then + READ(Line,*,IOSTAT=ErrStat2) m%LineTypeList(l)%name, m%LineTypeList(l)%d, & + m%LineTypeList(l)%w, tempString1, tempString2, tempString3, & + m%LineTypeList(l)%Cdn, m%LineTypeList(l)%Can, m%LineTypeList(l)%Cdt, m%LineTypeList(l)%Cat, & + m%LineTypeList(l)%Cl, m%LineTypeList(l)%dF, m%LineTypeList(l)%cF + else + CALL SetErrStat( ErrID_Fatal, ' Unable to parse Line type '//trim(Num2LStr(l))//' on row '//trim(Num2LStr(i))//' in input file. Row has wrong number of columns. Must have 10, 11, or 13 columns.', ErrStat, ErrMsg, RoutineName ) + CALL CleanUp() + RETURN + endif + IF ( ErrStat2 /= ErrID_None ) THEN CALL SetErrStat( ErrID_Fatal, 'Failed to process line type inputs of entry '//trim(Num2LStr(l))//'. Check formatting and correct number of columns.', ErrStat, ErrMsg, RoutineName ) @@ -721,10 +740,16 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er write(p%UnLog, '(A12,A)' ) " name: ", trim(m%LineTypeList(l)%name) write(p%UnLog, '(A12,f12.4)') " d : ", m%LineTypeList(l)%d write(p%UnLog, '(A12,f12.4)') " w : ", m%LineTypeList(l)%w + write(p%UnLog, '(A12,A)' ) " EA : ", tempString1 + write(p%UnLog, '(A12,A)' ) " EA_D: ", tempString2 + write(p%UnLog, '(A12,A)' ) " EI : ", tempString3 write(p%UnLog, '(A12,f12.4)') " Cdn : ", m%LineTypeList(l)%Cdn write(p%UnLog, '(A12,f12.4)') " Can : ", m%LineTypeList(l)%Can write(p%UnLog, '(A12,f12.4)') " Cdt : ", m%LineTypeList(l)%Cdt write(p%UnLog, '(A12,f12.4)') " Cat : ", m%LineTypeList(l)%Cat + write(p%UnLog, '(A12,f12.4)') " Cl : ", m%LineTypeList(l)%Cl + write(p%UnLog, '(A12,f12.4)') " dF : ", m%LineTypeList(l)%dF + write(p%UnLog, '(A12,f12.4)') " cF : ", m%LineTypeList(l)%cF end if IF ( ErrStat2 /= ErrID_None ) THEN @@ -1459,14 +1484,18 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! account for states of line m%LineStateIs1(l) = Nx + 1 + Nx = Nx + 6*m%LineList(l)%N - 6 ! normal case, just 6 states per internal node + if (m%LineTypeList(m%LineList(l)%PropsIdNum)%ElasticMod > 1) then ! todo add an error check here? or change to 2 or 3? - Nx = Nx + 7*m%LineList(l)%N - 6 ! if using viscoelastic model, need one more state per segment - m%LineStateIsN(l) = Nx - else - Nx = Nx + 6*m%LineList(l)%N - 6 ! normal case, just 6 states per internal node - m%LineStateIsN(l) = Nx + Nx = Nx + m%LineList(l)%N ! if using viscoelastic model, need one more state per segment end if + if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 1) then + Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). TODO: when we change to only internal nodes need to make N - 1. + endif + + m%LineStateIsN(l) = Nx + ! Process attachment identfiers and attach line ends ! First for the anchor (or end A)... @@ -1559,9 +1588,10 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er IF ( scan( LineOutString, 'U') > 0 ) m%LineList(l)%OutFlagList(4) = 1 IF ( scan( LineOutString, 'D') > 0 ) m%LineList(l)%OutFlagList(5) = 1 IF ( scan( LineOutString, 'b') > 0 ) m%LineList(l)%OutFlagList(6) = 1 ! seabed contact forces + IF ( scan( LineOutString, 'V') > 0 ) m%LineList(l)%OutFlagList(7) = 1 ! VIV forces ! per node 1 component - IF ( scan( LineOutString, 'W') > 0 ) m%LineList(l)%OutFlagList(7) = 1 ! node weight/buoyancy (positive up) - IF ( scan( LineOutString, 'K') > 0 ) m%LineList(l)%OutFlagList(8) = 1 ! curvature at node + IF ( scan( LineOutString, 'W') > 0 ) m%LineList(l)%OutFlagList(8) = 1 ! node weight/buoyancy (positive up) + IF ( scan( LineOutString, 'K') > 0 ) m%LineList(l)%OutFlagList(9) = 1 ! curvature at node ! per element 1 component IF ( scan( LineOutString, 't') > 0 ) m%LineList(l)%OutFlagList(10) = 1 ! segment tension force (just EA) IF ( scan( LineOutString, 'c') > 0 ) m%LineList(l)%OutFlagList(11) = 1 ! segment internal damping force @@ -1583,7 +1613,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er write(p%UnLog, '(A15,f12.4)') " Len : ", m%LineList(l)%UnstrLen write(p%UnLog, '(A15,A)' ) " Node A : ", " "//tempString2 write(p%UnLog, '(A15,A)' ) " Node B : ", " "//tempString3 - write(p%UnLog, '(A15,I2)' ) " NumSegs: ", m%LineList(l)%N + write(p%UnLog, '(A15,I4)' ) " NumSegs: ", m%LineList(l)%N end if ! check for sequential IdNums @@ -1606,6 +1636,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er RETURN END IF + IF (wordy > 0) print *, "Set up Line", l, "of type", m%LineList(l)%PropsIdNum + END DO ! l = 1,p%nLines !------------------------------------------------------------------------------------------- @@ -2244,7 +2276,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! prepare state vector etc. !------------------------------------------------------------------------------------ - ! the number of states is Nx and Nxtra includes additional states for potential line failures + ! the number of states is Nx and Nxtra includes additional states (points) for potential line failures m%Nx = Nx m%Nxtra = m%Nx + 6*2*p%nLines @@ -2656,6 +2688,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! boost drag coefficient of each line type <<<<<<<< does this actually do anything or do lines hold these coefficients??? 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 @@ -2708,10 +2741,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er !loop through line integration time steps DO J = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) - Call MD_Step(t, dtM, u_interp, u_array, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) - IF ( ErrStat2 /= ErrID_None ) THEN - CALL CheckError(ErrStat2, ErrMsg2) - END IF + Call MD_Step(t, dtM, u_interp, u_array, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) ; if(Failed()) return ! check for NaNs - is this a good place/way to do it? DO K = 1, m%Nx @@ -2816,6 +2846,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er ! UNboost drag coefficient of each line type <<< 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 @@ -2896,7 +2927,7 @@ END FUNCTION Failed SUBROUTINE CheckError(ErrID,Msg) - ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev + ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLe ! Passed arguments INTEGER(IntKi), INTENT(IN) :: ErrID ! The error identifier (ErrStat) @@ -2949,7 +2980,7 @@ CHARACTER(1024) function NextLine(i) NextLine="---" ! Set as a separator so we can escape some of the while loops else NextLine=trim(FileInfo_In%Lines(i)) - !TODO: add comment character recognition here? (discard any characters past a #) + ! # is comment character handled by file loading stuff earlier on (in NWTC routine ProcessComFile) endif end function NextLine @@ -3036,7 +3067,7 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er !loop through line integration time steps DO I = 1, NdtM ! for (double ts=t; ts<=t+ICdt-dts; ts+=dts) - Call MD_Step(t2, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) + Call MD_Step(t2, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat2, ErrMsg2) IF ( ErrStat2 /= ErrID_None ) THEN CALL CheckError(ErrStat2, ErrMsg2) END IF @@ -3045,13 +3076,13 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er DO J = 1, m%Nx IF (Is_NaN(x%states(J))) THEN ErrStat = ErrID_Fatal - ErrMsg = ' NaN state detected.' + ErrMsg = ' NaN state detected at time '//TRIM(Num2LStr(t2)) IF (wordy > 1) THEN print *, ". Here is the state vector: " print *, x%states END IF CALL CheckError(ErrStat, ErrMsg) - EXIT + RETURN END IF END DO @@ -3068,22 +3099,7 @@ SUBROUTINE MD_UpdateStates( t, n, u, t_array, p, x, xd, z, other, m, ErrStat, Er END IF ! CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MD_UpdateStates') - - ! check for NaNs - is this a good place/way to do it? - DO J = 1, m%Nx - IF (Is_NaN(x%states(J))) THEN - ErrStat = ErrID_Fatal - ErrMsg = ' NaN state detected.' - IF (wordy > 1) THEN - print *, ". Here is the state vector: " - print *, x%states - END IF - CALL CheckError(ErrStat, ErrMsg) - EXIT - END IF - END DO - - ! do we want to check failures here (at the coupling step level? Or at the dtM level?) + ! do we want to check failures here (at the coupling step level? Or at the dtM level?) TODO: move this to the dtM level ! --------------- check for line failures (detachments!) ---------------- DO l= 1,p%nFails @@ -3905,9 +3921,10 @@ SUBROUTINE MD_Step ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrSta ELSE ErrStat = ErrID_Fatal ErrMsg = ' Unrecognized tScheme option in MD_Step' - RETURN ENDIF + IF ( ErrStat >= AbortErrLev ) RETURN + END SUBROUTINE MD_Step ! RK2 integrator (part of what was in TimeStep) @@ -3940,9 +3957,11 @@ SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat ! step 1 CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) - ! k0 = m%xdTemp + if ( ErrStat >= AbortErrLev ) return + ! k0 = m%xdTemp DO J = 1, m%Nx m%xTemp%states(J) = x%states(J) + 0.5_DbKi*dtM*m%xdTemp%states(J) !x1 = x0 + dt*k0/2.0; END DO @@ -3950,8 +3969,10 @@ SUBROUTINE MD_RK2 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat ! step 2 CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t+0.5*dtM) - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) !called with updated states x2 and time = t + dt/2.0 + if ( ErrStat >= AbortErrLev ) return ! k1 = m%xdTemp DO J = 1, m%Nx x%states(J) = x%states(J) + dtM*m%xdTemp%states(J) @@ -3992,8 +4013,10 @@ SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat ! k0 CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t , ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( t, u_interp, p, x, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + if ( ErrStat >= AbortErrLev ) return ! k0 = m%xdTemp ! k1 @@ -4003,8 +4026,10 @@ SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat END DO CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + if ( ErrStat >= AbortErrLev ) return ! k1 = m%xdTemp ! k2 @@ -4014,8 +4039,10 @@ SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat END DO CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + 0.5_DbKi*dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) TODO: is this needed, it is already called for k1 - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( (t + 0.5_DbKi*dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + if ( ErrStat >= AbortErrLev ) return ! k2 = m%xdTemp ! k3 @@ -4025,8 +4052,10 @@ SUBROUTINE MD_RK4 ( t, dtM, u_interp, u, t_array, p, x, xd, z, other, m, ErrStat END DO CALL MD_Input_ExtrapInterp(u, t_array, u_interp, t + dtM, ErrStat, ErrMsg) ! interpolate input mesh to correct time (t) - + if ( ErrStat >= AbortErrLev ) return + CALL MD_CalcContStateDeriv( (t + dtM), u_interp, p, m%xTemp, xd, z, other, m, m%xdTemp, ErrStat, ErrMsg ) + if ( ErrStat >= AbortErrLev ) return ! k3 = m%xdTemp ! Apply diff --git a/modules/moordyn/src/MoorDyn_IO.f90 b/modules/moordyn/src/MoorDyn_IO.f90 index 4eb8a6e2a0..efe59ca038 100644 --- a/modules/moordyn/src/MoorDyn_IO.f90 +++ b/modules/moordyn/src/MoorDyn_IO.f90 @@ -785,8 +785,8 @@ SUBROUTINE MDIO_ProcessOutList(OutList, p, m, y, InitOut, ErrStat, ErrMsg ) ! calculate number of output entries (excluding time) to write for this line - LineNumOuts = 3*(m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(2:6)) & - + (m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(7:9)) & + LineNumOuts = 3*(m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(2:7)) & + + (m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(8:9)) & + m%LineList(I)%N*SUM(m%LineList(I)%OutFlagList(10:18)) ALLOCATE(m%LineList(I)%LineWrOutput( 1 + LineNumOuts), STAT = ErrStat) @@ -934,8 +934,8 @@ SUBROUTINE MDIO_OpenOutput( MD_ProgDesc, p, m, InitOut, ErrStat, ErrMsg ) ! calculate number of output entries (excluding time) to write for this line - LineNumOuts = 3*(m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(2:6)) & - + (m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(7:9)) & + LineNumOuts = 3*(m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(2:7)) & + + (m%LineList(I)%N + 1)*SUM(m%LineList(I)%OutFlagList(8:9)) & + m%LineList(I)%N*SUM(m%LineList(I)%OutFlagList(10:18)) if (wordy > 2) PRINT *, LineNumOuts, " output channels" @@ -966,12 +966,17 @@ SUBROUTINE MDIO_OpenOutput( MD_ProgDesc, p, m, InitOut, ErrStat, ErrMsg ) WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((3+3*m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, 'Node'//TRIM(Int2Lstr(J))//'bx', p%Delim, 'Node'//TRIM(Int2Lstr(J))//'by', p%Delim, 'Node'//TRIM(Int2Lstr(J))//'bz', J=0,(m%LineList(I)%N) ) END IF - + IF (m%LineList(I)%OutFlagList(7) == 1) THEN + WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((3+3*m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & + ( p%Delim, 'Node'//TRIM(Int2Lstr(J))//'Vx', p%Delim, 'Node'//TRIM(Int2Lstr(J))//'Vy', p%Delim, 'Node'//TRIM(Int2Lstr(J))//'Vz', J=0,(m%LineList(I)%N) ) ! TODO adjust these when force to internal nodes + END IF + + IF (m%LineList(I)%OutFlagList(8) == 1) THEN WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, 'Node'//TRIM(Int2Lstr(J))//'Wz', J=0,(m%LineList(I)%N) ) END IF - IF (m%LineList(I)%OutFlagList(8) == 1) THEN + IF (m%LineList(I)%OutFlagList(9) == 1) THEN WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, 'Node'//TRIM(Int2Lstr(J))//'Kurv', J=0,(m%LineList(I)%N) ) END IF @@ -1023,12 +1028,16 @@ SUBROUTINE MDIO_OpenOutput( MD_ProgDesc, p, m, InitOut, ErrStat, ErrMsg ) WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((3+3*m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, '(N)', p%Delim, '(N)', p%Delim, '(N)', J=0,(m%LineList(I)%N) ) END IF - IF (m%LineList(I)%OutFlagList(7) == 1) THEN + WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((3+3*m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & + ( p%Delim, '(N)', p%Delim, '(N)', p%Delim, '(N)', J=0,(m%LineList(I)%N) ) + END IF + + IF (m%LineList(I)%OutFlagList(8) == 1) THEN WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, '(Nup)', J=0,(m%LineList(I)%N) ) END IF - IF (m%LineList(I)%OutFlagList(8) == 1) THEN + IF (m%LineList(I)%OutFlagList(9) == 1) THEN WRITE(m%LineList(I)%LineUnOut,'('//TRIM(Int2LStr((m%LineList(I)%N)))//'(A1,A15))', advance='no', IOSTAT=ErrStat2) & ( p%Delim, '(1/m)', J=0,(m%LineList(I)%N) ) END IF @@ -1669,22 +1678,31 @@ SUBROUTINE MDIO_WriteOutputs( Time, p, m, y, ErrStat, ErrMsg ) END DO END IF + ! Node VIV force + IF (m%LineList(I)%OutFlagList(7) == 1) THEN + DO J = 0,m%LineList(I)%N + DO K = 1,3 + m%LineList(I)%LineWrOutput(L) = m%LineList(I)%Lf(K,J) + L = L+1 + END DO + END DO + END IF ! Node weights - IF (m%LineList(I)%OutFlagList(7) == 1) THEN + IF (m%LineList(I)%OutFlagList(8) == 1) THEN DO J = 0,m%LineList(I)%N m%LineList(I)%LineWrOutput(L) = m%LineList(I)%W(3,J) L = L+1 END DO END IF - ! ! Node curvatures - ! IF (m%LineList(I)%OutFlagList(8) == 1) THEN - ! DO J = 0,m%LineList(I)%N - ! m%LineList(I)%LineWrOutput(L) = m%LineList(I)%W(3,J) - ! L = L+1 - ! END DO - ! END IF + ! Node curvatures + IF (m%LineList(I)%OutFlagList(9) == 1) THEN + DO J = 0,m%LineList(I)%N + m%LineList(I)%LineWrOutput(L) = m%LineList(I)%Kurv(J) + L = L+1 + END DO + END IF ! Segment tension force (excludes damping term, just EA) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 20c90417b8..65f9bdf87c 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -29,7 +29,7 @@ MODULE MoorDyn_Line PRIVATE - INTEGER(IntKi), PARAMETER :: wordy = 0 ! verbosity level. >1 = more console output + INTEGER(IntKi), PARAMETER :: wordy = 3 ! verbosity level. >1 = more console output PUBLIC :: SetupLine PUBLIC :: Line_Initialize @@ -80,7 +80,12 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%Can = LineProp%Can Line%Cat = LineProp%Cat Line%Cdn = LineProp%Cdn - Line%Cdt = LineProp%Cdt + Line%Cdt = LineProp%Cdt + Line%Cl = LineProp%Cl + Line%dF = LineProp%dF + Line%cF = LineProp%cF + + Line%n_m = 500 ! hardcode n_m to 500 (for VIV) ! copy over elasticity data Line%ElasticMod = LineProp%ElasticMod @@ -120,7 +125,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! - highest axial vibration mode of a segment is wn = sqrt(k/m) = 2N/UnstrLen*sqrt(EA/w) Line%BA = -LineProp%BA * Line%UnstrLen / Line%N * SQRT(LineProp%EA * LineProp%w) if (p%writeLog > 0) then - write(p%UnLog,'(A)') "Based on -zeta of "//trim(Num2LStr(LineProp%BA))//", BA set to "//trim(Num2LStr(Line%BA)) + write(p%UnLog,'(A)') " "//"Based on -zeta of "//trim(Num2LStr(LineProp%BA))//", BA set to "//trim(Num2LStr(Line%BA)) ! extra space at front to make nice formatting in log file end if IF (wordy > 1) print *, 'Based on zeta, BA set to ', Line%BA @@ -162,6 +167,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! if using viscoelastic model, allocate additional state quantities if (Line%ElasticMod > 1) then + if (wordy > 1) print *, "Using the viscoelastic model" ALLOCATE ( Line%dl_1(N), STAT = ErrStat ) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating dl_1 array.' @@ -171,6 +177,22 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! initialize to zero Line%dl_1 = 0.0_DbKi end if + + ! if using viscoelastic model, allocate additional state quantities + if (Line%Cl > 0) then + if (wordy > 1) print *, "Using the VIV model" + ! allocate old acclerations [for VIV] + ALLOCATE ( Line%phi(0:N), Line%rdd_old(3,0:N), Line%yd_rms_old(0:N), Line%ydd_rms_old(0:N), STAT = ErrStat ) + IF ( ErrStat /= ErrID_None ) THEN + ErrMsg = ' Error allocating VIV arrays.' + !CALL CleanUp() + RETURN + END IF + ! initialize to unique values on the range 0-2pi + do I=0, Line%N + Line%phi(I) = (I/Line%N)*2*Pi + enddo + end if ! allocate node and segment tangent vectors ALLOCATE ( Line%q(3, 0:N), Line%qs(3, N), STAT = ErrStat ) @@ -213,8 +235,9 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) END IF ! allocate node force vectors - ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), Line%Ap(3, 0:N), & - Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), Line%Fnet(3, 0:N), STAT = ErrStat ) + ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), & + Line%Ap(3, 0:N), Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), & + Line%Lf(3, 0:N), Line%Fnet(3, 0:N), STAT = ErrStat ) ! TODO: LF only to internal nodes after testing IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating node force arrays.' !CALL CleanUp() @@ -996,6 +1019,7 @@ SUBROUTINE Line_SetState(Line, X, t) INTEGER(IntKi) :: i ! index of segments or nodes along line INTEGER(IntKi) :: J ! index + ! print*, "begining set state" ! store current time Line%time = t @@ -1013,9 +1037,20 @@ SUBROUTINE Line_SetState(Line, X, t) ! if using viscoelastic model, also set the static stiffness stretch if (Line%ElasticMod > 1) then do I=1,Line%N - Line%dl_1(I) = X( 6*Line%N-6 + I) ! these will be the last N entries in the state vector + Line%dl_1(I) = X( 6*Line%N-6 + I) ! these will be the N entries in the state vector passed the internal node states end do end if + + ! if using the viv mdodel, also set the lift force phase + if (Line%Cl > 0) then + do I=0, Line%N ! TODO: after checking change to internal + 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 + else ! if only VIV state, then N+1 entries after internal node states + Line%phi(I) = X( 6*Line%N-6 + I + 1) - (2 * Pi * floor(X( 6*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 + endif + enddo + endif END SUBROUTINE Line_SetState !-------------------------------------------------------------- @@ -1103,6 +1138,13 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Real(DbKi) :: FfT(3) ! total friction force in the transverse direction Real(DbKi) :: FfA(3) ! total friction force in the axial direction Real(DbKi) :: Ff(3) ! total friction force on the line node + ! VIV stuff + Real(DbKi) :: yd ! Crossflow velocity + Real(DbKi) :: ydd ! Crossflow acceleration + Real(DbKi) :: yd_rms ! Rolling rms of crossflow velocity + Real(DbKi) :: ydd_rms ! Rolling rms of crossflow acceleration + Real(DbKi) :: phi_dot ! frequency of lift force (rad/s) + Real(DbKi) :: f_hat ! non-dimensional frequency N = Line%N ! for convenience @@ -1119,7 +1161,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Line%rd(J,0) = m%PointList(Line%AnchPoint)%rd(J) ! END DO - + ! print*, "begining get state deriv" ! -------------------- calculate various kinematic quantities --------------------------- DO I = 1, N @@ -1300,9 +1342,24 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai else if (Line%ElasticMod > 1) then if (Line%ElasticMod == 3) then - if (Line%dl_1(I) >= 0.0) then + if (Line%dl_1(I) > 0.0) then ! Mean load dependent dynamic stiffness: from combining eqn. 2 and eqn. 10 from original MD viscoelastic paper, taking mean load = k1 delta_L1 / MBL, and solving for k_D using WolframAlpha with following conditions: k_D > k_s, (MBL,alpha,beta,unstrLen,delta_L1) > 0 EA_D = 0.5 * ((Line%alphaMBL) + (Line%vbeta*Line%dl_1(I)*(Line%EA / Line%l(I))) + Line%EA + sqrt((Line%alphaMBL * Line%alphaMBL) + (2*Line%alphaMBL*(Line%EA / Line%l(I)) * (Line%vbeta*Line%dl_1(I) - Line%l(I))) + ((Line%EA / Line%l(I))*(Line%EA / Line%l(I)) * (Line%vbeta*Line%dl_1(I) + Line%l(I))*(Line%vbeta*Line%dl_1(I) + Line%l(I))))) + + ! Double check none of the assumptions were violated (this should never happen) + IF (Line%alphaMBL <= 0 .OR. Line%vbeta <= 0 .OR. Line%l(I) <= 0 .OR. Line%dl_1(I) <= 0 .OR. EA_D < Line%EA) THEN + ErrStat = ErrID_Warn + ErrMsg = "Viscoelastic model: Assumption for mean laod dependent dynamic stiffness violated" + if (wordy > 2) then + print *, "Line%alphaMBL", Line%alphaMBL + print *, "Line%vbeta", Line%vbeta + print *, "Line%l(I)", Line%l(I) + print *, "Line%dl_1(I)", Line%dl_1(I) + print *, "EA_D", EA_D + print *, "Line%EA", Line%EA + endif + ENDIF + else EA_D = Line%alphaMBL ! mean load is considered to be 0 in this case. The second term in the above equation is not valid for delta_L1 < 0. endif @@ -1336,7 +1393,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai MagTd = Line%BA*ld_1 / Line%l(I) ! compute tension based on static portion (dynamic portion would give same). See eqn. 14 in paper - ! update state derivative for static stiffness stretch (last N entries in the state vector) + ! update state derivative for static stiffness stretch (last N entries in the line state vector if no VIV model, otherwise N entries past the 6*N-6 entries in this vector) Xd( 6*N-6 + I) = ld_1 end if @@ -1449,14 +1506,14 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! relative flow velocities DO J = 1, 3 - Vi(J) = Line%U(J,I) - Line%rd(J,I) ! relative flow velocity over node -- this is where wave velicites would be added + Vi(J) = Line%U(J,I) - Line%rd(J,I) ! relative flow velocity over node -- this is where wave velocities would be added END DO ! decomponse relative flow into components SumSqVp = 0.0_DbKi ! start sums of squares at zero SumSqVq = 0.0_DbKi DO J = 1, 3 - Vq(J) = DOT_PRODUCT( Vi , Line%q(:,I) ) * Line%q(J,I); ! tangential relative flow component + Vq(J) = DOT_PRODUCT( Vi , Line%q(:,I) ) * Line%q(J,I) ! tangential relative flow component Vp(J) = Vi(J) - Vq(J) ! transverse relative flow component SumSqVq = SumSqVq + Vq(J)*Vq(J) SumSqVp = SumSqVp + Vp(J)*Vp(J) @@ -1476,9 +1533,126 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%Dq(:,I) = 0.25*p%rhoW*Line%Cdt* Pi*d*(Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) * MagVq * Vq END IF + ! print *, "Made it to VIV model" + ! Vortex Induced Vibration (VIV) cross-flow lift force + IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen then VIV to be calculated + + ! ----- The Synchronization Model ------ + ! Crossflow velocity and acceleration. rd component in the crossflow direction + + yd = dot_product(Line%rd(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) + ydd = dot_product(Line%rdd_old(0:2,I), cross_product(Line%q(:,I), normalize(Vp) ) ) + + ! ! for checking rdd_old fix + ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. (I > 95 .or. I <5)) then + ! print*, "rdd_old at t = ", Line%time + ! print*, "I =", I, "rdd_old(0:2,I) =", Line%rdd_old(0:2,I) + ! endif + + ! print*, "crossflow stuff done" + + ! Rolling RMS calculation + yd_rms = sqrt((((Line%n_m-1) * Line%yd_rms_old(I) * Line%yd_rms_old(I)) + (yd * yd)) / Line%n_m) ! RMS approximation from Thorsen + ydd_rms = sqrt((((Line%n_m-1) * Line%ydd_rms_old(I) * Line%ydd_rms_old(I)) + (ydd * ydd)) / Line%n_m) + + ! print*, "end RMS calc" + + IF ((Line%time >= Line%t_old + p%dtM0) .OR. (Line%time == 0.0)) THEN ! Update the stormed RMS vaues + ! update back indexing one moordyn time step (regardless of time integration scheme or coupling step size). T_old is handled at end of getStateDeriv when rdd_old is updated. + Line%yd_rms_old(I) = yd_rms ! for rms back indexing (one moordyn timestep back) + Line%ydd_rms_old(I) = ydd_rms ! for rms back indexing (one moordyn timestep back) + ENDIF + + IF ((yd_rms==0.0) .OR. (ydd_rms == 0.0)) THEN + Line%phi_yd = atan2(-ydd, yd) ! To avoid divide by zero + ELSE + Line%phi_yd = atan2(-ydd/ydd_rms, yd/yd_rms) + ENDIF + + IF (Line%phi_yd < 0) THEN + Line%phi_yd = 2*Pi + Line%phi_yd ! atan2 to 0-2Pi range + ENDIF + + ! print*, "end phase stuff" + + ! Note: amplitude calculations and states commented out. Would be needed if a Cl vs A lookup table was ever implemented + + ! const real A_int = Misc(I)[1]; + ! const real As = Misc(I)[2]; + + ! non-dimensional frequency + f_hat = Line%cF + Line%dF *sin(Line%phi_yd - Line%phi(I)) ! phi is integrated from state deriv phi_dot + ! frequency of lift force (rad/s) + phi_dot = 2*Pi*f_hat*MagVp / d ! to be added to state + + ! ----- The rest of the model ----- + + ! ! Oscillation amplitude + ! const real A_int_dot = abs(yd); + ! ! Note: Check if this actually measures zero crossings + ! if ((yd * yd_old[i]) < 0) { ! if sign changed, i.e. a zero crossing + ! Amp[i] = A_int-A_int_old[i]; ! amplitude calculation since last zero crossing + ! A_int_old[i] = A_int; ! stores amplitude of previous zero crossing for finding Amp + ! } + ! ! Careful with integrating smoothed amplitude, as 0.1 was a calibarated value based on a very simple integration method + ! const real As_dot = (0.1/dtm)*(Amp[i]-As); ! As to be variable integrated from the state. stands for amplitude smoothed + + ! ! Lift coefficient from lookup table + ! const real C_l = cl_lookup(x = As/d); ! create function in Line.hpp that uses lookup table + + ! The Lift force + IF (I==0) THEN ! TODO: drop for only internal nodes + Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(1)*Line%l(1) + ELSE IF (I==N) THEN ! TODO: drop for only internal nodes + Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(N)*Line%l(N) + ELSE + Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * (Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) + END IF + + if (wordy > 0) then + if (Is_NaN(norm2(Line%Lf(:,I)))) print*, "Lf nan at node", I, "for line", Line%IdNum + endif + + ! finding nans + if (Line%time == 2*p%dtM0 .and. I==N-2) then + print *, "-------" + print *, "Line%time", Line%time + print *, "I", I + print *, "Line%Lf(:,I)", Line%Lf(:,I) + print *, "MagVp", MagVp + print *, "Line%U(:,I)", Line%U(:,I) + print *, "Line%phi(I)", Line%phi(I) + print *, "phi_dot", phi_dot + print *, "f_hat", f_hat + print *, "Line%phi_yd", Line%phi_yd + print *, "yd", yd + print *, "ydd", ydd + print *, "Line%rd(:,I)", Line%rd(:,I) + print *, "Line%rdd_old(0:2,I)", Line%rdd_old(0:2,I) + print *, "Line%q(:,I)", Line%q(:,I) + print *, "Vp", Vp + endif + + + ! print*, "end force" + ! update state derivative with lift force frequency + + if (Line%ElasticMod > 1) then ! if both additional states are included then N+1 entries after internal node states and visco segment states + Xd( 7*Line%N-6 + I + 1) = phi_dot + else ! if only VIV state, then N+1 entries after internal node states + Xd( 6*Line%N-6 + I + 1) = phi_dot + endif + + ! print*, "end update state deriv" + ! Miscd(I)[2] = As_dot; ! unused state that could be used for future amplitude calculations + + ENDIF + + ! print*, "made it past VIV model" + ! ------ fluid acceleration components for current node (from MD-C) ------ DO J = 1, 3 - aq(J) = DOT_PRODUCT( Line%Ud(:,I) , Line%q(:,I) ) * Line%q(J,I); ! tangential fluid acceleration component + aq(J) = DOT_PRODUCT( Line%Ud(:,I) , Line%q(:,I) ) * Line%q(J,I) ! tangential fluid acceleration component ap(J) = Line%Ud(J,I) - aq(J) ! transverse fluid acceleration component ENDDO @@ -1570,14 +1744,16 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! total forces IF (I==0) THEN - Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ! TODO: remove LF here ELSE IF (I==N) THEN - Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ! TODO: remove LF here ELSE - Line%Fnet(:,I) = Line%T(:,I+1) - Line%T(:,I) + Line%Td(:,I+1) - Line%Td(:,I) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Fnet(:,I) = Line%T(:,I+1) - Line%T(:,I) + Line%Td(:,I+1) - Line%Td(:,I) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) END IF END DO ! I - done looping through nodes + + ! print *, "Done looping through nodes. Time to find acceleration" ! loop through internal nodes and update their states <<< should/could convert to matrix operations instead of all these loops DO I=1, N-1 @@ -1593,9 +1769,33 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Xd(3*N-3 + 3*I-3 + J) = Line%rd(J,I); ! dxdt = V (velocities) Xd( 3*I-3 + J) = Sum1 ! dVdt = RHS * A (accelerations) + IF (Line%Cl > 0) THEN + ! ! for checking rdd_old fix + ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. I > 95) then + ! print*, "I in the rdd_old loop", I + ! print*, "Sum1 for above I at J =", J, "is", Sum1 + ! endif + Line%rdd_old(J-1,I-1) = Sum1 ! saving the acceleration for VIV RMS calculation. WARNING: I-1 is intential to match the incorrect approach in MD-C (map internal node accel to the first N-2 elements of rdd_old.) After testing REMOVE the I-1 and make the VIV model only apply to internal nodes. J-1 here to shift things correctly. Not sure why we need this but it doesnt work otherwise ... + ENDIF + END DO ! J END DO ! I + if ((Line%time >= Line%t_old + p%dtM0) .OR. (Line%time == 0.0)) then ! update back indexing one moordyn time step (regardless of time integration scheme) + Line%t_old = Line%time ! for updating back indexing if statements + endif + + ! ! for checking rdd_old fix + ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen) then + ! print*, "rdd_old at t = ", Line%time + ! DO I = 0, 4 + ! print*, "I =", I, "rdd_old =", Line%rdd_old(0,I), Line%rdd_old(1,I), Line%rdd_old(2,I) + ! enddo + ! print*, "..." + ! DO I = N-4, N + ! print*, "I =", I, "rdd_old =", Line%rdd_old(0,I), Line%rdd_old(1,I), Line%rdd_old(2,I) + ! enddo + ! endif ! check for NaNs DO J = 1, 6*(N-1) @@ -1646,6 +1846,27 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! END DO ! END DO + CONTAINS + + FUNCTION normalize(vector) RESULT(normalized) + ! function to normalize a vector. Returns the input vector if the magnitude is sufficiently small + + REAL(DbKi), DIMENSION(3), INTENT(IN ) :: vector ! The input array + REAL(DbKi), DIMENSION(3) :: normalized ! The normalized array + REAL(DbKi) :: mag ! the magnitude + + mag = norm2(vector) + + if (mag > 0.0_DbKi) then + DO J=1,3 + normalized(J) = vector(J) / mag + ENDDO + else + normalized = vector + endif + + END FUNCTION normalize + END SUBROUTINE Line_GetStateDeriv !===================================================================== diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 289dc4ac13..95adb6df1d 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -64,6 +64,9 @@ typedef ^ ^ DbKi Can - typedef ^ ^ DbKi Cat - - - "tangential added mass coefficient" typedef ^ ^ DbKi Cdn - - - "transverse drag coefficient" typedef ^ ^ DbKi Cdt - - - "tangential drag coefficient" +typedef ^ ^ DbKi Cl - - - "VIV lift coefficient. If 0, VIV turned off" +typedef ^ ^ DbKi dF - - - "+- range of VIV synchronization in non-dimensional frequency" +typedef ^ ^ DbKi cF - - - "Center VIV synchronization in non-dimensional frequency" typedef ^ ^ IntKi ElasticMod - - - "Which elasticity model to use: {1 basic, 2 viscoelastic, 3 viscoelastic+meanload} " - typedef ^ ^ IntKi nEApoints - - - "number of values in stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi stiffXs {30} - - "x array for stress-strain lookup table (up to nCoef)" @@ -219,9 +222,9 @@ typedef ^ ^ DbKi Bquad {2} # this is the Line type, which holds data for each line object typedef ^ MD_Line IntKi IdNum - - - "integer identifier of this Line" -#typedef ^ ^ CHARACTER(10) type - - - "type of line. should match one of LineProp names" +#typedef ^ ^ CHARACTER(10) type - - - "type of line. should match one of LineProp names" typedef ^ ^ IntKi PropsIdNum - - - "the IdNum of the associated line properties" - -typedef ^ ^ IntKi ElasticMod - - - "Which elasticity model to use: {0 basic, 1 viscoelastic, 2 future SYCOM} " - +typedef ^ ^ IntKi ElasticMod - - - "Which elasticity model to use: {1 basic, 2 viscoelastic, 3 viscoelastic+meanload} " - typedef ^ ^ IntKi OutFlagList {20} - - "array specifying what line quantities should be output (1 vs 0)" - typedef ^ ^ IntKi CtrlChan - - - "index of control channel that will drive line active tensioning (0 for none)" - typedef ^ ^ IntKi FairPoint - - - "IdNum of Point at fairlead" @@ -243,6 +246,9 @@ typedef ^ ^ DbKi Can - typedef ^ ^ DbKi Cat - - - "" "[-]" typedef ^ ^ DbKi Cdn - - - "" "[-]" typedef ^ ^ DbKi Cdt - - - "" "[-]" +typedef ^ ^ DbKi Cl - - - "VIV lift coefficient. If 0, VIV turned off" +typedef ^ ^ DbKi dF - - - "+- range of VIV synchronization in non-dimensional frequency" +typedef ^ ^ DbKi cF - - - "Center VIV synchronization in non-dimensional frequency" typedef ^ ^ IntKi nEApoints - - - "number of values in stress-strain lookup table (0 means using constant E)" typedef ^ ^ DbKi stiffXs {30} - - "x array for stress-strain lookup table (up to nCoef)" typedef ^ ^ DbKi stiffYs {30} - - "y array for stress-strain lookup table" @@ -278,6 +284,7 @@ typedef ^ ^ DbKi Ap {:}{:} typedef ^ ^ DbKi Aq {:}{:} - - "node added mass forcing (axial)" "[N]" typedef ^ ^ DbKi B {:}{:} - - "node bottom contact force" "[N]" typedef ^ ^ DbKi Bs {:}{:} - - "node force due to bending moments" "[N]" +typedef ^ ^ DbKi Lf {:}{:} - - "node crossflow viv lift force. TODO this will have to be allocated for only internal after fix" "[N]" typedef ^ ^ DbKi Fnet {:}{:} - - "total force on node" "[N]" typedef ^ ^ DbKi S {:}{:}{:} - - "node inverse mass matrix" "[kg]" typedef ^ ^ DbKi M {:}{:}{:} - - "node mass matrix" "[kg]" @@ -285,6 +292,16 @@ typedef ^ ^ DbKi EndMomentA {3} typedef ^ ^ DbKi EndMomentB {3} - - "vector of end moments due to bending at line end B" "[N-m]" typedef ^ ^ IntKi LineUnOut - - - "unit number of line output file" 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" +typedef ^ ^ DbKi yd_rms_old {:} - - "node old cf vel rms" "m/s" +typedef ^ ^ DbKi ydd_rms_old {:} - - "node old cf accel rms" "m/s^2" +typedef ^ ^ DbKi rdd_old {:}{:} - - "node accelerations previous iteration" "m/s^2" + # this is the ExtLd type, which holds data for each external load specification typedef ^ MD_ExtLd IntKi IdNum - - - "integer identifier of this external load entry" diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index cba62a494a..63eba5f6ba 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -77,6 +77,9 @@ MODULE MoorDyn_Types REAL(DbKi) :: Cat = 0.0_R8Ki !< tangential added mass coefficient [-] REAL(DbKi) :: Cdn = 0.0_R8Ki !< transverse drag coefficient [-] REAL(DbKi) :: Cdt = 0.0_R8Ki !< tangential drag coefficient [-] + REAL(DbKi) :: Cl = 0.0_R8Ki !< VIV lift coefficient. If 0, VIV turned off [-] + REAL(DbKi) :: dF = 0.0_R8Ki !< +- range of VIV synchronization in non-dimensional frequency [-] + REAL(DbKi) :: cF = 0.0_R8Ki !< Center VIV synchronization in non-dimensional frequency [-] INTEGER(IntKi) :: ElasticMod = 0_IntKi !< Which elasticity model to use: {1 basic, 2 viscoelastic, 3 viscoelastic+meanload} [-] INTEGER(IntKi) :: nEApoints = 0_IntKi !< number of values in stress-strain lookup table (0 means using constant E) [-] REAL(DbKi) , DIMENSION(1:30) :: stiffXs = 0.0_R8Ki !< x array for stress-strain lookup table (up to nCoef) [-] @@ -241,7 +244,7 @@ MODULE MoorDyn_Types TYPE, PUBLIC :: MD_Line INTEGER(IntKi) :: IdNum = 0_IntKi !< integer identifier of this Line [-] INTEGER(IntKi) :: PropsIdNum = 0_IntKi !< the IdNum of the associated line properties [-] - INTEGER(IntKi) :: ElasticMod = 0_IntKi !< Which elasticity model to use: {0 basic, 1 viscoelastic, 2 future SYCOM} [-] + INTEGER(IntKi) :: ElasticMod = 0_IntKi !< Which elasticity model to use: {1 basic, 2 viscoelastic, 3 viscoelastic+meanload} [-] INTEGER(IntKi) , DIMENSION(1:20) :: OutFlagList = 0_IntKi !< array specifying what line quantities should be output (1 vs 0) [-] INTEGER(IntKi) :: CtrlChan = 0_IntKi !< index of control channel that will drive line active tensioning (0 for none) [-] INTEGER(IntKi) :: FairPoint = 0_IntKi !< IdNum of Point at fairlead [-] @@ -263,6 +266,9 @@ MODULE MoorDyn_Types REAL(DbKi) :: Cat = 0.0_R8Ki !< [[-]] REAL(DbKi) :: Cdn = 0.0_R8Ki !< [[-]] REAL(DbKi) :: Cdt = 0.0_R8Ki !< [[-]] + REAL(DbKi) :: Cl = 0.0_R8Ki !< VIV lift coefficient. If 0, VIV turned off [-] + REAL(DbKi) :: dF = 0.0_R8Ki !< +- range of VIV synchronization in non-dimensional frequency [-] + REAL(DbKi) :: cF = 0.0_R8Ki !< Center VIV synchronization in non-dimensional frequency [-] INTEGER(IntKi) :: nEApoints = 0_IntKi !< number of values in stress-strain lookup table (0 means using constant E) [-] REAL(DbKi) , DIMENSION(1:30) :: stiffXs = 0.0_R8Ki !< x array for stress-strain lookup table (up to nCoef) [-] REAL(DbKi) , DIMENSION(1:30) :: stiffYs = 0.0_R8Ki !< y array for stress-strain lookup table [-] @@ -298,6 +304,7 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Aq !< node added mass forcing (axial) [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: B !< node bottom contact force [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Bs !< node force due to bending moments [[N]] + REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Lf !< node crossflow viv lift force. TODO this will have to be allocated for only internal after fix [[N]] REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: Fnet !< total force on node [[N]] REAL(DbKi) , DIMENSION(:,:,:), ALLOCATABLE :: S !< node inverse mass matrix [[kg]] REAL(DbKi) , DIMENSION(:,:,:), ALLOCATABLE :: M !< node mass matrix [[kg]] @@ -305,6 +312,14 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(1:3) :: EndMomentB = 0.0_R8Ki !< vector of end moments due to bending at line end B [[N-m]] 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] + REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: yd_rms_old !< node old cf vel rms [m/s] + REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: ydd_rms_old !< node old cf accel rms [m/s^2] + REAL(DbKi) , DIMENSION(:,:), ALLOCATABLE :: rdd_old !< node accelerations previous iteration [m/s^2] END TYPE MD_Line ! ======================= ! ========= MD_ExtLd ======= @@ -725,6 +740,9 @@ subroutine MD_CopyLineProp(SrcLinePropData, DstLinePropData, CtrlCode, ErrStat, DstLinePropData%Cat = SrcLinePropData%Cat DstLinePropData%Cdn = SrcLinePropData%Cdn DstLinePropData%Cdt = SrcLinePropData%Cdt + DstLinePropData%Cl = SrcLinePropData%Cl + DstLinePropData%dF = SrcLinePropData%dF + DstLinePropData%cF = SrcLinePropData%cF DstLinePropData%ElasticMod = SrcLinePropData%ElasticMod DstLinePropData%nEApoints = SrcLinePropData%nEApoints DstLinePropData%stiffXs = SrcLinePropData%stiffXs @@ -766,6 +784,9 @@ subroutine MD_PackLineProp(RF, Indata) call RegPack(RF, InData%Cat) call RegPack(RF, InData%Cdn) call RegPack(RF, InData%Cdt) + call RegPack(RF, InData%Cl) + call RegPack(RF, InData%dF) + call RegPack(RF, InData%cF) call RegPack(RF, InData%ElasticMod) call RegPack(RF, InData%nEApoints) call RegPack(RF, InData%stiffXs) @@ -799,6 +820,9 @@ subroutine MD_UnPackLineProp(RF, OutData) call RegUnpack(RF, OutData%Cat); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Cdn); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Cdt); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Cl); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%dF); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%cF); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%ElasticMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nEApoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%stiffXs); if (RegCheckErr(RF, RoutineName)) return @@ -1700,6 +1724,9 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) DstLineData%Cat = SrcLineData%Cat DstLineData%Cdn = SrcLineData%Cdn DstLineData%Cdt = SrcLineData%Cdt + DstLineData%Cl = SrcLineData%Cl + DstLineData%dF = SrcLineData%dF + DstLineData%cF = SrcLineData%cF DstLineData%nEApoints = SrcLineData%nEApoints DstLineData%stiffXs = SrcLineData%stiffXs DstLineData%stiffYs = SrcLineData%stiffYs @@ -2010,6 +2037,18 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) end if DstLineData%Bs = SrcLineData%Bs end if + if (allocated(SrcLineData%Lf)) then + LB(1:2) = lbound(SrcLineData%Lf) + UB(1:2) = ubound(SrcLineData%Lf) + if (.not. allocated(DstLineData%Lf)) then + allocate(DstLineData%Lf(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%Lf.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%Lf = SrcLineData%Lf + end if if (allocated(SrcLineData%Fnet)) then LB(1:2) = lbound(SrcLineData%Fnet) UB(1:2) = ubound(SrcLineData%Fnet) @@ -2061,6 +2100,58 @@ subroutine MD_CopyLine(SrcLineData, DstLineData, CtrlCode, ErrStat, ErrMsg) end if DstLineData%LineWrOutput = SrcLineData%LineWrOutput end if + if (allocated(SrcLineData%phi)) then + LB(1:1) = lbound(SrcLineData%phi) + UB(1:1) = ubound(SrcLineData%phi) + if (.not. allocated(DstLineData%phi)) then + allocate(DstLineData%phi(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%phi.', ErrStat, ErrMsg, RoutineName) + return + end if + 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 + if (allocated(SrcLineData%yd_rms_old)) then + LB(1:1) = lbound(SrcLineData%yd_rms_old) + UB(1:1) = ubound(SrcLineData%yd_rms_old) + if (.not. allocated(DstLineData%yd_rms_old)) then + allocate(DstLineData%yd_rms_old(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%yd_rms_old.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%yd_rms_old = SrcLineData%yd_rms_old + end if + if (allocated(SrcLineData%ydd_rms_old)) then + LB(1:1) = lbound(SrcLineData%ydd_rms_old) + UB(1:1) = ubound(SrcLineData%ydd_rms_old) + if (.not. allocated(DstLineData%ydd_rms_old)) then + allocate(DstLineData%ydd_rms_old(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%ydd_rms_old.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%ydd_rms_old = SrcLineData%ydd_rms_old + end if + if (allocated(SrcLineData%rdd_old)) then + LB(1:2) = lbound(SrcLineData%rdd_old) + UB(1:2) = ubound(SrcLineData%rdd_old) + if (.not. allocated(DstLineData%rdd_old)) then + allocate(DstLineData%rdd_old(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstLineData%rdd_old.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstLineData%rdd_old = SrcLineData%rdd_old + end if end subroutine subroutine MD_DestroyLine(LineData, ErrStat, ErrMsg) @@ -2145,6 +2236,9 @@ subroutine MD_DestroyLine(LineData, ErrStat, ErrMsg) if (allocated(LineData%Bs)) then deallocate(LineData%Bs) end if + if (allocated(LineData%Lf)) then + deallocate(LineData%Lf) + end if if (allocated(LineData%Fnet)) then deallocate(LineData%Fnet) end if @@ -2157,6 +2251,18 @@ subroutine MD_DestroyLine(LineData, ErrStat, ErrMsg) if (allocated(LineData%LineWrOutput)) then deallocate(LineData%LineWrOutput) end if + if (allocated(LineData%phi)) then + deallocate(LineData%phi) + end if + if (allocated(LineData%yd_rms_old)) then + deallocate(LineData%yd_rms_old) + end if + if (allocated(LineData%ydd_rms_old)) then + deallocate(LineData%ydd_rms_old) + end if + if (allocated(LineData%rdd_old)) then + deallocate(LineData%rdd_old) + end if end subroutine subroutine MD_PackLine(RF, Indata) @@ -2188,6 +2294,9 @@ subroutine MD_PackLine(RF, Indata) call RegPack(RF, InData%Cat) call RegPack(RF, InData%Cdn) call RegPack(RF, InData%Cdt) + call RegPack(RF, InData%Cl) + call RegPack(RF, InData%dF) + call RegPack(RF, InData%cF) call RegPack(RF, InData%nEApoints) call RegPack(RF, InData%stiffXs) call RegPack(RF, InData%stiffYs) @@ -2223,6 +2332,7 @@ subroutine MD_PackLine(RF, Indata) call RegPackAlloc(RF, InData%Aq) call RegPackAlloc(RF, InData%B) call RegPackAlloc(RF, InData%Bs) + call RegPackAlloc(RF, InData%Lf) call RegPackAlloc(RF, InData%Fnet) call RegPackAlloc(RF, InData%S) call RegPackAlloc(RF, InData%M) @@ -2230,6 +2340,14 @@ subroutine MD_PackLine(RF, Indata) call RegPack(RF, InData%EndMomentB) 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) + call RegPackAlloc(RF, InData%yd_rms_old) + call RegPackAlloc(RF, InData%ydd_rms_old) + call RegPackAlloc(RF, InData%rdd_old) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -2265,6 +2383,9 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpack(RF, OutData%Cat); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Cdn); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Cdt); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Cl); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%dF); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%cF); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%nEApoints); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%stiffXs); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%stiffYs); if (RegCheckErr(RF, RoutineName)) return @@ -2300,6 +2421,7 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpackAlloc(RF, OutData%Aq); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%B); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Bs); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%Lf); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%Fnet); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%S); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%M); if (RegCheckErr(RF, RoutineName)) return @@ -2307,6 +2429,14 @@ subroutine MD_UnPackLine(RF, OutData) call RegUnpack(RF, OutData%EndMomentB); if (RegCheckErr(RF, RoutineName)) return 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 + call RegUnpackAlloc(RF, OutData%yd_rms_old); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%ydd_rms_old); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%rdd_old); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine MD_CopyExtLd(SrcExtLdData, DstExtLdData, CtrlCode, ErrStat, ErrMsg) From fd0b36b7a1491636719d9a7c3652fe7ead9fc684 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Mon, 24 Feb 2025 14:04:00 -0700 Subject: [PATCH 05/10] MD: VIV working, verified against MD-C --- modules/moordyn/src/MoorDyn.f90 | 5 +- modules/moordyn/src/MoorDyn_Line.f90 | 76 ++++++++++++++-------------- 2 files changed, 42 insertions(+), 39 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 4fc13f5b30..1593eff2db 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -36,7 +36,7 @@ MODULE MoorDyn TYPE(ProgDesc), PARAMETER :: MD_ProgDesc = ProgDesc( 'MoorDyn', 'v2.2.2', '2024-01-16' ) - INTEGER(IntKi), PARAMETER :: wordy = 1 ! verbosity level. >1 = more console output + INTEGER(IntKi), PARAMETER :: wordy = 0 ! verbosity level. >1 = more console output PUBLIC :: MD_Init PUBLIC :: MD_UpdateStates @@ -1490,8 +1490,9 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er Nx = Nx + m%LineList(l)%N ! if using viscoelastic model, need one more state per segment end if - if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 1) then + if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 0) then Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). TODO: when we change to only internal nodes need to make N - 1. + print *, "Added state for VIV. Nx is now", Nx endif m%LineStateIsN(l) = Nx diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 65f9bdf87c..62aad4518e 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -29,7 +29,7 @@ MODULE MoorDyn_Line PRIVATE - INTEGER(IntKi), PARAMETER :: wordy = 3 ! verbosity level. >1 = more console output + INTEGER(IntKi), PARAMETER :: wordy = 0 ! verbosity level. >1 = more console output PUBLIC :: SetupLine PUBLIC :: Line_Initialize @@ -178,7 +178,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%dl_1 = 0.0_DbKi end if - ! if using viscoelastic model, allocate additional state quantities + ! if using VIV model, allocate additional state quantities. TODO: interal nodes only if (Line%Cl > 0) then if (wordy > 1) print *, "Using the VIV model" ! allocate old acclerations [for VIV] @@ -189,8 +189,8 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) RETURN END IF ! initialize to unique values on the range 0-2pi - do I=0, Line%N - Line%phi(I) = (I/Line%N)*2*Pi + do I=0, Line%N ! TODO: internal ndoes only + Line%phi(I) = (REAL(I)/Line%N)*2*Pi enddo end if @@ -1042,12 +1042,13 @@ SUBROUTINE Line_SetState(Line, X, t) end if ! if using the viv mdodel, also set the lift force phase - if (Line%Cl > 0) then + 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 (t < 0.0002) print *, "Line State" , X do I=0, Line%N ! TODO: after checking change to internal - 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 + 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 else ! if only VIV state, then N+1 entries after internal node states - Line%phi(I) = X( 6*Line%N-6 + I + 1) - (2 * Pi * floor(X( 6*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 + Line%phi(I) = X( 6*Line%N-6 + I+1) - (2 * Pi * floor(X( 6*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 endif enddo endif @@ -1535,13 +1536,13 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! print *, "Made it to VIV model" ! Vortex Induced Vibration (VIV) cross-flow lift force - IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen then VIV to be calculated + IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen ! Ignore the following: and internal node then VIV to be calculated .AND. (I /= 0) .AND. (I /= N) ! ----- The Synchronization Model ------ ! Crossflow velocity and acceleration. rd component in the crossflow direction yd = dot_product(Line%rd(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) - ydd = dot_product(Line%rdd_old(0:2,I), cross_product(Line%q(:,I), normalize(Vp) ) ) + ydd = dot_product(Line%rdd_old(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) ! ! for checking rdd_old fix ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. (I > 95 .or. I <5)) then @@ -1602,37 +1603,17 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! The Lift force IF (I==0) THEN ! TODO: drop for only internal nodes - Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(1)*Line%l(1) + Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(1)*Line%l(1) ELSE IF (I==N) THEN ! TODO: drop for only internal nodes - Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(N)*Line%l(N) + Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(N)*Line%l(N) ELSE - Line%Lf(:,I) = 0.5 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * (Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) + Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * (Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) END IF if (wordy > 0) then if (Is_NaN(norm2(Line%Lf(:,I)))) print*, "Lf nan at node", I, "for line", Line%IdNum endif - ! finding nans - if (Line%time == 2*p%dtM0 .and. I==N-2) then - print *, "-------" - print *, "Line%time", Line%time - print *, "I", I - print *, "Line%Lf(:,I)", Line%Lf(:,I) - print *, "MagVp", MagVp - print *, "Line%U(:,I)", Line%U(:,I) - print *, "Line%phi(I)", Line%phi(I) - print *, "phi_dot", phi_dot - print *, "f_hat", f_hat - print *, "Line%phi_yd", Line%phi_yd - print *, "yd", yd - print *, "ydd", ydd - print *, "Line%rd(:,I)", Line%rd(:,I) - print *, "Line%rdd_old(0:2,I)", Line%rdd_old(0:2,I) - print *, "Line%q(:,I)", Line%q(:,I) - print *, "Vp", Vp - endif - ! print*, "end force" ! update state derivative with lift force frequency @@ -1643,6 +1624,27 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Xd( 6*Line%N-6 + I + 1) = phi_dot endif + ! ! debugging + ! if (Line%time < 2*p%dtM0 .AND. I == N/2) then + ! print *, "Line%time", Line%time + ! print *, "I", I + ! print *, "Line%Lf(:,I)", Line%Lf(:,I) + ! print *, "MagVp", MagVp + ! print *, "Line%U(:,I)", Line%U(:,I) + ! print *, "Line%phi(I)", Line%phi(I) + ! print *, "phi_dot", phi_dot + ! print *, "f_hat", f_hat + ! ! print *, "Line%phi_yd", Line%phi_yd + ! ! print *, "yd", yd + ! ! print *, "ydd", ydd + ! ! print *, "Line%rd(:,I)", Line%rd(:,I) + ! ! print *, "Line%rdd_old(0:2,I)", Line%rdd_old(:,I) + ! print *, "Line%q(:,I)", Line%q(:,I) + ! print *, "Vp", Vp + ! print*, "State Deriv", Xd + ! print *, "-------" + ! endif + ! print*, "end update state deriv" ! Miscd(I)[2] = As_dot; ! unused state that could be used for future amplitude calculations @@ -1771,11 +1773,11 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai IF (Line%Cl > 0) THEN ! ! for checking rdd_old fix - ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. I > 95) then + ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. I > N-5) then ! print*, "I in the rdd_old loop", I ! print*, "Sum1 for above I at J =", J, "is", Sum1 ! endif - Line%rdd_old(J-1,I-1) = Sum1 ! saving the acceleration for VIV RMS calculation. WARNING: I-1 is intential to match the incorrect approach in MD-C (map internal node accel to the first N-2 elements of rdd_old.) After testing REMOVE the I-1 and make the VIV model only apply to internal nodes. J-1 here to shift things correctly. Not sure why we need this but it doesnt work otherwise ... + Line%rdd_old(J,I-1) = Sum1 ! saving the acceleration for VIV RMS calculation. WARNING: I-1 is intential to match the incorrect approach in MD-C (map internal node accel to the first N-2 elements of rdd_old.) After testing REMOVE the I-1 and make the VIV model only apply to internal nodes. J-1 here to shift things correctly. Not sure why we need this but it doesnt work otherwise ... ENDIF END DO ! J @@ -1789,11 +1791,11 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen) then ! print*, "rdd_old at t = ", Line%time ! DO I = 0, 4 - ! print*, "I =", I, "rdd_old =", Line%rdd_old(0,I), Line%rdd_old(1,I), Line%rdd_old(2,I) + ! print*, "I =", I, "rdd_old =", Line%rdd_old(:,I) ! enddo ! print*, "..." ! DO I = N-4, N - ! print*, "I =", I, "rdd_old =", Line%rdd_old(0,I), Line%rdd_old(1,I), Line%rdd_old(2,I) + ! print*, "I =", I, "rdd_old =", Line%rdd_old(:,I) ! enddo ! endif From fdd47fc8bd2dd4223975ae11a5debfac1bb64856 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Tue, 25 Feb 2025 11:25:32 -0700 Subject: [PATCH 06/10] MD: VIV for internal nodes only, removing print statements --- modules/moordyn/src/MoorDyn.f90 | 2 +- modules/moordyn/src/MoorDyn_Line.f90 | 78 ++++++---------------------- 2 files changed, 16 insertions(+), 64 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 1593eff2db..c43b6e21bf 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -1491,7 +1491,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er end if if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 0) then - Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). TODO: when we change to only internal nodes need to make N - 1. + Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). print *, "Added state for VIV. Nx is now", Nx endif diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 62aad4518e..527e7687fb 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -178,7 +178,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) Line%dl_1 = 0.0_DbKi end if - ! if using VIV model, allocate additional state quantities. TODO: interal nodes only + ! if using VIV model, allocate additional state quantities. if (Line%Cl > 0) then if (wordy > 1) print *, "Using the VIV model" ! allocate old acclerations [for VIV] @@ -189,7 +189,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) RETURN END IF ! initialize to unique values on the range 0-2pi - do I=0, Line%N ! TODO: internal ndoes only + do I=0, Line%N Line%phi(I) = (REAL(I)/Line%N)*2*Pi enddo end if @@ -237,7 +237,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) ! allocate node force vectors ALLOCATE ( Line%W(3, 0:N), Line%Dp(3, 0:N), Line%Dq(3, 0:N), & Line%Ap(3, 0:N), Line%Aq(3, 0:N), Line%B(3, 0:N), Line%Bs(3, 0:N), & - Line%Lf(3, 0:N), Line%Fnet(3, 0:N), STAT = ErrStat ) ! TODO: LF only to internal nodes after testing + Line%Lf(3, 0:N), Line%Fnet(3, 0:N), STAT = ErrStat ) IF ( ErrStat /= ErrID_None ) THEN ErrMsg = ' Error allocating node force arrays.' !CALL CleanUp() @@ -393,7 +393,6 @@ SUBROUTINE Line_Initialize (Line, LineProp, p, ErrStat, ErrMsg) IF (wordy == 1) THEN CALL WrScr(' Message from catenary solver: '//ErrMsg2) ENDIF - ! print *, "Node positions: " DO J = 0,N ! Loop through all nodes per line where the line position and tension can be output Line%r(1,J) = Line%r(1,0) + (Line%r(1,N) - Line%r(1,0))*REAL(J, DbKi)/REAL(N, DbKi) @@ -1018,8 +1017,6 @@ SUBROUTINE Line_SetState(Line, X, t) INTEGER(IntKi) :: i ! index of segments or nodes along line INTEGER(IntKi) :: J ! index - - ! print*, "begining set state" ! store current time Line%time = t @@ -1043,8 +1040,7 @@ SUBROUTINE Line_SetState(Line, X, t) ! 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 (t < 0.0002) print *, "Line State" , X - do I=0, Line%N ! TODO: after checking change to internal + 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 else ! if only VIV state, then N+1 entries after internal node states @@ -1162,8 +1158,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Line%rd(J,0) = m%PointList(Line%AnchPoint)%rd(J) ! END DO - ! print*, "begining get state deriv" - ! -------------------- calculate various kinematic quantities --------------------------- DO I = 1, N @@ -1534,7 +1528,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%Dq(:,I) = 0.25*p%rhoW*Line%Cdt* Pi*d*(Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) * MagVq * Vq END IF - ! print *, "Made it to VIV model" ! Vortex Induced Vibration (VIV) cross-flow lift force IF ((Line%Cl > 0.0) .AND. (.NOT. Line%IC_gen)) THEN ! If non-zero lift coefficient and not during IC_gen ! Ignore the following: and internal node then VIV to be calculated .AND. (I /= 0) .AND. (I /= N) @@ -1542,22 +1535,12 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! Crossflow velocity and acceleration. rd component in the crossflow direction yd = dot_product(Line%rd(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) - ydd = dot_product(Line%rdd_old(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) + ydd = dot_product(Line%rdd_old(:,I), cross_product(Line%q(:,I), normalize(Vp) ) ) ! note: rdd_old initializes as 0's. End nodes don't ever get updated, thus stay at zero - ! ! for checking rdd_old fix - ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. (I > 95 .or. I <5)) then - ! print*, "rdd_old at t = ", Line%time - ! print*, "I =", I, "rdd_old(0:2,I) =", Line%rdd_old(0:2,I) - ! endif - - ! print*, "crossflow stuff done" - ! Rolling RMS calculation yd_rms = sqrt((((Line%n_m-1) * Line%yd_rms_old(I) * Line%yd_rms_old(I)) + (yd * yd)) / Line%n_m) ! RMS approximation from Thorsen ydd_rms = sqrt((((Line%n_m-1) * Line%ydd_rms_old(I) * Line%ydd_rms_old(I)) + (ydd * ydd)) / Line%n_m) - ! print*, "end RMS calc" - IF ((Line%time >= Line%t_old + p%dtM0) .OR. (Line%time == 0.0)) THEN ! Update the stormed RMS vaues ! update back indexing one moordyn time step (regardless of time integration scheme or coupling step size). T_old is handled at end of getStateDeriv when rdd_old is updated. Line%yd_rms_old(I) = yd_rms ! for rms back indexing (one moordyn timestep back) @@ -1574,8 +1557,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%phi_yd = 2*Pi + Line%phi_yd ! atan2 to 0-2Pi range ENDIF - ! print*, "end phase stuff" - ! Note: amplitude calculations and states commented out. Would be needed if a Cl vs A lookup table was ever implemented ! const real A_int = Misc(I)[1]; @@ -1602,10 +1583,12 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! const real C_l = cl_lookup(x = As/d); ! create function in Line.hpp that uses lookup table ! The Lift force - IF (I==0) THEN ! TODO: drop for only internal nodes - Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(1)*Line%l(1) - ELSE IF (I==N) THEN ! TODO: drop for only internal nodes - Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(N)*Line%l(N) + IF (I==0) THEN ! Disable for end nodes for now, node acceleration needed for synch model + Line%Lf(:,0) = 0.0_DbKi + ! Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(1)*Line%l(1) + ELSE IF (I==N) THEN ! Disable for end nodes for now, node acceleration needed for synch model + Line%Lf(:,I) = 0.0_DbKi + ! Line%Lf(:,N) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * Line%F(N)*Line%l(N) ELSE Line%Lf(:,I) = 0.25 * p%rhoW * d * MagVp * Line%Cl * cos(Line%phi(I)) * cross_product(Line%q(:,I), Vp) * (Line%F(I)*Line%l(I) + Line%F(I+1)*Line%l(I+1)) END IF @@ -1614,8 +1597,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai if (Is_NaN(norm2(Line%Lf(:,I)))) print*, "Lf nan at node", I, "for line", Line%IdNum endif - - ! print*, "end force" ! update state derivative with lift force frequency if (Line%ElasticMod > 1) then ! if both additional states are included then N+1 entries after internal node states and visco segment states @@ -1624,34 +1605,10 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Xd( 6*Line%N-6 + I + 1) = phi_dot endif - ! ! debugging - ! if (Line%time < 2*p%dtM0 .AND. I == N/2) then - ! print *, "Line%time", Line%time - ! print *, "I", I - ! print *, "Line%Lf(:,I)", Line%Lf(:,I) - ! print *, "MagVp", MagVp - ! print *, "Line%U(:,I)", Line%U(:,I) - ! print *, "Line%phi(I)", Line%phi(I) - ! print *, "phi_dot", phi_dot - ! print *, "f_hat", f_hat - ! ! print *, "Line%phi_yd", Line%phi_yd - ! ! print *, "yd", yd - ! ! print *, "ydd", ydd - ! ! print *, "Line%rd(:,I)", Line%rd(:,I) - ! ! print *, "Line%rdd_old(0:2,I)", Line%rdd_old(:,I) - ! print *, "Line%q(:,I)", Line%q(:,I) - ! print *, "Vp", Vp - ! print*, "State Deriv", Xd - ! print *, "-------" - ! endif - - ! print*, "end update state deriv" ! Miscd(I)[2] = As_dot; ! unused state that could be used for future amplitude calculations ENDIF - ! print*, "made it past VIV model" - ! ------ fluid acceleration components for current node (from MD-C) ------ DO J = 1, 3 aq(J) = DOT_PRODUCT( Line%Ud(:,I) , Line%q(:,I) ) * Line%q(J,I) ! tangential fluid acceleration component @@ -1746,9 +1703,9 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai ! total forces IF (I==0) THEN - Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ! TODO: remove LF here + Line%Fnet(:,I) = Line%T(:,1) + Line%Td(:,1) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ELSE IF (I==N) THEN - Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ! TODO: remove LF here + Line%Fnet(:,I) = -Line%T(:,N) - Line%Td(:,N) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) ELSE Line%Fnet(:,I) = Line%T(:,I+1) - Line%T(:,I) + Line%Td(:,I+1) - Line%Td(:,I) + Line%W(:,I) + Line%Dp(:,I) + Line%Dq(:,I) + Line%Ap(:,I) + Line%Aq(:,I) + Line%B(:,I) + Line%Bs(:,I) + Line%Lf(:,I) END IF @@ -1772,12 +1729,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Xd( 3*I-3 + J) = Sum1 ! dVdt = RHS * A (accelerations) IF (Line%Cl > 0) THEN - ! ! for checking rdd_old fix - ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen .and. I > N-5) then - ! print*, "I in the rdd_old loop", I - ! print*, "Sum1 for above I at J =", J, "is", Sum1 - ! endif - Line%rdd_old(J,I-1) = Sum1 ! saving the acceleration for VIV RMS calculation. WARNING: I-1 is intential to match the incorrect approach in MD-C (map internal node accel to the first N-2 elements of rdd_old.) After testing REMOVE the I-1 and make the VIV model only apply to internal nodes. J-1 here to shift things correctly. Not sure why we need this but it doesnt work otherwise ... + Line%rdd_old(J,I) = Sum1 ! saving the acceleration for VIV RMS calculation. End nodes are left at zero, VIV disabled for end nodes ENDIF END DO ! J @@ -1787,7 +1739,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai Line%t_old = Line%time ! for updating back indexing if statements endif - ! ! for checking rdd_old fix + ! ! for checking rdd_old ! if (Line%time <0.5+p%dtM0 .and. Line%time >0.5-p%dtM0 .and. .not. Line%IC_gen) then ! print*, "rdd_old at t = ", Line%time ! DO I = 0, 4 From 20bbb03fe55794f2086fc28faaf2b43bbdd88fb1 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Tue, 25 Feb 2025 12:39:48 -0700 Subject: [PATCH 07/10] MD: update OpenfAST_io for VIV --- openfast_io/openfast_io/FAST_reader.py | 19 +++++++++++++++++-- openfast_io/openfast_io/FAST_writer.py | 13 ++++++++++--- 2 files changed, 27 insertions(+), 5 deletions(-) diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index 8bd0d3438e..28df0aaf95 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -2909,8 +2909,11 @@ def read_MoorDyn(self, moordyn_file): self.fst_vt['MoorDyn']['Ca'] = [] self.fst_vt['MoorDyn']['CdAx'] = [] self.fst_vt['MoorDyn']['CaAx'] = [] + self.fst_vt['MoorDyn']['Cl'] = [] + self.fst_vt['MoorDyn']['dF'] = [] + self.fst_vt['MoorDyn']['cF'] = [] data_line = f.readline().strip().split() - while data_line[0] and data_line[0][:3] != '---': # OpenFAST searches for ---, so we'll do the same + while data_line[0] and data_line[0][:3] != '---': # OpenFAST searches for ---, so we'll do the same. TODO: make compatible with MD viscoelastic model self.fst_vt['MoorDyn']['Name'].append(str(data_line[0])) self.fst_vt['MoorDyn']['Diam'].append(float(data_line[1])) self.fst_vt['MoorDyn']['MassDen'].append(float(data_line[2])) @@ -2921,6 +2924,18 @@ def read_MoorDyn(self, moordyn_file): self.fst_vt['MoorDyn']['Ca'].append(float(data_line[7])) self.fst_vt['MoorDyn']['CdAx'].append(float(data_line[8])) self.fst_vt['MoorDyn']['CaAx'].append(float(data_line[9])) + if len(data_line) == 10: + self.fst_vt['MoorDyn']['Cl'].append(None) + self.fst_vt['MoorDyn']['dF'].append(None) + self.fst_vt['MoorDyn']['cF'].append(None) + elif (len(data_line) == 11): + self.fst_vt['MoorDyn']['Cl'].append(float(data_line[10])) + self.fst_vt['MoorDyn']['dF'].append(None) + self.fst_vt['MoorDyn']['cF'].append(None) + elif (len(data_line) == 13): + self.fst_vt['MoorDyn']['Cl'].append(float(data_line[10])) + self.fst_vt['MoorDyn']['dF'].append(float(data_line[11])) + self.fst_vt['MoorDyn']['cF'].append(float(data_line[12])) data_line = f.readline().strip().split() data_line = ''.join(data_line) # re-join @@ -3071,7 +3086,7 @@ def read_MoorDyn(self, moordyn_file): data_line = f.readline().strip().split() data_line = ''.join(data_line) # re-join for reading next section uniformly - elif 'control' in data_line.lower(): + elif 'control' in data_line.lower(): # TODO: add in failures section, add in user specified loads and damping section f.readline() f.readline() diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index 4046357cd5..805128e57e 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -2327,7 +2327,7 @@ def write_MAP(self): os.fsync(f) f.close() - def write_MoorDyn(self): + def write_MoorDyn(self): # TODO: see TODO's in FAST_reader/read_MoorDyn.py and make corresponding changes here self.fst_vt['Fst']['MooringFile'] = self.FAST_namingOut + '_MoorDyn.dat' moordyn_file = os.path.join(self.FAST_runDirectory, self.fst_vt['Fst']['MooringFile']) @@ -2337,8 +2337,8 @@ def write_MoorDyn(self): f.write('Generated with OpenFAST_IO\n') f.write('{!s:<22} {:<11} {:}'.format(self.fst_vt['MoorDyn']['Echo'], 'Echo', '- echo the input file data (flag)\n')) f.write('----------------------- LINE TYPES ------------------------------------------\n') - f.write(" ".join(['{:^11s}'.format(i) for i in ['Name', 'Diam', 'MassDen', 'EA', 'BA/-zeta', 'EI', 'Cd', 'Ca', 'CdAx', 'CaAx']])+'\n') - f.write(" ".join(['{:^11s}'.format(i) for i in ['(-)', '(m)', '(kg/m)', '(N)', '(N-s/-)', '(-)', '(-)', '(-)', '(-)', '(-)']])+'\n') + f.write(" ".join(['{:^11s}'.format(i) for i in ['Name', 'Diam', 'MassDen', 'EA', 'BA/-zeta', 'EI', 'Cd', 'Ca', 'CdAx', 'CaAx', 'Cl (optional)', 'dF (optional)', 'cF (optional)']])+'\n') + f.write(" ".join(['{:^11s}'.format(i) for i in ['(-)', '(m)', '(kg/m)', '(N)', '(N-s/-)', '(N-m^2)', '(-)', '(-)', '(-)', '(-)', '(-)', '(-)', '(-)']])+'\n') for i in range(len(self.fst_vt['MoorDyn']['Name'])): ln = [] ln.append('{:^11}'.format(self.fst_vt['MoorDyn']['Name'][i])) @@ -2351,7 +2351,14 @@ def write_MoorDyn(self): ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['Ca'][i])) ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['CdAx'][i])) ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['CaAx'][i])) + if self.fst_vt['MoorDyn']['Cl'][i] != None: + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['Cl'][i])) + if self.fst_vt['MoorDyn']['dF'][i] != None: + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['dF'][i])) + if self.fst_vt['MoorDyn']['cF'][i] != None: + ln.append('{:^11.4f}'.format(self.fst_vt['MoorDyn']['cF'][i])) f.write(" ".join(ln) + '\n') + if 'Rod_Name' in self.fst_vt['MoorDyn'] and self.fst_vt['MoorDyn']['Rod_Name']: f.write('----------------------- ROD TYPES ------------------------------------------\n') f.write(" ".join(['{:^11s}'.format(i) for i in ['TypeName', 'Diam', 'Mass/m', 'Cd', 'Ca', 'CdEnd', 'CaEnd']])+'\n') From 1c883b541da4f91c84cd4d072dd8cc8926852ad9 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Tue, 25 Feb 2025 12:51:18 -0700 Subject: [PATCH 08/10] MD: Code clean up --- modules/moordyn/src/MoorDyn.f90 | 3 +-- modules/moordyn/src/MoorDyn_Line.f90 | 2 -- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index c43b6e21bf..74b6ad7f72 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -1492,7 +1492,6 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 0) then Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). - print *, "Added state for VIV. Nx is now", Nx endif m%LineStateIsN(l) = Nx @@ -2928,7 +2927,7 @@ END FUNCTION Failed SUBROUTINE CheckError(ErrID,Msg) - ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLe + ! This subroutine sets the error message and level and cleans up if the error is >= AbortErrLev ! Passed arguments INTEGER(IntKi), INTENT(IN) :: ErrID ! The error identifier (ErrStat) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index 527e7687fb..cca56348bd 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -1712,8 +1712,6 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai END DO ! I - done looping through nodes - ! print *, "Done looping through nodes. Time to find acceleration" - ! loop through internal nodes and update their states <<< should/could convert to matrix operations instead of all these loops DO I=1, N-1 DO J=1,3 From 1f3e89a54ea86fde634c8763d1c1dd6ea5126bb0 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Thu, 27 Feb 2025 15:02:39 -0700 Subject: [PATCH 09/10] MD: Error check to prevent viv or visco with linearized MD --- modules/moordyn/src/MoorDyn.f90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/modules/moordyn/src/MoorDyn.f90 b/modules/moordyn/src/MoorDyn.f90 index 74b6ad7f72..f6c2eac93f 100644 --- a/modules/moordyn/src/MoorDyn.f90 +++ b/modules/moordyn/src/MoorDyn.f90 @@ -98,6 +98,8 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er REAL(DbKi) :: dtM ! actual mooring dynamics time step INTEGER(IntKi) :: NdtM ! number of time steps to integrate through with RK2 ! INTEGER(IntKi) :: ntWave ! number of time steps of wave data + LOGICAL :: compVIV = .FALSE. ! flag to check if simulating VIV + LOGICAL :: compVisco = .FALSE. ! flag to check if simulating viscoelastic line TYPE(MD_InputType) :: u_array(1) ! a size-one array for u to make call to TimeStep happy REAL(DbKi) :: t_array(1) ! a size-one array saying time is 0 to make call to TimeStep happy @@ -1488,10 +1490,12 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er if (m%LineTypeList(m%LineList(l)%PropsIdNum)%ElasticMod > 1) then ! todo add an error check here? or change to 2 or 3? Nx = Nx + m%LineList(l)%N ! if using viscoelastic model, need one more state per segment + compVisco = .TRUE. ! turn of flag for linearization error checking end if if (m%LineTypeList(m%LineList(l)%PropsIdNum)%Cl > 0) then Nx = Nx + m%LineList(l)%N+1 ! if using VIV model, need one more state per node (note here N is the num sgemnts, so N+1 is number of nodes). + compVIV = .TRUE. ! turn of flag for linearization error checking endif m%LineStateIsN(l) = Nx @@ -2875,7 +2879,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er z%dummy = 0 if (InitInp%Linearize) then - call MD_Init_Jacobian(InitInp, p, u, y, m, InitOut, ErrStat2, ErrMsg2); if(Failed()) return + IF ((compVIV) .OR. (compVisco)) THEN + ErrStat2 = ErrID_Fatal + ErrMsg2 = "Linearization cannot be used with the VIV or Viscoelastic model in MoorDyn" + CALL CheckError( ErrStat2, ErrMsg2 ) + RETURN + ELSE + call MD_Init_Jacobian(InitInp, p, u, y, m, InitOut, ErrStat2, ErrMsg2); if(Failed()) return + ENDIF endif CALL WrScr(' MoorDyn initialization completed.') From 02f07ec540a14041c1af08940294eb27cf31e3a1 Mon Sep 17 00:00:00 2001 From: RyanDavies19 Date: Thu, 27 Feb 2025 16:07:32 -0700 Subject: [PATCH 10/10] MD: initialize unused VIV things to zero --- modules/moordyn/src/MoorDyn_Line.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/modules/moordyn/src/MoorDyn_Line.f90 b/modules/moordyn/src/MoorDyn_Line.f90 index cca56348bd..a001d1597e 100644 --- a/modules/moordyn/src/MoorDyn_Line.f90 +++ b/modules/moordyn/src/MoorDyn_Line.f90 @@ -192,6 +192,11 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg) do I=0, Line%N Line%phi(I) = (REAL(I)/Line%N)*2*Pi enddo + + ! initialize other things to 0 + Line%rdd_old = 0.0_DbKi + Line%yd_rms_old = 0.0_DbKi + Line%yd_rms_old = 0.0_DbKi end if ! allocate node and segment tangent vectors @@ -1529,6 +1534,7 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, Fai END IF ! 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 ! Ignore the following: and internal node then VIV to be calculated .AND. (I /= 0) .AND. (I /= N) ! ----- The Synchronization Model ------