Skip to content

Minor numerical improvements #2975

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 11 commits into from
Sep 2, 2025
Merged

Minor numerical improvements #2975

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
10182e4
MD driver: fix cleanup/end routine
bjonkman Aug 12, 2025
9d4da52
Bug fix: NWTC Library: don't use optional argument if not present
bjonkman Aug 13, 2025
e5549b1
MD: remove unnecessary array
bjonkman Aug 18, 2025
930d29e
ModMesh: fix LogMap/Euler angle usage
bjonkman Aug 18, 2025
edab463
NWTC Lib: numeric fixes
bjonkman Aug 19, 2025
cc15cc9
Merge branch 'dev' into b/MD_updates
bjonkman Aug 19, 2025
9bd2512
HD/Morison: Remove double operator (extension to standard)
bjonkman Aug 22, 2025
4126a64
HD: adjust spacing in registry txt file
bjonkman Aug 22, 2025
00b8679
Update modules/nwtc-library/src/NWTC_Num.f90
bjonkman Sep 2, 2025
8a5b02f
Update modules/nwtc-library/src/NWTC_Num.f90
bjonkman Sep 2, 2025
eab1548
Update modules/nwtc-library/src/NWTC_Num.f90
bjonkman Sep 2, 2025
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
7 changes: 1 addition & 6 deletions modules/aerodyn/src/AirfoilInfo.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1745,12 +1745,7 @@ subroutine AFI_ComputeAirfoilCoefs1D( AOA, p, AFI_interp, errStat, errMsg, Table


! Spline interpolation of lower table based on requested AOA

IntAFCoefs(1:s1) = CubicSplineInterpM( Alpha &
, p%Table(iTab)%Alpha &
, p%Table(iTab)%Coefs &
, p%Table(iTab)%SplineCoefs &
, ErrStat, ErrMsg )
CALL CubicSplineInterpM( Alpha, p%Table(iTab)%Alpha, p%Table(iTab)%Coefs, p%Table(iTab)%SplineCoefs, IntAFCoefs(1:s1) )
end if

AFI_interp%Cl = IntAFCoefs(p%ColCl)
Expand Down
2 changes: 1 addition & 1 deletion modules/hydrodyn/src/HydroDyn.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ typedef HydroDyn/HydroDyn InitInputType CHARACTER(1
typedef ^ ^ LOGICAL UseInputFile - .TRUE. - "Supplied by Driver: .TRUE. if using a input file, .FALSE. if all inputs are being passed in by the caller" -
typedef ^ ^ FileInfoType PassedFileData - - - "If we don't use the input file, pass everything through this" -
typedef ^ ^ CHARACTER(1024) OutRootName - - - "Supplied by Driver: The name of the root file (without extension) including the full path" -
typedef ^ ^ Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." -
typedef ^ ^ Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." -
typedef ^ ^ ReKi Gravity - - - "Supplied by Driver: Gravitational acceleration" "(m/s^2)"
typedef ^ ^ DbKi TMax - - - "Supplied by Driver: The total simulation time" "(sec)"
typedef ^ ^ logical VisMeshes - .false. - "Output visualization meshes" -
Expand Down
3 changes: 1 addition & 2 deletions modules/hydrodyn/src/Morison.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4988,8 +4988,7 @@ SUBROUTINE GetSectionHstLds_Rec( origin, pos0, k_hat, x_hat, y_hat, Sa, Sb, dSad
s2 = 0.0_DbKi;
end if

dFdl(1:3) = dFdl(1:3) + &
-n(:,i) * ( z1*(s2-s1) + 0.5_DbKi*(z2-z1)/s*(s2*s2-s1*s1) )
dFdl(1:3) = dFdl(1:3) -n(:,i) * ( z1*(s2-s1) + 0.5_DbKi*(z2-z1)/s*(s2*s2-s1*s1) )
C(1) = (z2-z1)*(x2-x1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(x1-x0)/s+(x2-x1)*z1/s)*(s2*s2-s1*s1) + z1*(x1-x0)*(s2-s1)
C(2) = (z2-z1)*(y2-y1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(y1-y0)/s+(y2-y1)*z1/s)*(s2*s2-s1*s1) + z1*(y1-y0)*(s2-s1)
C(3) = (z2-z1)*(z2-z1)/3.0_DbKi/(s*s)*(s2**3-s1**3) + 0.5_DbKi*((z2-z1)*(z1-z0)/s+(z2-z1)*z1/s)*(s2*s2-s1*s1) + z1*(z1-z0)*(s2-s1)
Expand Down
114 changes: 57 additions & 57 deletions modules/moordyn/src/MoorDyn_Driver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,15 +78,14 @@ PROGRAM MoorDyn_Driver
! SeaState types
TYPE(SeaSt_InitInputType) :: InitInData_SeaSt ! Input data for initialization
TYPE(SeaSt_InitOutputType) :: InitOutData_SeaSt ! Output data from initialization
type(SeaSt_ContinuousStateType) :: x_SeaSt ! Continuous states
type(SeaSt_DiscreteStateType) :: xd_SeaSt ! Discrete states
type(SeaSt_ConstraintStateType) :: z_SeaSt ! Constraint states
type(SeaSt_OtherStateType) :: OtherState_SeaSt ! Other states
type(SeaSt_MiscVarType) :: m_SeaSt ! Misc/optimization variables
type(SeaSt_ParameterType) :: p_SeaSt ! Parameters
type(SeaSt_InputType) :: u_SeaSt(1) ! System inputs
type(SeaSt_OutputType) :: y_SeaSt ! System outputs
LOGICAL :: SeaState_Initialized = .FALSE.
type(SeaSt_ContinuousStateType) :: x_SeaSt ! Continuous states
type(SeaSt_DiscreteStateType) :: xd_SeaSt ! Discrete states
type(SeaSt_ConstraintStateType) :: z_SeaSt ! Constraint states
type(SeaSt_OtherStateType) :: OtherState_SeaSt ! Other states
type(SeaSt_MiscVarType) :: m_SeaSt ! Misc/optimization variables
type(SeaSt_ParameterType) :: p_SeaSt ! Parameters
type(SeaSt_InputType) :: u_SeaSt(1) ! System inputs
type(SeaSt_OutputType) :: y_SeaSt ! System outputs

! Motion file parsing
type(FileInfoType) :: FileInfo_PrescribeMtn !< The derived type for holding the prescribed forces input file for parsing -- we may pass this in the future
Expand Down Expand Up @@ -138,6 +137,7 @@ PROGRAM MoorDyn_Driver
ErrStat = ErrID_None
UnEcho=-1 ! set to -1 as echo is no longer used by MD
UnIn =-1


! TODO: Sort out error handling (two sets of flags currently used)

Expand Down Expand Up @@ -210,8 +210,6 @@ PROGRAM MoorDyn_Driver

! allocate Input and Output arrays; used for interpolation and extrapolation
Allocate(MD_uTimes(MD_interp_order + 1))

! @bonnie : This is in the FAST developers glue code example, but it's probably not needed here.
Allocate(MD_u(MD_interp_order + 1))


Expand All @@ -237,8 +235,8 @@ PROGRAM MoorDyn_Driver
InitInData_SeaSt%TMax = MD_InitInp%TMax
InitInData_SeaSt%Linearize = MD_InitInp%Linearize

CALL SeaSt_Init( InitInData_SeaSt, u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, dtC, InitOutData_SeaSt, ErrStat2, ErrMsg2 ); call AbortIfFailed()
SeaState_Initialized = .TRUE.
CALL SeaSt_Init( InitInData_SeaSt, u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, dtC, InitOutData_SeaSt, ErrStat2, ErrMsg2 )
call AbortIfFailed()

IF ( dtC /= drvrInitInp%dtC) THEN
ErrMsg2 = 'The SeaState Module attempted to change the coupling timestep, but this is not allowed. The SeaState Module must use the Driver coupling timestep.'
Expand All @@ -252,13 +250,10 @@ PROGRAM MoorDyn_Driver
END IF

! call the initialization routine
CALL MD_Init( MD_InitInp, MD_u(1), MD_p, MD_x , MD_xd, MD_xc, MD_xo, MD_y, MD_m, dtC, MD_InitOut, ErrStat2, ErrMsg2 ); call AbortIfFailed()
CALL MD_Init( MD_InitInp, MD_u(1), MD_p, MD_x , MD_xd, MD_xc, MD_xo, MD_y, MD_m, dtC, MD_InitOut, ErrStat2, ErrMsg2 )
call AbortIfFailed()

CALL MD_DestroyInitInput ( MD_InitInp , ErrStat2, ErrMsg2 ); call AbortIfFailed()
CALL MD_DestroyInitOutput ( MD_InitOut , ErrStat2, ErrMsg2 ); call AbortIfFailed()

CALL DispNVD( MD_InitOut%Ver )

CALL DispNVD( MD_InitOut%Ver )

! determine number of input channels expected from driver input file time series (DOFs including active tensioning channels)
if (allocated(MD_u(1)%DeltaL)) then
Expand Down Expand Up @@ -693,50 +688,47 @@ PROGRAM MoorDyn_Driver
CALL RunTimes( ProgStrtTime, ProgStrtCPU, SimStrtTime, SimStrtCPU, t )

! Destroy all objects
IF (SeaState_Initialized) THEN
CALL SeaSt_End( u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, ErrStat2, ErrMsg2); call AbortIfFailed()
ENDIF
CALL MD_End( MD_u(1), MD_p, MD_x, MD_xd, MD_xc , MD_xo, MD_y, MD_m, ErrStat2, ErrMsg2 ); call AbortIfFailed()

do j = 2,MD_interp_order+1
call MD_DestroyInput( MD_u(j), ErrStat2, ErrMsg2)
end do

if ( ErrStat /= ErrID_None ) THEN ! Display all errors
CALL WrScr1( "Errors: " )
CALL WrScr( trim(GetErrStr(ErrStat))//': '//trim(ErrMsg) )
endif

!close (un)
call CleanUp()
call EndAndCleanUp()
CALL NormStop()


CONTAINS

!-------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE AbortIfFailed()

call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver')

if (ErrStat >= AbortErrLev) then
if (SeaState_Initialized) then
call SeaSt_End( u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, ErrStat2, ErrMsg2)
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )
end if

CALL SeaSt_DestroyInitOutput( InitOutData_SeaSt, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )
CALL SeaSt_DestroyInitInput( InitInData_SeaSt, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )

call CleanUp()
Call ProgAbort(trim(ErrMsg))
elseif ( ErrStat2 /= ErrID_None ) THEN
CALL WrScr1( trim(GetErrStr(ErrStat2))//': '//trim(ErrMsg2)//NewLine)
end if

if (ErrStat >= AbortErrLev .OR. ErrStat2 >= AbortErrLev) then
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver')

call EndAndCleanUp()
Call ProgAbort(trim(ErrMsg))
elseif ( ErrStat2 /= ErrID_None ) THEN ! print messages as we get them (but don't call SetErrStat or they will be printed 2x)
CALL WrScr1( trim(GetErrStr(ErrStat2))//': '//trim(ErrMsg2)//NewLine)
end if

END SUBROUTINE AbortIfFailed
!-------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE EndAndCleanUp()
CALL SeaSt_DestroyInitOutput( InitOutData_SeaSt, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )
CALL SeaSt_DestroyInitInput( InitInData_SeaSt, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )

CALL MD_DestroyInitOutput( MD_InitOut, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )
CALL MD_DestroyInitInput( MD_InitInp, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )

call SeaSt_End( u_SeaSt(1), p_SeaSt, x_SeaSt, xd_SeaSt, z_SeaSt, OtherState_SeaSt, y_SeaSt, m_SeaSt, ErrStat2, ErrMsg2)
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )

CALL MD_End( MD_u(1), MD_p, MD_x, MD_xd, MD_xc , MD_xo, MD_y, MD_m, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )

do j = 2,MD_interp_order+1
call MD_DestroyInput( MD_u(j), ErrStat2, ErrMsg2)
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'MoorDyn_Driver' )
end do

SUBROUTINE CleanUp()
if(UnEcho >0) CLOSE( UnEcho )
if(UnIn >0) CLOSE( UnIn )

Expand All @@ -749,7 +741,15 @@ SUBROUTINE CleanUp()
IF (ALLOCATED(rd_in2 )) DEALLOCATE(rd_in2 )
IF (ALLOCATED(rdd_in )) DEALLOCATE(rdd_in )
IF (ALLOCATED(rdd_in2 )) DEALLOCATE(rdd_in2 )
END SUBROUTINE CleanUp
IF (ALLOCATED(TmpRe )) DEALLOCATE(TmpRe )



if ( ErrStat /= ErrID_None ) THEN ! Display all errors
CALL WrScr1( "Errors: " )
CALL WrScr( trim(GetErrStr(ErrStat))//': '//trim(ErrMsg) )
endif
END SUBROUTINE EndAndCleanUp

!-------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE ReadDriverInputFile( inputFile, InitInp)
Expand Down
19 changes: 5 additions & 14 deletions modules/moordyn/src/MoorDyn_Misc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1071,7 +1071,6 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
REAL(SiKi) :: WaveDOmega ! frequency step
REAL(SiKi), ALLOCATABLE :: SinWaveDir(:) ! SIN( WaveDirArr(I) ) -- Each wave frequency has a unique wave direction.
REAL(SiKi), ALLOCATABLE :: CosWaveDir(:) ! COS( WaveDirArr(I) ) -- Each wave frequency has a unique wave direction.
LOGICAL :: WaveMultiDir = .FALSE. ! Indicates the waves are multidirectional -- set by WaveField pointer if enabled

REAL(SiKi), ALLOCATABLE :: TmpFFTWaveElev(:) ! Data for the FFT calculation
TYPE(FFT_DataType) :: FFT_Data ! the instance of the FFT module we're using
Expand All @@ -1083,7 +1082,6 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
COMPLEX(SiKi), PARAMETER :: ImagNmbr = (0.0,1.0) ! The imaginary number, SQRT(-1.0)
COMPLEX(SiKi) :: ImagOmega ! = ImagNmbr*Omega (rad/s)
REAL(DbKi), ALLOCATABLE :: WaveNmbr(:) ! wave number for frequency array
REAL(SiKi), ALLOCATABLE :: WaveDirArr(:) ! Wave direction array. Each frequency has a unique direction of WaveNDir > 1 (degrees). 0's for WaveKin = 1 or if disabled in SeaState.
REAL(SiKi), ALLOCATABLE :: WaveElevC0(:,:) ! Discrete Fourier transform of the instantaneous elevation of incident waves at the ref point (meters)
COMPLEX(SiKi), ALLOCATABLE :: WaveElevC( :) ! Discrete Fourier transform of the instantaneous elevation of incident waves at the ref point (meters)
COMPLEX(SiKi), ALLOCATABLE :: WaveAccCHx(:) ! Discrete Fourier transform of the instantaneous horizontal acceleration in x-direction of incident waves before applying stretching at the zi-coordinates for points (m/s^2)
Expand Down Expand Up @@ -1412,7 +1410,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
END IF

! Warning check to make sure SeaState and MoorDyn have the same wave dir. For now, no wave spreading. This can be updated
IF (p%WaveField%WaveDir /= WaveDir) THEN
IF (p%WaveField%WaveDir /= WaveDir) THEN !bjj: the local WaveDir doesn't appear to be used when WaveKin = 2 and p%WaveField%WaveDirArr is true, so I don't think this error message is completely valid.
IF (p%writeLog > 0) THEN
WRITE(p%UnLog, '(A)' ) " WARNING SeaState WaveDir does not match MoorDyn WaveDir. Using MoorDyn values for interpolating SeaState data to MoorDyn grid."
ENDIF
Expand All @@ -1422,7 +1420,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
! Info check for if MoorDyn dtWave is non-zero. Users may set this accidentially, or could be left over in an input file
IF (p%dtWave > 0) THEN
IF (p%writeLog > 0) THEN
WRITE(p%UnLog, '(A)' ) " MoorDyn dtWave is ignored when using WaveKinMod = 2 becasue wave frequency information is supplied by SeaState"
WRITE(p%UnLog, '(A)' ) " MoorDyn dtWave is ignored when using WaveKinMod = 2 because wave frequency information is supplied by SeaState"
ENDIF
END IF

Expand All @@ -1445,15 +1443,8 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
NStepWave = p%WaveField%NStepWave

! Pull some other things out of the WaveField pointer
WaveMultiDir = p%WaveField%WaveMultiDir
p%ntWave = NStepWave ! set ntWave to NStepWave

! Set wave spreading array if enabled in SeaState, otherwise set to zero
If (WaveMultiDir) THEN
! Note: allocations not needed here because they are already allocated in SeaState
WaveDirArr = p%WaveField%WaveDirArr
ENDIF

ELSEIF (p%WaveKin == 1) THEN ! must be a filepath therefore read wave elevations from timeseries

! NOTE: there is a decent ammount of code duplication (intentional for now) with what is in SeaState that eventually
Expand Down Expand Up @@ -1659,9 +1650,9 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg)
ALLOCATE ( SinWaveDir(0:NStepWave2), STAT=ErrStat2); ErrMsg2 = 'Cannot allocate SinWaveDir.'; IF (Failed0()) RETURN

! Set the CosWaveDir and SinWaveDir values.
IF (WaveMultiDir) THEN ! This is only possible with WaveKinMod = 2
CosWaveDir=COS(D2R*WaveDirArr)
SinWaveDir=SIN(D2R*WaveDirArr)
IF (p%WaveKin == 2 .and. p%WaveField%WaveMultiDir) THEN ! This is only possible with WaveKinMod = 2
CosWaveDir=COS(D2R*p%WaveField%WaveDirArr)
SinWaveDir=SIN(D2R*p%WaveField%WaveDirArr)
ELSE
CosWaveDir=COS(D2R*WaveDir)
SinWaveDir=SIN(D2R*WaveDir)
Expand Down
2 changes: 1 addition & 1 deletion modules/nwtc-library/src/ModMesh.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2928,7 +2928,7 @@ SUBROUTINE PackMotionMesh(M, Ary, indx_first, FieldMask, TrimOP)

if (PackForTrimSolution) then
do i=1,M%NNodes
call DCM_logMap(M%Orientation(:,:,i), logmap, ErrStat2, ErrMsg2) !NOTE: we cannot use GetSmllRotAngs because we CANNOT assume that all DCMs in the code are small.
logmap = EulerExtract(M%Orientation(:,:,i)) !NOTE: we cannot use GetSmllRotAngs because we CANNOT assume that all DCMs in the code are small.
do k=1,3
Ary(indx_first) = logmap(k)
indx_first = indx_first + 1
Expand Down
Loading