From 4f0535e9c7ca1c29839a2ea1b2e9540e21ef2781 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 13 Jun 2023 12:59:10 -0600 Subject: [PATCH 1/4] Added a switch to the SeaState WaveField functions that forces the input point to be treated as in water. This is needed when querying the wave kinematics at the estimated intersection between the free surface and a Morison member to ensure nonzero wave kinematics are returned. --- modules/hydrodyn/src/Morison.f90 | 8 ++++---- modules/seastate/src/SeaSt_WaveField.f90 | 12 ++++++++---- modules/seastate/src/SeaState.f90 | 2 +- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index a88200f2ac..c1fee0283f 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -2596,7 +2596,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, END IF ! Get the wave elevation and wave kinematics at each node - CALL WaveField_GetWaveKin( p%WaveField, Time, pos1, m%nodeInWater(j), m%WaveElev1(j), m%WaveElev2(j), m%WaveElev(j), FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) + CALL WaveField_GetWaveKin( p%WaveField, Time, pos1, .FALSE., m%nodeInWater(j), m%WaveElev1(j), m%WaveElev2(j), m%WaveElev(j), FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) m%FDynP(j) = REAL(FDynP,ReKi) m%FV(:, j) = REAL(FV, ReKi) @@ -3048,8 +3048,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !----------------------------------------------------------------------------------------------------! ! Compute the distributed loads at the point of intersection between the member and the free surface ! !----------------------------------------------------------------------------------------------------! - ! Get wave dynamics at the free surface intersection - CALL WaveField_GetWaveKin( p%WaveField, Time, FSInt, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) + ! Get wave kinematics at the free-surface intersection. Set forceNodeInWater=.TRUE. to guarantee the free-surface intersection is in water. + CALL WaveField_GetWaveKin( p%WaveField, Time, FSInt, .TRUE., nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FDynPFSInt = REAL(FDynP,ReKi) FVFSInt = REAL(FV, ReKi) @@ -4229,7 +4229,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat END IF ! Get fluid velocity at the joint - CALL WaveField_GetWaveVel( p%WaveField, Time, pos, nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL WaveField_GetWaveVel( p%WaveField, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FV = REAL(FVTmp, ReKi) vrel = ( FV - u%Mesh%TranslationVel(:,J) ) * nodeInWater diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 4d6efea5ca..c7f8769772 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -136,10 +136,11 @@ SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg END SUBROUTINE WaveField_GetWaveNormal !-------------------- Subroutine for full wave field kinematics --------------------! -SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) + LOGICAL, INTENT( IN ) :: forceNodeInWater REAL(SiKi), INTENT( OUT ) :: WaveElev1 REAL(SiKi), INTENT( OUT ) :: WaveElev2 REAL(SiKi), INTENT( OUT ) :: WaveElev @@ -200,7 +201,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W ELSE ! Wave stretching enabled - IF ( pos(3) <= WaveElev ) THEN ! Node is submerged + IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged nodeInWater = 1_IntKi @@ -259,6 +260,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W ! Map the node z-position linearly from [-EffWtrDpth,m%WaveElev(j)] to [-EffWtrDpth,0] posPrime = pos posPrime(3) = WaveField%EffWtrDpth*(WaveField%EffWtrDpth+pos(3))/(WaveField%EffWtrDpth+WaveElev)-WaveField%EffWtrDpth + posPrime(3) = MIN( posPrime(3), 0.0_ReKi) ! Clamp z-position to zero. Needed when forceNodeInWater=.TRUE. ! Obtain the wave-field variables by interpolation with the mapped position. CALL SeaSt_Interp_Setup( Time, posPrime, WaveField%seast_interp_p, seast_interp_m, ErrStat2, ErrMsg2 ) @@ -290,10 +292,11 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, nodeInWater, WaveElev1, W END SUBROUTINE WaveField_GetWaveKin !-------------------- Subroutine for wave field velocity only --------------------! -SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) + LOGICAL, INTENT( IN ) :: forceNodeInWater INTEGER(IntKi), INTENT( OUT ) :: nodeInWater REAL(SiKi), INTENT( OUT ) :: FV(3) INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation @@ -333,7 +336,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ELSE ! Wave stretching enabled - IF ( pos(3) <= WaveElev ) THEN ! Node is submerged + IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged nodeInWater = 1_IntKi @@ -368,6 +371,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, nodeInWater, FV, ErrStat, ! Map the node z-position linearly from [-EffWtrDpth,m%WaveElev(j)] to [-EffWtrDpth,0] posPrime = pos posPrime(3) = WaveField%EffWtrDpth*(WaveField%EffWtrDpth+pos(3))/(WaveField%EffWtrDpth+WaveElev)-WaveField%EffWtrDpth + posPrime(3) = MIN( posPrime(3), 0.0_ReKi) ! Clamp z-position to zero. Needed when forceNodeInWater=.TRUE. ! Obtain the wave-field variables by interpolation with the mapped position. CALL SeaSt_Interp_Setup( Time, posPrime, WaveField%seast_interp_p, seast_interp_m, ErrStat2, ErrMsg2 ) diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90 index 91dd0e054c..6845379989 100644 --- a/modules/seastate/src/SeaState.f90 +++ b/modules/seastate/src/SeaState.f90 @@ -857,7 +857,7 @@ SUBROUTINE SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, Er DO i = 1, p%NWaveKin positionXYZ = (/p%WaveKinxi(i),p%WaveKinyi(i),p%WaveKinzi(i)/) - CALL WaveField_GetWaveKin( p%WaveField, Time, positionXYZ, nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) + CALL WaveField_GetWaveKin( p%WaveField, Time, positionXYZ, .FALSE., nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END DO From 594e48d23b904c8391fab97e33ab2f216d6ca62f Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 13 Jun 2023 16:05:43 -0600 Subject: [PATCH 2/4] Added a function to the SeaState WaveField module that computes the wave kinematics at a list of node positions. Update Morison_CalcOutput to use this new WaveField function. --- modules/hydrodyn/src/Morison.f90 | 75 ++++++------- modules/hydrodyn/src/Morison.txt | 4 +- modules/hydrodyn/src/Morison_Types.f90 | 132 +++++++++++++++++++++++ modules/seastate/src/SeaSt_WaveField.f90 | 118 +++++++++++++------- modules/seastate/src/SeaState.f90 | 6 +- 5 files changed, 257 insertions(+), 78 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index c1fee0283f..0c9e8451ba 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -2335,7 +2335,9 @@ SUBROUTINE AllocateNodeLoadVariables(InitInp, p, m, NNodes, errStat, errMsg ) ! Initialize errStat errStat = ErrID_None errMsg = "" - + + call AllocAry( m%DispNodePosHdn, 3, NNodes , 'm%DispNodePosHdn', errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) + call AllocAry( m%DispNodePosHst, 3, NNodes , 'm%DispNodePosHst', errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%nodeInWater , NNodes , 'm%nodeInWater' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%vrel , 3, NNodes , 'm%vrel' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) call AllocAry( m%FV , 3, NNodes , 'm%FV' , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, routineName) @@ -2366,6 +2368,8 @@ SUBROUTINE AllocateNodeLoadVariables(InitInp, p, m, NNodes, errStat, errMsg ) if (errStat >= AbortErrLev) return + m%DispNodePosHdn = 0.0_ReKi + m%DispNodePosHst = 0.0_ReKi m%nodeInWater = 0 m%vrel = 0.0_ReKi m%FV = 0.0_ReKi @@ -2571,42 +2575,21 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, errMsg = "" Imat = 0.0_ReKi g = p%Gravity - WtrDpth = p%WtrDpth + p%MSL2SWL ! Water depth measured from the free surface + WtrDpth = p%WtrDpth + p%MSL2SWL ! Water depth measured from the still water level + + !=============================================================================================== + ! Get displaced positions of the hydrodynamic nodes + CALL GetDisplacedNodePosition( .FALSE., m%DispNodePosHdn ) ! For hydrodynamic loads; depends on WaveDisp and WaveStMod + CALL GetDisplacedNodePosition( .TRUE. , m%DispNodePosHst ) ! For hydrostatic loads; always use actual displaced position - !InterpolationSlope = GetInterpolationSlope(Time, p, m, IntWrapIndx) - !=============================================================================================== ! Calculate the fluid kinematics at all mesh nodes and store for use in the equations below - + CALL WaveField_GetWaveKin( p%WaveField, Time, m%DispNodePosHdn, .FALSE., m%nodeInWater, m%WaveElev1, m%WaveElev2, m%WaveElev, m%FDynP, m%FV, m%FA, m%FAMCF, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + ! Compute fluid velocity relative to the structure DO j = 1, p%NNodes - IF (p%WaveDisp == 0 ) THEN - ! use the initial X,Y location - pos1(1) = u%Mesh%Position(1,j) - pos1(2) = u%Mesh%Position(2,j) - ELSE - ! Use current X,Y location - pos1(1) = u%Mesh%TranslationDisp(1,j) + u%Mesh%Position(1,j) - pos1(2) = u%Mesh%TranslationDisp(2,j) + u%Mesh%Position(2,j) - END IF - - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Wave stretching enabled - pos1(3) = u%Mesh%Position(3,j) + u%Mesh%TranslationDisp(3,j) - p%MSL2SWL ! Use the current Z location. - ELSE ! Wave stretching disabled - pos1(3) = u%Mesh%Position(3,j) - p%MSL2SWL ! We are intentionally using the undisplaced Z position of the node. - END IF - - ! Get the wave elevation and wave kinematics at each node - CALL WaveField_GetWaveKin( p%WaveField, Time, pos1, .FALSE., m%nodeInWater(j), m%WaveElev1(j), m%WaveElev2(j), m%WaveElev(j), FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - m%FDynP(j) = REAL(FDynP,ReKi) - m%FV(:, j) = REAL(FV, ReKi) - m%FA(:, j) = REAL(FA, ReKi) - IF (ALLOCATED(m%FAMCF)) THEN - m%FAMCF(:,j) = REAL(FAMCF,ReKi) - END IF m%vrel(:,j) = ( m%FV(:,j) - u%Mesh%TranslationVel(:,j) ) * m%nodeInWater(j) - - END DO ! j = 1, p%NNodes + END DO ! ============================================================================================== ! Calculate instantaneous loads on each member except for the hydrodynamic loads on member ends. @@ -3049,7 +3032,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ! Compute the distributed loads at the point of intersection between the member and the free surface ! !----------------------------------------------------------------------------------------------------! ! Get wave kinematics at the free-surface intersection. Set forceNodeInWater=.TRUE. to guarantee the free-surface intersection is in water. - CALL WaveField_GetWaveKin( p%WaveField, Time, FSInt, .TRUE., nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveKin( p%WaveField, Time, FSInt, .TRUE., nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FDynPFSInt = REAL(FDynP,ReKi) FVFSInt = REAL(FV, ReKi) @@ -3580,6 +3563,26 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, CONTAINS + + SUBROUTINE GetDisplacedNodePosition( forceDisplaced, pos ) + LOGICAL, INTENT( IN ) :: forceDisplaced ! Set to true to return the exact displaced position no matter WaveDisp or WaveStMod + REAL(ReKi), INTENT( OUT ) :: pos(:,:) ! Displaced node positions + + ! Undisplaced node position + pos = u%Mesh%Position + pos(3,:) = pos(3,:) - p%MSL2SWL ! Z position measured from the SWL + IF ( (p%WaveDisp /= 0) .OR. forceDisplaced ) THEN + ! Use displaced X and Y position + pos(1,:) = pos(1,:) + u%Mesh%TranslationDisp(1,:) + pos(2,:) = pos(2,:) + u%Mesh%TranslationDisp(2,:) + IF ( (p%WaveStMod > 0) .OR. forceDisplaced ) THEN + ! Use displaced Z position only when wave stretching is enabled + pos(3,:) = pos(3,:) + u%Mesh%TranslationDisp(3,:) + END IF + END IF + + END SUBROUTINE GetDisplacedNodePosition + SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. @@ -3592,7 +3595,7 @@ SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = "" - Zeta = WaveField_GetTotalWaveElev( p%WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta = WaveField_GetNodeTotalWaveElev( p%WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetTotalWaveElev @@ -3610,7 +3613,7 @@ SUBROUTINE GetFreeSurfaceNormal( Time, pos, r, n, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - CALL WaveField_GetWaveNormal( p%WaveField, Time, pos, r, n, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveNormal( p%WaveField, Time, pos, r, n, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetFreeSurfaceNormal @@ -4229,7 +4232,7 @@ SUBROUTINE Morison_UpdateDiscState( Time, u, p, x, xd, z, OtherState, m, errStat END IF ! Get fluid velocity at the joint - CALL WaveField_GetWaveVel( p%WaveField, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveVel( p%WaveField, Time, pos, .FALSE., nodeInWater, FVTmp, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FV = REAL(FVTmp, ReKi) vrel = ( FV - u%Mesh%TranslationVel(:,J) ) * nodeInWater diff --git a/modules/hydrodyn/src/Morison.txt b/modules/hydrodyn/src/Morison.txt index a12a4a0cde..41be828b76 100644 --- a/modules/hydrodyn/src/Morison.txt +++ b/modules/hydrodyn/src/Morison.txt @@ -306,7 +306,9 @@ typedef ^ OtherStateType IntKi # ..... Misc/Optimization variables................................................................................................. # Define any data that are used only for efficiency purposes (these variables are not associated with time): # e.g. indices for searching in an array, large arrays that are local variables in any routine called multiple times, etc. -typedef ^ MiscVarType ReKi FV {:}{:} - - "Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times" - +typedef ^ MiscVarType ReKi DispNodePosHdn {:}{:} - - "Instantaneous displaced position of the line element nodes at time t for hydrodynamic load calculation" (m) +typedef ^ ^ ReKi DispNodePosHst {:}{:} - - "Instantaneous displaced position of the line element nodes at time t for hydrostatic and other load calcuation" (m) +typedef ^ ^ ReKi FV {:}{:} - - "Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FA {:}{:} - - "Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FAMCF {:}{:} - - "Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times" - typedef ^ ^ ReKi FDynP {:} - - "Fluid dynamic pressure at line element node at time t, which may not correspond to the WaveTime array of times" - diff --git a/modules/hydrodyn/src/Morison_Types.f90 b/modules/hydrodyn/src/Morison_Types.f90 index ac212b7386..165080b7af 100644 --- a/modules/hydrodyn/src/Morison_Types.f90 +++ b/modules/hydrodyn/src/Morison_Types.f90 @@ -369,6 +369,8 @@ MODULE Morison_Types ! ======================= ! ========= Morison_MiscVarType ======= TYPE, PUBLIC :: Morison_MiscVarType + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: DispNodePosHdn !< Instantaneous displaced position of the line element nodes at time t for hydrodynamic load calculation [(m)] + REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: DispNodePosHst !< Instantaneous displaced position of the line element nodes at time t for hydrostatic and other load calcuation [(m)] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FV !< Fluid velocity at line element node at time t, which may not correspond to the WaveTime array of times [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FA !< Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: FAMCF !< Fluid acceleration at line element node at time t, which may not correspond to the WaveTime array of times [-] @@ -9017,6 +9019,34 @@ SUBROUTINE Morison_CopyMisc( SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg ! ErrStat = ErrID_None ErrMsg = "" +IF (ALLOCATED(SrcMiscData%DispNodePosHdn)) THEN + i1_l = LBOUND(SrcMiscData%DispNodePosHdn,1) + i1_u = UBOUND(SrcMiscData%DispNodePosHdn,1) + i2_l = LBOUND(SrcMiscData%DispNodePosHdn,2) + i2_u = UBOUND(SrcMiscData%DispNodePosHdn,2) + IF (.NOT. ALLOCATED(DstMiscData%DispNodePosHdn)) THEN + ALLOCATE(DstMiscData%DispNodePosHdn(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%DispNodePosHdn.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstMiscData%DispNodePosHdn = SrcMiscData%DispNodePosHdn +ENDIF +IF (ALLOCATED(SrcMiscData%DispNodePosHst)) THEN + i1_l = LBOUND(SrcMiscData%DispNodePosHst,1) + i1_u = UBOUND(SrcMiscData%DispNodePosHst,1) + i2_l = LBOUND(SrcMiscData%DispNodePosHst,2) + i2_u = UBOUND(SrcMiscData%DispNodePosHst,2) + IF (.NOT. ALLOCATED(DstMiscData%DispNodePosHst)) THEN + ALLOCATE(DstMiscData%DispNodePosHst(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%DispNodePosHst.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstMiscData%DispNodePosHst = SrcMiscData%DispNodePosHst +ENDIF IF (ALLOCATED(SrcMiscData%FV)) THEN i1_l = LBOUND(SrcMiscData%FV,1) i1_u = UBOUND(SrcMiscData%FV,1) @@ -9273,6 +9303,12 @@ SUBROUTINE Morison_DestroyMisc( MiscData, ErrStat, ErrMsg ) ErrStat = ErrID_None ErrMsg = "" +IF (ALLOCATED(MiscData%DispNodePosHdn)) THEN + DEALLOCATE(MiscData%DispNodePosHdn) +ENDIF +IF (ALLOCATED(MiscData%DispNodePosHst)) THEN + DEALLOCATE(MiscData%DispNodePosHst) +ENDIF IF (ALLOCATED(MiscData%FV)) THEN DEALLOCATE(MiscData%FV) ENDIF @@ -9368,6 +9404,16 @@ SUBROUTINE Morison_PackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Re_BufSz = 0 Db_BufSz = 0 Int_BufSz = 0 + Int_BufSz = Int_BufSz + 1 ! DispNodePosHdn allocated yes/no + IF ( ALLOCATED(InData%DispNodePosHdn) ) THEN + Int_BufSz = Int_BufSz + 2*2 ! DispNodePosHdn upper/lower bounds for each dimension + Re_BufSz = Re_BufSz + SIZE(InData%DispNodePosHdn) ! DispNodePosHdn + END IF + Int_BufSz = Int_BufSz + 1 ! DispNodePosHst allocated yes/no + IF ( ALLOCATED(InData%DispNodePosHst) ) THEN + Int_BufSz = Int_BufSz + 2*2 ! DispNodePosHst upper/lower bounds for each dimension + Re_BufSz = Re_BufSz + SIZE(InData%DispNodePosHst) ! DispNodePosHst + END IF Int_BufSz = Int_BufSz + 1 ! FV allocated yes/no IF ( ALLOCATED(InData%FV) ) THEN Int_BufSz = Int_BufSz + 2*2 ! FV upper/lower bounds for each dimension @@ -9505,6 +9551,46 @@ SUBROUTINE Morison_PackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Db_Xferred = 1 Int_Xferred = 1 + IF ( .NOT. ALLOCATED(InData%DispNodePosHdn) ) THEN + IntKiBuf( Int_Xferred ) = 0 + Int_Xferred = Int_Xferred + 1 + ELSE + IntKiBuf( Int_Xferred ) = 1 + Int_Xferred = Int_Xferred + 1 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHdn,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHdn,1) + Int_Xferred = Int_Xferred + 2 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHdn,2) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHdn,2) + Int_Xferred = Int_Xferred + 2 + + DO i2 = LBOUND(InData%DispNodePosHdn,2), UBOUND(InData%DispNodePosHdn,2) + DO i1 = LBOUND(InData%DispNodePosHdn,1), UBOUND(InData%DispNodePosHdn,1) + ReKiBuf(Re_Xferred) = InData%DispNodePosHdn(i1,i2) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF + IF ( .NOT. ALLOCATED(InData%DispNodePosHst) ) THEN + IntKiBuf( Int_Xferred ) = 0 + Int_Xferred = Int_Xferred + 1 + ELSE + IntKiBuf( Int_Xferred ) = 1 + Int_Xferred = Int_Xferred + 1 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHst,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHst,1) + Int_Xferred = Int_Xferred + 2 + IntKiBuf( Int_Xferred ) = LBOUND(InData%DispNodePosHst,2) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%DispNodePosHst,2) + Int_Xferred = Int_Xferred + 2 + + DO i2 = LBOUND(InData%DispNodePosHst,2), UBOUND(InData%DispNodePosHst,2) + DO i1 = LBOUND(InData%DispNodePosHst,1), UBOUND(InData%DispNodePosHst,1) + ReKiBuf(Re_Xferred) = InData%DispNodePosHst(i1,i2) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF IF ( .NOT. ALLOCATED(InData%FV) ) THEN IntKiBuf( Int_Xferred ) = 0 Int_Xferred = Int_Xferred + 1 @@ -9883,6 +9969,52 @@ SUBROUTINE Morison_UnPackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err Re_Xferred = 1 Db_Xferred = 1 Int_Xferred = 1 + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! DispNodePosHdn not allocated + Int_Xferred = Int_Xferred + 1 + ELSE + Int_Xferred = Int_Xferred + 1 + i1_l = IntKiBuf( Int_Xferred ) + i1_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + i2_l = IntKiBuf( Int_Xferred ) + i2_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%DispNodePosHdn)) DEALLOCATE(OutData%DispNodePosHdn) + ALLOCATE(OutData%DispNodePosHdn(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%DispNodePosHdn.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i2 = LBOUND(OutData%DispNodePosHdn,2), UBOUND(OutData%DispNodePosHdn,2) + DO i1 = LBOUND(OutData%DispNodePosHdn,1), UBOUND(OutData%DispNodePosHdn,1) + OutData%DispNodePosHdn(i1,i2) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! DispNodePosHst not allocated + Int_Xferred = Int_Xferred + 1 + ELSE + Int_Xferred = Int_Xferred + 1 + i1_l = IntKiBuf( Int_Xferred ) + i1_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + i2_l = IntKiBuf( Int_Xferred ) + i2_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%DispNodePosHst)) DEALLOCATE(OutData%DispNodePosHst) + ALLOCATE(OutData%DispNodePosHst(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%DispNodePosHst.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i2 = LBOUND(OutData%DispNodePosHst,2), UBOUND(OutData%DispNodePosHst,2) + DO i1 = LBOUND(OutData%DispNodePosHst,1), UBOUND(OutData%DispNodePosHst,1) + OutData%DispNodePosHst(i1,i2) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO + END DO + END IF IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! FV not allocated Int_Xferred = Int_Xferred + 1 ELSE diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index c7f8769772..a01caaa831 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -8,27 +8,29 @@ MODULE SeaSt_WaveField PRIVATE ! Public functions and subroutines -PUBLIC WaveField_GetWaveElev1 -PUBLIC WaveField_GetWaveElev2 -PUBLIC WaveField_GetTotalWaveElev -PUBLIC WaveField_GetWaveNormal +PUBLIC WaveField_GetNodeWaveElev1 +PUBLIC WaveField_GetNodeWaveElev2 +PUBLIC WaveField_GetNodeTotalWaveElev +PUBLIC WaveField_GetNodeWaveNormal +PUBLIC WaveField_GetNodeWaveKin +PUBLIC WaveField_GetNodeWaveVel + PUBLIC WaveField_GetWaveKin -PUBLIC WaveField_GetWaveVel CONTAINS !-------------------- Subroutine for wave elevation ------------------! -FUNCTION WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetWaveElev1 + REAL(SiKi) :: WaveField_GetNodeWaveElev1 REAL(SiKi) :: Zeta LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveElev1' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveElev1' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -42,21 +44,21 @@ FUNCTION WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat, ErrMsg ) Zeta = 0.0_SiKi END IF - WaveField_GetWaveElev1 = Zeta + WaveField_GetNodeWaveElev1 = Zeta -END FUNCTION WaveField_GetWaveElev1 +END FUNCTION WaveField_GetNodeWaveElev1 -FUNCTION WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetWaveElev2 + REAL(SiKi) :: WaveField_GetNodeWaveElev2 REAL(SiKi) :: Zeta LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveElev2' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveElev2' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -70,36 +72,36 @@ FUNCTION WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat, ErrMsg ) Zeta = 0.0_SiKi END IF - WaveField_GetWaveElev2 = Zeta + WaveField_GetNodeWaveElev2 = Zeta -END FUNCTION WaveField_GetWaveElev2 +END FUNCTION WaveField_GetNodeWaveElev2 -FUNCTION WaveField_GetTotalWaveElev( WaveField, Time, pos, ErrStat, ErrMsg ) +FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, Time, pos, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(*) ! Position at which free-surface elevation is to be calculated. Third entry ignored if present. INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None - REAL(SiKi) :: WaveField_GetTotalWaveElev + REAL(SiKi) :: WaveField_GetNodeTotalWaveElev REAL(SiKi) :: Zeta1, Zeta2 LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetTotalWaveElev' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeTotalWaveElev' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None ErrMsg = "" - Zeta1 = WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta1 = WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - Zeta2 = WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + Zeta2 = WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveField_GetTotalWaveElev = Zeta1 + Zeta2 + WaveField_GetNodeTotalWaveElev = Zeta1 + Zeta2 -END FUNCTION WaveField_GetTotalWaveElev +END FUNCTION WaveField_GetNodeTotalWaveElev -SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time @@ -110,7 +112,7 @@ SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None REAL(SiKi) :: ZetaP,ZetaM REAL(ReKi) :: r1,dZetadx,dZetady - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetFreeSurfaceNormal' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveNormal' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 ErrStat = ErrID_None @@ -118,25 +120,25 @@ SUBROUTINE WaveField_GetWaveNormal( WaveField, Time, pos, r, n, ErrStat, ErrMsg r1 = MAX(r,1.0e-6) ! In case r is zero - ZetaP = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ) + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ZetaM = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ) + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) dZetadx = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) - ZetaP = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1),pos(2)+r1/), ErrStat2, ErrMsg2 ) + ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1),pos(2)+r1/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ZetaM = WaveField_GetTotalWaveElev( WaveField, Time, (/pos(1),pos(2)-r1/), ErrStat2, ErrMsg2 ) + ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, Time, (/pos(1),pos(2)-r1/), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) dZetady = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) n = (/-dZetadx,-dZetady,1.0_ReKi/) n = n / SQRT(Dot_Product(n,n)) -END SUBROUTINE WaveField_GetWaveNormal +END SUBROUTINE WaveField_GetNodeWaveNormal !-------------------- Subroutine for full wave field kinematics --------------------! -SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) @@ -156,7 +158,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInW REAL(ReKi) :: posXY(2), posPrime(3), posXY0(3) TYPE(SeaSt_Interp_MiscVarType) :: SeaSt_Interp_m LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveKin' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveKin' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -168,9 +170,9 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInW FAMCF(:) = 0.0 ! Wave elevation - WaveElev1 = WaveField_GetWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev1 = WaveField_GetNodeWaveElev1( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveElev2 = WaveField_GetWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev2 = WaveField_GetNodeWaveElev2( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) WaveElev = WaveElev1 + WaveElev2 @@ -289,10 +291,10 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInW END IF ! If wave stretching is on or off -END SUBROUTINE WaveField_GetWaveKin +END SUBROUTINE WaveField_GetNodeWaveKin !-------------------- Subroutine for wave field velocity only --------------------! -SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField REAL(DbKi), INTENT( IN ) :: Time REAL(ReKi), INTENT( IN ) :: pos(3) @@ -306,7 +308,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInW REAL(ReKi) :: posXY(2), posPrime(3), posXY0(3) TYPE(SeaSt_Interp_MiscVarType) :: SeaSt_Interp_m LOGICAL :: FirstWarn_Clamp - CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveVel' + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetNodeWaveVel' INTEGER(IntKi) :: errStat2 CHARACTER(ErrMsgLen) :: errMsg2 @@ -317,7 +319,7 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInW posXY0 = (/pos(1),pos(2),0.0_ReKi/) ! Wave elevation - WaveElev = WaveField_GetTotalWaveElev( WaveField, Time, pos, ErrStat2, ErrMsg2 ) + WaveElev = WaveField_GetNodeTotalWaveElev( WaveField, Time, pos, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) IF (WaveField%WaveStMod == 0) THEN ! No wave stretching @@ -390,6 +392,46 @@ SUBROUTINE WaveField_GetWaveVel( WaveField, Time, pos, forceNodeInWater, nodeInW END IF ! If wave stretching is on or off -END SUBROUTINE WaveField_GetWaveVel +END SUBROUTINE WaveField_GetNodeWaveVel + +SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) + TYPE(SeaSt_WaveFieldType), INTENT( IN ) :: WaveField + REAL(DbKi), INTENT( IN ) :: Time + REAL(ReKi), INTENT( IN ) :: pos(:,:) + LOGICAL, INTENT( IN ) :: forceNodeInWater + REAL(SiKi), INTENT( OUT ) :: WaveElev1(:) + REAL(SiKi), INTENT( OUT ) :: WaveElev2(:) + REAL(SiKi), INTENT( OUT ) :: WaveElev(:) + REAL(ReKi), INTENT( OUT ) :: FV(:,:) + REAL(ReKi), INTENT( OUT ) :: FA(:,:) + REAL(ReKi), ALLOCATABLE, INTENT( OUT ) :: FAMCF(:,:) + REAL(ReKi), INTENT( OUT ) :: FDynP(:) + INTEGER(IntKi), INTENT( OUT ) :: nodeInWater(:) + INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation + CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None + + CHARACTER(*), PARAMETER :: RoutineName = 'WaveField_GetWaveKin' + INTEGER(IntKi) :: errStat2 + CHARACTER(ErrMsgLen) :: errMsg2 + + INTEGER(IntKi) :: NumPoints, i + REAL(SiKi) :: FDynP_node, FV_node(3), FA_node(3), FAMCF_node(3) + + ErrStat = ErrID_None + ErrMsg = "" + + NumPoints = size(pos, dim=2) + DO i = 1, NumPoints + CALL WaveField_GetNodeWaveKin( WaveField, Time, pos(:,i), forceNodeInWater, nodeInWater(i), WaveElev1(i), WaveElev2(i), WaveElev(i), FDynP_node, FV_node, FA_node, FAMCF_node, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + FDynP(i) = REAL(FDynP_node,ReKi) + FV(:, i) = REAL(FV_node, ReKi) + FA(:, i) = REAL(FA_node, ReKi) + IF (ALLOCATED(FAMCF)) THEN + FAMCF(:,i) = REAL(FAMCF_node,ReKi) + END IF + END DO + +END SUBROUTINE WaveField_GetWaveKin END MODULE SeaSt_WaveField diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90 index 6845379989..0855e36a45 100644 --- a/modules/seastate/src/SeaState.f90 +++ b/modules/seastate/src/SeaState.f90 @@ -857,7 +857,7 @@ SUBROUTINE SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, Er DO i = 1, p%NWaveKin positionXYZ = (/p%WaveKinxi(i),p%WaveKinyi(i),p%WaveKinzi(i)/) - CALL WaveField_GetWaveKin( p%WaveField, Time, positionXYZ, .FALSE., nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveKin( p%WaveField, Time, positionXYZ, .FALSE., nodeInWater, zeta1, zeta2, zeta, WaveDynP(i), WaveVel(:,i), WaveAcc(:,i), WaveAccMCF(:,i), ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END DO @@ -865,9 +865,9 @@ SUBROUTINE SeaSt_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, ErrStat, Er DO i = 1, p%NWaveElev positionXY = (/p%WaveElevxi(i),p%WaveElevyi(i)/) - WaveElev1(i) = WaveField_GetWaveElev1( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) + WaveElev1(i) = WaveField_GetNodeWaveElev1( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - WaveElev2(i) = WaveField_GetWaveElev2( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) + WaveElev2(i) = WaveField_GetNodeWaveElev2( p%WaveField, Time, positionXY, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) WaveElev(i) = WaveElev1(i) + WaveElev2(i) END DO From b7b5c681fe9d7b9b58c164c02f4b866a5d5a1082 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 13 Jun 2023 16:46:24 -0600 Subject: [PATCH 3/4] Fix segfault with MCF --- modules/seastate/src/SeaSt_WaveField.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index a01caaa831..c33c053aaf 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -404,7 +404,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInW REAL(SiKi), INTENT( OUT ) :: WaveElev(:) REAL(ReKi), INTENT( OUT ) :: FV(:,:) REAL(ReKi), INTENT( OUT ) :: FA(:,:) - REAL(ReKi), ALLOCATABLE, INTENT( OUT ) :: FAMCF(:,:) + REAL(ReKi), INTENT( OUT ) :: FAMCF(:,:) REAL(ReKi), INTENT( OUT ) :: FDynP(:) INTEGER(IntKi), INTENT( OUT ) :: nodeInWater(:) INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation @@ -427,7 +427,7 @@ SUBROUTINE WaveField_GetWaveKin( WaveField, Time, pos, forceNodeInWater, nodeInW FDynP(i) = REAL(FDynP_node,ReKi) FV(:, i) = REAL(FV_node, ReKi) FA(:, i) = REAL(FA_node, ReKi) - IF (ALLOCATED(FAMCF)) THEN + IF (ALLOCATED(WaveField%WaveAccMCF)) THEN FAMCF(:,i) = REAL(FAMCF_node,ReKi) END IF END DO From b6a4fe8d6089aacd3b4eb52dbc3ae3845ca992b1 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Wed, 14 Jun 2023 13:35:25 -0600 Subject: [PATCH 4/4] Streamline Morison_CalcOutput by using precomputed displaced node positions --- modules/hydrodyn/src/Morison.f90 | 87 ++++++++------------------------ 1 file changed, 22 insertions(+), 65 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index 0c9e8451ba..1192036db1 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -2651,12 +2651,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed ! calculate instantaneous incline angle and heading, and related trig values - ! the first and last NodeIndx values point to the corresponding Joint nodes idices which are at the start of the Mesh - - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh + pos1 = m%DispNodePosHst(:, mem%NodeIndx(i )) + pos2 = m%DispNodePosHst(:, mem%NodeIndx(i+1)) call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -2785,12 +2782,9 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = max(mem%i_floor,1), N ! loop through member elements that are not completely buried in the seabed ! calculate instantaneous incline angle and heading, and related trig values - ! the first and last NodeIndx values point to the corresponding Joint nodes idices which are at the start of the Mesh - - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! the first and last NodeIndx values point to the corresponding Joint nodes indices which are at the start of the Mesh + pos1 = m%DispNodePosHst(:,mem%NodeIndx(i )) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(i+1)) call GetOrientationAngles( pos1, pos2, phi, sinPhi, cosPhi, tanPhi, sinBeta, cosBeta, k_hat, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -2912,16 +2906,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, DO i = mem%i_floor+1,N ! loop through member nodes starting from the first node above seabed, but skip the last node which should not be submerged anyways ! Get positions of node i and i+1 - IF (p%WaveDisp /= 0) THEN ! Use current position - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i)) + u%Mesh%Position(:, mem%NodeIndx(i)) - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(i+1)) + u%Mesh%Position(:, mem%NodeIndx(i+1)) - ELSE ! Use initial position - pos1 = u%Mesh%Position(:, mem%NodeIndx(i)) - pos2 = u%Mesh%Position(:, mem%NodeIndx(i+1)) - END if - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1(3) = pos1(3) - p%MSL2SWL - pos2(3) = pos2(3) - p%MSL2SWL + pos1 = m%DispNodePosHdn(:,mem%NodeIndx(i )) + pos2 = m%DispNodePosHdn(:,mem%NodeIndx(i+1)) ! Free surface elevation above or below node i and i+1 Zeta1 = m%WaveElev(mem%NodeIndx(i)) @@ -3176,15 +3162,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ELSE IF ( MemSubStat .NE. 3_IntKi) THEN ! Skip members with centerline completely out of water !----------------------------No load smoothing----------------------------! DO i = mem%i_floor+1,N+1 ! loop through member nodes starting from the first node above seabed - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - ! Using the initial z-position to be consistent with the evaluation of wave kinematics - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN - ! Use current z-position - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i)) - p%MSL2SWL - ELSE - ! Use initial z-position - z1 = u%Mesh%Position(3, mem%NodeIndx(i)) - p%MSL2SWL - END IF + z1 = m%DispNodePosHdn(3, mem%NodeIndx(i)) !---------------------------------------------Compute deltal and h_c------------------------------------------! ! Cannot make any assumption about WaveStMod and member orientation IF ( m%NodeInWater(mem%NodeIndx(i)) .EQ. 0_IntKi ) THEN ! Node is out of water @@ -3202,11 +3180,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( m%NodeInWater(mem%NodeIndx(i-1)) .EQ. 1_IntKi ) THEN ! Node to the left is submerged deltalLeft = 0.5_ReKi * mem%dl ELSE ! Element i-1 crosses the free surface - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Use current z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i-1)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i-1)) - p%MSL2SWL - ELSE ! Use initial z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i-1)) - p%MSL2SWL - END IF + z2 = m%DispNodePosHdn(3, mem%NodeIndx(i-1)) IF ( p%WaveStMod > 0_IntKi ) THEN ! Wave stretching enabled zeta1 = m%WaveElev(mem%NodeIndx(i )) zeta2 = m%WaveElev(mem%NodeIndx(i-1)) @@ -3225,11 +3199,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF ( m%NodeInWater(mem%NodeIndx(i+1)) .EQ. 1_IntKi ) THEN ! Node to the right is submerged deltalRight = 0.5_ReKi * mem%dl ELSE ! Element i crosses the free surface - IF (p%WaveStMod > 0 .AND. p%WaveDisp /= 0) THEN ! Use current z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) + u%Mesh%TranslationDisp(3, mem%NodeIndx(i+1)) - p%MSL2SWL - ELSE ! Use initial z-position - z2 = u%Mesh%Position(3, mem%NodeIndx(i+1)) - p%MSL2SWL - END IF + z2 = m%DispNodePosHdn(3, mem%NodeIndx(i+1)) IF ( p%WaveStMod > 0_IntKi ) THEN ! Wave stretching enabled zeta1 = m%WaveElev(mem%NodeIndx(i )) zeta2 = m%WaveElev(mem%NodeIndx(i+1)) @@ -3331,12 +3301,8 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !-----------------------------------------------------------------------------------------------------! ! reassign convenience variables to correspond to member ends ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(1)) + u%Mesh%Position(:, mem%NodeIndx(1)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(2)) + u%Mesh%Position(:, mem%NodeIndx(2)) - pos2(3) = pos2(3) - p%MSL2SWL - z1 = pos1(3) - + pos1 = m%DispNodePosHst(:,mem%NodeIndx(1)) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(2)) call GetOrientationAngles( pos1, pos2, phi1, sinPhi1, cosPhi1, tanPhi, sinBeta1, cosBeta1, k_hat1, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) if ( N == 1 ) then ! Only one element in member @@ -3347,18 +3313,14 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, k_hat2 = k_hat1 else ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N)) + u%Mesh%Position(:, mem%NodeIndx(N)) - pos1(3) = pos1(3) - p%MSL2SWL - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL + pos1 = m%DispNodePosHst(:, mem%NodeIndx(N )) + pos2 = m%DispNodePosHst(:, mem%NodeIndx(N+1)) call GetOrientationAngles( pos1, pos2, phi2, sinPhi2, cosPhi2, tanPhi, sinBeta2, cosBeta2, k_hat2, errStat2, errMsg2 ) call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) end if - - ! We need to subtract the MSL2SWL offset to place this in the SWL reference system - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL - z2 = pos2(3) + ! z-coordinates of the two ends of the member + z1 = m%DispNodePosHst(3,mem%NodeIndx( 1)) + z2 = m%DispNodePosHst(3,mem%NodeIndx(N+1)) !----------------------------------- filled buoyancy loads: starts -----------------------------------! !TODO: Do the equations below still work if z1 > z2 ? @@ -3403,13 +3365,10 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !---------------------------------- external buoyancy loads: starts ----------------------------------! if ( (.not. mem%PropPot) .AND. (mem%MHstLMod /= 0) ) then - - ! Get positions of member end nodes - pos1 = u%Mesh%TranslationDisp(:, mem%NodeIndx(1 )) + u%Mesh%Position(:, mem%NodeIndx(1 )) - pos1(3) = pos1(3) - p%MSL2SWL - r1 = mem%RMGB(1 ) - pos2 = u%Mesh%TranslationDisp(:, mem%NodeIndx(N+1)) + u%Mesh%Position(:, mem%NodeIndx(N+1)) - pos2(3) = pos2(3) - p%MSL2SWL + ! Get positions and scaled radii of member end nodes + pos1 = m%DispNodePosHst(:,mem%NodeIndx( 1)) + pos2 = m%DispNodePosHst(:,mem%NodeIndx(N+1)) + r1 = mem%RMGB( 1) r2 = mem%RMGB(N+1) if (mem%i_floor == 0) then ! both ends above or at seabed ! Compute loads on the end plate of node 1 @@ -3559,8 +3518,6 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, !---------------------------------------------------------------------------------------------------------------! ! Map calculated results into the y%WriteOutput Array CALL MrsnOut_MapOutputs(y, p, u, m) - - CONTAINS