Skip to content

Add back free-surface normal computation based on finite differencing #3503

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
andrew-platt merged 1 commit into from
Nov 3, 2026
Merged

Add back free-surface normal computation based on finite differencing #3503

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 16 additions & 8 deletions modules/hydrodyn/src/Morison.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3503,7 +3503,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
REAL(ReKi) :: sinBeta, sinBeta1, sinBeta2
REAL(ReKi) :: cosBeta, cosBeta1, cosBeta2
REAL(ReKi) :: CMatrix(3,3), CMatrix1(3,3), CMatrix2(3,3), CTrans(3,3) ! Direction cosine matrix for element, and its transpose
REAL(ReKi) :: l, z1, z2, zMid, r1, r2, r1b, r2b, r1In, r2In, rMidIn, z_hi, zFillGroup
REAL(ReKi) :: l, z1, z2, zMid, r1, r2, r1b, r2b, r1In, r2In, rMidIn, rn, z_hi, zFillGroup
REAL(ReKi) :: Sa1, Sa2, Sa1b, Sa2b, SaMidb, Sa1In, Sa2In, SaMidIn
REAL(ReKi) :: Sb1, Sb2, Sb1b, Sb2b, SbMidb, Sb1In, Sb2In, SbMidIn
REAL(ReKi) :: dRdl_mg, dSadl_mg, dSbdl_mg ! shorthand for taper including marine growth of element i
Expand Down Expand Up @@ -3822,7 +3822,12 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
IF (p%HstMod > 0_IntKi) THEN
CALL GetTotalWaveElev(p, m, Time, posMid, ZetaMid, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL GetFreeSurfaceNormal( p, m, Time, posMid, n_hat, ErrStat2, ErrMsg2 )
IF (mem%MSecGeom == MSecGeom_Cyl) THEN
rn = 0.5 * (r1b + r2b)
ELSE
rn = 0.25 * MAX(Sa1b+Sa2b, Sb1b+Sb2b)
END IF
CALL GetFreeSurfaceNormal( p, m, Time, posMid, rn, n_hat, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
FSPt = (/posMid(1),posMid(2),ZetaMid/) ! Reference point on the free surface
ELSE
Expand Down Expand Up @@ -4623,13 +4628,15 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
Sa2 = mem%SaMGB(N+1)
Sb1 = mem%SbMGB( 1)
Sb2 = mem%SbMGB(N+1)
r1 = 0.5 * MAX(Sa1,Sb1) ! Only used for computing the free-surface normal with finite differencing
r2 = 0.5 * MAX(Sa2,Sb2) ! Only used for computing the free-surface normal with finite differencing
end if
if (mem%i_floor == 0) then ! both ends above or at seabed
! Compute loads on the end plate of node 1
IF (p%HstMod > 0_IntKi) THEN
CALL GetTotalWaveElev(p, m, Time, pos1, Zeta1, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL GetFreeSurfaceNormal(p, m, Time, pos1, n_hat, ErrStat2, ErrMsg2)
CALL GetFreeSurfaceNormal(p, m, Time, pos1, r1, n_hat, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
FSPt = (/pos1(1),pos1(2),Zeta1/) ! Reference point on the free surface
ELSE
Expand Down Expand Up @@ -4657,7 +4664,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
IF (p%HstMod > 0_IntKi) THEN
CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL GetFreeSurfaceNormal(p, m, Time, pos2, n_hat, ErrStat2, ErrMsg2)
CALL GetFreeSurfaceNormal(p, m, Time, pos2, r2, n_hat, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface
ELSE
Expand Down Expand Up @@ -4686,7 +4693,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat,
IF (p%HstMod > 0_IntKi) THEN
CALL GetTotalWaveElev(p, m, Time, pos2, Zeta2, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL GetFreeSurfaceNormal(p, m, Time, pos2, n_hat, ErrStat2, ErrMsg2)
CALL GetFreeSurfaceNormal(p, m, Time, pos2, r2, n_hat, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface
ELSE
Expand Down Expand Up @@ -4835,11 +4842,12 @@ SUBROUTINE GetTotalWaveElev(p, m, Time, pos, Zeta, ErrStat, ErrMsg )

END SUBROUTINE GetTotalWaveElev

SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, n, ErrStat, ErrMsg)
SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, r, n, ErrStat, ErrMsg)
TYPE(Morison_ParameterType), INTENT( IN ) :: p
TYPE(Morison_MiscVarType), INTENT( INOUT ) :: m
REAL(DbKi), INTENT( In ) :: Time
REAL(ReKi), INTENT( In ) :: pos(:) ! Position at which free-surface normal is to be calculated. Third entry ignored if present.
REAL(ReKi), INTENT( In ) :: r ! Characteristic section half-width used to filter the free-surface slope
REAL(ReKi), INTENT( OUT ) :: n(3) ! Free-surface normal vector
INTEGER(IntKi), INTENT( OUT ) :: ErrStat ! Error status of the operation
CHARACTER(*), INTENT( OUT ) :: ErrMsg ! Error message if errStat /= ErrID_None
Expand All @@ -4849,7 +4857,7 @@ SUBROUTINE GetFreeSurfaceNormal(p, m, Time, pos, n, ErrStat, ErrMsg)
ErrStat = ErrID_None
ErrMsg = ""

CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, n, ErrStat2, ErrMsg2 )
CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, n, ErrStat2, ErrMsg2, r )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

END SUBROUTINE GetFreeSurfaceNormal
Expand Down Expand Up @@ -5974,7 +5982,7 @@ SUBROUTINE getElementHstLds_Mod1(p, m, mem, Time, pos1, pos2, Zeta1, Zeta2, k_ha
rh = r1 + h0*dRdl
! Estimate the free-surface normal at the free-surface intersection, n_hat
IF (p%HstMod > 0_IntKi) THEN ! If wave stretching is enabled, compute free surface normal
CALL GetFreeSurfaceNormal(p, m, Time, FSInt, n_hat, ErrStat2, ErrMsg2 )
CALL GetFreeSurfaceNormal(p, m, Time, FSInt, rh, n_hat, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
ELSE ! Without wave stretching, use the normal of the SWL
n_hat = (/0.0_ReKi,0.0_ReKi,1.0_ReKi/)
Expand Down
32 changes: 27 additions & 5 deletions modules/seastate/src/SeaSt_WaveField.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrS
ErrStat = ErrID_None

IF (ALLOCATED(WaveField%WaveElev1) .or. ALLOCATED(WaveField%WaveElev2)) then
CALL WaveField_Interp_Setup3D(Time, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2)
CALL WaveField_Interp_Setup3D(Time+WaveField%WaveTimeShift, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2)
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
end if
Expand Down Expand Up @@ -99,23 +99,45 @@ subroutine WaveField_GetMinMaxWaveElevEstimate( WaveField, MinElev, MaxElev, Err
endif
end subroutine WaveField_GetMinMaxWaveElevEstimate

SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, ErrStat, ErrMsg )
SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, ErrStat, ErrMsg, r )
type(SeaSt_WaveFieldType), intent(in ) :: WaveField
type(GridInterp_MiscVarType), intent(inout) :: WaveField_m
real(DbKi), intent(in ) :: Time
real(ReKi), intent(in ) :: pos(:) ! Position at which free-surface normal is to be calculated. Third entry ignored if present.
real(ReKi), intent( out) :: n(3) ! Free-surface normal vector
integer(IntKi), intent( out) :: ErrStat ! Error status of the operation
character(*), intent( out) :: ErrMsg ! Error message if errStat /= ErrID_None
real(ReKi), optional, intent(in ) :: r ! Section half-width; if > 0, slope is central-differenced over +/- r (filters waves shorter than the section)

real(SiKi) :: slope(2)
real(SiKi) :: ZetaP, ZetaM
real(DbKi) :: TimeShifted
character(*), parameter :: RoutineName = 'WaveField_GetNodeWaveNormal'
integer(IntKi) :: errStat2
character(ErrMsgLen) :: errMsg2

ErrStat = ErrID_None
ErrMsg = ""

if (present(r)) then
if (r > 0.0_ReKi) then
ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)+r,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return;
ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)-r,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return;
slope(1) = (ZetaP-ZetaM) / real(2.0_ReKi*r,SiKi)
ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)+r/), ErrStat2, ErrMsg2 ); if (Failed()) return;
ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)-r/), ErrStat2, ErrMsg2 ); if (Failed()) return;
slope(2) = (ZetaP-ZetaM) / real(2.0_ReKi*r,SiKi)
n = Real( (/-slope(1),-slope(2),1.0_SiKi/), ReKi)
n = n / TwoNorm(n)
return
end if
end if

call GridInterpSetupN( (/Real(Time+WaveField%WaveTimeShift,ReKi),pos(1),pos(2)/), WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 )
n = (/0.0_ReKi,0.0_ReKi,1.0_ReKi/)
if (.not. ALLOCATED(WaveField%WaveElev1)) return

TimeShifted = Time + WaveField%WaveTimeShift
call GridInterpSetupN( (/Real(TimeShifted,ReKi),pos(1),pos(2)/), WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
slope = GridInterpS( WaveField%WaveElev1, WaveField%SrfGridParams, WaveField_m )
if (ALLOCATED(WaveField%WaveElev2)) then
slope = slope + GridInterpS( WaveField%WaveElev2, WaveField%SrfGridParams, WaveField_m )
Expand Down Expand Up @@ -460,7 +482,7 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod
ELSE ! Node is above SWL - need wave stretching

! Vertical wave stretching
CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m )

! Extrapolated wave stretching
Expand Down Expand Up @@ -585,7 +607,7 @@ SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, force
ELSE ! Node is above SWL - need wave stretching

! Vertical wave stretching
CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return;
FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m )
FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m )

Expand Down
Loading