From 7c38ef86eb8898e3ce335e4a61bcde3a2680e671 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Wed, 17 Sep 2025 19:55:15 -0600 Subject: [PATCH 1/9] Added a set of generic uniform grid interpolation tools * Currently only supports 3D and 4D interpolation. Can add 2D interpolation in the future if needed. * Supports interpolation of scalar data, 3-by-1 vector data, and 6-by-1 vector data. * Will use C1 continuous cubic interpolation if possible. Toward grid boundary, quadratic interpolation is used. Linear interpolation is used if only two grid points are present. * If one of the grid dimensions has only 1 grid point, this dimension is effectively ignored. * Each dimension can be indicated as either periodic or not. * Periodic dimensions will always use cubic interpolation because the grid is effectively infinite in this direction. * For non-periodic dimensions, a query point will be snapped back to the grid boundary if out of bound. A warning will be given at the first occurrence. * Note that for periodic dimensions, the subroutines assume the first and last grid points are NOT duplicates of each other. --- modules/nwtc-library/CMakeLists.txt | 4 + modules/nwtc-library/src/GridInterp.f90 | 657 ++++++++++++++++++ modules/nwtc-library/src/GridInterp.txt | 11 + modules/nwtc-library/src/GridInterp_Types.f90 | 148 ++++ 4 files changed, 820 insertions(+) create mode 100644 modules/nwtc-library/src/GridInterp.f90 create mode 100644 modules/nwtc-library/src/GridInterp.txt create mode 100644 modules/nwtc-library/src/GridInterp_Types.f90 diff --git a/modules/nwtc-library/CMakeLists.txt b/modules/nwtc-library/CMakeLists.txt index d3aced33d0..71b5a0fd9a 100644 --- a/modules/nwtc-library/CMakeLists.txt +++ b/modules/nwtc-library/CMakeLists.txt @@ -17,6 +17,7 @@ if (GENERATE_TYPES) generate_f90_types(src/Registry_NWTC_Library_base.txt ${CMAKE_CURRENT_SOURCE_DIR}/src/NWTC_Library_Types.f90 -noextrap) generate_f90_types(src/Registry_NWTC_Library_mesh.txt ${CMAKE_CURRENT_SOURCE_DIR}/src/NWTC_Library_IncSubs.f90 -incsubs -noextrap) + generate_f90_types(src/GridInterp.txt ${CMAKE_CURRENT_SOURCE_DIR}/src/GridInterp_Types.f90 -noextrap) # Generate Registry_NWTC_Library.txt by concatenating _base.txt and _mesh.txt set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS src/Registry_NWTC_Library_mesh.txt @@ -85,6 +86,9 @@ set(NWTCLIBS_SOURCES src/YAML.f90 src/JSON.f90 + src/GridInterp.f90 + src/GridInterp_Types.f90 + # RanLux sources src/ranlux/RANLUX.f90 diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 new file mode 100644 index 0000000000..294e541672 --- /dev/null +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -0,0 +1,657 @@ +MODULE GridInterp + +USE GridInterp_Types + +IMPLICIT NONE + +PRIVATE SetIndex +PRIVATE GetN1D + +PUBLIC GridInterp_SetParams +PUBLIC GridInterpSetup3D +PUBLIC GridInterpSetup4D + +INTERFACE GridInterp3D + MODULE PROCEDURE GridInterp3DR4 + MODULE PROCEDURE GridInterp3DR8 +END INTERFACE + +INTERFACE GridInterp3DVec + MODULE PROCEDURE GridInterp3DVecR4 + MODULE PROCEDURE GridInterp3DVecR8 +END INTERFACE + +INTERFACE GridInterp3DVec6 + MODULE PROCEDURE GridInterp3DVec6R4 + MODULE PROCEDURE GridInterp3DVec6R8 +END INTERFACE + +INTERFACE GridInterp4D + MODULE PROCEDURE GridInterp4DR4 + MODULE PROCEDURE GridInterp4DR8 +END INTERFACE + +INTERFACE GridInterp4DVec + MODULE PROCEDURE GridInterp4DVecR4 + MODULE PROCEDURE GridInterp4DVecR8 +END INTERFACE + +INTERFACE GridInterp4DVec6 + MODULE PROCEDURE GridInterp4DVec6R4 + MODULE PROCEDURE GridInterp4DVec6R8 +END INTERFACE + + +CONTAINS + +Subroutine GridInterp_SetParams(dim, n, delta, pZero, IsPeriodic, p, ErrStat, ErrMsg) + + Integer(IntKi), intent(in ) :: dim + Integer(IntKi), intent(in ) :: n(:) + Real(ReKi), intent(in ) :: delta(:) + Real(ReKi), intent(in ) :: pZero(:) + Logical, intent(in ) :: IsPeriodic(:) + Type(GridInterp_ParameterType), intent(inout) :: p + + Integer(IntKi), INTENT(OUT) :: ErrStat + Character(*), INTENT(OUT) :: ErrMsg + + Integer(IntKi) :: i + + ErrStat = ErrID_None + ErrMsg = "" + + if (dim/=3 .and. dim/=4) then + ErrStat = ErrID_Fatal + ErrMsg = 'GridInterp_Init: dim must be 3 or 4' + return + end if + + do i = 1,dim + p%n(i) = n(i) + p%delta(i) = delta(i) + p%pZero(i) = pZero(i) + p%IsPeriodic(i) = IsPeriodic(i) + end do + +End Subroutine GridInterp_SetParams + +Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn,ErrStat,ErrMsg) + + Real(ReKi), intent(in ) :: pIn + Real(ReKi), intent(in ) :: pZero + Real(ReKi), intent(in ) :: delta + Integer(IntKi), intent(in ) :: nMax + Logical, intent(in ) :: IsPeriodic + Integer(IntKi), intent(inout) :: Indx(:) + Real(ReKi), intent(inout) :: isopc + Integer(IntKi), intent(inout) :: Support ! = 0 for linear interpolation, = 1 for quadratic with one point to the left, = 2 for quadratic with one point to the right, = 3 for cubic + Logical, intent(inout) :: FirstWarn + Integer(IntKi), intent( out) :: ErrStat + Character(*), intent( out) :: ErrMsg + + Real(ReKi) :: p, pMax + Integer(IntKi) :: i + + ErrStat = ErrID_None + ErrMsg = "" + + isopc = 0.0_ReKi + Indx = 0_IntKi + + if ( nMax .EQ. 1_IntKi ) then ! Only one grid point, effectively ignore this dimension + ! Construct a dummy linear interpolation for now + Indx(1) = -1_IntKi + Indx(2) = 1_IntKi + Indx(3) = 1_IntKi + Indx(4) = -1_IntKi + isopc = 0.5_ReKi + Support = 0 + return + end if + + ! Compute normalized coordinate + p = (pIn-pZero) / delta + + if (isPeriodic) then + + ! Calculate normalized coordinate between 0 and 1 + isopc = p - floor(p,ReKi) + + ! Get the normalized coordinates of the two nearest nodes to the left and to the right + Indx = floor( p, IntKi ) + (/-1,0,1,2/) + + ! Make sure the coordinates are not out of bound using periodicity + do i = 1,4 + if ( Indx(i) < 0 ) then + Indx(i) = Indx(i) - nMax*floor( Real(Indx(i),ReKi) / Real(nMax,ReKi) ) + end if + Indx(i) = mod(Indx(i),nMax) + end do + + else + + pMax = Real(nMax-1_IntKi,ReKi) + + if (p<0) then + + p = 0.0_ReKi + + if (FirstWarn) then + call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetIndex') + FirstWarn = .false. + end if + + else if (p>pMax) then + + p = pMax + + if (FirstWarn) then ! don't warn if we are exactly at the boundary + call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetIndex') + FirstWarn = .false. + end if + + end if + + if ( EqualRealNos( p, pMax ) ) then + ! Calculate normalized coordinate between 0 and 1 + isopc = 1.0_ReKi + ! Get the normalized coordinates of the two nearest nodes to the left and to the right + Indx = nMax + (/-3,-2,-1,0/) + else + ! Calculate normalized coordinate between 0 and 1 + isopc = p - floor(p,ReKi) + ! Get the normalized coordinates of the two nearest nodes to the left and to the right + Indx = floor( p, IntKi ) + (/-1,0,1,2/) + end if + + end if + + ! Supported interpolation method + if ( Indx(1) < 0_IntKi ) then + if (Indx(4) > (nMax-1_IntKi) ) then + support = 0 ! Linear interpolation + else + support = 1 ! Quadratic interpolation with only one node to the left + end if + else if ( Indx(4) > (nMax-1_IntKi) ) then + support = 2 ! Quadratic interpolation with only one node to the right + else + support = 3 ! Cubic interpolation + end if + +End Subroutine SetIndex + +Subroutine GetN1D(isopc, support, N1D) + + real(ReKi), intent(in ) :: isopc ! isoparametric coordinates + integer(IntKi), intent(in ) :: support + real(ReKi), intent(inout) :: N1D(4) + real(ReKi) :: isopc2,isopc3 + + if ( Support == 3 ) then ! Cubic interpolation + + isopc2 = isopc*isopc + isopc3 = isopc*isopc2 + + N1D(1) = -0.5*isopc3 + isopc2 - 0.5*isopc + N1D(2) = 1.5*isopc3 - 2.5*isopc2 + 1.0 + N1D(3) = -1.5*isopc3 + 2.0*isopc2 + 0.5*isopc + N1D(4) = 0.5*isopc3 - 0.5*isopc2 + + else if ( Support == 1 ) then ! Quadratic interpolation with only one node to the left + + isopc2 = isopc*isopc + + N1D(1) = 0.0 + N1D(2) = 0.5*isopc2 - 1.5*isopc + 1.0 + N1D(3) = -1.0*isopc2 + 2.0*isopc + N1D(4) = 0.5*isopc2 - 0.5*isopc + + else if ( Support == 2 ) then ! Quadratic interpolation with only one node to the right + + isopc2 = isopc*isopc + + N1D(1) = 0.5*isopc2 - 0.5*isopc + N1D(2) = -1.0*isopc2 + 1.0 + N1D(3) = 0.5*isopc2 + 0.5*isopc + N1D(4) = 0.0 + + else ! if ( Support == 0 ) then ! Linear interpolation + + N1D(1) = 0.0 + N1D(2) = -isopc + 1.0 + N1D(3) = isopc + N1D(4) = 0.0 + + end if + +End Subroutine GetN1D + +Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) + + real(ReKi), intent(in ) :: Position(3) !< Array of 3 coordinates + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(inout) :: m !< MiscVars + integer(IntKi), intent( out) :: ErrStat !< Error status + character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + character(*), parameter :: RoutineName = 'GridInterpSetup3D' + integer(IntKi) :: dim,i,j,k + integer(IntKi) :: support + real(ReKi) :: N1D(4,3) + real(ReKi) :: isopc ! isoparametric coordinates + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + + ErrStat = ErrID_None + ErrMsg = "" + + do dim = 1,3 + call SetIndex(Position(dim), p%pZero(dim), p%delta(dim), p%n(dim), p%IsPeriodic(dim), m%Indx(:,dim), isopc, Support, m%FirstWarn_Clamp, ErrStat2, ErrMsg2) + if (Failed()) return; + call GetN1D(isopc, Support, N1D(:,dim)) + end do + + do k = 1,4 + do j = 1,4 + do i = 1,4 + m%N3D(i,j,k) = N1D(i,1)*N1D(j,2)*N1D(k,3) + end do + end do + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + +End Subroutine GridInterpSetup3D + +Subroutine GridInterpSetup4D( position, p, m, ErrStat, ErrMsg ) + + real(ReKi), intent(in ) :: Position(4) !< Array of 4 coordinates + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(inout) :: m !< MiscVars + integer(IntKi), intent( out) :: ErrStat !< Error status + character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + character(*), parameter :: RoutineName = 'GridInterpSetup4D' + integer(IntKi) :: dim,i,j,k,l + integer(IntKi) :: support + real(ReKi) :: N1D(4,4) + real(ReKi) :: isopc ! isoparametric coordinates + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + + ErrStat = ErrID_None + ErrMsg = "" + + do dim = 1,4 + call SetIndex(Position(dim), p%pZero(dim), p%delta(dim), p%n(dim), p%IsPeriodic(dim), m%Indx(:,dim), isopc, Support, m%FirstWarn_Clamp, ErrStat2, ErrMsg2) + if (Failed()) return; + call GetN1D(isopc, Support, N1D(:,dim)) + end do + + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + m%N4D(i,j,k,l) = N1D(i,1)*N1D(j,2)*N1D(k,3)*N1D(l,4) + end do + end do + end do + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + +End Subroutine GridInterpSetup4D + +!============================================================================================================= +! INTERFACE GridInterp3D +! - GridInterp3DR4 +! - GridInterp3DR8 +!============================================================================================================= +function GridInterp3DR4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DR4' + real(SiKi) :: GridInterp3DR4 + integer(IntKi) :: i,j,k + + ! interpolate + GridInterp3DR4 = 0.0_SiKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + GridInterp3DR4 = GridInterp3DR4 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DR4 + +function GridInterp3DR8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DR8' + real(DbKi) :: GridInterp3DR8 + integer(IntKi) :: i,j,k + + ! interpolate + GridInterp3DR8 = 0.0_DbKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + GridInterp3DR8 = GridInterp3DR8 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DR8 + +!============================================================================================================= +! INTERFACE GridInterp3DVec +! - GridInterp3DVecR4 +! - GridInterp3DVecR8 +!============================================================================================================= +function GridInterp3DVecR4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:,:) !< 3D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DVecR4' + integer(IntKi), parameter :: vDim = 3 + integer(IntKi) :: i,j,k,vi + real(SiKi) :: GridInterp3DVecR4(vDim) + + ! interpolate + GridInterp3DVecR4 = 0.0_SiKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + do vi = 1,vDim + GridInterp3DVecR4(vi) = GridInterp3DVecR4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DVecR4 + +function GridInterp3DVecR8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:,:) !< 3D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DVecR8' + integer(IntKi), parameter :: vDim = 3 + integer(IntKi) :: i,j,k,vi + real(DbKi) :: GridInterp3DVecR8(vDim) + + ! interpolate + GridInterp3DVecR8 = 0.0_DbKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + do vi = 1,vDim + GridInterp3DVecR8(vi) = GridInterp3DVecR8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DVecR8 + +!============================================================================================================= +! INTERFACE GridInterp3DVec6 +! - GridInterp3DVec6R4 +! - GridInterp3DVec6R8 +!============================================================================================================= +function GridInterp3DVec6R4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:,:) !< 3D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DVec6R4' + integer(IntKi), parameter :: vDim = 6 + integer(IntKi) :: i,j,k,vi + real(SiKi) :: GridInterp3DVec6R4(vDim) + + ! interpolate + GridInterp3DVec6R4 = 0.0_SiKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + do vi = 1,vDim + GridInterp3DVec6R4(vi) = GridInterp3DVec6R4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DVec6R4 + +function GridInterp3DVec6R8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:,:) !< 3D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp3DVec6R8' + integer(IntKi), parameter :: vDim = 6 + integer(IntKi) :: i,j,k,vi + real(DbKi) :: GridInterp3DVec6R8(vDim) + + ! interpolate + GridInterp3DVec6R8 = 0.0_DbKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then + do vi = 1,vDim + GridInterp3DVec6R8(vi) = GridInterp3DVec6R8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + +end function GridInterp3DVec6R8 + +!============================================================================================================= +! INTERFACE GridInterp4D +! - GridInterp4DR4 +! - GridInterp4DR8 +!============================================================================================================= +function GridInterp4DR4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:,0:) !< 4D grid of scalar data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DR4' + real(SiKi) :: GridInterp4DR4 + integer(IntKi) :: i,j,k,l + + ! interpolate + GridInterp4DR4 = 0.0_SiKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + GridInterp4DR4 = GridInterp4DR4 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DR4 + +function GridInterp4DR8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:,0:) !< 4D grid of scalar data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DR8' + real(DbKi) :: GridInterp4DR8 + integer(IntKi) :: i,j,k,l + + ! interpolate + GridInterp4DR8 = 0.0_DbKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + GridInterp4DR8 = GridInterp4DR8 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DR8 + +!============================================================================================================= +! INTERFACE GridInterp4DVec +! - GridInterp4DVecR4 +! - GridInterp4DVecR8 +!============================================================================================================= +function GridInterp4DVecR4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:,0:,:) !< 4D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DVecR4' + integer(IntKi), parameter :: vDim = 3 + integer(IntKi) :: i,j,k,l,vi + real(SiKi) :: GridInterp4DVecR4(vDim) + + ! interpolate + GridInterp4DVecR4 = 0.0_SiKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + do vi = 1,vDim + GridInterp4DVecR4(vi) = GridInterp4DVecR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DVecR4 + +function GridInterp4DVecR8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:,0:,:) !< 4D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DVecR8' + integer(IntKi), parameter :: vDim = 3 + integer(IntKi) :: i,j,k,l,vi + real(DbKi) :: GridInterp4DVecR8(vDim) + + ! interpolate + GridInterp4DVecR8 = 0.0_DbKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + do vi = 1,vDim + GridInterp4DVecR8(vi) = GridInterp4DVecR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DVecR8 + +!============================================================================================================= +! INTERFACE GridInterp4DVec6 +! - GridInterp4DVec6R4 +! - GridInterp4DVec6R8 +!============================================================================================================= +function GridInterp4DVec6R4( data, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:,0:,:) !< 4D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DVec6R4' + integer(IntKi), parameter :: vDim = 6 + integer(IntKi) :: i,j,k,l,vi + real(SiKi) :: GridInterp4DVec6R4(vDim) + + ! interpolate + GridInterp4DVec6R4 = 0.0_SiKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + do vi = 1,vDim + GridInterp4DVec6R4(vi) = GridInterp4DVec6R4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DVec6R4 + +function GridInterp4DVec6R8( data, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:,0:,:) !< 4D grid of vector data + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterp4DVec6R8' + integer(IntKi), parameter :: vDim = 6 + integer(IntKi) :: i,j,k,l,vi + real(DbKi) :: GridInterp4DVec6R8(vDim) + + ! interpolate + GridInterp4DVec6R8 = 0.0_DbKi + do l = 1,4 + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then + do vi = 1,vDim + GridInterp4DVec6R8(vi) = GridInterp4DVec6R8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do + ! else - Indx out of bound + end if + end do + end do + end do + end do + +end function GridInterp4DVec6R8 + +END MODULE GridInterp \ No newline at end of file diff --git a/modules/nwtc-library/src/GridInterp.txt b/modules/nwtc-library/src/GridInterp.txt new file mode 100644 index 0000000000..cf562f7cb7 --- /dev/null +++ b/modules/nwtc-library/src/GridInterp.txt @@ -0,0 +1,11 @@ +#--------------------------------------------------------------------------------------------------------------------------------------------------------- +# +#--------------------------------------------------------------------------------------------------------------------------------------------------------- +typedef GridInterp ParameterType IntKi n {4} - - "number of evenly-spaced grid points in each dimension" - +typedef ^ ParameterType ReKi delta {4} - - "size between 2 consecutive grid points in each grid direction" - +typedef ^ ParameterType ReKi pZero {4} - - "fixed position of the grid starting corner (i.e., XYZW coordinates of m%V(0,0,0,0,:))" - +typedef ^ ParameterType logical IsPeriodic {4} .false. - "flag to indicate whether this dimension should be treated as periodic" - +typedef ^ MiscVarType ReKi N3D {4}{4}{4} - - "this is the weights for 3D grid values" - +typedef ^ MiscVarType ReKi N4D {4}{4}{4}{4} - - "this is the weights for 4D grid values" - +typedef ^ MiscVarType integer Indx {4}{4} - - "this is the neighboring node index into the grid" - +typedef ^ MiscVarType logical FirstWarn_Clamp - .true. - "used to avoid too many 'Position has been clamped to the grid boundary' warning messages " - diff --git a/modules/nwtc-library/src/GridInterp_Types.f90 b/modules/nwtc-library/src/GridInterp_Types.f90 new file mode 100644 index 0000000000..f8a89950f6 --- /dev/null +++ b/modules/nwtc-library/src/GridInterp_Types.f90 @@ -0,0 +1,148 @@ +!STARTOFREGISTRYGENERATEDFILE 'GridInterp_Types.f90' +! +! WARNING This file is generated automatically by the FAST registry. +! Do not edit. Your changes to this file will be lost. +! +! FAST Registry +!********************************************************************************************************************************* +! GridInterp_Types +!................................................................................................................................. +! This file is part of GridInterp. +! +! Copyright (C) 2012-2016 National Renewable Energy Laboratory +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +! +! W A R N I N G : This file was automatically generated from the FAST registry. Changes made to this file may be lost. +! +!********************************************************************************************************************************* +!> This module contains the user-defined types needed in GridInterp. It also contains copy, destroy, pack, and +!! unpack routines associated with each defined data type. This code is automatically generated by the FAST Registry. +MODULE GridInterp_Types +!--------------------------------------------------------------------------------------------------------------------------------- +USE NWTC_Library +IMPLICIT NONE +! ========= GridInterp_ParameterType ======= + TYPE, PUBLIC :: GridInterp_ParameterType + INTEGER(IntKi) , DIMENSION(1:4) :: n = 0_IntKi !< number of evenly-spaced grid points in each dimension [-] + REAL(ReKi) , DIMENSION(1:4) :: delta = 0.0_ReKi !< size between 2 consecutive grid points in each grid direction [-] + REAL(ReKi) , DIMENSION(1:4) :: pZero = 0.0_ReKi !< fixed position of the grid starting corner (i.e., XYZW coordinates of m%V(0,0,0,0,:)) [-] + LOGICAL , DIMENSION(1:4) :: IsPeriodic = .false. !< flag to indicate whether this dimension should be treated as periodic [-] + END TYPE GridInterp_ParameterType +! ======================= +! ========= GridInterp_MiscVarType ======= + TYPE, PUBLIC :: GridInterp_MiscVarType + REAL(ReKi) , DIMENSION(1:4,1:4,1:4) :: N3D = 0.0_ReKi !< this is the weights for 3D grid values [-] + REAL(ReKi) , DIMENSION(1:4,1:4,1:4,1:4) :: N4D = 0.0_ReKi !< this is the weights for 4D grid values [-] + INTEGER(IntKi) , DIMENSION(1:4,1:4) :: Indx = 0_IntKi !< this is the neighboring node index into the grid [-] + LOGICAL :: FirstWarn_Clamp = .true. !< used to avoid too many 'Position has been clamped to the grid boundary' warning messages [-] + END TYPE GridInterp_MiscVarType +! ======================= +CONTAINS + +subroutine GridInterp_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) + type(GridInterp_ParameterType), intent(in) :: SrcParamData + type(GridInterp_ParameterType), intent(inout) :: DstParamData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'GridInterp_CopyParam' + ErrStat = ErrID_None + ErrMsg = '' + DstParamData%n = SrcParamData%n + DstParamData%delta = SrcParamData%delta + DstParamData%pZero = SrcParamData%pZero + DstParamData%IsPeriodic = SrcParamData%IsPeriodic +end subroutine + +subroutine GridInterp_DestroyParam(ParamData, ErrStat, ErrMsg) + type(GridInterp_ParameterType), intent(inout) :: ParamData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'GridInterp_DestroyParam' + ErrStat = ErrID_None + ErrMsg = '' +end subroutine + +subroutine GridInterp_PackParam(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GridInterp_ParameterType), intent(in) :: InData + character(*), parameter :: RoutineName = 'GridInterp_PackParam' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%n) + call RegPack(RF, InData%delta) + call RegPack(RF, InData%pZero) + call RegPack(RF, InData%IsPeriodic) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine GridInterp_UnPackParam(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GridInterp_ParameterType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'GridInterp_UnPackParam' + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%n); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%delta); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%pZero); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%IsPeriodic); if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine GridInterp_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) + type(GridInterp_MiscVarType), intent(in) :: SrcMiscData + type(GridInterp_MiscVarType), intent(inout) :: DstMiscData + integer(IntKi), intent(in ) :: CtrlCode + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'GridInterp_CopyMisc' + ErrStat = ErrID_None + ErrMsg = '' + DstMiscData%N3D = SrcMiscData%N3D + DstMiscData%N4D = SrcMiscData%N4D + DstMiscData%Indx = SrcMiscData%Indx + DstMiscData%FirstWarn_Clamp = SrcMiscData%FirstWarn_Clamp +end subroutine + +subroutine GridInterp_DestroyMisc(MiscData, ErrStat, ErrMsg) + type(GridInterp_MiscVarType), intent(inout) :: MiscData + integer(IntKi), intent( out) :: ErrStat + character(*), intent( out) :: ErrMsg + character(*), parameter :: RoutineName = 'GridInterp_DestroyMisc' + ErrStat = ErrID_None + ErrMsg = '' +end subroutine + +subroutine GridInterp_PackMisc(RF, Indata) + type(RegFile), intent(inout) :: RF + type(GridInterp_MiscVarType), intent(in) :: InData + character(*), parameter :: RoutineName = 'GridInterp_PackMisc' + if (RF%ErrStat >= AbortErrLev) return + call RegPack(RF, InData%N3D) + call RegPack(RF, InData%N4D) + call RegPack(RF, InData%Indx) + call RegPack(RF, InData%FirstWarn_Clamp) + if (RegCheckErr(RF, RoutineName)) return +end subroutine + +subroutine GridInterp_UnPackMisc(RF, OutData) + type(RegFile), intent(inout) :: RF + type(GridInterp_MiscVarType), intent(inout) :: OutData + character(*), parameter :: RoutineName = 'GridInterp_UnPackMisc' + if (RF%ErrStat /= ErrID_None) return + call RegUnpack(RF, OutData%N3D); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%N4D); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Indx); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%FirstWarn_Clamp); if (RegCheckErr(RF, RoutineName)) return +end subroutine +END MODULE GridInterp_Types +!ENDOFREGISTRYGENERATEDFILE From e9cc4ac2c1718b98467fa2a35bd6379e972e809c Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Mon, 22 Sep 2025 15:32:39 -0600 Subject: [PATCH 2/9] Bug fix in the new grid interpolation SetIndex subroutine --- modules/nwtc-library/src/GridInterp.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index 294e541672..9f30e56573 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -102,8 +102,8 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn if ( nMax .EQ. 1_IntKi ) then ! Only one grid point, effectively ignore this dimension ! Construct a dummy linear interpolation for now Indx(1) = -1_IntKi - Indx(2) = 1_IntKi - Indx(3) = 1_IntKi + Indx(2) = 0_IntKi + Indx(3) = 0_IntKi Indx(4) = -1_IntKi isopc = 0.5_ReKi Support = 0 @@ -179,7 +179,7 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn else support = 3 ! Cubic interpolation end if - +support = 0 End Subroutine SetIndex Subroutine GetN1D(isopc, support, N1D) From 02f3f66202f0d4e18bda87885e0dd7c80ed71e8b Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Mon, 22 Sep 2025 16:47:35 -0600 Subject: [PATCH 3/9] Adopted the new grid interpolation module in SeaState, HydroDyn, and MoorDyn --- modules/hydrodyn/src/Morison.txt | 2 +- modules/hydrodyn/src/Morison_Types.f90 | 10 +- modules/hydrodyn/src/SS_Excitation.txt | 2 +- modules/hydrodyn/src/SS_Excitation_Types.f90 | 10 +- modules/hydrodyn/src/WAMIT.f90 | 3 +- modules/hydrodyn/src/WAMIT.txt | 5 +- modules/hydrodyn/src/WAMIT2.f90 | 3 +- modules/hydrodyn/src/WAMIT2.txt | 5 +- modules/hydrodyn/src/WAMIT2_Types.f90 | 20 +- modules/hydrodyn/src/WAMIT_Interp.f90 | 65 +-- modules/hydrodyn/src/WAMIT_Types.f90 | 20 +- modules/moordyn/src/MoorDyn_Misc.f90 | 2 +- modules/moordyn/src/MoorDyn_Registry.txt | 2 +- modules/moordyn/src/MoorDyn_Types.f90 | 10 +- modules/seastate/src/SeaSt_WaveField.f90 | 539 +++--------------- modules/seastate/src/SeaSt_WaveField.txt | 5 +- .../seastate/src/SeaSt_WaveField_Types.f90 | 23 +- modules/seastate/src/SeaState.f90 | 43 +- modules/seastate/src/SeaState.txt | 2 +- modules/seastate/src/SeaState_Output.f90 | 2 +- modules/seastate/src/SeaState_Types.f90 | 10 +- 21 files changed, 189 insertions(+), 594 deletions(-) diff --git a/modules/hydrodyn/src/Morison.txt b/modules/hydrodyn/src/Morison.txt index 41fb1fbdc8..fe912f470c 100644 --- a/modules/hydrodyn/src/Morison.txt +++ b/modules/hydrodyn/src/Morison.txt @@ -462,7 +462,7 @@ typedef ^ ^ ReKi typedef ^ ^ ReKi V_rel_n_HiPass {:} - - "High-pass filtered normal relative flow velocity at joints" m/s typedef ^ ^ ReKi zFillGroup {:} - - "Instantaneous highest point of each filled group" m typedef ^ ^ MeshMapType VisMeshMap - - - "Mesh mapping for visualization mesh" - -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - # ..... Parameters ................................................................................................................ # Define parameters here: diff --git a/modules/hydrodyn/src/Morison_Types.f90 b/modules/hydrodyn/src/Morison_Types.f90 index 44f27a819d..166dc8f5dd 100644 --- a/modules/hydrodyn/src/Morison_Types.f90 +++ b/modules/hydrodyn/src/Morison_Types.f90 @@ -530,7 +530,7 @@ MODULE Morison_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: V_rel_n_HiPass !< High-pass filtered normal relative flow velocity at joints [m/s] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: zFillGroup !< Instantaneous highest point of each filled group [m] TYPE(MeshMapType) :: VisMeshMap !< Mesh mapping for visualization mesh [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] END TYPE Morison_MiscVarType ! ======================= ! ========= Morison_ParameterType ======= @@ -4693,7 +4693,7 @@ subroutine Morison_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) call NWTC_Library_CopyMeshMapType(SrcMiscData%VisMeshMap, DstMiscData%VisMeshMap, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -4780,7 +4780,7 @@ subroutine Morison_DestroyMisc(MiscData, ErrStat, ErrMsg) end if call NWTC_Library_DestroyMeshMapType(MiscData%VisMeshMap, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -4821,7 +4821,7 @@ subroutine Morison_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%V_rel_n_HiPass) call RegPackAlloc(RF, InData%zFillGroup) call NWTC_Library_PackMeshMapType(RF, InData%VisMeshMap) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4868,7 +4868,7 @@ subroutine Morison_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%V_rel_n_HiPass); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%zFillGroup); if (RegCheckErr(RF, RoutineName)) return call NWTC_Library_UnpackMeshMapType(RF, OutData%VisMeshMap) ! VisMeshMap - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m end subroutine subroutine Morison_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) diff --git a/modules/hydrodyn/src/SS_Excitation.txt b/modules/hydrodyn/src/SS_Excitation.txt index 1372e9a823..01a573b589 100644 --- a/modules/hydrodyn/src/SS_Excitation.txt +++ b/modules/hydrodyn/src/SS_Excitation.txt @@ -43,7 +43,7 @@ typedef ^ ^ SS_Exc_ContinuousStateType # 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 INTEGER LastIndWave - 1 - "last used index in the WaveTime array" - -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - # ..... Parameters ......................... diff --git a/modules/hydrodyn/src/SS_Excitation_Types.f90 b/modules/hydrodyn/src/SS_Excitation_Types.f90 index 22ba2cac7b..aa0d25a528 100644 --- a/modules/hydrodyn/src/SS_Excitation_Types.f90 +++ b/modules/hydrodyn/src/SS_Excitation_Types.f90 @@ -73,7 +73,7 @@ MODULE SS_Excitation_Types ! ========= SS_Exc_MiscVarType ======= TYPE, PUBLIC :: SS_Exc_MiscVarType INTEGER(IntKi) :: LastIndWave = 1 !< last used index in the WaveTime array [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] END TYPE SS_Exc_MiscVarType ! ======================= ! ========= SS_Exc_ParameterType ======= @@ -494,7 +494,7 @@ subroutine SS_Exc_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = '' DstMiscData%LastIndWave = SrcMiscData%LastIndWave - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -508,7 +508,7 @@ subroutine SS_Exc_DestroyMisc(MiscData, ErrStat, ErrMsg) character(*), parameter :: RoutineName = 'SS_Exc_DestroyMisc' ErrStat = ErrID_None ErrMsg = '' - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -518,7 +518,7 @@ subroutine SS_Exc_PackMisc(RF, Indata) character(*), parameter :: RoutineName = 'SS_Exc_PackMisc' if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%LastIndWave) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -528,7 +528,7 @@ subroutine SS_Exc_UnPackMisc(RF, OutData) character(*), parameter :: RoutineName = 'SS_Exc_UnPackMisc' if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%LastIndWave); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m end subroutine subroutine SS_Exc_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) diff --git a/modules/hydrodyn/src/WAMIT.f90 b/modules/hydrodyn/src/WAMIT.f90 index bf1260f4ca..e12174a110 100644 --- a/modules/hydrodyn/src/WAMIT.f90 +++ b/modules/hydrodyn/src/WAMIT.f90 @@ -241,7 +241,7 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS ! Set up wave excitation grid - Can no longer use the WaveField parameters due to different headings ! Copy WaveField grid parameters - call SeaSt_WaveField_CopyParam(p%WaveField%GridParams, p%ExctnGridParams, 0, ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call GridInterp_CopyParam(p%WaveField%VolGridParams, p%ExctnGridParams, 0, ErrStat2, ErrMsg2); CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if ( p%ExctnDisp == 0 ) then p%ExctnGridParams%n(2:3) = 1_IntKi p%ExctnGridParams%delta(2:3) = 0.0_SiKi @@ -257,7 +257,6 @@ SUBROUTINE WAMIT_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, ErrS p%ExctnGridParams%pZero(4) = -Pi end if p%ExctnGridParams%n(4) = p%NExctnHdg+1 - p%ExctnGridParams%Z_depth = -1.0 ! Set to Z_depth to a negative value to indicate uniform "z" grid for platform heading ! This module's implementation requires that if NBodyMod = 2 or 3, then there is one instance of a WAMIT module for each body, therefore, HydroDyn may have NBody > 1, but this WAMIT module will have NBody = 1 if ( (p%NBodyMod > 1) .and. (p%NBody > 1) ) then diff --git a/modules/hydrodyn/src/WAMIT.txt b/modules/hydrodyn/src/WAMIT.txt index 0d07cfe542..f826ab565a 100644 --- a/modules/hydrodyn/src/WAMIT.txt +++ b/modules/hydrodyn/src/WAMIT.txt @@ -17,6 +17,7 @@ usefrom Conv_Radiation.txt usefrom SS_Radiation.txt usefrom SS_Excitation.txt usefrom SeaSt_WaveField.txt +usefrom GridInterp.txt typedef WAMIT/WAMIT InitInputType INTEGER NBody - - - "[>=1; only used when PotMod=1. If NBodyMod=1, the WAMIT data contains a vector of size 6*NBody x 1 and matrices of size 6*NBody x 6*NBody; if NBodyMod>1, there are NBody sets of WAMIT data each with a vector of size 6 x 1 and matrices of size 6 x 6]" - typedef ^ ^ INTEGER NBodyMod - - - "Body coupling model {1: include coupling terms between each body and NBody in HydroDyn equals NBODY in WAMIT, 2: neglect coupling terms between each body and NBODY=1 with XBODY=0 in WAMIT, 3: Neglect coupling terms between each body and NBODY=1 with XBODY=/0 in WAMIT} (switch) [only used when PotMod=1]" - @@ -95,7 +96,7 @@ typedef ^ ^ SS_Exc_Outp typedef ^ ^ Conv_Rdtn_MiscVarType Conv_Rdtn - - - "" - typedef ^ ^ Conv_Rdtn_InputType Conv_Rdtn_u - - - "" - typedef ^ ^ Conv_Rdtn_OutputType Conv_Rdtn_y - - - "" - -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - # ..... Parameters ................................................................................................................ # Define parameters here: @@ -120,7 +121,7 @@ typedef ^ ^ SS_Exc_Para typedef ^ ^ DbKi DT - - - "" - typedef ^ ^ SeaSt_WaveFieldType *WaveField - - - "Pointer to wave field" typedef ^ ^ INTEGER PtfmYMod - - - "Large yaw model" - -typedef ^ ^ SeaSt_WaveField_ParameterType ExctnGridParams - - - "Parameters of WaveExctnGrid" - +typedef ^ ^ GridInterp_ParameterType ExctnGridParams - - - "Parameters of WaveExctnGrid" - # # # ..... Inputs .................................................................................................................... diff --git a/modules/hydrodyn/src/WAMIT2.f90 b/modules/hydrodyn/src/WAMIT2.f90 index 3c38462510..ef92fab0fa 100644 --- a/modules/hydrodyn/src/WAMIT2.f90 +++ b/modules/hydrodyn/src/WAMIT2.f90 @@ -2646,7 +2646,7 @@ SUBROUTINE CheckInitInput( InitInp, p, MnDriftData, NewmanAppData, DiffQTFData, ! Set up 2nd-order wave excitation grid ! Copy WaveField grid parameters - call SeaSt_WaveField_CopyParam(InitInp%WaveField%GridParams, p%Exctn2GridParams, 0, ErrStatTmp, ErrMsgTmp); CALL SetErrStat( ErrStatTmp, ErrMsgTmp, ErrStat, ErrMsg, RoutineName) + call GridInterp_CopyParam(InitInp%WaveField%VolGridParams, p%Exctn2GridParams, 0, ErrStatTmp, ErrMsgTmp); CALL SetErrStat( ErrStatTmp, ErrMsgTmp, ErrStat, ErrMsg, RoutineName) ! x and y grids are currently not used for second-order wave excitation p%Exctn2GridParams%n(2:3) = 1_IntKi p%Exctn2GridParams%delta(2:3) = 0.0_SiKi @@ -2662,7 +2662,6 @@ SUBROUTINE CheckInitInput( InitInp, p, MnDriftData, NewmanAppData, DiffQTFData, p%Exctn2GridParams%pZero(4) = -Pi end if p%Exctn2GridParams%n(4) = p%NExctnHdg+1 - p%Exctn2GridParams%Z_depth = -1.0 ! Set to Z_depth to a negative value to indicate uniform "z" grid for platform heading !> 1. Check that we only specified one of MnDrift, NewmanApp, or DiffQTF diff --git a/modules/hydrodyn/src/WAMIT2.txt b/modules/hydrodyn/src/WAMIT2.txt index ad3ec0d6f3..9fa2e20306 100644 --- a/modules/hydrodyn/src/WAMIT2.txt +++ b/modules/hydrodyn/src/WAMIT2.txt @@ -14,6 +14,7 @@ # make sure that the file name does not have any trailing white spaces! include Registry_NWTC_Library.txt usefrom SeaSt_WaveField.txt +usefrom GridInterp.txt param WAMIT2/WAMIT2 unused INTEGER MaxWAMIT2Outputs - 6 - "" - @@ -50,7 +51,7 @@ typedef ^ ^ LOGICAL SumQTFF # e.g. indices for searching in an array, large arrays that are local variables in any routine called multiple times, etc. typedef ^ MiscVarType INTEGER LastIndWave : - - "Index for last interpolation step of 2nd order forces" - typedef ^ ^ ReKi F_Waves2 {:} - - "2nd order force from this timestep" - -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - # ..... Parameters ................................................................................................................ # Define parameters here: @@ -61,7 +62,7 @@ typedef ^ ^ INTEGER NBodyMod #The 2nd order force time series grid typedef ^ ^ SiKi WaveExctn2Grid {:}{:}{:}{:}{:} - - "Grid of time series of the resulting 2nd order force (Index 1: Time, Index 2: x, Index 3: y, Index 4: platform heading, and Index 5: load component)" (N) -typedef ^ ^ SeaSt_WaveField_ParameterType Exctn2GridParams - - - "Parameters of WaveExctn2Grid" - +typedef ^ ^ GridInterp_ParameterType Exctn2GridParams - - - "Parameters of WaveExctn2Grid" - #Flags set for dimensions to use with each method (MnDrift, NewmanApp, etc). These are stored by method because .8 files that can be used in MnDrift or NewmanApp don't have some of the dimensions. typedef ^ ^ LOGICAL MnDriftDims {6} - - "Flags for which dimensions to calculate in MnDrift calculations" - diff --git a/modules/hydrodyn/src/WAMIT2_Types.f90 b/modules/hydrodyn/src/WAMIT2_Types.f90 index 5f928caa56..ba8b8c4490 100644 --- a/modules/hydrodyn/src/WAMIT2_Types.f90 +++ b/modules/hydrodyn/src/WAMIT2_Types.f90 @@ -65,7 +65,7 @@ MODULE WAMIT2_Types TYPE, PUBLIC :: WAMIT2_MiscVarType INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: LastIndWave !< Index for last interpolation step of 2nd order forces [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: F_Waves2 !< 2nd order force from this timestep [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] END TYPE WAMIT2_MiscVarType ! ======================= ! ========= WAMIT2_ParameterType ======= @@ -73,7 +73,7 @@ MODULE WAMIT2_Types INTEGER(IntKi) :: NBody = 0_IntKi !< [>=1; only used when PotMod=1. If NBodyMod=1, the WAMIT data contains a vector of size 6*NBody x 1 and matrices of size 6*NBody x 6*NBody; if NBodyMod>1, there are NBody sets of WAMIT data each with a vector of size 6 x 1 and matrices of size 6 x 6] [-] INTEGER(IntKi) :: NBodyMod = 0_IntKi !< Body coupling model {1: include coupling terms between each body and NBody in HydroDyn equals NBODY in WAMIT, 2: neglect coupling terms between each body and NBODY=1 with XBODY=0 in WAMIT, 3: Neglect coupling terms between each body and NBODY=1 with XBODY=/0 in WAMIT} (switch) [only used when PotMod=1] [-] REAL(SiKi) , DIMENSION(:,:,:,:,:), ALLOCATABLE :: WaveExctn2Grid !< Grid of time series of the resulting 2nd order force (Index 1: Time, Index 2: x, Index 3: y, Index 4: platform heading, and Index 5: load component) [(N)] - TYPE(SeaSt_WaveField_ParameterType) :: Exctn2GridParams !< Parameters of WaveExctn2Grid [-] + TYPE(GridInterp_ParameterType) :: Exctn2GridParams !< Parameters of WaveExctn2Grid [-] LOGICAL , DIMENSION(1:6) :: MnDriftDims = .false. !< Flags for which dimensions to calculate in MnDrift calculations [-] LOGICAL , DIMENSION(1:6) :: NewmanAppDims = .false. !< Flags for which dimensions to calculate in NewmanApp calculations [-] LOGICAL , DIMENSION(1:6) :: DiffQTFDims = .false. !< Flags for which dimensions to calculate in DiffQTF calculations [-] @@ -321,7 +321,7 @@ subroutine WAMIT2_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%F_Waves2 = SrcMiscData%F_Waves2 end if - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -341,7 +341,7 @@ subroutine WAMIT2_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%F_Waves2)) then deallocate(MiscData%F_Waves2) end if - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -352,7 +352,7 @@ subroutine WAMIT2_PackMisc(RF, Indata) if (RF%ErrStat >= AbortErrLev) return call RegPackAlloc(RF, InData%LastIndWave) call RegPackAlloc(RF, InData%F_Waves2) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -366,7 +366,7 @@ subroutine WAMIT2_UnPackMisc(RF, OutData) if (RF%ErrStat /= ErrID_None) return call RegUnpackAlloc(RF, OutData%LastIndWave); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%F_Waves2); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m end subroutine subroutine WAMIT2_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) @@ -395,7 +395,7 @@ subroutine WAMIT2_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMs end if DstParamData%WaveExctn2Grid = SrcParamData%WaveExctn2Grid end if - call SeaSt_WaveField_CopyParam(SrcParamData%Exctn2GridParams, DstParamData%Exctn2GridParams, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyParam(SrcParamData%Exctn2GridParams, DstParamData%Exctn2GridParams, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstParamData%MnDriftDims = SrcParamData%MnDriftDims @@ -422,7 +422,7 @@ subroutine WAMIT2_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%WaveExctn2Grid)) then deallocate(ParamData%WaveExctn2Grid) end if - call SeaSt_WaveField_DestroyParam(ParamData%Exctn2GridParams, ErrStat2, ErrMsg2) + call GridInterp_DestroyParam(ParamData%Exctn2GridParams, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -434,7 +434,7 @@ subroutine WAMIT2_PackParam(RF, Indata) call RegPack(RF, InData%NBody) call RegPack(RF, InData%NBodyMod) call RegPackAlloc(RF, InData%WaveExctn2Grid) - call SeaSt_WaveField_PackParam(RF, InData%Exctn2GridParams) + call GridInterp_PackParam(RF, InData%Exctn2GridParams) call RegPack(RF, InData%MnDriftDims) call RegPack(RF, InData%NewmanAppDims) call RegPack(RF, InData%DiffQTFDims) @@ -459,7 +459,7 @@ subroutine WAMIT2_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%NBody); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NBodyMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%WaveExctn2Grid); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackParam(RF, OutData%Exctn2GridParams) ! Exctn2GridParams + call GridInterp_UnpackParam(RF, OutData%Exctn2GridParams) ! Exctn2GridParams call RegUnpack(RF, OutData%MnDriftDims); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NewmanAppDims); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DiffQTFDims); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/hydrodyn/src/WAMIT_Interp.f90 b/modules/hydrodyn/src/WAMIT_Interp.f90 index 90a95e3432..a41191bcf4 100644 --- a/modules/hydrodyn/src/WAMIT_Interp.f90 +++ b/modules/hydrodyn/src/WAMIT_Interp.f90 @@ -29,8 +29,9 @@ MODULE WAMIT_Interp USE NWTC_Library - use SeaSt_WaveField_Types, only: SeaSt_WaveField_ParameterType, SeaSt_WaveField_MiscVarType - use SeaSt_WaveField, only: WaveField_Interp_Setup3D, WaveField_Interp_Setup4D + use GridInterp_Types, only: GridInterp_ParameterType, GridInterp_MiscVarType + use GridInterp, only: GridInterpSetup3D, GridInterpSetup4D, GridInterp3DVec6, GridInterp4DVec6 + IMPLICIT NONE PRIVATE @@ -654,30 +655,16 @@ function WAMIT_ForceWaves_Interp_3D_vec6(Time, pos, pKinXX, WF_p, WF_m, ErrStat3 real(DbKi), intent(in ) :: Time real(ReKi), intent(in ) :: pos(2) !< position real(SiKi), intent(in ) :: pKinXX(0:,:,:,:) !< 3D Wave excitation data (SiKi for storage space reasons) - type(SeaSt_WaveField_ParameterType), intent(in ) :: WF_p !< wavefield parameters - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WF_m !< wavefield misc/optimization variables + type(GridInterp_ParameterType), intent(in ) :: WF_p !< Wave excitation grid parameters + type(GridInterp_MiscVarType), intent(inout) :: WF_m !< GridInterp misc/optimization variables integer(IntKi), intent( out) :: ErrStat3 character(*), intent( out) :: ErrMsg3 real(SiKi) :: WAMIT_ForceWaves_Interp_3D_vec6(6) - real(SiKi) :: u(8) - integer(IntKi) :: i - - ! get the bounding indices from the WaveField info (same indexing used in WAMIT) - call WaveField_Interp_Setup3D( Time, pos, WF_p, WF_m, ErrStat3, ErrMsg3 ) - - ! interpolate - do i = 1,6 - u(1) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), i ) - u(2) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), i ) - u(3) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), i ) - u(4) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), i ) - u(5) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), i ) - u(6) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), i ) - u(7) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), i ) - u(8) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), i ) - WAMIT_ForceWaves_Interp_3D_vec6(i) = SUM ( WF_m%N3D * u ) - end do + + call GridInterpSetup3D( (/Real(Time,ReKi),pos(1),pos(2)/), WF_p, WF_m, ErrStat3, ErrMsg3 ) + WAMIT_ForceWaves_Interp_3D_vec6 = GridInterp3DVec6( pKinXX, WF_m ) + end function @@ -687,38 +674,16 @@ function WAMIT_ForceWaves_Interp_4D_vec6(Time, pos, pKinXX, WF_p, WF_m, ErrStat3 real(DbKi), intent(in ) :: Time real(ReKi), intent(in ) :: pos(3) !< position real(SiKi), intent(in ) :: pKinXX(0:,:,:,:,:) !< 4D Wave excitation data (SiKi for storage space reasons) - type(SeaSt_WaveField_ParameterType), intent(in ) :: WF_p !< wavefield parameters - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WF_m !< wavefield misc/optimization variables + type(GridInterp_ParameterType), intent(in ) :: WF_p !< Wave excitation grid parameters + type(GridInterp_MiscVarType), intent(inout) :: WF_m !< GridInterp misc/optimization variables integer(IntKi), intent( out) :: ErrStat3 character(*), intent( out) :: ErrMsg3 real(SiKi) :: WAMIT_ForceWaves_Interp_4D_vec6(6) - real(SiKi) :: u(16) - integer(IntKi) :: i - - ! get the bounding indices from the WaveField info (same indexing used in WAMIT) - call WaveField_Interp_Setup4D( Time, pos, WF_p, WF_m, ErrStat3, ErrMsg3 ) - - ! interpolate - do i = 1,6 - u( 1) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), WF_m%Indx_Lo(4), i ) - u( 2) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), WF_m%Indx_Hi(4), i ) - u( 3) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), WF_m%Indx_Lo(4), i ) - u( 4) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), WF_m%Indx_Hi(4), i ) - u( 5) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), WF_m%Indx_Lo(4), i ) - u( 6) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), WF_m%Indx_Hi(4), i ) - u( 7) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), WF_m%Indx_Lo(4), i ) - u( 8) = pKinXX( WF_m%Indx_Lo(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), WF_m%Indx_Hi(4), i ) - u( 9) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), WF_m%Indx_Lo(4), i ) - u(10) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Lo(3), WF_m%Indx_Hi(4), i ) - u(11) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), WF_m%Indx_Lo(4), i ) - u(12) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Lo(2), WF_m%Indx_Hi(3), WF_m%Indx_Hi(4), i ) - u(13) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), WF_m%Indx_Lo(4), i ) - u(14) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Lo(3), WF_m%Indx_Hi(4), i ) - u(15) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), WF_m%Indx_Lo(4), i ) - u(16) = pKinXX( WF_m%Indx_Hi(1), WF_m%Indx_Hi(2), WF_m%Indx_Hi(3), WF_m%Indx_Hi(4), i ) - WAMIT_ForceWaves_Interp_4D_vec6(i) = SUM ( WF_m%N4D * u ) - end do + + call GridInterpSetup4D( (/Real(Time,ReKi),pos(1),pos(2),pos(3)/), WF_p, WF_m, ErrStat3, ErrMsg3 ) + WAMIT_ForceWaves_Interp_4D_vec6 = GridInterp4DVec6( pKinXX, WF_m ) + end function diff --git a/modules/hydrodyn/src/WAMIT_Types.f90 b/modules/hydrodyn/src/WAMIT_Types.f90 index cfa511e855..f74679ba14 100644 --- a/modules/hydrodyn/src/WAMIT_Types.f90 +++ b/modules/hydrodyn/src/WAMIT_Types.f90 @@ -109,7 +109,7 @@ MODULE WAMIT_Types TYPE(Conv_Rdtn_MiscVarType) :: Conv_Rdtn !< [-] TYPE(Conv_Rdtn_InputType) :: Conv_Rdtn_u !< [-] TYPE(Conv_Rdtn_OutputType) :: Conv_Rdtn_y !< [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] END TYPE WAMIT_MiscVarType ! ======================= ! ========= WAMIT_ParameterType ======= @@ -133,7 +133,7 @@ MODULE WAMIT_Types REAL(DbKi) :: DT = 0.0_R8Ki !< [-] TYPE(SeaSt_WaveFieldType) , POINTER :: WaveField => NULL() !< Pointer to wave field [-] INTEGER(IntKi) :: PtfmYMod = 0_IntKi !< Large yaw model [-] - TYPE(SeaSt_WaveField_ParameterType) :: ExctnGridParams !< Parameters of WaveExctnGrid [-] + TYPE(GridInterp_ParameterType) :: ExctnGridParams !< Parameters of WaveExctnGrid [-] END TYPE WAMIT_ParameterType ! ======================= ! ========= WAMIT_InputType ======= @@ -742,7 +742,7 @@ subroutine WAMIT_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) call Conv_Rdtn_CopyOutput(SrcMiscData%Conv_Rdtn_y, DstMiscData%Conv_Rdtn_y, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -786,7 +786,7 @@ subroutine WAMIT_DestroyMisc(MiscData, ErrStat, ErrMsg) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call Conv_Rdtn_DestroyOutput(MiscData%Conv_Rdtn_y, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -809,7 +809,7 @@ subroutine WAMIT_PackMisc(RF, Indata) call Conv_Rdtn_PackMisc(RF, InData%Conv_Rdtn) call Conv_Rdtn_PackInput(RF, InData%Conv_Rdtn_u) call Conv_Rdtn_PackOutput(RF, InData%Conv_Rdtn_y) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -835,7 +835,7 @@ subroutine WAMIT_UnPackMisc(RF, OutData) call Conv_Rdtn_UnpackMisc(RF, OutData%Conv_Rdtn) ! Conv_Rdtn call Conv_Rdtn_UnpackInput(RF, OutData%Conv_Rdtn_u) ! Conv_Rdtn_u call Conv_Rdtn_UnpackOutput(RF, OutData%Conv_Rdtn_y) ! Conv_Rdtn_y - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m end subroutine subroutine WAMIT_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) @@ -930,7 +930,7 @@ subroutine WAMIT_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg DstParamData%DT = SrcParamData%DT DstParamData%WaveField => SrcParamData%WaveField DstParamData%PtfmYMod = SrcParamData%PtfmYMod - call SeaSt_WaveField_CopyParam(SrcParamData%ExctnGridParams, DstParamData%ExctnGridParams, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyParam(SrcParamData%ExctnGridParams, DstParamData%ExctnGridParams, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -966,7 +966,7 @@ subroutine WAMIT_DestroyParam(ParamData, ErrStat, ErrMsg) call SS_Exc_DestroyParam(ParamData%SS_Exctn, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) nullify(ParamData%WaveField) - call SeaSt_WaveField_DestroyParam(ParamData%ExctnGridParams, ErrStat2, ErrMsg2) + call GridInterp_DestroyParam(ParamData%ExctnGridParams, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -1001,7 +1001,7 @@ subroutine WAMIT_PackParam(RF, Indata) end if end if call RegPack(RF, InData%PtfmYMod) - call SeaSt_WaveField_PackParam(RF, InData%ExctnGridParams) + call GridInterp_PackParam(RF, InData%ExctnGridParams) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1051,7 +1051,7 @@ subroutine WAMIT_UnPackParam(RF, OutData) OutData%WaveField => null() end if call RegUnpack(RF, OutData%PtfmYMod); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackParam(RF, OutData%ExctnGridParams) ! ExctnGridParams + call GridInterp_UnpackParam(RF, OutData%ExctnGridParams) ! ExctnGridParams end subroutine subroutine WAMIT_CopyInput(SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg) diff --git a/modules/moordyn/src/MoorDyn_Misc.f90 b/modules/moordyn/src/MoorDyn_Misc.f90 index a1db2de91a..79e557fc9c 100644 --- a/modules/moordyn/src/MoorDyn_Misc.f90 +++ b/modules/moordyn/src/MoorDyn_Misc.f90 @@ -1158,7 +1158,7 @@ SUBROUTINE setupWaterKin(WaterKinString, p, Tmax, ErrStat, ErrMsg) END IF ! Check for if SeaState grid does not match water depth - IF (p%WaveField%GridParams%Z_Depth /= p%WtrDpth) THEN + IF (p%WaveField%GridDepth /= p%WtrDpth) THEN IF (p%writeLog > 0) THEN WRITE(p%UnLog, '(A)' ) " INFO SeaState grid depth does not match MoorDyn water depth." ENDIF diff --git a/modules/moordyn/src/MoorDyn_Registry.txt b/modules/moordyn/src/MoorDyn_Registry.txt index 4543387108..c4bb4fc8fa 100644 --- a/modules/moordyn/src/MoorDyn_Registry.txt +++ b/modules/moordyn/src/MoorDyn_Registry.txt @@ -397,7 +397,7 @@ typedef ^ ^ DbKi BathymetryGrid {:}{:} typedef ^ ^ DbKi BathGrid_Xs {:} - - "array of x-coordinates in the bathymetry grid" typedef ^ ^ DbKi BathGrid_Ys {:} - - "array of y-coordinates in the bathymetry grid" typedef ^ ^ IntKi BathGrid_npoints {:} - - "number of grid points to describe the bathymetry grid" -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - typedef ^ ^ LOGICAL IC_gen - .FALSE. - "boolean to indicate dynamic relaxation occuring" "-" diff --git a/modules/moordyn/src/MoorDyn_Types.f90 b/modules/moordyn/src/MoorDyn_Types.f90 index 407f009b8f..6e3f83f7ee 100644 --- a/modules/moordyn/src/MoorDyn_Types.f90 +++ b/modules/moordyn/src/MoorDyn_Types.f90 @@ -436,7 +436,7 @@ MODULE MoorDyn_Types REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_Xs !< array of x-coordinates in the bathymetry grid [-] REAL(DbKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_Ys !< array of y-coordinates in the bathymetry grid [-] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: BathGrid_npoints !< number of grid points to describe the bathymetry grid [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] LOGICAL :: IC_gen = .FALSE. !< boolean to indicate dynamic relaxation occuring [-] END TYPE MD_MiscVarType ! ======================= @@ -3507,7 +3507,7 @@ subroutine MD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%BathGrid_npoints = SrcMiscData%BathGrid_npoints end if - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstMiscData%IC_gen = SrcMiscData%IC_gen @@ -3661,7 +3661,7 @@ subroutine MD_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%BathGrid_npoints)) then deallocate(MiscData%BathGrid_npoints) end if - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -3773,7 +3773,7 @@ subroutine MD_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%BathGrid_Xs) call RegPackAlloc(RF, InData%BathGrid_Ys) call RegPackAlloc(RF, InData%BathGrid_npoints) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) call RegPack(RF, InData%IC_gen) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -3920,7 +3920,7 @@ subroutine MD_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%BathGrid_Xs); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BathGrid_Ys); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%BathGrid_npoints); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m call RegUnpack(RF, OutData%IC_gen); if (RegCheckErr(RF, RoutineName)) return end subroutine diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 0c1fb951e2..4af7168ac8 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -1,6 +1,8 @@ MODULE SeaSt_WaveField +USE GridInterp USE SeaSt_WaveField_Types +USE GridInterp_Types IMPLICIT NONE @@ -13,17 +15,14 @@ MODULE SeaSt_WaveField PUBLIC WaveField_GetNodeWaveNormal PUBLIC WaveField_GetNodeWaveKin PUBLIC WaveField_GetNodeWaveVel - PUBLIC WaveField_GetWaveKin -public WaveField_Interp_Setup3D, WaveField_Interp_Setup4D - CONTAINS !-------------------- Subroutine for wave elevation ------------------! function WaveField_GetNodeWaveElev1( WaveField, WaveField_m, Time, pos, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m 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 @@ -39,9 +38,9 @@ function WaveField_GetNodeWaveElev1( WaveField, WaveField_m, Time, pos, ErrStat, ErrMsg = "" IF (ALLOCATED(WaveField%WaveElev1)) THEN - CALL WaveField_Interp_Setup3D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ) + CALL WaveField_Interp_Setup3D( Time, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - Zeta = WaveField_Interp_3D( WaveField%WaveElev1, WaveField_m ) + Zeta = GridInterp3D(WaveField%WaveElev1,WaveField_m) ELSE Zeta = 0.0_SiKi END IF @@ -53,7 +52,7 @@ end function WaveField_GetNodeWaveElev1 function WaveField_GetNodeWaveElev2( WaveField, WaveField_m, Time, pos, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m 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 @@ -69,9 +68,9 @@ function WaveField_GetNodeWaveElev2( WaveField, WaveField_m, Time, pos, ErrStat, ErrMsg = "" IF (ALLOCATED(WaveField%WaveElev2)) THEN - CALL WaveField_Interp_Setup3D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ) + CALL WaveField_Interp_Setup3D( Time, pos, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - Zeta = WaveField_Interp_3D( WaveField%WaveElev2, WaveField_m ) + Zeta = GridInterp3D(WaveField%WaveElev2,WaveField_m) ELSE Zeta = 0.0_SiKi END IF @@ -83,7 +82,7 @@ end function WaveField_GetNodeWaveElev2 FUNCTION WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, pos, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m 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 @@ -112,7 +111,7 @@ END FUNCTION WaveField_GetNodeTotalWaveElev SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + 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(in ) :: r ! Distance for central differencing @@ -153,7 +152,7 @@ END SUBROUTINE WaveField_GetNodeWaveNormal !-------------------- Subroutine for full wave field kinematics --------------------! SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m real(DbKi), intent(in ) :: Time real(ReKi), intent(in ) :: pos(3) logical, intent(in ) :: forceNodeInWater @@ -190,12 +189,12 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL nodeInWater = 1_IntKi ! Use location to obtain interpolated values of kinematics - CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) - FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) - FDynP = WaveField_Interp_4D ( WaveField%WaveDynP, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) + FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m ) + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) IF ( ALLOCATED(WaveField%WaveAccMCF) ) THEN - FAMCF(:) = WaveField_Interp_4D_Vec( WaveField%WaveAccMCF, WaveField_m ) + FAMCF(:) = GridInterp4DVec( WaveField%WaveAccMCF, WaveField_m ) END IF ELSE ! Node is above the SWL nodeInWater = 0_IntKi @@ -216,33 +215,33 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod IF ( pos(3) <= 0.0_SiKi) THEN ! Node is below the SWL - evaluate wave dynamics as usual ! Use location to obtain interpolated values of kinematics - CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) - FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) - FDynP = WaveField_Interp_4D ( WaveField%WaveDynP, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) + FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m ) + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) IF ( ALLOCATED(WaveField%WaveAccMCF) ) THEN - FAMCF(:) = WaveField_Interp_4D_Vec( WaveField%WaveAccMCF, WaveField_m ) + FAMCF(:) = GridInterp4DVec( WaveField%WaveAccMCF, WaveField_m ) END IF ELSE ! Node is above SWL - need wave stretching ! Vertical wave stretching - CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_vec( WaveField%WaveVel, WaveField_m ) - FA(:) = WaveField_Interp_4D_vec( WaveField%WaveAcc, WaveField_m ) - FDynP = WaveField_Interp_4D ( WaveField%WaveDynP, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) + FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m ) + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) IF ( ALLOCATED(WaveField%WaveAccMCF) ) THEN - FAMCF(:) = WaveField_Interp_4D_vec( WaveField%WaveAccMCF, WaveField_m ) + FAMCF(:) = GridInterp4DVec( WaveField%WaveAccMCF, WaveField_m ) END IF ! Extrapoled wave stretching IF (WaveField%WaveStMod == 2) THEN - CALL WaveField_Interp_Setup3D( Time, posXY, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = FV(:) + WaveField_Interp_3D_vec( WaveField%PWaveVel0, WaveField_m ) * pos(3) - FA(:) = FA(:) + WaveField_Interp_3D_vec( WaveField%PWaveAcc0, WaveField_m ) * pos(3) - FDynP = FDynP + WaveField_Interp_3D ( WaveField%PWaveDynP0, WaveField_m ) * pos(3) + CALL WaveField_Interp_Setup3D( Time, posXY, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = FV(:) + GridInterp3DVec( WaveField%PWaveVel0, WaveField_m ) * pos(3) + FA(:) = FA(:) + GridInterp3DVec( WaveField%PWaveAcc0, WaveField_m ) * pos(3) + FDynP = FDynP + GridInterp3D ( WaveField%PWaveDynP0, WaveField_m ) * pos(3) IF ( ALLOCATED(WaveField%WaveAccMCF) ) THEN - FAMCF(:) = FAMCF(:) + WaveField_Interp_3D_vec( WaveField%PWaveAccMCF0, WaveField_m ) * pos(3) + FAMCF(:) = FAMCF(:) + GridInterp3DVec( WaveField%PWaveAccMCF0, WaveField_m ) * pos(3) END IF END IF @@ -256,12 +255,12 @@ SUBROUTINE WaveField_GetNodeWaveKin( WaveField, WaveField_m, Time, pos, forceNod 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 WaveField_Interp_Setup4D( Time, posPrime, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) - FA(:) = WaveField_Interp_4D_Vec( WaveField%WaveAcc, WaveField_m ) - FDynP = WaveField_Interp_4D ( WaveField%WaveDynP, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, posPrime, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) + FA(:) = GridInterp4DVec( WaveField%WaveAcc, WaveField_m ) + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) IF ( ALLOCATED(WaveField%WaveAccMCF) ) THEN - FAMCF(:) = WaveField_Interp_4D_Vec( WaveField%WaveAccMCF, WaveField_m ) + FAMCF(:) = GridInterp4DVec( WaveField%WaveAccMCF, WaveField_m ) END IF END IF @@ -288,7 +287,7 @@ END SUBROUTINE WaveField_GetNodeWaveKin !-------------------- Subroutine for wave field velocity only --------------------! SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FV, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m real(DbKi), intent(in ) :: Time real(ReKi), intent(in ) :: pos(3) logical, intent(in ) :: forceNodeInWater @@ -317,8 +316,8 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod IF ( pos(3) <= 0.0_ReKi) THEN ! Node is at or below the SWL nodeInWater = 1_IntKi ! Use location to obtain interpolated values of kinematics - CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) ELSE ! Node is above the SWL nodeInWater = 0_IntKi FV(:) = 0.0 @@ -335,19 +334,19 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod IF ( pos(3) <= 0.0_SiKi) THEN ! Node is below the SWL - evaluate wave dynamics as usual ! Use location to obtain interpolated values of kinematics - CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) ELSE ! Node is above SWL - need wave stretching ! Vertical wave stretching - CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_vec( WaveField%WaveVel, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) ! Extrapoled wave stretching IF (WaveField%WaveStMod == 2) THEN - CALL WaveField_Interp_Setup3D( Time, posXY, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = FV(:) + WaveField_Interp_3D_vec( WaveField%PWaveVel0, WaveField_m ) * pos(3) + CALL WaveField_Interp_Setup3D( Time, posXY, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = FV(:) + GridInterp3DVec( WaveField%PWaveVel0, WaveField_m ) * pos(3) END IF END IF ! Node is submerged @@ -360,8 +359,8 @@ SUBROUTINE WaveField_GetNodeWaveVel( WaveField, WaveField_m, Time, pos, forceNod 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 WaveField_Interp_Setup4D( Time, posPrime, WaveField%GridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; - FV(:) = WaveField_Interp_4D_Vec( WaveField%WaveVel, WaveField_m ) + CALL WaveField_Interp_Setup4D( Time, posPrime, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FV(:) = GridInterp4DVec( WaveField%WaveVel, WaveField_m ) END IF @@ -384,7 +383,7 @@ END SUBROUTINE WaveField_GetNodeWaveVel SUBROUTINE WaveField_GetWaveKin( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, WaveElev1, WaveElev2, WaveElev, FDynP, FV, FA, FAMCF, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField - type(SeaSt_WaveField_MiscVarType), intent(inout) :: WaveField_m + type(GridInterp_MiscVarType), intent(inout) :: WaveField_m real(DbKi), intent(in ) :: Time real(ReKi), intent(in ) :: pos(:,:) logical, intent(in ) :: forceNodeInWater @@ -433,460 +432,62 @@ end subroutine WaveField_GetWaveKin ! Interpolation related functions !---------------------------------------------------------------------------------------------------- -subroutine SetCartesianXYIndex(p, pZero, delta, nMax, Indx_Lo, Indx_Hi, isopc, FirstWarn, ErrStat, ErrMsg) - REAL(ReKi), intent(in ) :: p - REAL(ReKi), intent(in ) :: pZero - REAL(ReKi), intent(in ) :: delta - INTEGER(IntKi), intent(in ) :: nMax - INTEGER(IntKi), intent(inout) :: Indx_Lo - INTEGER(IntKi), intent(inout) :: Indx_Hi - real(SiKi), intent(inout) :: isopc - logical, intent(inout) :: FirstWarn - INTEGER(IntKi), intent( out) :: ErrStat - CHARACTER(*), intent( out) :: ErrMsg - - real(ReKi) :: Tmp - - ErrStat = ErrID_None - ErrMsg = "" - - isopc = -1.0 - Indx_Lo = 0 - Indx_Hi = 0 - - if ( nMax .EQ. 1_IntKi ) then ! Only one grid point - Indx_Lo = 1_IntKi - Indx_Hi = 1_IntKi - isopc = 0_SiKi - return - end if - - Tmp = (p-pZero) / delta - Indx_Lo = INT( Tmp ) + 1 ! convert REAL to INTEGER, then add one since our grid indices start at 1, not 0 - isopc = 2.0_ReKi * (Tmp - REAL(Indx_Lo - 1, ReKi)) - 1.0_ReKi ! convert to value between -1 and 1 - - if ( Indx_Lo < 1 ) then - Indx_Lo = 1 - isopc = -1.0 - if (FirstWarn) then - call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetCartesianXYIndex') !error out if time is outside the lower bounds - FirstWarn = .false. - end if - end if - - Indx_Hi = min( Indx_Lo + 1, nMax ) ! make sure it's a valid index, zero-based - - if ( Indx_Lo >= Indx_Hi ) then - ! Need to clamp to grid boundary - if (FirstWarn .and. Indx_Lo /= Indx_Hi) then ! don't warn if we are exactly at the boundary - call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetCartesianXYIndex') !error out if time is outside the lower bounds - FirstWarn = .false. - end if - Indx_Lo = max(Indx_Hi - 1, 1) - isopc = 1.0 - end if - - !------------------------------------------------------------------------------------------------- - ! to verify that we don't extrapolate, make sure isopc is bound between -1 and 1 (effectively nearest neighbor) - !------------------------------------------------------------------------------------------------- - isopc = min( 1.0_SiKi, isopc ) - isopc = max(-1.0_SiKi, isopc ) - -end subroutine SetCartesianXYIndex - - -subroutine SetCartesianZIndex(p, z_depth, delta, nMax, Indx_Lo, Indx_Hi, isopc, FirstWarn, ErrStat, ErrMsg) - real(ReKi), intent(in ) :: p - real(ReKi), intent(in ) :: z_depth - real(ReKi), intent(in ) :: delta - integer(IntKi), intent(in ) :: nMax - integer(IntKi), intent(inout) :: Indx_Lo - integer(IntKi), intent(inout) :: Indx_Hi - real(SiKi), intent(inout) :: isopc - logical, intent(inout) :: FirstWarn - integer(IntKi), intent( out) :: ErrStat - character(*), intent( out) :: ErrMsg - - real(ReKi) :: Tmp - - ErrStat = ErrID_None - ErrMsg = "" - - isopc = -1.0 - Indx_Lo = 0 - Indx_Hi = 0 - - - !Tmp = acos(-p / z_depth) / delta - Tmp = acos( max(-1.0_ReKi, min(1.0_ReKi, 1+(p / z_depth)) ) ) / delta - Tmp = nmax - 1 - Tmp - Indx_Lo = INT( Tmp ) + 1 ! convert REAL to INTEGER, then add one since our grid indices start at 1, not 0 - isopc = 2.0_ReKi * (Tmp - REAL(Indx_Lo - 1, ReKi)) - 1.0_ReKi ! convert to value between -1 and 1 - - if ( Indx_Lo < 1 ) then - Indx_Lo = 1 - isopc = -1.0 - if (FirstWarn) then - call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetCartesianZIndex') !error out if z is outside the lower bounds - FirstWarn = .false. - end if - end if - - Indx_Hi = min( Indx_Lo + 1, nMax ) ! make sure it's a valid index, one-based - - if ( Indx_Lo >= Indx_Hi ) then - ! Need to clamp to grid boundary - if (FirstWarn .and. Indx_Lo /= Indx_Hi) then ! don't warn if we are exactly at the boundary - call SetErrStat(ErrID_Warn,'Position has been clamped to the grid boundary. Warning will not be repeated though condition may persist.',ErrStat,ErrMsg,'SetCartesianZIndex') !error out if z is outside the upper bounds - FirstWarn = .false. - end if - Indx_Lo = max(Indx_Hi - 1, 1) - isopc = 1.0 - end if - - !------------------------------------------------------------------------------------------------- - ! to verify that we don't extrapolate, make sure isopc is bound between -1 and 1 (effectively nearest neighbor) - !------------------------------------------------------------------------------------------------- - isopc = min( 1.0_SiKi, isopc ) - isopc = max(-1.0_SiKi, isopc ) - -end subroutine SetCartesianZIndex - - -subroutine SetTimeIndex(Time, deltaT, nMax, Indx_Lo, Indx_Hi, isopc, ErrStat, ErrMsg) - real(DbKi), intent(in ) :: Time !< time from the start of the simulation - real(ReKi), intent(in ) :: deltaT - integer(IntKi), intent(in ) :: nMax - integer(IntKi), intent(inout) :: Indx_Lo - integer(IntKi), intent(inout) :: Indx_Hi - real(SiKi), intent(inout) :: isopc - integer(IntKi), intent( out) :: ErrStat - character(*), intent( out) :: ErrMsg - - real(ReKi) :: Tmp - - ErrStat = ErrID_None - ErrMsg = "" - - isopc = -1.0 - Indx_Lo = 0 - Indx_Hi = 0 - if ( Time < 0.0_DbKi ) then - CALL SetErrStat(ErrID_Fatal,'Time value must be greater than or equal to zero!',ErrStat,ErrMsg,'SetTimeIndex') !error out if time is outside the lower bounds - RETURN - end if - - ! if there are no timesteps, don't proceed - if (EqualRealNos(deltaT,0.0_ReKi) .or. deltaT < 0.0_ReKi) return; - -! NOTE: nMax is the total number of time values in the grid, since this is zero-based indexing, the max index is nMax-1 -! for example: in a time grid with 11 grid points, the indices run from 0,1,2,3,4,5,6,7,8,9,10 -! for the repeating waves feature, index 10 is the same as index 0, so if Indx_Lo = 10 then we want to -! wrap it back to index 0, if Indx_Lo = 11 we want to wrap back to index 1. - - Tmp = real( (Time/ real(deltaT,DbKi)) ,ReKi) - Tmp = MOD(Tmp,real((nMax), ReKi)) - Indx_Lo = INT( Tmp ) ! convert REAL to INTEGER - - isopc = 2.0_ReKi * (Tmp - REAL(Indx_Lo , ReKi)) - 1.0_ReKi ! convert to value between -1 and 1 - - !------------------------------------------------------------------------------------------------- - ! to verify that we don't extrapolate, make sure isopc is bound between -1 and 1 (effectively nearest neighbor) - !------------------------------------------------------------------------------------------------- - isopc = min( 1.0_SiKi, isopc ) - isopc = max(-1.0_SiKi, isopc ) - - Indx_Hi = min( Indx_Lo + 1, nMax ) ! make sure it's a valid index, zero-based - -end subroutine SetTimeIndex - - !==================================================================================================== !> This routine sets up interpolation of a 3-d or 4-d dataset. !! This method is described here: http://rjwagner49.com/Mathematics/Interpolation.pdf -subroutine WaveField_Interp_Setup4D( Time, Position, p, m, ErrStat, ErrMsg ) - real(DbKi), intent(in ) :: Time !< time from the start of the simulation - real(ReKi), intent(in ) :: Position(3) !< Array of XYZ coordinates, 3 - type(SeaSt_WaveField_ParameterType), intent(in ) :: p !< Parameters - type(SeaSt_WaveField_MiscVarType), intent(inout) :: m !< MiscVars +subroutine WaveField_Interp_Setup3D( Time, Position, p, m, ErrStat, ErrMsg ) + real(DbKi), intent(in ) :: Time !< Time from the start of the simulation + real(ReKi), intent(in ) :: Position(2) !< Array of XY coordinates, 2 + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(inout) :: m !< MiscVars integer(IntKi), intent( out) :: ErrStat !< Error status character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None - - character(*), parameter :: RoutineName = 'WaveField_Interp_Setup4D' - integer(IntKi) :: i - real(SiKi) :: isopc(4) ! isoparametric coordinates + character(*), parameter :: RoutineName = 'WaveField_Interp_Setup3D' integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None ErrMsg = "" - ! Find the bounding indices for time - call SetTimeIndex(Time, p%delta(1), p%n(1), m%Indx_Lo(1), m%Indx_Hi(1), isopc(1), ErrStat2, ErrMsg2) - if (Failed()) return; - - ! Find the bounding indices for XY position - do i=2,3 ! x and y components - call SetCartesianXYIndex(Position(i-1), p%pZero(i), p%delta(i), p%n(i), m%Indx_Lo(i), m%Indx_Hi(i), isopc(i), m%FirstWarn_Clamp, ErrStat2, ErrMsg2) - if (Failed()) return; - enddo - - ! Find the bounding indices for Z position - i=4 ! z component - if (p%Z_Depth>0) then - call SetCartesianZIndex(Position(i-1), p%Z_Depth, p%delta(i), p%n(i), m%Indx_Lo(i), m%Indx_Hi(i), isopc(i), m%FirstWarn_Clamp, ErrStat2, ErrMsg2) - if (Failed()) return; - else ! Regular z-grid - call SetCartesianXYIndex(Position(i-1), p%pZero(i), p%delta(i), p%n(i), m%Indx_Lo(i), m%Indx_Hi(i), isopc(i), m%FirstWarn_Clamp, ErrStat2, ErrMsg2) - if (Failed()) return; - end if - - ! compute weighting factors - m%N4D( 1) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D( 2) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D( 3) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D( 4) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D( 5) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D( 6) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D( 7) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D( 8) = ( 1.0_SiKi - isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D( 9) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D(10) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D(11) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D(12) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi - isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D(13) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D(14) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi - isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D(15) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi - isopc(4) ) - m%N4D(16) = ( 1.0_SiKi + isopc(1) ) * ( 1.0_SiKi + isopc(2) ) * ( 1.0_SiKi + isopc(3) ) * ( 1.0_SiKi + isopc(4) ) - m%N4D = m%N4D / REAL( SIZE(m%N4D), SiKi ) ! normalize + CALL GridInterpSetup3D((/Real(Time,ReKi),Position(1),Position(2)/), p, m, ErrStat2, ErrMsg2 ) + if (Failed()) return; contains logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev end function -END Subroutine WaveField_Interp_Setup4D - +END Subroutine WaveField_Interp_Setup3D -subroutine WaveField_Interp_Setup3D( Time, Position, p, m, ErrStat, ErrMsg ) - real(DbKi), intent(in ) :: Time !< time from the start of the simulation - real(ReKi), intent(in ) :: Position(2) !< Array of XYZ coordinates, 3 - type(SeaSt_WaveField_ParameterType), intent(in ) :: p !< Parameters - type(SeaSt_WaveField_MiscVarType), intent(inout) :: m !< MiscVars +subroutine WaveField_Interp_Setup4D( Time, Position, GridDepth, p, m, ErrStat, ErrMsg ) + real(DbKi), intent(in ) :: Time !< Time from the start of the simulation + real(ReKi), intent(in ) :: Position(3) !< Array of XYZ coordinates, 3 + real(SiKi), intent(in ) :: GridDepth !< Depth (>0) of the wave grid below SWL + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(inout) :: m !< MiscVars integer(IntKi), intent( out) :: ErrStat !< Error status character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None - character(*), parameter :: RoutineName = 'WaveField_Interp_Setup3D' - integer(IntKi) :: i - real(SiKi) :: isopc(4) ! isoparametric coordinates + real(ReKi) :: kz + + character(*), parameter :: RoutineName = 'WaveField_Interp_Setup4D' integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 ErrStat = ErrID_None ErrMsg = "" - ! Find the bounding indices for time - call SetTimeIndex(Time, p%delta(1), p%n(1), m%Indx_Lo(1), m%Indx_Hi(1), isopc(1), ErrStat2, ErrMsg2) - if (Failed()) return; - - ! Find the bounding indices for XY position - do i=2,3 ! x and y components - call SetCartesianXYIndex(Position(i-1), p%pZero(i), p%delta(i), p%n(i), m%Indx_Lo(i), m%Indx_Hi(i), isopc(i), m%FirstWarn_Clamp, ErrStat2, ErrMsg2) - if (Failed()) return; - enddo - - ! compute weighting factors - m%N3D(1) = ( 1.0_ReKi + isopc(1) )*( 1.0_ReKi - isopc(2) )*( 1.0_ReKi - isopc(3) ) - m%N3D(2) = ( 1.0_ReKi + isopc(1) )*( 1.0_ReKi + isopc(2) )*( 1.0_ReKi - isopc(3) ) - m%N3D(3) = ( 1.0_ReKi - isopc(1) )*( 1.0_ReKi + isopc(2) )*( 1.0_ReKi - isopc(3) ) - m%N3D(4) = ( 1.0_ReKi - isopc(1) )*( 1.0_ReKi - isopc(2) )*( 1.0_ReKi - isopc(3) ) - m%N3D(5) = ( 1.0_ReKi + isopc(1) )*( 1.0_ReKi - isopc(2) )*( 1.0_ReKi + isopc(3) ) - m%N3D(6) = ( 1.0_ReKi + isopc(1) )*( 1.0_ReKi + isopc(2) )*( 1.0_ReKi + isopc(3) ) - m%N3D(7) = ( 1.0_ReKi - isopc(1) )*( 1.0_ReKi + isopc(2) )*( 1.0_ReKi + isopc(3) ) - m%N3D(8) = ( 1.0_ReKi - isopc(1) )*( 1.0_ReKi - isopc(2) )*( 1.0_ReKi + isopc(3) ) - m%N3D = m%N3D / REAL( SIZE(m%N3D), ReKi ) ! normalize + ! Map physical z-coordinate to grid index space + kz = 0.5_ReKi*Pi - acos( max( -1.0_ReKi, min( 1.0_ReKi, 1.0_ReKi + (Position(3) / GridDepth) ) ) ) + call GridInterpSetup4D( (/Real(Time,ReKi),Position(1),Position(2),kz/), p, m, ErrStat, ErrMsg ) contains logical function Failed() call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev end function -END Subroutine WaveField_Interp_Setup3D - - -!==================================================================================================== -!> This routine interpolates a 4-d dataset. -!! This method is described here: http://rjwagner49.com/Mathematics/WaveFieldolation.pdf -function WaveField_Interp_4D( pKinXX, m ) - real(SiKi), intent(in ) :: pKinXX(0:,:,:,:) - type(SeaSt_WaveField_MiscVarType), intent(in ) :: m - - real(SiKi) :: WaveField_Interp_4D - real(SiKi) :: u(16) ! size 2^n - - ! interpolate - u( 1) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4) ) - u( 2) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4) ) - u( 3) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4) ) - u( 4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4) ) - u( 5) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4) ) - u( 6) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4) ) - u( 7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4) ) - u( 8) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4) ) - u( 9) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4) ) - u(10) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4) ) - u(11) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4) ) - u(12) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4) ) - u(13) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4) ) - u(14) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4) ) - u(15) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4) ) - u(16) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4) ) - WaveField_Interp_4D = SUM ( m%N4D * u ) -end function WaveField_Interp_4D - - -!==================================================================================================== -!> This routine interpolates a 4-d dataset. -!! This method is described here: http://rjwagner49.com/Mathematics/Interpolation.pdf -function WaveField_Interp_4D_Vec( pKinXX, m) - real(SiKi), intent(in ) :: pKinXX(0:,:,:,:,:) - type(SeaSt_WaveField_MiscVarType), intent(in ) :: m !< misc vars for interpolation - - real(SiKi) :: WaveField_Interp_4D_Vec(3) - real(SiKi) :: u(16) ! size 2^n - integer(IntKi) :: iDir - - ! interpolate - do iDir = 1,3 - u( 1) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u( 2) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u( 3) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u( 4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u( 5) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u( 6) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u( 7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u( 8) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u( 9) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u(10) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u(11) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u(12) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u(13) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u(14) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u(15) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u(16) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - WaveField_Interp_4D_Vec(iDir) = SUM ( m%N4D * u ) - end do -END FUNCTION WaveField_Interp_4D_Vec - - -!==================================================================================================== -!> This routine interpolates a 4-d dataset. -!! This method is described here: http://rjwagner49.com/Mathematics/Interpolation.pdf -function WaveField_Interp_4D_Vec6( pKinXX, m) - real(SiKi), intent(in ) :: pKinXX(0:,:,:,:,:) - type(SeaSt_WaveField_MiscVarType), intent(in ) :: m !< misc vars for interpolation - - real(SiKi) :: WaveField_Interp_4D_Vec6(6) - real(SiKi) :: u(16) ! size 2^n - integer(IntKi) :: iDir - - ! interpolate - do iDir = 1,6 - u( 1) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u( 2) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u( 3) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u( 4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u( 5) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u( 6) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u( 7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u( 8) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u( 9) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u(10) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u(11) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u(12) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - u(13) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Lo(4), iDir ) - u(14) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), m%Indx_Hi(4), iDir ) - u(15) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Lo(4), iDir ) - u(16) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), m%Indx_Hi(4), iDir ) - WaveField_Interp_4D_Vec6(iDir) = SUM ( m%N4D * u ) - end do -END FUNCTION WaveField_Interp_4D_Vec6 - - -!==================================================================================================== -!> This routine interpolates a 3-d dataset with index 1 = time (zero-based indexing), 2 = x-coordinate (1-based indexing), 3 = y-coordinate (1-based indexing) -!! This method is described here: http://rjwagner49.com/Mathematics/Interpolation.pdf -!FIXME: do like the above and call the WaveField_Interp_Setup3D routine ahead -function WaveField_Interp_3D( pKinXX, m ) - real(SiKi), intent(in ) :: pKinXX(0:,:,:) !< 3D Wave elevation data (SiKi for storage space reasons) - type(SeaSt_WaveField_MiscVarType), intent(inout) :: m !< MiscVars - - character(*), parameter :: RoutineName = 'WaveField_Interp_3D' - real(SiKi) :: WaveField_Interp_3D - real(SiKi) :: u(8) - integer(IntKi) :: i - - ! interpolate - u(1) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3) ) - u(2) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3) ) - u(3) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3) ) - u(4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3) ) - u(5) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3) ) - u(6) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3) ) - u(7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3) ) - u(8) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3) ) - WaveField_Interp_3D = SUM ( m%N3D * u ) -end function WaveField_Interp_3D - - -FUNCTION WaveField_Interp_3D_VEC( pKinXX, m ) - real(SiKi), intent(in ) :: pKinXX(0:,:,:,:) !< 3D Wave excitation data (SiKi for storage space reasons) - type(SeaSt_WaveField_MiscVarType), intent(inout) :: m !< MiscVars - - character(*), parameter :: RoutineName = 'WaveField_Interp_3D_VEC' - real(SiKi) :: WaveField_Interp_3D_VEC(3) - real(SiKi) :: u(8) - integer(IntKi) :: i - - ! interpolate - do i = 1,3 - u(1) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), i ) - u(2) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), i ) - u(3) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), i ) - u(4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), i ) - u(5) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), i ) - u(6) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), i ) - u(7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), i ) - u(8) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), i ) - WaveField_Interp_3D_VEC(i) = SUM ( m%N3D * u ) - end do -end function WaveField_Interp_3D_VEC - - -function Wavefield_Interp_3D_VEC6( pKinXX, m ) - real(SiKi), intent(in ) :: pKinXX(0:,:,:,:) !< 3D Wave excitation data (SiKi for storage space reasons) - type(SeaSt_WaveField_MiscVarType), intent(inout) :: m !< Miscvars - - character(*), parameter :: RoutineName = 'Wavefield_Interp_3D_VEC6' - real(SiKi) :: Wavefield_Interp_3D_VEC6(6) - real(SiKi) :: u(8) - integer(IntKi) :: i - - ! interpolate - do i = 1,6 - u(1) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Lo(3), i ) - u(2) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Lo(3), i ) - u(3) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Lo(3), i ) - u(4) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Lo(3), i ) - u(5) = pKinXX( m%Indx_Hi(1), m%Indx_Lo(2), m%Indx_Hi(3), i ) - u(6) = pKinXX( m%Indx_Hi(1), m%Indx_Hi(2), m%Indx_Hi(3), i ) - u(7) = pKinXX( m%Indx_Lo(1), m%Indx_Hi(2), m%Indx_Hi(3), i ) - u(8) = pKinXX( m%Indx_Lo(1), m%Indx_Lo(2), m%Indx_Hi(3), i ) - Wavefield_Interp_3D_VEC6(i) = SUM ( m%N3D * u ) - end do -end function Wavefield_Interp_3D_VEC6 - +END Subroutine WaveField_Interp_Setup4D END MODULE SeaSt_WaveField diff --git a/modules/seastate/src/SeaSt_WaveField.txt b/modules/seastate/src/SeaSt_WaveField.txt index 4dc9ba0732..380de0f028 100644 --- a/modules/seastate/src/SeaSt_WaveField.txt +++ b/modules/seastate/src/SeaSt_WaveField.txt @@ -1,5 +1,6 @@ # ...... Include files ..... usefrom Current.txt +usefrom GridInterp.txt #--------------------------------------------------------------------------------------------------------------------------------------------------------- # Data structures for representing wave fields. # @@ -47,7 +48,8 @@ typedef ^ ^ SiKi PWaveVel0 typedef ^ ^ SiKi WaveElev0 {:} - - "Instantaneous elevation time-series of incident waves at the platform reference point (NOTE THAT THIS CAN GET MODIFIED IN WAMIT)" (m) typedef ^ ^ SiKi WaveElev1 {:}{:}{:} - - "First order wave elevation" (m) typedef ^ ^ SiKi WaveElev2 {:}{:}{:} - - "Second order wave elevation" (m) -typedef ^ ^ SeaSt_WaveField_ParameterType GridParams - - - "Parameters for grid spacing" (-) +typedef ^ ^ GridInterp_ParameterType SrfGridParams - - - "Parameters of the wave free surface grid needed for interpolation" - +typedef ^ ^ GridInterp_ParameterType VolGridParams - - - "Parameters of the wave field volume grid needed for interpolation" - typedef ^ ^ IntKi WaveStMod - - - "Wave stretching model" typedef ^ ^ ReKi EffWtrDpth - - - "Water depth" (-) typedef ^ ^ ReKi MSL2SWL - - - "Vertical distance from mean sea level to still water level" (m) @@ -74,4 +76,5 @@ typedef ^ ^ INTEGER WaveMod typedef ^ ^ INTEGER NStepWave - - - "Total number of frequency components = total number of time steps in the incident wave" - typedef ^ ^ INTEGER NStepWave2 - - - "NStepWave / 2" - +typedef ^ ^ SiKi GridDepth - - - "Depth (>0) of wave grid below SWL" m typedef ^ ^ Current_InitInputType Current_InitInput - - - "InitInputs in the Current Module. For coupling with MD." - diff --git a/modules/seastate/src/SeaSt_WaveField_Types.f90 b/modules/seastate/src/SeaSt_WaveField_Types.f90 index 6db5590929..de8385a0bc 100644 --- a/modules/seastate/src/SeaSt_WaveField_Types.f90 +++ b/modules/seastate/src/SeaSt_WaveField_Types.f90 @@ -32,6 +32,7 @@ MODULE SeaSt_WaveField_Types !--------------------------------------------------------------------------------------------------------------------------------- USE Current_Types +USE GridInterp_Types USE NWTC_Library IMPLICIT NONE INTEGER(IntKi), PUBLIC, PARAMETER :: WaveDirMod_None = 0 ! WaveDirMod = 0 [Directional spreading function is NONE] [-] @@ -79,7 +80,8 @@ MODULE SeaSt_WaveField_Types REAL(SiKi) , DIMENSION(:), ALLOCATABLE :: WaveElev0 !< Instantaneous elevation time-series of incident waves at the platform reference point (NOTE THAT THIS CAN GET MODIFIED IN WAMIT) [(m)] REAL(SiKi) , DIMENSION(:,:,:), ALLOCATABLE :: WaveElev1 !< First order wave elevation [(m)] REAL(SiKi) , DIMENSION(:,:,:), ALLOCATABLE :: WaveElev2 !< Second order wave elevation [(m)] - TYPE(SeaSt_WaveField_ParameterType) :: GridParams !< Parameters for grid spacing [(-)] + TYPE(GridInterp_ParameterType) :: SrfGridParams !< Parameters of the wave free surface grid needed for interpolation [-] + TYPE(GridInterp_ParameterType) :: VolGridParams !< Parameters of the wave field volume grid needed for interpolation [-] INTEGER(IntKi) :: WaveStMod = 0_IntKi !< Wave stretching model [-] REAL(ReKi) :: EffWtrDpth = 0.0_ReKi !< Water depth [(-)] REAL(ReKi) :: MSL2SWL = 0.0_ReKi !< Vertical distance from mean sea level to still water level [(m)] @@ -104,6 +106,7 @@ MODULE SeaSt_WaveField_Types INTEGER(IntKi) :: WaveMod = 0_IntKi !< Incident wave kinematics model: See valid values in SeaSt_WaveField module parameters. [-] INTEGER(IntKi) :: NStepWave = 0_IntKi !< Total number of frequency components = total number of time steps in the incident wave [-] INTEGER(IntKi) :: NStepWave2 = 0_IntKi !< NStepWave / 2 [-] + REAL(SiKi) :: GridDepth = 0.0_R4Ki !< Depth (>0) of wave grid below SWL [m] TYPE(Current_InitInputType) :: Current_InitInput !< InitInputs in the Current Module. For coupling with MD. [-] END TYPE SeaSt_WaveFieldType ! ======================= @@ -362,7 +365,10 @@ subroutine SeaSt_WaveField_CopySeaSt_WaveFieldType(SrcSeaSt_WaveFieldTypeData, D end if DstSeaSt_WaveFieldTypeData%WaveElev2 = SrcSeaSt_WaveFieldTypeData%WaveElev2 end if - call SeaSt_WaveField_CopyParam(SrcSeaSt_WaveFieldTypeData%GridParams, DstSeaSt_WaveFieldTypeData%GridParams, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyParam(SrcSeaSt_WaveFieldTypeData%SrfGridParams, DstSeaSt_WaveFieldTypeData%SrfGridParams, CtrlCode, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (ErrStat >= AbortErrLev) return + call GridInterp_CopyParam(SrcSeaSt_WaveFieldTypeData%VolGridParams, DstSeaSt_WaveFieldTypeData%VolGridParams, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstSeaSt_WaveFieldTypeData%WaveStMod = SrcSeaSt_WaveFieldTypeData%WaveStMod @@ -422,6 +428,7 @@ subroutine SeaSt_WaveField_CopySeaSt_WaveFieldType(SrcSeaSt_WaveFieldTypeData, D DstSeaSt_WaveFieldTypeData%WaveMod = SrcSeaSt_WaveFieldTypeData%WaveMod DstSeaSt_WaveFieldTypeData%NStepWave = SrcSeaSt_WaveFieldTypeData%NStepWave DstSeaSt_WaveFieldTypeData%NStepWave2 = SrcSeaSt_WaveFieldTypeData%NStepWave2 + DstSeaSt_WaveFieldTypeData%GridDepth = SrcSeaSt_WaveFieldTypeData%GridDepth call Current_CopyInitInput(SrcSeaSt_WaveFieldTypeData%Current_InitInput, DstSeaSt_WaveFieldTypeData%Current_InitInput, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -472,7 +479,9 @@ subroutine SeaSt_WaveField_DestroySeaSt_WaveFieldType(SeaSt_WaveFieldTypeData, E if (allocated(SeaSt_WaveFieldTypeData%WaveElev2)) then deallocate(SeaSt_WaveFieldTypeData%WaveElev2) end if - call SeaSt_WaveField_DestroyParam(SeaSt_WaveFieldTypeData%GridParams, ErrStat2, ErrMsg2) + call GridInterp_DestroyParam(SeaSt_WaveFieldTypeData%SrfGridParams, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call GridInterp_DestroyParam(SeaSt_WaveFieldTypeData%VolGridParams, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (allocated(SeaSt_WaveFieldTypeData%WaveElevC)) then deallocate(SeaSt_WaveFieldTypeData%WaveElevC) @@ -504,7 +513,8 @@ subroutine SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, Indata) call RegPackAlloc(RF, InData%WaveElev0) call RegPackAlloc(RF, InData%WaveElev1) call RegPackAlloc(RF, InData%WaveElev2) - call SeaSt_WaveField_PackParam(RF, InData%GridParams) + call GridInterp_PackParam(RF, InData%SrfGridParams) + call GridInterp_PackParam(RF, InData%VolGridParams) call RegPack(RF, InData%WaveStMod) call RegPack(RF, InData%EffWtrDpth) call RegPack(RF, InData%MSL2SWL) @@ -529,6 +539,7 @@ subroutine SeaSt_WaveField_PackSeaSt_WaveFieldType(RF, Indata) call RegPack(RF, InData%WaveMod) call RegPack(RF, InData%NStepWave) call RegPack(RF, InData%NStepWave2) + call RegPack(RF, InData%GridDepth) call Current_PackInitInput(RF, InData%Current_InitInput) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -553,7 +564,8 @@ subroutine SeaSt_WaveField_UnPackSeaSt_WaveFieldType(RF, OutData) call RegUnpackAlloc(RF, OutData%WaveElev0); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%WaveElev1); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%WaveElev2); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackParam(RF, OutData%GridParams) ! GridParams + call GridInterp_UnpackParam(RF, OutData%SrfGridParams) ! SrfGridParams + call GridInterp_UnpackParam(RF, OutData%VolGridParams) ! VolGridParams call RegUnpack(RF, OutData%WaveStMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%EffWtrDpth); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MSL2SWL); if (RegCheckErr(RF, RoutineName)) return @@ -578,6 +590,7 @@ subroutine SeaSt_WaveField_UnPackSeaSt_WaveFieldType(RF, OutData) call RegUnpack(RF, OutData%WaveMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NStepWave); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NStepWave2); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%GridDepth); if (RegCheckErr(RF, RoutineName)) return call Current_UnpackInitInput(RF, OutData%Current_InitInput) ! Current_InitInput end subroutine END MODULE SeaSt_WaveField_Types diff --git a/modules/seastate/src/SeaState.f90 b/modules/seastate/src/SeaState.f90 index af5d0b15bf..f3dae11151 100644 --- a/modules/seastate/src/SeaState.f90 +++ b/modules/seastate/src/SeaState.f90 @@ -32,6 +32,7 @@ MODULE SeaState USE SeaState_Output USE Current USE Waves2 + USE GridInterp IMPLICIT NONE PRIVATE @@ -260,16 +261,28 @@ SUBROUTINE SeaSt_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, Init CALL SeaStOut_WrSummaryFile(InitInp, InputFileData, p, ErrStat2, ErrMsg2); if(Failed()) return; - - ! Setup the 4D grid information for the Interpolation Module - p%WaveField%GridParams%n = (/p%WaveField%NStepWave,p%nGrid(1),p%nGrid(2),p%nGrid(3)/) - p%WaveField%GridParams%delta = (/real(p%WaveDT,ReKi),p%deltaGrid(1),p%deltaGrid(2),p%deltaGrid(3)/) - p%WaveField%GridParams%pZero(1) = 0.0 !Time - p%WaveField%GridParams%pZero(2) = -InputFileData%X_HalfWidth - p%WaveField%GridParams%pZero(3) = -InputFileData%Y_HalfWidth - p%WaveField%GridParams%pZero(4) = -InputFileData%Z_Depth ! zi - p%WaveField%GridParams%Z_Depth = InputFileData%Z_Depth + ! Setup the 3D and 4D grid information for the Interpolation Module + p%WaveField%GridDepth = InputFileData%Z_Depth + + ! Set the parameters of the 3D free surface grid + CALL GridInterp_SetParams(3_IntKi, & ! dimension + (/p%WaveField%NStepWave,p%nGrid(1),p%nGrid(2)/), & ! n + (/real(p%WaveDT,ReKi),p%deltaGrid(1),p%deltaGrid(2)/), & ! delta + (/0.0,-InputFileData%X_HalfWidth,-InputFileData%Y_HalfWidth/), & ! pZero + (/.true.,.false.,.false./), & ! periodicity + p%WaveField%SrfGridParams, ErrStat2, ErrMsg2 ) + if(Failed()) return; + + ! Set the parameters of the 4D volume grid + CALL GridInterp_SetParams(4_IntKi, & ! dimension + (/p%WaveField%NStepWave,p%nGrid(1),p%nGrid(2),p%nGrid(3)/), & ! n + (/real(p%WaveDT,ReKi),p%deltaGrid(1),p%deltaGrid(2),p%deltaGrid(3)/), & ! delta + (/0.0,-InputFileData%X_HalfWidth,-InputFileData%Y_HalfWidth,0.0/), & ! pZero + (/.true.,.false.,.false.,.false./), & ! periodicity + p%WaveField%VolGridParams, ErrStat2, ErrMsg2 ) + if(Failed()) return; + IF ( p%OutSwtch == 1 ) THEN ! Only SeaSt-level output writing ! HACK WE can tell FAST not to write any SeaState outputs by simply deallocating the WriteOutputHdr array! @@ -390,15 +403,15 @@ subroutine SurfaceVisGenerate(ErrStat3, ErrMsg3) ErrMsg3 = "" ! Grid half width from the WaveField - HWidX = (real(p%WaveField%GridParams%n(2)-1,SiKi)) / 2.0_SiKi * p%WaveField%GridParams%delta(2) - HWidY = (real(p%WaveField%GridParams%n(3)-1,SiKi)) / 2.0_SiKi * p%WaveField%GridParams%delta(3) + HWidX = (real(p%WaveField%SrfGridParams%n(2)-1,SiKi)) / 2.0_SiKi * p%WaveField%SrfGridParams%delta(2) + HWidY = (real(p%WaveField%SrfGridParams%n(3)-1,SiKi)) / 2.0_SiKi * p%WaveField%SrfGridParams%delta(3) if ((InitInp%SurfaceVisNx <= 0) .or. (InitInp%SurfaceVisNy <= 0))then ! use the SeaState points exactly ! Set number of points to the number of seastate grid points in each direction - Nx = p%WaveField%GridParams%n(2) - Ny = p%WaveField%GridParams%n(3) - dx = p%WaveField%GridParams%delta(2) - dy = p%WaveField%GridParams%delta(3) + Nx = p%WaveField%SrfGridParams%n(2) + Ny = p%WaveField%SrfGridParams%n(3) + dx = p%WaveField%SrfGridParams%delta(2) + dy = p%WaveField%SrfGridParams%delta(3) call SetErrStat(ErrID_Info,"Setting wavefield visualization grid to "//trim(Num2LStr(Nx))//" x "//trim(Num2LStr(Ny))//"points",ErrStat3,ErrMsg3,RoutineName) elseif ((InitInp%SurfaceVisNx < 3) .or. (InitInp%SurfaceVisNx < 3)) then ! Set to 3 for minimum Nx = 3 diff --git a/modules/seastate/src/SeaState.txt b/modules/seastate/src/SeaState.txt index f38dfdf231..8b2635bfed 100644 --- a/modules/seastate/src/SeaState.txt +++ b/modules/seastate/src/SeaState.txt @@ -126,7 +126,7 @@ typedef ^ OtherStateType R8Ki Unu typedef ^ MiscVarType INTEGER Decimate - - - "The output decimation counter" - typedef ^ ^ DbKi LastOutTime - - - "Last time step which was written to the output file (sec)" - typedef ^ ^ INTEGER LastIndWave - - - "The last index used in the wave kinematics arrays, used to optimize interpolation" - -typedef ^ ^ SeaSt_WaveField_MiscVarType WaveField_m - - - "misc var information from the SeaState Interpolation module" - +typedef ^ ^ GridInterp_MiscVarType WaveField_m - - - "misc var information from the Grid Interpolation module" - # .... Linearization params ....................................................................................................... # NOTE: This is overkill given how limited linearization is. For completeness and similarity to other modules, keeping all this here. Also note some diff --git a/modules/seastate/src/SeaState_Output.f90 b/modules/seastate/src/SeaState_Output.f90 index 7dacda0be8..4f56ea9339 100644 --- a/modules/seastate/src/SeaState_Output.f90 +++ b/modules/seastate/src/SeaState_Output.f90 @@ -272,7 +272,7 @@ SUBROUTINE SeaStOut_WriteWvKinFiles( Rootname, SeaSt_Prog, WaveField, WaveDT, X_ y_gridPts(i+1) = -Y_HalfWidth + deltaGrid(2)*i end do do i = 0, NGrid(3)-1 - z_gridPts(i+1) = - ( 1.0 - cos( real((NGrid(3) - 1) - i, ReKi) * deltaGrid(3) ) ) * WaveField%GridParams%Z_Depth + z_gridPts(i+1) = - ( 1.0 - cos( real((NGrid(3) - 1) - i, ReKi) * deltaGrid(3) ) ) * WaveField%GridDepth end do ! Write the increments from [0, NStepWave] even though for OpenFAST data, NStepWave = 0, but for arbitrary user data this may not be true. diff --git a/modules/seastate/src/SeaState_Types.f90 b/modules/seastate/src/SeaState_Types.f90 index bdd95ed4de..03017a6f11 100644 --- a/modules/seastate/src/SeaState_Types.f90 +++ b/modules/seastate/src/SeaState_Types.f90 @@ -146,7 +146,7 @@ MODULE SeaState_Types INTEGER(IntKi) :: Decimate = 0_IntKi !< The output decimation counter [-] REAL(DbKi) :: LastOutTime = 0.0_R8Ki !< Last time step which was written to the output file (sec) [-] INTEGER(IntKi) :: LastIndWave = 0_IntKi !< The last index used in the wave kinematics arrays, used to optimize interpolation [-] - TYPE(SeaSt_WaveField_MiscVarType) :: WaveField_m !< misc var information from the SeaState Interpolation module [-] + TYPE(GridInterp_MiscVarType) :: WaveField_m !< misc var information from the Grid Interpolation module [-] END TYPE SeaSt_MiscVarType ! ======================= ! ========= Jac_u_idxStarts ======= @@ -986,7 +986,7 @@ subroutine SeaSt_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) DstMiscData%Decimate = SrcMiscData%Decimate DstMiscData%LastOutTime = SrcMiscData%LastOutTime DstMiscData%LastIndWave = SrcMiscData%LastIndWave - call SeaSt_WaveField_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) + call GridInterp_CopyMisc(SrcMiscData%WaveField_m, DstMiscData%WaveField_m, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return end subroutine @@ -1000,7 +1000,7 @@ subroutine SeaSt_DestroyMisc(MiscData, ErrStat, ErrMsg) character(*), parameter :: RoutineName = 'SeaSt_DestroyMisc' ErrStat = ErrID_None ErrMsg = '' - call SeaSt_WaveField_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) + call GridInterp_DestroyMisc(MiscData%WaveField_m, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end subroutine @@ -1012,7 +1012,7 @@ subroutine SeaSt_PackMisc(RF, Indata) call RegPack(RF, InData%Decimate) call RegPack(RF, InData%LastOutTime) call RegPack(RF, InData%LastIndWave) - call SeaSt_WaveField_PackMisc(RF, InData%WaveField_m) + call GridInterp_PackMisc(RF, InData%WaveField_m) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1024,7 +1024,7 @@ subroutine SeaSt_UnPackMisc(RF, OutData) call RegUnpack(RF, OutData%Decimate); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%LastOutTime); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%LastIndWave); if (RegCheckErr(RF, RoutineName)) return - call SeaSt_WaveField_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m + call GridInterp_UnpackMisc(RF, OutData%WaveField_m) ! WaveField_m end subroutine subroutine SeaSt_CopyJac_u_idxStarts(SrcJac_u_idxStartsData, DstJac_u_idxStartsData, CtrlCode, ErrStat, ErrMsg) From 5db120c689e6bdb0dadac980eeee987276b650ca Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 23 Sep 2025 16:35:02 -0600 Subject: [PATCH 4/9] Updated the method for computing the free-surface normal --- modules/nwtc-library/src/GridInterp.f90 | 230 ++++++++++++++++++++++- modules/seastate/src/SeaSt_WaveField.f90 | 21 +-- 2 files changed, 237 insertions(+), 14 deletions(-) diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index 9f30e56573..6df26b738d 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -6,10 +6,12 @@ MODULE GridInterp PRIVATE SetIndex PRIVATE GetN1D +PRIVATE GetN1Ddx PUBLIC GridInterp_SetParams PUBLIC GridInterpSetup3D PUBLIC GridInterpSetup4D +PUBLIC GridInterpSetupN INTERFACE GridInterp3D MODULE PROCEDURE GridInterp3DR4 @@ -41,6 +43,15 @@ MODULE GridInterp MODULE PROCEDURE GridInterp4DVec6R8 END INTERFACE +INTERFACE GridInterpN + MODULE PROCEDURE GridInterpNR4 + MODULE PROCEDURE GridInterpNR8 +END INTERFACE + +INTERFACE GridInterpS + MODULE PROCEDURE GridInterpSR4 + MODULE PROCEDURE GridInterpSR8 +END INTERFACE CONTAINS @@ -179,7 +190,7 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn else support = 3 ! Cubic interpolation end if -support = 0 + End Subroutine SetIndex Subroutine GetN1D(isopc, support, N1D) @@ -228,6 +239,47 @@ Subroutine GetN1D(isopc, support, N1D) End Subroutine GetN1D +Subroutine GetN1Ddx(isopc, support, N1D) + + real(ReKi), intent(in ) :: isopc ! isoparametric coordinates + integer(IntKi), intent(in ) :: support + real(ReKi), intent(inout) :: N1D(4) + real(ReKi) :: isopc2 + + if ( Support == 3 ) then ! Cubic interpolation + + isopc2 = isopc*isopc + + N1D(1) = -1.5*isopc2 + 2.0*isopc - 0.5 + N1D(2) = 4.5*isopc2 - 5.0*isopc + N1D(3) = -4.5*isopc2 + 4.0*isopc + 0.5 + N1D(4) = 1.5*isopc2 - isopc + + else if ( Support == 1 ) then ! Quadratic interpolation with only one node to the left + + N1D(1) = 0.0 + N1D(2) = isopc - 1.5 + N1D(3) = -2.0*isopc + 2.0 + N1D(4) = isopc - 0.5 + + else if ( Support == 2 ) then ! Quadratic interpolation with only one node to the right + + N1D(1) = isopc - 0.5 + N1D(2) = -2.0*isopc + N1D(3) = isopc + 0.5 + N1D(4) = 0.0 + + else ! if ( Support == 0 ) then ! Linear interpolation + + N1D(1) = 0.0 + N1D(2) = -1.0 + N1D(3) = 1.0 + N1D(4) = 0.0 + + end if + +End Subroutine GetN1Ddx + Subroutine GridInterpSetup3D( position, p, m, ErrStat, ErrMsg ) real(ReKi), intent(in ) :: Position(3) !< Array of 3 coordinates @@ -312,6 +364,54 @@ logical function Failed() End Subroutine GridInterpSetup4D +Subroutine GridInterpSetupN( position, p, m, ErrStat, ErrMsg ) + + real(ReKi), intent(in ) :: Position(3) !< Array of 3 coordinates + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(inout) :: m !< MiscVars + integer(IntKi), intent( out) :: ErrStat !< Error status + character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + character(*), parameter :: RoutineName = 'GridInterpSetupN' + integer(IntKi) :: dim,i,j,k + integer(IntKi) :: support + real(ReKi) :: N1D(4,3) + real(ReKi) :: N1Ddx(4,2:3) + real(ReKi) :: isopc ! isoparametric coordinates + integer(IntKi) :: ErrStat2 + character(ErrMsgLen) :: ErrMsg2 + + ErrStat = ErrID_None + ErrMsg = "" + + do dim = 1,3 + call SetIndex(Position(dim), p%pZero(dim), p%delta(dim), p%n(dim), p%IsPeriodic(dim), m%Indx(:,dim), isopc, Support, m%FirstWarn_Clamp, ErrStat2, ErrMsg2) + if (Failed()) return; + call GetN1D(isopc, Support, N1D(:,dim)) + if (dim>1) then + call GetN1Ddx(isopc, Support, N1Ddx(:,dim)) + end if + end do + + ! Need two sets of weights for d(.)/dx and d(.)/dy. Borrow m%N4D for this. + do k = 1,4 + do j = 1,4 + do i = 1,4 + m%N4D(i,j,k,1) = N1D(i,1)*N1Ddx(j,2)*N1D (k,3) + m%N4D(i,j,k,2) = N1D(i,1)*N1D (j,2)*N1Ddx(k,3) + end do + end do + end do + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function + +End Subroutine GridInterpSetupN + + !============================================================================================================= ! INTERFACE GridInterp3D ! - GridInterp3DR4 @@ -654,4 +754,132 @@ function GridInterp4DVec6R8( data, m ) end function GridInterp4DVec6R8 +!============================================================================================================= +! INTERFACE GridInterpN +! - GridInterpNR4 +! - GridInterpNR8 +!============================================================================================================= +function GridInterpNR4( data, p, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterpNR4' + real(SiKi) :: GridInterpNR4(3) + real(SiKi) :: dZetadx, dZetady + integer(IntKi) :: i,j,k + + ! interpolate slope + dZetadx = 0.0_SiKi + dZetady = 0.0_SiKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,1),0.0_ReKi)) then + dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + if (.not.EqualRealNos(m%N4D(i,j,k,2),0.0_ReKi)) then + dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + end do + end do + end do + dZetadx = dZetadx / p%delta(2) + dZetady = dZetady / p%delta(3) + + GridInterpNR4 = (/-dZetadx,-dZetady,1.0_SiKi/) + GridInterpNR4 = GridInterpNR4 / TwoNorm(GridInterpNR4) + +end function GridInterpNR4 + +function GridInterpNR8( data, p, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterpNR8' + real(DbKi) :: GridInterpNR8(3) + real(DbKi) :: dZetadx, dZetady + integer(IntKi) :: i,j,k + + ! interpolate slope + dZetadx = 0.0_DbKi + dZetady = 0.0_DbKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + if (.not.EqualRealNos(m%N4D(i,j,k,1),0.0_ReKi)) then + dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + if (.not.EqualRealNos(m%N4D(i,j,k,2),0.0_ReKi)) then + dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + end do + end do + end do + dZetadx = dZetadx / p%delta(2) + dZetady = dZetady / p%delta(3) + + GridInterpNR8 = (/-dZetadx,-dZetady,1.0_DbKi/) + GridInterpNR8 = GridInterpNR8 / TwoNorm(GridInterpNR8) + +end function GridInterpNR8 + +!============================================================================================================= +! INTERFACE GridInterpS +! - GridInterpSR4 +! - GridInterpSR8 +!============================================================================================================= +function GridInterpSR4( data, p, m ) + real(SiKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterpSR4' + real(SiKi) :: GridInterpSR4(2) + integer(IntKi) :: i,j,k,dir + + ! interpolate slope + GridInterpSR4 = 0.0_SiKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + do dir = 1,2 + if (.not.EqualRealNos(m%N4D(i,j,k,dir),0.0_ReKi)) then + GridInterpSR4(dir) = GridInterpSR4(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + end do + end do + end do + end do + GridInterpSR4 = GridInterpSR4 / p%delta(2:3) + +end function GridInterpSR4 + +function GridInterpSR8( data, p, m ) + real(DbKi), intent(in ) :: data(0:,0:,0:) !< 3D grid of scalar data + type(GridInterp_ParameterType), intent(in ) :: p !< Parameters + type(GridInterp_MiscVarType), intent(in ) :: m !< MiscVars + + character(*), parameter :: RoutineName = 'GridInterpSR8' + real(DbKi) :: GridInterpSR8(2) + integer(IntKi) :: i,j,k,dir + + ! interpolate slope + GridInterpSR8 = 0.0_DbKi + do k = 1,4 + do j = 1,4 + do i = 1,4 + do dir = 1,2 + if (.not.EqualRealNos(m%N4D(i,j,k,dir),0.0_ReKi)) then + GridInterpSR8(dir) = GridInterpSR8(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + end if + end do + end do + end do + end do + GridInterpSR8 = GridInterpSR8 / p%delta(2:3) + +end function GridInterpSR8 + END MODULE GridInterp \ No newline at end of file diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 4af7168ac8..6c39bf8535 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -119,8 +119,7 @@ SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, integer(IntKi), intent( out) :: ErrStat ! Error status of the operation character(*), intent( out) :: ErrMsg ! Error message if errStat /= ErrID_None - real(SiKi) :: ZetaP,ZetaM - real(ReKi) :: r1,dZetadx,dZetady + real(SiKi) :: slope(2) character(*), parameter :: RoutineName = 'WaveField_GetNodeWaveNormal' integer(IntKi) :: errStat2 character(ErrMsgLen) :: errMsg2 @@ -128,18 +127,14 @@ SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, ErrStat = ErrID_None ErrMsg = "" - r1 = MAX(r,real(1.0e-6,ReKi)) ! In case r is zero + call GridInterpSetupN( (/Real(Time,ReKi),pos(1),pos(2)/), WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ) + slope = GridInterpS( WaveField%WaveElev1, WaveField%SrfGridParams, WaveField_m ) + if (ALLOCATED(WaveField%WaveElev2)) then + slope = slope + GridInterpS( WaveField%WaveElev2, WaveField%SrfGridParams, WaveField_m ) + end if - ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)+r1,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; - ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1)-r1,pos(2)/), ErrStat2, ErrMsg2 ); if (Failed()) return; - dZetadx = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) - - ZetaP = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)+r1/), ErrStat2, ErrMsg2 ); if (Failed()) return; - ZetaM = WaveField_GetNodeTotalWaveElev( WaveField, WaveField_m, Time, (/pos(1),pos(2)-r1/), ErrStat2, ErrMsg2 ); if (Failed()) return; - dZetady = REAL(ZetaP-ZetaM,ReKi)/(2.0_ReKi*r1) - - n = (/-dZetadx,-dZetady,1.0_ReKi/) - n = n / SQRT(Dot_Product(n,n)) + n = Real( (/-slope(1),-slope(2),1.0_SiKi/), ReKi) + n = n / TwoNorm(n) contains logical function Failed() From 1d0ee6c2c4df01a84437e85e8a65f36e68e5b258 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Tue, 23 Sep 2025 17:13:01 -0600 Subject: [PATCH 5/9] Cleanup --- modules/hydrodyn/src/Morison.f90 | 27 +++++++----------------- modules/seastate/src/SeaSt_WaveField.f90 | 3 +-- 2 files changed, 9 insertions(+), 21 deletions(-) diff --git a/modules/hydrodyn/src/Morison.f90 b/modules/hydrodyn/src/Morison.f90 index fb7003e210..008e0ec3dc 100644 --- a/modules/hydrodyn/src/Morison.f90 +++ b/modules/hydrodyn/src/Morison.f90 @@ -3480,7 +3480,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, rn, rn1, rn2, z_hi, zFillGroup + REAL(ReKi) :: l, z1, z2, zMid, r1, r2, r1b, r2b, r1In, r2In, rMidIn, 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 @@ -3777,16 +3777,10 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, ELSE IF (mem%MHstLMod == 2) THEN ! Alternative hydrostatic load calculation ! Get free surface elevation and normal at the element midpoint (both assumed constant over the element) posMid = 0.5 * (pos1+pos2) - ! rn is only used to estimate free surface normal numerically - IF (mem%MSecGeom == MSecGeom_Cyl) THEN - rn = 0.5 * (r1b +r2b ) - ELSE IF (mem%MSecGeom == MSecGeom_Rec) THEN - rn = MAX( 0.5*(Sa1b+Sa2b), 0.5*(Sb1b+Sb2b) ) - END IF IF (p%WaveField%WaveStMod > 0) THEN CALL GetTotalWaveElev( Time, posMid, ZetaMid, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, posMid, rn, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( Time, posMid, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/posMid(1),posMid(2),ZetaMid/) ! Reference point on the free surface ELSE @@ -4561,22 +4555,18 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, if (mem%MSecGeom==MSecGeom_Cyl) then r1 = mem%RMGB( 1) r2 = mem%RMGB(N+1) - rn1 = r1 - rn2 = r2 else if (mem%MSecGeom==MSecGeom_Rec) then Sa1 = mem%SaMGB( 1) Sa2 = mem%SaMGB(N+1) Sb1 = mem%SbMGB( 1) Sb2 = mem%SbMGB(N+1) - rn1 = MAX(Sa1,Sb1) - rn2 = MAX(Sa2,Sb2) 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%WaveField%WaveStMod > 0) THEN CALL GetTotalWaveElev( Time, pos1, Zeta1, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos1, rn1, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( Time, pos1, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos1(1),pos1(2),Zeta1/) ! Reference point on the free surface ELSE @@ -4604,7 +4594,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF (p%WaveField%WaveStMod > 0) THEN CALL GetTotalWaveElev( Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( Time, pos2, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4633,7 +4623,7 @@ SUBROUTINE Morison_CalcOutput( Time, u, p, x, xd, z, OtherState, y, m, errStat, IF (p%WaveField%WaveStMod > 0) THEN CALL GetTotalWaveElev( Time, pos2, Zeta2, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL GetFreeSurfaceNormal( Time, pos2, rn2, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( Time, pos2, n_hat, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) FSPt = (/pos2(1),pos2(2),Zeta2/) ! Reference point on the free surface ELSE @@ -4767,10 +4757,9 @@ SUBROUTINE GetTotalWaveElev( Time, pos, Zeta, ErrStat, ErrMsg ) END SUBROUTINE GetTotalWaveElev - SUBROUTINE GetFreeSurfaceNormal( Time, pos, r, n, ErrStat, ErrMsg) + SUBROUTINE GetFreeSurfaceNormal( Time, pos, n, ErrStat, ErrMsg) 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 ! Distance for central differencing 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 @@ -4780,7 +4769,7 @@ SUBROUTINE GetFreeSurfaceNormal( Time, pos, r, n, ErrStat, ErrMsg) ErrStat = ErrID_None ErrMsg = "" - CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, r, n, ErrStat2, ErrMsg2 ) + CALL WaveField_GetNodeWaveNormal( p%WaveField, m%WaveField_m, Time, pos, n, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) END SUBROUTINE GetFreeSurfaceNormal @@ -5714,7 +5703,7 @@ SUBROUTINE getElementHstLds_Mod1( Time, pos1, pos2, Zeta1, Zeta2, k_hat, r1, r2, rh = r1 + h0*dRdl ! Estimate the free-surface normal at the free-surface intersection, n_hat IF ( p%WaveField%WaveStMod > 0_IntKi ) THEN ! If wave stretching is enabled, compute free surface normal - CALL GetFreeSurfaceNormal( Time, FSInt, rh, n_hat, ErrStat2, ErrMsg2 ) + CALL GetFreeSurfaceNormal( Time, FSInt, 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/) diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 6c39bf8535..970a7c21f3 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -109,12 +109,11 @@ logical function Failed() END FUNCTION WaveField_GetNodeTotalWaveElev -SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, r, n, ErrStat, ErrMsg ) +SUBROUTINE WaveField_GetNodeWaveNormal( WaveField, WaveField_m, Time, pos, n, ErrStat, ErrMsg ) 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(in ) :: r ! Distance for central differencing 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 From 90ae7c1c540546c6f1e48055679b164995832717 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Mon, 29 Sep 2025 17:12:45 -0600 Subject: [PATCH 6/9] Update r-test pointer --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index c8f61832d2..bc4986fab9 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit c8f61832d2854eac9ad0eeec3a0e2c2f4938e13f +Subproject commit bc4986fab93d00a8cc497039e8fc1360253d90d4 From 31474e24d92c2f8925cb64709205891bb126a488 Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Mon, 29 Sep 2025 17:50:17 -0600 Subject: [PATCH 7/9] Update r-test pointer --- reg_tests/r-test | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reg_tests/r-test b/reg_tests/r-test index bc4986fab9..966e0a432b 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit bc4986fab93d00a8cc497039e8fc1360253d90d4 +Subproject commit 966e0a432b1f1acf039065e2b6c8282829bab09d From 383159911fde6a6c27006f2675964c1311d770d0 Mon Sep 17 00:00:00 2001 From: andrew-platt Date: Mon, 24 Nov 2025 17:31:44 -0700 Subject: [PATCH 8/9] GridInterp: update VS projects - FASTlib.vfproj - HydroDynDriver.vfproj - HydroDyn_c_binding.vfproj - MoorDynDriver.vfproj - MoorDyn_c_binding.vfproj - RunRegistry.bat - SeaStateDriver.vfproj - SeaState_c_binding.vfproj --- vs-build/FASTlib/FASTlib.vfproj | 29 ++++ vs-build/HydroDyn/HydroDynDriver.vfproj | 23 ++- .../HydroDyn_c_binding.vfproj | 13 +- vs-build/MoorDyn/MoorDynDriver.vfproj | 23 ++- .../MoorDyn_c_binding.vfproj | 146 +++--------------- vs-build/RunRegistry.bat | 6 + vs-build/SeaState/SeaStateDriver.vfproj | 21 +++ .../SeaState_c_binding.vfproj | 37 ++++- 8 files changed, 169 insertions(+), 129 deletions(-) diff --git a/vs-build/FASTlib/FASTlib.vfproj b/vs-build/FASTlib/FASTlib.vfproj index e34ecf918e..f4cacd3b7c 100644 --- a/vs-build/FASTlib/FASTlib.vfproj +++ b/vs-build/FASTlib/FASTlib.vfproj @@ -1943,6 +1943,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1968,6 +1996,7 @@ + diff --git a/vs-build/HydroDyn/HydroDynDriver.vfproj b/vs-build/HydroDyn/HydroDynDriver.vfproj index 6b4bd4fe78..1b60204924 100644 --- a/vs-build/HydroDyn/HydroDynDriver.vfproj +++ b/vs-build/HydroDyn/HydroDynDriver.vfproj @@ -290,6 +290,26 @@ + + + + + + + + + + + + + + + + + + + + @@ -308,7 +328,7 @@ - + @@ -317,6 +337,7 @@ + diff --git a/vs-build/HydroDyn_c_binding/HydroDyn_c_binding.vfproj b/vs-build/HydroDyn_c_binding/HydroDyn_c_binding.vfproj index bc9116ccf7..801edf151e 100644 --- a/vs-build/HydroDyn_c_binding/HydroDyn_c_binding.vfproj +++ b/vs-build/HydroDyn_c_binding/HydroDyn_c_binding.vfproj @@ -224,7 +224,6 @@ - @@ -241,10 +240,19 @@ + + + + + + + + + - + @@ -253,6 +261,7 @@ + diff --git a/vs-build/MoorDyn/MoorDynDriver.vfproj b/vs-build/MoorDyn/MoorDynDriver.vfproj index 9f7ac9fdfc..9d936b59cc 100644 --- a/vs-build/MoorDyn/MoorDynDriver.vfproj +++ b/vs-build/MoorDyn/MoorDynDriver.vfproj @@ -109,6 +109,26 @@ + + + + + + + + + + + + + + + + + + + + @@ -123,13 +143,14 @@ - + + diff --git a/vs-build/MoorDyn_c_binding/MoorDyn_c_binding.vfproj b/vs-build/MoorDyn_c_binding/MoorDyn_c_binding.vfproj index 45710d681e..a9e9da955d 100644 --- a/vs-build/MoorDyn_c_binding/MoorDyn_c_binding.vfproj +++ b/vs-build/MoorDyn_c_binding/MoorDyn_c_binding.vfproj @@ -129,7 +129,7 @@ - + @@ -142,30 +142,14 @@ - - + - + - + - - - - - - - - - - - - - - - - - + + @@ -178,112 +162,28 @@ - - - - - + + - + - + - - - + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - + + diff --git a/vs-build/RunRegistry.bat b/vs-build/RunRegistry.bat index bce124d4ec..76558bfe37 100644 --- a/vs-build/RunRegistry.bat +++ b/vs-build/RunRegistry.bat @@ -71,6 +71,12 @@ SET Output_Loc=%CURR_LOC% %REGISTRY% "%CURR_LOC%\Registry_NWTC_Library_base.txt" -I "%NWTC_Lib_Loc%" -I "%CURR_LOC%" -O "%Output_Loc%" -noextrap GOTO checkError +:GridInterp +SET CURR_LOC=%NWTC_Lib_Loc% +SET Output_Loc=%CURR_LOC% +%REGISTRY% "%CURR_LOC%\GridInterp.txt " -I "%NWTC_Lib_Loc%" -I "%CURR_LOC%" -O "%Output_Loc%" -noextrap +GOTO checkError + :MAP SET CURR_LOC=%MAP_Loc% SET Output_Loc=%CURR_LOC% diff --git a/vs-build/SeaState/SeaStateDriver.vfproj b/vs-build/SeaState/SeaStateDriver.vfproj index f0f29c7847..f7825846a6 100644 --- a/vs-build/SeaState/SeaStateDriver.vfproj +++ b/vs-build/SeaState/SeaStateDriver.vfproj @@ -86,6 +86,26 @@ + + + + + + + + + + + + + + + + + + + + @@ -107,6 +127,7 @@ + diff --git a/vs-build/SeaState_c_binding/SeaState_c_binding.vfproj b/vs-build/SeaState_c_binding/SeaState_c_binding.vfproj index 7ea7b9851b..817bff1514 100644 --- a/vs-build/SeaState_c_binding/SeaState_c_binding.vfproj +++ b/vs-build/SeaState_c_binding/SeaState_c_binding.vfproj @@ -88,7 +88,6 @@ - @@ -105,8 +104,41 @@ + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -117,6 +149,7 @@ + From 7ce01497b6cac8b6513e3ec1b467d4f753a15b3a Mon Sep 17 00:00:00 2001 From: Lu Wang Date: Wed, 26 Nov 2025 14:17:53 -0700 Subject: [PATCH 9/9] Update GridInterp based on review comments --- modules/nwtc-library/src/GridInterp.f90 | 240 ++++++++++-------------- 1 file changed, 100 insertions(+), 140 deletions(-) diff --git a/modules/nwtc-library/src/GridInterp.f90 b/modules/nwtc-library/src/GridInterp.f90 index 6df26b738d..3398ca5269 100644 --- a/modules/nwtc-library/src/GridInterp.f90 +++ b/modules/nwtc-library/src/GridInterp.f90 @@ -112,10 +112,10 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn if ( nMax .EQ. 1_IntKi ) then ! Only one grid point, effectively ignore this dimension ! Construct a dummy linear interpolation for now - Indx(1) = -1_IntKi + Indx(1) = 0_IntKi Indx(2) = 0_IntKi Indx(3) = 0_IntKi - Indx(4) = -1_IntKi + Indx(4) = 0_IntKi isopc = 0.5_ReKi Support = 0 return @@ -130,7 +130,7 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn isopc = p - floor(p,ReKi) ! Get the normalized coordinates of the two nearest nodes to the left and to the right - Indx = floor( p, IntKi ) + (/-1,0,1,2/) + Indx = floor( p, IntKi ) + [-1,0,1,2] ! Make sure the coordinates are not out of bound using periodicity do i = 1,4 @@ -140,6 +140,9 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn Indx(i) = mod(Indx(i),nMax) end do + ! Always use cubic interpolation for periodic dimensions + support = 3 ! Cubic interpolation + else pMax = Real(nMax-1_IntKi,ReKi) @@ -168,27 +171,30 @@ Subroutine SetIndex(pIn,pZero,delta,nMax,IsPeriodic,Indx,isopc,Support,FirstWarn ! Calculate normalized coordinate between 0 and 1 isopc = 1.0_ReKi ! Get the normalized coordinates of the two nearest nodes to the left and to the right - Indx = nMax + (/-3,-2,-1,0/) + Indx = nMax + [-3,-2,-1,0] else ! Calculate normalized coordinate between 0 and 1 isopc = p - floor(p,ReKi) ! Get the normalized coordinates of the two nearest nodes to the left and to the right - Indx = floor( p, IntKi ) + (/-1,0,1,2/) + Indx = floor( p, IntKi ) + [-1,0,1,2] end if - end if - - ! Supported interpolation method - if ( Indx(1) < 0_IntKi ) then - if (Indx(4) > (nMax-1_IntKi) ) then - support = 0 ! Linear interpolation + ! Supported interpolation method + if ( Indx(1) < 0_IntKi ) then + Indx(1) = 0_IntKi + if (Indx(4) > (nMax-1_IntKi) ) then + support = 0 ! Linear interpolation + Indx(4) = nMax-1_IntKi + else + support = 1 ! Quadratic interpolation with only one node to the left + end if + else if ( Indx(4) > (nMax-1_IntKi) ) then + support = 2 ! Quadratic interpolation with only one node to the right + Indx(4) = nMax-1_IntKi else - support = 1 ! Quadratic interpolation with only one node to the left + support = 3 ! Cubic interpolation end if - else if ( Indx(4) > (nMax-1_IntKi) ) then - support = 2 ! Quadratic interpolation with only one node to the right - else - support = 3 ! Cubic interpolation + end if End Subroutine SetIndex @@ -200,42 +206,43 @@ Subroutine GetN1D(isopc, support, N1D) real(ReKi), intent(inout) :: N1D(4) real(ReKi) :: isopc2,isopc3 - if ( Support == 3 ) then ! Cubic interpolation + select case ( Support ) + case ( 3 ) ! Cubic interpolation isopc2 = isopc*isopc isopc3 = isopc*isopc2 - N1D(1) = -0.5*isopc3 + isopc2 - 0.5*isopc - N1D(2) = 1.5*isopc3 - 2.5*isopc2 + 1.0 - N1D(3) = -1.5*isopc3 + 2.0*isopc2 + 0.5*isopc - N1D(4) = 0.5*isopc3 - 0.5*isopc2 + N1D(1) = -0.5_ReKi*isopc3 + isopc2 - 0.5_ReKi*isopc + N1D(2) = 1.5_ReKi*isopc3 - 2.5_ReKi*isopc2 + 1.0_ReKi + N1D(3) = -1.5_ReKi*isopc3 + 2.0_ReKi*isopc2 + 0.5_ReKi*isopc + N1D(4) = 0.5_ReKi*isopc3 - 0.5_ReKi*isopc2 - else if ( Support == 1 ) then ! Quadratic interpolation with only one node to the left + case ( 1 ) ! Quadratic interpolation with only one node to the left isopc2 = isopc*isopc - N1D(1) = 0.0 - N1D(2) = 0.5*isopc2 - 1.5*isopc + 1.0 - N1D(3) = -1.0*isopc2 + 2.0*isopc - N1D(4) = 0.5*isopc2 - 0.5*isopc + N1D(1) = 0.0_ReKi + N1D(2) = 0.5_ReKi*isopc2 - 1.5_ReKi*isopc + 1.0_ReKi + N1D(3) = -1.0_ReKi*isopc2 + 2.0_ReKi*isopc + N1D(4) = 0.5_ReKi*isopc2 - 0.5_ReKi*isopc - else if ( Support == 2 ) then ! Quadratic interpolation with only one node to the right + case ( 2 ) ! Quadratic interpolation with only one node to the right isopc2 = isopc*isopc - N1D(1) = 0.5*isopc2 - 0.5*isopc - N1D(2) = -1.0*isopc2 + 1.0 - N1D(3) = 0.5*isopc2 + 0.5*isopc - N1D(4) = 0.0 + N1D(1) = 0.5_ReKi*isopc2 - 0.5_ReKi*isopc + N1D(2) = -1.0_ReKi*isopc2 + 1.0_ReKi + N1D(3) = 0.5_ReKi*isopc2 + 0.5_ReKi*isopc + N1D(4) = 0.0_ReKi - else ! if ( Support == 0 ) then ! Linear interpolation + case default ! Support == 0 Linear interpolation - N1D(1) = 0.0 - N1D(2) = -isopc + 1.0 + N1D(1) = 0.0_ReKi + N1D(2) = -isopc + 1.0_ReKi N1D(3) = isopc - N1D(4) = 0.0 + N1D(4) = 0.0_ReKi - end if + end select End Subroutine GetN1D @@ -246,37 +253,38 @@ Subroutine GetN1Ddx(isopc, support, N1D) real(ReKi), intent(inout) :: N1D(4) real(ReKi) :: isopc2 - if ( Support == 3 ) then ! Cubic interpolation + select case ( Support ) + case ( 3 ) ! Cubic interpolation isopc2 = isopc*isopc - N1D(1) = -1.5*isopc2 + 2.0*isopc - 0.5 - N1D(2) = 4.5*isopc2 - 5.0*isopc - N1D(3) = -4.5*isopc2 + 4.0*isopc + 0.5 - N1D(4) = 1.5*isopc2 - isopc + N1D(1) = -1.5_ReKi*isopc2 + 2.0_ReKi*isopc - 0.5_ReKi + N1D(2) = 4.5_ReKi*isopc2 - 5.0_ReKi*isopc + N1D(3) = -4.5_ReKi*isopc2 + 4.0_ReKi*isopc + 0.5_ReKi + N1D(4) = 1.5_ReKi*isopc2 - isopc - else if ( Support == 1 ) then ! Quadratic interpolation with only one node to the left + case ( 1 ) ! Quadratic interpolation with only one node to the left - N1D(1) = 0.0 - N1D(2) = isopc - 1.5 - N1D(3) = -2.0*isopc + 2.0 - N1D(4) = isopc - 0.5 + N1D(1) = 0.0_ReKi + N1D(2) = isopc - 1.5_ReKi + N1D(3) = -2.0_ReKi*isopc + 2.0_ReKi + N1D(4) = isopc - 0.5_ReKi - else if ( Support == 2 ) then ! Quadratic interpolation with only one node to the right + case ( 2 ) ! Quadratic interpolation with only one node to the right - N1D(1) = isopc - 0.5 - N1D(2) = -2.0*isopc - N1D(3) = isopc + 0.5 - N1D(4) = 0.0 + N1D(1) = isopc - 0.5_ReKi + N1D(2) = -2.0_ReKi*isopc + N1D(3) = isopc + 0.5_ReKi + N1D(4) = 0.0_ReKi - else ! if ( Support == 0 ) then ! Linear interpolation + case default ! Support == 0 Linear interpolation - N1D(1) = 0.0 - N1D(2) = -1.0 - N1D(3) = 1.0 - N1D(4) = 0.0 + N1D(1) = 0.0_ReKi + N1D(2) = -1.0_ReKi + N1D(3) = 1.0_ReKi + N1D(4) = 0.0_ReKi - end if + end select End Subroutine GetN1Ddx @@ -430,10 +438,7 @@ function GridInterp3DR4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - GridInterp3DR4 = GridInterp3DR4 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - ! else - Indx out of bound - end if + GridInterp3DR4 = GridInterp3DR4 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do @@ -453,10 +458,7 @@ function GridInterp3DR8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - GridInterp3DR8 = GridInterp3DR8 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - ! else - Indx out of bound - end if + GridInterp3DR8 = GridInterp3DR8 + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do @@ -482,12 +484,9 @@ function GridInterp3DVecR4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - do vi = 1,vDim - GridInterp3DVecR4(vi) = GridInterp3DVecR4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp3DVecR4(vi) = GridInterp3DVecR4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do end do end do end do @@ -508,12 +507,9 @@ function GridInterp3DVecR8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - do vi = 1,vDim - GridInterp3DVecR8(vi) = GridInterp3DVecR8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp3DVecR8(vi) = GridInterp3DVecR8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do end do end do end do @@ -539,12 +535,9 @@ function GridInterp3DVec6R4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - do vi = 1,vDim - GridInterp3DVec6R4(vi) = GridInterp3DVec6R4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp3DVec6R4(vi) = GridInterp3DVec6R4(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do end do end do end do @@ -565,12 +558,9 @@ function GridInterp3DVec6R8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N3D(i,j,k),0.0_ReKi)) then - do vi = 1,vDim - GridInterp3DVec6R8(vi) = GridInterp3DVec6R8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp3DVec6R8(vi) = GridInterp3DVec6R8(vi) + m%N3D(i,j,k) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), vi ) + end do end do end do end do @@ -596,10 +586,7 @@ function GridInterp4DR4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - GridInterp4DR4 = GridInterp4DR4 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) - ! else - Indx out of bound - end if + GridInterp4DR4 = GridInterp4DR4 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) end do end do end do @@ -621,10 +608,7 @@ function GridInterp4DR8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - GridInterp4DR8 = GridInterp4DR8 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) - ! else - Indx out of bound - end if + GridInterp4DR8 = GridInterp4DR8 + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4) ) end do end do end do @@ -652,12 +636,9 @@ function GridInterp4DVecR4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - do vi = 1,vDim - GridInterp4DVecR4(vi) = GridInterp4DVecR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp4DVecR4(vi) = GridInterp4DVecR4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do end do end do end do @@ -680,12 +661,9 @@ function GridInterp4DVecR8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - do vi = 1,vDim - GridInterp4DVecR8(vi) = GridInterp4DVecR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp4DVecR8(vi) = GridInterp4DVecR8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do end do end do end do @@ -713,12 +691,9 @@ function GridInterp4DVec6R4( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - do vi = 1,vDim - GridInterp4DVec6R4(vi) = GridInterp4DVec6R4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp4DVec6R4(vi) = GridInterp4DVec6R4(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do end do end do end do @@ -741,12 +716,9 @@ function GridInterp4DVec6R8( data, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,l),0.0_ReKi)) then - do vi = 1,vDim - GridInterp4DVec6R8(vi) = GridInterp4DVec6R8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) - end do - ! else - Indx out of bound - end if + do vi = 1,vDim + GridInterp4DVec6R8(vi) = GridInterp4DVec6R8(vi) + m%N4D(i,j,k,l) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3), m%Indx(l,4), vi ) + end do end do end do end do @@ -775,19 +747,15 @@ function GridInterpNR4( data, p, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,1),0.0_ReKi)) then - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if - if (.not.EqualRealNos(m%N4D(i,j,k,2),0.0_ReKi)) then - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if + dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do dZetadx = dZetadx / p%delta(2) dZetady = dZetady / p%delta(3) - GridInterpNR4 = (/-dZetadx,-dZetady,1.0_SiKi/) + GridInterpNR4 = [-dZetadx,-dZetady,1.0_SiKi] GridInterpNR4 = GridInterpNR4 / TwoNorm(GridInterpNR4) end function GridInterpNR4 @@ -808,12 +776,8 @@ function GridInterpNR8( data, p, m ) do k = 1,4 do j = 1,4 do i = 1,4 - if (.not.EqualRealNos(m%N4D(i,j,k,1),0.0_ReKi)) then - dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if - if (.not.EqualRealNos(m%N4D(i,j,k,2),0.0_ReKi)) then - dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if + dZetadx = dZetadx + m%N4D(i,j,k,1) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) + dZetady = dZetady + m%N4D(i,j,k,2) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do @@ -845,9 +809,7 @@ function GridInterpSR4( data, p, m ) do j = 1,4 do i = 1,4 do dir = 1,2 - if (.not.EqualRealNos(m%N4D(i,j,k,dir),0.0_ReKi)) then - GridInterpSR4(dir) = GridInterpSR4(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if + GridInterpSR4(dir) = GridInterpSR4(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do @@ -871,9 +833,7 @@ function GridInterpSR8( data, p, m ) do j = 1,4 do i = 1,4 do dir = 1,2 - if (.not.EqualRealNos(m%N4D(i,j,k,dir),0.0_ReKi)) then - GridInterpSR8(dir) = GridInterpSR8(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) - end if + GridInterpSR8(dir) = GridInterpSR8(dir) + m%N4D(i,j,k,dir) * data( m%Indx(i,1), m%Indx(j,2), m%Indx(k,3) ) end do end do end do