From fd8e84eaf267eef93717d2c120be101f1de2daca Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Wed, 26 Jul 2017 18:49:03 -0600 Subject: [PATCH 1/5] bjonkman/f/Envision-BeamDyn: the BD input file specified in the driver was defined relative to the working directory instead of the directory where the driver input file is located --- modules-local/beamdyn/src/Driver_Beam_Subs.f90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 index 9363b6326e..c937990f48 100644 --- a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 +++ b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 @@ -83,6 +83,7 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& INTEGER(IntKi) :: UnEc CHARACTER(1024) :: FTitle ! "File Title": the 2nd line of the input file, which contains a description of its contents + CHARACTER(1024) :: PriPath ! Path name of the primary file INTEGER(IntKi) :: i !------------------------------------------------------------------------------------ @@ -98,6 +99,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) IF (ErrStat >= AbortErrLev) RETURN + CALL GetPath( DvrInputFile, PriPath ) ! Input files will be relative to the path where the primary input file is located. + !-------------------------- HEADER --------------------------------------------- CALL ReadCom(UnIn,DvrInputFile,'File Header: Module Version (line 1)',ErrStat2,ErrMsg2,UnEc) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) @@ -270,6 +273,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,& !---------------------- BEAM SECTIONAL PARAMETER ---------------------------------------- CALL ReadVar ( UnIn, DvrInputFile, InitInputData%InputFile, 'InputFile', 'Name of the primary input file', ErrStat2,ErrMsg2, UnEc ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF ( PathIsRelative( InitInputData%InputFile ) ) InitInputData%InputFile = TRIM(PriPath)//TRIM(InitInputData%InputFile) + if (ErrStat >= AbortErrLev) then call cleanup() return From 062941b0ae1eb60b6e9b69dcc5f0ec5fd9f1126a Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Wed, 26 Jul 2017 18:50:31 -0600 Subject: [PATCH 2/5] bjonkman/f/Envision-BeamDyn: rotational velocity vector is now double precision --- .../beamdyn/src/Driver_Beam_Subs.f90 | 38 +++++++++---------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 index c937990f48..7229e56caa 100644 --- a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 +++ b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 @@ -27,26 +27,22 @@ module BeamDyn_driver_subs TYPE , PUBLIC :: BD_DriverInternalType REAL(ReKi) , DIMENSION(1:6) :: DistrLoad !< Constant distributed load along beam axis, 3 forces and 3 moments [-] REAL(ReKi) , DIMENSION(1:6) :: TipLoad !< Constant point load applied at tip, 3 forces and 3 moments [-] - INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments [-] - REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis (index 1= Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-] -! INTEGER(IntKi) :: NumPointLoads !< Number of multi-point loads in the driver input file -! REAL(DbKi) :: MultiPointLoads{:}{:} !< The array of multipoint loads Index 1: [1, NumPointLoads]; - !< Index 2: [1,7] (index of Loads 1 = Relative position along blade span; - !! 2-7 = Fx, Fy, Fz, Mx, My, Mz ) - - TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations - TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads - TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion - TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad - TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh - TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point - - TYPE(MeshType) :: RotationCenter - TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion + INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments, from the driver input file [-] + REAL(BDKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis, size (NumPointLoads,7); (dimension 2: index 1=Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-] + + TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations + TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads + TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion + TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad + TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh + TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point + + TYPE(MeshType) :: RotationCenter + TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion - REAL(DbKi) :: t_initial - REAL(DbKi) :: t_final - REAL(ReKi) :: w ! magnitude of rotational velocity vector + REAL(DbKi) :: t_initial + REAL(DbKi) :: t_final + REAL(R8Ki) :: w ! magnitude of rotational velocity vector END TYPE @@ -591,8 +587,8 @@ SUBROUTINE Init_RotationCenterMesh(DvrData, InitInputData, RootMotionMesh, ErrSt DvrData%w = TwoNorm( InitInputData%RootVel(4:6) ) - if (EqualRealNos(DvrData%w,0.0_ReKi)) then - DvrData%w = 0.0_ReKi + if (EqualRealNos(DvrData%w,0.0_R8Ki)) then + DvrData%w = 0.0_R8Ki ! the beam is not rotating, so pick an orientation call eye(orientation, ErrStat2, ErrMsg2) else From d76d4b49621b8ea418e460adf3deff4049d725c0 Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Thu, 27 Jul 2017 14:27:01 -0600 Subject: [PATCH 3/5] use R8Ki in BeamDyn_Driver for support with single precision This forces the use of double over float when compiling in single precision --- modules-local/beamdyn/src/Driver_Beam_Subs.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 index 7229e56caa..19927d800f 100644 --- a/modules-local/beamdyn/src/Driver_Beam_Subs.f90 +++ b/modules-local/beamdyn/src/Driver_Beam_Subs.f90 @@ -594,13 +594,13 @@ SUBROUTINE Init_RotationCenterMesh(DvrData, InitInputData, RootMotionMesh, ErrSt else z_hat = InitInputData%RootVel(4:6) / DvrData%w - if ( EqualRealNos( z_hat(3), 1.0_ReKi ) ) then + if ( EqualRealNos( z_hat(3), 1.0_R8Ki ) ) then call eye(orientation, ErrStat2, ErrMsg2) - elseif ( EqualRealNos( z_hat(3), -1.0_ReKi ) ) then + elseif ( EqualRealNos( z_hat(3), -1.0_R8Ki ) ) then orientation = 0.0_ReKi - orientation(1,1) = -1.0_ReKi - orientation(2,2) = 1.0_ReKi - orientation(3,3) = -1.0_ReKi + orientation(1,1) = -1.0_R8Ki + orientation(2,2) = 1.0_R8Ki + orientation(3,3) = -1.0_R8Ki else Z_unit = (/0, 0, 1/) From 6fc9086ae15f91d0f2f21cc4ca96c03d6c83d1ae Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Thu, 27 Jul 2017 14:28:40 -0600 Subject: [PATCH 4/5] bjonkman/f/Envision-BeamDyn: BD: updated BD_CrvExtractCrv for more robust arithmetic 4 sets of equations should give same results, assuming the "pivot" is not zero. New implementation picks set of equations with largest pivot value, not just one that is slightly larger than 0. (The values are all bounded by 0 and 4 [or less].) --- modules-local/beamdyn/src/BeamDyn_Subs.f90 | 100 ++++++++++----------- 1 file changed, 46 insertions(+), 54 deletions(-) diff --git a/modules-local/beamdyn/src/BeamDyn_Subs.f90 b/modules-local/beamdyn/src/BeamDyn_Subs.f90 index 89dd81e459..6cc143f969 100644 --- a/modules-local/beamdyn/src/BeamDyn_Subs.f90 +++ b/modules-local/beamdyn/src/BeamDyn_Subs.f90 @@ -339,15 +339,10 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg) CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None !Local variables - REAL(BDKi) :: pivot(4) ! Trace of the rotation matrix and diagonal elements - REAL(BDKi) :: sm0 - REAL(BDKi) :: sm1 - REAL(BDKi) :: sm2 - REAL(BDKi) :: sm3 - REAL(BDKi) :: em - REAL(BDKi) :: temp + integer :: ipivot(1) + REAL(BDKi) :: sm(0:3) + REAL(BDKi) :: em REAL(BDKi) :: Rr(3,3) ! (possibly) corrected rotation matrix - INTEGER :: i ! case indicator INTEGER(IntKi) :: ErrStat2 ! Temporary Error status CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary Error message @@ -376,54 +371,51 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg) ! mjs--determine whether R is a valid rotation matrix and correct it if not call BD_CheckRotMat(R, Rr, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - if (ErrStat >= AbortErrLev) return + if (ErrStat >= AbortErrLev) return - ! mjs--find max value of T := Tr(Rr) and diagonal elements of Rr - pivot = (/ Rr(1, 1), Rr(2, 2), Rr(3, 3), Rr(1,1) + Rr(2,2) + Rr(3,3) /) - i = maxloc(pivot, 1) - - !> - Condition 1: \f$ \underline{\underline{R(3,3)}} > T \f$ - !! This implies that \f$ c_3^2 > c_0^2 \f$ - IF (i == 3) THEN - sm0 = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0} - sm1 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} - sm2 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} - sm3 = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3^2 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm3)), sm0 ) - - !> - Condition 2: \f$ \underline{\underline{R(2,2)}} > T \f$ - !! This implies that \f$ c_2^2 > c_0^2 \f$ - ELSEIF (i == 2) THEN - sm0 = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0} - sm1 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} - sm2 = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2^2 t_{r0} - sm3 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm2)), sm0 ) - - !> - Condition 3: \f$ \underline{\underline{R(1,1)}} > T \f$ - !! This implies that \f$ c_1^2 > c_0^2 \f$ - ELSEIF (i == 1) THEN - sm0 = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0} - sm1 = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1^2 t_{r0} - sm2 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} - sm3 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm1)), sm0 ) - - !> - Condition 4: all diagonal terms are less than the trace - ELSE + ! Note that 0 <= c_0 <= 2 since 0 <= dot_product(cc,cc) <= 16 + + ! Let's first find out which parameter is largest so we can divide by it later + sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0} + sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0} + sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0} + sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0} + ipivot = maxloc(sm) - 1 ! our sm array starts at 0, so we need to subtract 1 here to get the correct index + + select case (ipivot(1)) + case (3) + !! We need the condition that c_3 is not zero. + sm(0) = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0} + sm(1) = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0} + sm(2) = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0} + !sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0} + + case (2) + !! We need the condition that c_2 is not zero. + sm(0) = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0} + sm(1) = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0} + !sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0} + sm(3) = Rr(2,3) + Rr(3,2) ! 4 c_3 c_2 t_{r0} + + case (1) + !! We need the condition that c_1 is not zero. + sm(0) = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0} + !sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0} + sm(2) = Rr(1,2) + Rr(2,1) ! 4 c_2 c_1 t_{r0} + sm(3) = Rr(1,3) + Rr(3,1) ! 4 c_3 c_1 t_{r0} + + case (0) + !! We need the condition that c_0 is not zero. !ipivot = 0 - sm0 = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0^2 t_{r0} - sm1 = Rr(3,2) - Rr(2,3) ! -4 c_0 c_1 t_{r0} - sm2 = Rr(1,3) - Rr(3,1) ! -4 c_0 c_2 t_{r0} - sm3 = Rr(2,1) - Rr(1,2) ! -4 c_0 c_3 t_{r0} - temp = SIGN( 2.0_BDKi*SQRT(ABS(sm0)), sm0 ) - ENDIF - - em = sm0 + temp - em = 4.0_BDKi/em - cc(1) = em*sm1 - cc(2) = em*sm2 - cc(3) = em*sm3 + !sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0} + sm(1) = Rr(3,2) - Rr(2,3) ! 4 c_1 c_0 t_{r0} + sm(2) = Rr(1,3) - Rr(3,1) ! 4 c_2 c_0 t_{r0} + sm(3) = Rr(2,1) - Rr(1,2) ! 4 c_3 c_0 t_{r0} + end select + + em = sm(0) + SIGN( 2.0_BDKi*SQRT(sm(ipivot(1))), sm(0) ) + em = 4.0_BDKi/em ! 1 / ( 4 t_{r0} c_{ipivot} ), assuming 0 <= c_0 < 4 and c_{ipivot} > 0 + cc = em*sm(1:3) END SUBROUTINE BD_CrvExtractCrv !------------------------------------------------------------------------------ From f65ba70debcaf1e7a202d16aa7a814673f49aabd Mon Sep 17 00:00:00 2001 From: Rafael M Mudafort Date: Thu, 27 Jul 2017 22:41:12 -0600 Subject: [PATCH 5/5] use ReKi in OrcaDriver_Subs for support with single precision --- .../src/OrcaDriver_Subs.f90 | 43 +++++++++---------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 b/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 index e1306b5b48..e07b6105c1 100644 --- a/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 +++ b/modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90 @@ -1907,12 +1907,13 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName ! Temporary local variables INTEGER(IntKi) :: ErrStatTmp !< Temporary variable for the status of error message - CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite' + CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite' CHARACTER(2048) :: ErrMsgTmp !< Temporary variable for the error message INTEGER(IntKi) :: LenErrMsgTmp !< Length of ErrMsgTmp (for getting WindGrid info) REAL(ReKi) :: rotdisp(3) !< Rotational displacement (euler angles) INTEGER(IntKi) :: I !< Generic counter - + REAL(ReKi) :: outputArray(13) + CHARACTER(47) :: PointsOutputFmt !< Format specifier for the output file for wave elevation series CHARACTER(3) :: AngleUnit !< Units for the angle @@ -1972,32 +1973,30 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' ENDIF - rotdisp = GetSmllRotAngs ( u%PtfmMesh%Orientation(:,:,1), ErrStatTmp, ErrMsgTmp ) CALL SetErrStat( ErrStatTmp, ErrMsgTmp, ErrStat, ErrMsg, RoutineName ) - IF ( AnglesInDegrees ) THEN - CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), & - u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), & - rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, & - u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & - u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /), & - '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' + outputArray = (/ REAL(Time, ReKi), & + REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), & + rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, & + u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & + u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /) ELSE - CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), & - u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), & - rotdisp(1), rotdisp(2), rotdisp(3), & - u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & - u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /), & - '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) - WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' + outputArray = (/ REAL(Time, ReKi), & + REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), & + REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), & + rotdisp(1), rotdisp(2), rotdisp(3), & + u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), & + u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /) ENDIF - - + + CALL WrNumAryFileNR( OutUnit, outputArray, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) + CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp ) + WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) '' END SUBROUTINE PointsForce_OutputWrite