diff --git a/modules/seastate/src/SeaSt_WaveField.f90 b/modules/seastate/src/SeaSt_WaveField.f90 index 04d66270d..392576750 100644 --- a/modules/seastate/src/SeaSt_WaveField.f90 +++ b/modules/seastate/src/SeaSt_WaveField.f90 @@ -17,6 +17,7 @@ MODULE SeaSt_WaveField PUBLIC WaveField_GetNodeWaveKin PUBLIC WaveField_GetNodeWaveVelAcc PUBLIC WaveField_GetWaveKin +PUBLIC WaveField_GetDynP CONTAINS @@ -310,6 +311,103 @@ logical function Failed() END SUBROUTINE WaveField_GetNodeWaveKin +!-------------------- Subroutine for dynamic pressure --------------------! +!NOTE: this is a subset of WaveField_GetNodeWaveKin +SUBROUTINE WaveField_GetDynP( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FDynP, 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(3) + logical, intent(in ) :: forceNodeInWater + real(SiKi), intent( out) :: FDynP + integer(IntKi), intent( out) :: nodeInWater + integer(IntKi), intent( out) :: ErrStat ! Error status of the operation + character(*), intent( out) :: ErrMsg ! Error message if errStat /= ErrID_None + + real(ReKi) :: posXY(2), posPrime(3), posXY0(3) + character(*), parameter :: RoutineName = 'WaveField_GetDynP' + integer(IntKi) :: errStat2 + character(ErrMsgLen) :: errMsg2 + + ! Temporary vars not kept + real(SiKi) :: WaveElev1 + real(SiKi) :: WaveElev2 + real(SiKi) :: WaveElev + + ErrStat = ErrID_None + ErrMsg = "" + + posXY = pos(1:2) + posXY0 = (/pos(1),pos(2),0.0_ReKi/) + + ! Wave elevation + WaveElev1 = WaveField_GetNodeWaveElev1( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return; + WaveElev2 = WaveField_GetNodeWaveElev2( WaveField, WaveField_m, Time, pos, ErrStat2, ErrMsg2 ); if (Failed()) return; + WaveElev = WaveElev1 + WaveElev2 + + IF (WaveField%WaveStMod == 0) THEN ! No wave stretching + + 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+WaveField%WaveTimeShift, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) + ELSE ! Node is above the SWL + nodeInWater = 0_IntKi + FDynP = 0.0 + END IF + + ELSE ! Wave stretching enabled + IF ( (pos(3) <= WaveElev) .OR. forceNodeInWater ) THEN ! Node is submerged + nodeInWater = 1_IntKi + IF ( WaveField%WaveStMod < 3 ) THEN ! Vertical or extrapolated wave stretching + + 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+WaveField%WaveTimeShift, pos, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) + ELSE ! Node is above SWL - need wave stretching + + ! Vertical wave stretching + CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posXY0, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) + + ! Extrapoled wave stretching + IF (WaveField%WaveStMod == 2) THEN + CALL WaveField_Interp_Setup3D( Time+WaveField%WaveTimeShift, posXY, WaveField%SrfGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FDynP = FDynP + GridInterp3D ( WaveField%PWaveDynP0, WaveField_m ) * pos(3) + END IF + + END IF ! Node is submerged + + ELSE ! Wheeler stretching - no need to check whether the node is above or below SWL + + ! Map the node z-position linearly from [-EffWtrDpth,m%WaveElev(j)] to [-EffWtrDpth,0] + posPrime = pos + posPrime(3) = WaveField%EffWtrDpth*(WaveField%EffWtrDpth+pos(3))/(WaveField%EffWtrDpth+WaveElev)-WaveField%EffWtrDpth + posPrime(3) = MIN( posPrime(3), 0.0_ReKi) ! Clamp z-position to zero. Needed when forceNodeInWater=.TRUE. + + ! Obtain the wave-field variables by interpolation with the mapped position. + CALL WaveField_Interp_Setup4D( Time+WaveField%WaveTimeShift, posPrime, WaveField%GridDepth, WaveField%VolGridParams, WaveField_m, ErrStat2, ErrMsg2 ); if (Failed()) return; + FDynP = GridInterp4D ( WaveField%WaveDynP, WaveField_m ) + END IF + + ELSE ! Node is out of water - zero-out all wave dynamics + + nodeInWater = 0_IntKi + FDynP = 0.0 + + END IF ! If node is in or out of water + END IF ! If wave stretching is on or off + +contains + logical function Failed() + call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + Failed = ErrStat >= AbortErrLev + end function +END SUBROUTINE WaveField_GetDynP + + !-------------------- Subroutine for wave field velocity only --------------------! SUBROUTINE WaveField_GetNodeWaveVelAcc( WaveField, WaveField_m, Time, pos, forceNodeInWater, nodeInWater, FV, FA, ErrStat, ErrMsg ) type(SeaSt_WaveFieldType), intent(in ) :: WaveField diff --git a/modules/seastate/src/SeaState_C_Binding.f90 b/modules/seastate/src/SeaState_C_Binding.f90 index 1a7a6d81d..747b786d7 100644 --- a/modules/seastate/src/SeaState_C_Binding.f90 +++ b/modules/seastate/src/SeaState_C_Binding.f90 @@ -44,6 +44,7 @@ MODULE SeaState_C_Binding PUBLIC :: SeaSt_C_GetDens PUBLIC :: SeaSt_C_GetDpth PUBLIC :: SeaSt_C_GetMSL2SWL + PUBLIC :: SeaSt_C_GetDynPressure !------------------------------------------------------------------------------------ ! Debugging: DebugLevel -- passed at PreInit @@ -898,6 +899,73 @@ subroutine CheckWaveFieldPtr(callingRoutine,valid,ErrStat,ErrMsg) end subroutine +!> return the surface elevation and normal vector at a point. +subroutine SeaSt_C_GetDynPressure(Time_C, Pos_C, DynP_C, ErrStat_C,ErrMsg_C) BIND (C, NAME='SeaSt_C_GetDynPressure') +#ifndef IMPLICIT_DLLEXPORT +!DEC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetDynPressure +!GCC$ ATTRIBUTES DLLEXPORT :: SeaSt_C_GetDynPressure +#endif + real(c_double), intent(in ) :: Time_C + real(c_float), intent(in ) :: Pos_c(3) + real(c_float), intent( out) :: DynP_C + integer(c_int), intent( out) :: ErrStat_C + character(kind=c_char), intent( out) :: ErrMsg_C(ErrMsgLen_C) + + real(DbKi) :: Time + real(ReKi) :: Pos(3) + real(SiKi) :: FDynP + integer :: ErrStat + character(ErrMsgLen) :: ErrMsg + character(*), parameter :: RoutineName = 'SeaSt_C_GetDynPressure' + logical :: valid + ! Temporary vars for calling GetNodeWaveKin -- remove if new interface developed + integer(IntKi) :: nodeInWater + + if (DebugLevel > 1) call ShowPassedData() + + ! convert position and time to fortran types + Time = real(Time_C, DbKi) + Pos = real(Pos_C(1:3), ReKi) + + ! verify there is actually wavefield data + call CheckWaveFieldPtr(RoutineName, valid, ErrStat, ErrMsg) + if (.not. valid) then + call Cleanup() + return + endif + + ! get wave elevation (total combined first and second order) + ! Notes: + ! - if position is outside the wave field boundary, it will simply return boundary edge value + ! - time must be positive or a fatal error occurs + ! - forceNodeInWater = .false. -- perhaps update this later to allow throug interface? + call WaveField_GetDynP( p%WaveField, m%WaveField_m, Time, pos, .false., nodeInWater, FDynP, ErrStat, ErrMsg ) + + + ! Store resulting elevation as C type + DynP_C = real(FDynP,c_float) + + if (DebugLevel > 1) call ShowReturnData() + call Cleanup() + return +contains + subroutine Cleanup() ! NOTE: we are ignoring any error reporting from here + CALL SetErrStat_F2C(ErrStat,ErrMsg,ErrStat_C,ErrMsg_C) + end subroutine Cleanup + subroutine ShowPassedData() + call WrScr("-----------------------------------------------------------") + call WrScr("Interface debugging: SeaSt_C_GetDynPressure") + call WrScr(" --------------------------------------------------------") + call WrScr(" Time_C -> "//trim(Num2LStr(Time_C))) + call WrScr(" Pos_C -> ("//trim(Num2LStr(Pos_C(1)))//","//trim(Num2LStr(Pos_C(2)))//","//trim(Num2LStr(Pos_C(3)))//")") + end subroutine ShowPassedData + subroutine ShowReturnData() + call WrScr(" DynP_C <- "//trim(Num2LStr(DynP_C))) + call WrScr("-----------------------------------------------------------") + end subroutine ShowReturnData +end subroutine SeaSt_C_GetDynPressure + + !FIXME: the following visualization writer should be merged into the library vtk.f90 ! this is a modified duplicate of the routine from FAST_Subs by the same name. !FIXME: this routine currently only grabs the closest timestep and does not do any interpolation