BeamDyn performance improvements

modules/beamdyn/src/BeamDyn_IO.f90

@@ -2472,20 +2472,24 @@ SUBROUTINE Perturb_x( p, fieldIndx, node, dof, perturb_sign, x, dx )
 
  REAL(R8Ki) :: orientation(3,3)
  REAL(R8Ki) :: rotation(3,3)
+ REAL(R8Ki) :: CrvPerturb(3), CrvBase(3)
 
  dx = p%dx(dof)
 
  if (fieldIndx==1) then
  if (dof < 4) then ! translational displacement
  x%q( dof, node ) = x%q( dof, node ) + dx * perturb_sign
  else ! w-m parameters
- call BD_CrvMatrixR( x%q( 4:6, node ), rotation ) ! returns the rotation matrix (transpose of DCM) that was stored in the state as a w-m parameter
- orientation = transpose(rotation)
 
- CALL PerturbOrientationMatrix( orientation, dx * perturb_sign, dof-3 ) ! NOTE: call not using DCM_logmap
-
- rotation = transpose(orientation)
- call BD_CrvExtractCrv( rotation, x%q( 4:6, node ), ErrStat2, ErrMsg2 ) ! return the w-m parameters of the new orientation
+ ! Calculate perturbation in WM parameters
+ CrvPerturb = 0.0_R8Ki
+ CrvPerturb(dof-3) = 4.0_R8Ki * tan(dx * perturb_sign / 4.0_R8Ki)
+
+ ! Get base rotation in WM parameters
+ CrvBase = x%q(4:6, node)
+
+ ! Compose pertubation and base rotation and store in state
+ call BD_CrvCompose(x%q(4:6, node), CrvPerturb, CrvBase, FLAG_R1R2)
  end if
  else
  x%dqdt( dof, node ) = x%dqdt( dof, node ) + dx * perturb_sign

modules/beamdyn/src/BeamDyn_Subs.f90

@@ -242,32 +242,31 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag)
  REAL(BDKi), INTENT( OUT):: rr(3) !< Composed rotation
  REAL(BDKi), INTENT(IN ):: pp(3) !< Input rotation 1
  REAL(BDKi), INTENT(IN ):: qq(3) !< Input rotation 2
- INTEGER ,INTENT(IN ):: flag !< Option flag
+ INTEGER, INTENT(IN ):: flag !< Option flag
 
- REAL(BDKi) :: pp0
- REAL(BDKi) :: p(3)
- REAL(BDKi) :: qq0
- REAL(BDKi) :: q(3)
+ REAL(BDKi) :: pp0, p(3)
+ REAL(BDKi) :: qq0, q(3)
  REAL(BDKi) :: tr1
- REAL(BDKi) :: Delta1
- REAL(BDKi) :: Delta2
+ REAL(BDKi) :: Delta1, Delta2
  REAL(BDKi) :: dd1
  REAL(BDKi) :: dd2
 
 
- ! Set the local values pp and qq allowing for the transpose
-
- IF(flag==FLAG_R1TR2 .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of first rotation parameter
- p = -pp
- ELSE
- p = pp
- ENDIF
-
- IF(flag==FLAG_R1R2T .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of second rotation parameter
- q = -qq
- ELSE
- q = qq
- ENDIF
+ ! Set the local values pp (R1) and qq (R2) and apply transpose based on flag value
+ select case (flag)
+ case (FLAG_R1R2)
+ p = pp ! R1
+ q = qq ! R2
+ case (FLAG_R1R2T)
+ p = pp ! R1
+ q = -qq ! R2^T
+ case (FLAG_R1TR2)
+ p = -pp ! R1^T
+ q = qq ! R2
+ case (FLAG_R1TR2T)
+ p = -pp ! R1^T
+ q = -qq ! R2^T
+ end select
 
  !> ## Composing the resulting Wiener-Milenkovic parameter
  !!
@@ -289,27 +288,23 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag)
  !!
  !!
 
-
  ! Calculate pp0 and qq0. See Bauchau for the mathematics here (equations 8 to 9 and interviening text)
 
  pp0 = 2.0_BDKi - dot_product(p,p) / 8.0_BDKi ! p_0
  qq0 = 2.0_BDKi - dot_product(q,q) / 8.0_BDKi ! q_0
 
  Delta1 = (4.0_BDKi - pp0) * (4.0_BDKi - qq0) ! Delta_1 in Bauchau
  Delta2 = pp0 * qq0 - dot_product(p,q) ! Delta_2 in Bauchau
- dd1 = Delta1 + Delta2 ! Denomimator term for \Delta_2 >= 0
- dd2 = Delta1 - Delta2 ! Denomimator term for \Delta_2 < 0
 
-  ! Rescaling to remove singularities at +/- 2 \pi
-  ! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation.
- IF ( Delta2 >= 0.0_BDKi ) THEN
- tr1 = 4.0_BDKi / dd1
+ ! Rescaling to remove singularities at +/- 2 \pi
+ ! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation.
+ IF (Delta2 >= 0.0_BDKi) THEN
+ tr1 = 4.0_BDKi / (Delta1 + Delta2)
  ELSE
- tr1 = -4.0_BDKi / dd2
+ tr1 = -4.0_BDKi / (Delta1 - Delta2)
  ENDIF
 
- rr = tr1 * (qq0*p + pp0*q + cross_product(p,q))
-
+ rr = tr1 * (qq0*p + pp0*q + Cross_Product(p,q))
 
 END SUBROUTINE BD_CrvCompose
 

modules/beamdyn/tests/test_tools.F90

@@ -173,8 +173,8 @@ type(BD_MiscVarType) function simpleMiscVarType(nqp, dof_node, elem_total, nodes
  call AllocAry(m%qp%RR0mEta, 3, nqp, elem_total, 'qp_RR0mEta', ErrStat, ErrMsg)
  call AllocAry(m%DistrLoad_QP, 6, nqp, elem_total, 'DistrLoad_QP', ErrStat, ErrMsg)
 
- CALL AllocAry(m%qp%uuu, dof_node ,nqp,elem_total, 'm%qp%uuu displacement at quadrature point',ErrStat,ErrMsg)
- CALL AllocAry(m%qp%uup, dof_node/2,nqp,elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat,ErrMsg)
+ call AllocAry(m%qp%uuu, dof_node, nqp, elem_total, 'm%qp%uuu displacement at quadrature point', ErrStat, ErrMsg)
+ call AllocAry(m%qp%uup, dof_node, nqp, elem_total, 'm%qp%uup displacement prime at quadrature point', ErrStat, ErrMsg)
 
  ! E1, kappa -- used in force calculations
  CALL AllocAry(m%qp%E1, dof_node/2,nqp,elem_total, 'm%qp%E1 at quadrature point',ErrStat,ErrMsg)

modules/nwtc-library/src/ModMesh.f90

@@ -2020,33 +2020,73 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
 
  IF (.NOT. SrcMesh%Initialized) RETURN !bjj: maybe we should first CALL MeshDestroy(DestMesh,ErrStat, ErrMess)
 
- IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_SIBLING .OR. CtrlCode .EQ. MESH_COUSIN ) THEN
-
- IF (CtrlCode .EQ. MESH_NEWCOPY) THEN
- IOS_l = SrcMesh%IOS
- Force_l = SrcMesh%FieldMask(MASKID_FORCE) 
- Moment_l = SrcMesh%FieldMask(MASKID_MOMENT) 
- Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION) 
- TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP) 
- TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL) 
- RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL) 
- TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC) 
- RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC) 
- nScalars_l = SrcMesh%nScalars 
- ELSE ! Sibling or cousin
- IOS_l = SrcMesh%IOS ; IF ( PRESENT(IOS) ) IOS_l = IOS
- Force_l = .FALSE. ; IF ( PRESENT(Force) ) Force_l = Force
- Moment_l = .FALSE. ; IF ( PRESENT(Moment) ) Moment_l = Moment
- Orientation_l = .FALSE. ; IF ( PRESENT(Orientation) ) Orientation_l = Orientation
- TranslationDisp_l = .FALSE. ; IF ( PRESENT(TranslationDisp) ) TranslationDisp_l = TranslationDisp
- TranslationVel_l = .FALSE. ; IF ( PRESENT(TranslationVel) ) TranslationVel_l = TranslationVel
- RotationVel_l = .FALSE. ; IF ( PRESENT(RotationVel) ) RotationVel_l = RotationVel
- TranslationAcc_l = .FALSE. ; IF ( PRESENT(TranslationAcc) ) TranslationAcc_l = TranslationAcc
- RotationAcc_l = .FALSE. ; IF ( PRESENT(RotationAcc) ) RotationAcc_l = RotationAcc
- nScalars_l = 0 ; IF ( PRESENT(nScalars) ) nScalars_l = nScalars
- END IF
-
- IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_COUSIN ) THEN
+ select case (CtrlCode)
+ case (MESH_NEWCOPY)
+ IOS_l = SrcMesh%IOS
+ Force_l = SrcMesh%FieldMask(MASKID_FORCE) 
+ Moment_l = SrcMesh%FieldMask(MASKID_MOMENT) 
+ Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION) 
+ TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP) 
+ TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL) 
+ RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL) 
+ TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC) 
+ RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC) 
+ nScalars_l = SrcMesh%nScalars 
+ case (MESH_SIBLING, MESH_COUSIN)
+ IF ( PRESENT(IOS) ) then
+ IOS_l = IOS
+ else
+ IOS_l = SrcMesh%IOS
+ end if
+ IF ( PRESENT(Force) ) then
+ Force_l = Force
+ else
+ Force_l = .FALSE.
+ end if
+ IF ( PRESENT(Moment) ) then
+ Moment_l = Moment
+ else
+ Moment_l = .FALSE.
+ end if
+ IF ( PRESENT(Orientation) ) then
+ Orientation_l = Orientation
+ else
+ Orientation_l = .FALSE.
+ end if
+ IF ( PRESENT(TranslationDisp) ) then
+ TranslationDisp_l = TranslationDisp
+ else
+ TranslationDisp_l = .FALSE.
+ end if
+ IF ( PRESENT(TranslationVel) ) then
+ TranslationVel_l = TranslationVel
+ else
+ TranslationVel_l = .FALSE.
+ end if
+ IF ( PRESENT(RotationVel) ) then
+ RotationVel_l = RotationVel
+ else
+ RotationVel_l = .FALSE.
+ end if
+ IF ( PRESENT(TranslationAcc) ) then
+ TranslationAcc_l = TranslationAcc
+ else
+ TranslationAcc_l = .FALSE.
+ end if
+ IF ( PRESENT(RotationAcc) ) then
+ RotationAcc_l = RotationAcc
+ else
+ RotationAcc_l = .FALSE.
+ end if
+ IF ( PRESENT(nScalars) ) then
+ nScalars_l = nScalars
+ else
+ nScalars_l = 0 
+ end if
+ end select
+
+ select case (CtrlCode) 
+ case (MESH_NEWCOPY, MESH_COUSIN)
 
  CALL MeshCreate( DestMesh, IOS=IOS_l, Nnodes=SrcMesh%Nnodes, ErrStat=ErrStat, ErrMess=ErrMess &
  ,Force=Force_l &
@@ -2123,7 +2163,7 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
 
  DestMesh%RemapFlag = SrcMesh%RemapFlag
 
-  ELSE IF ( CtrlCode .EQ. MESH_SIBLING ) THEN
+ case (MESH_SIBLING)
 !bjj: we should make sure the mesh has been committed, otherwise the element lists haven't been created, yet (and thus not shared)
  IF ( ASSOCIATED(SrcMesh%SiblingMesh) ) THEN
  ErrStat = ErrID_Fatal
@@ -2165,25 +2205,15 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
  DestMesh%maxelemlist = SrcMesh%maxelemlist
  DestMesh%nextelem = SrcMesh%nextelem
 
-
- ENDIF
-
- DO i = 1, NELEMKINDS
- IF ( ASSOCIATED(SrcMesh%ElemTable) ) THEN
- ENDIF
- IF ( ASSOCIATED(DestMesh%ElemTable) ) THEN
- ENDIF
- ENDDO
-
- ELSE IF ( CtrlCode .EQ. MESH_UPDATECOPY ) THEN
+ case (MESH_UPDATECOPY)
 
  IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN
  ErrStat = ErrID_Fatal
  ErrMess = "MeshCopy: MESH_UPDATECOPY of meshes with different numbers of nodes."
  RETURN
  ENDIF
 
- ELSE IF ( CtrlCode .EQ. MESH_UPDATEREFERENCE ) THEN
+ case (MESH_UPDATEREFERENCE)
 
  IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN
  ErrStat = ErrID_Fatal
@@ -2195,11 +2225,11 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
  DestMesh%RefOrientation = SrcMesh%RefOrientation
  DestMesh%RemapFlag = SrcMesh%RemapFlag 
 
- ELSE
+ case default
  ErrStat = ErrID_Fatal
  ErrMess = 'MeshCopy: Invalid CtrlCode.'
  RETURN
- ENDIF
+ end select
 
  ! These aren't shared between siblings, so they get copied, no matter what the CtrlCode:
 

modules/nwtc-library/src/ModMesh_Mapping.f90

@@ -3262,38 +3262,28 @@ FUNCTION GetLoadsScaleFactor( Src )
 
  ! LOCAL:
  INTEGER :: I, j
- REAL(ReKi) :: MaxLoad
+ REAL(ReKi) :: MaxLoad, MaxForce, MaxMoment
 
+ IF ( Src%FIELDMASK( MASKID_FORCE ) ) then
+ MaxForce = maxval(abs(src%Force))
+ else
+ MaxForce = 0.0_ReKi
+ end if
 
- GetLoadsScaleFactor = 1.0
- MaxLoad = 0.0
-
- IF ( Src%FIELDMASK( MASKID_FORCE ) ) THEN
-
- DO I=1,Src%Nnodes
- DO J=1,3
- MaxLoad = MAX(MaxLoad, ABS(Src%Force(j,I) ) )
- END DO 
- END DO
-
- END IF
-
+ IF ( Src%FIELDMASK( MASKID_MOMENT ) ) then
+ MaxMoment = maxval(abs(src%Moment))
+ else
+ MaxMoment = 0.0_ReKi
+ end if
+
+ MaxLoad = max(MaxForce, MaxMoment)
 
- IF ( Src%FIELDMASK( MASKID_MOMENT ) ) THEN
-
- DO I=1,Src%Nnodes
- DO J=1,3
- MaxLoad = MAX(MaxLoad, ABS(Src%Moment(j,I) ) )
- END DO 
- END DO
-
- END IF
-
  IF ( MaxLoad > 10. ) THEN
  GetLoadsScaleFactor = 10**MIN( NINT(log10(MaxLoad)), 15 ) ! Let's not get carried away and cause overflow; 10E15 is as far as we'll go
+ else
+ GetLoadsScaleFactor = 1.0_ReKi
  END IF
 
-
 END FUNCTION GetLoadsScaleFactor
 !----------------------------------------------------------------------------------------------------------------------------------
 SUBROUTINE CreateLoadMap_P_to_P( Src, Dest, MeshMap, ErrStat, ErrMsg )

modules/nwtc-library/src/NWTC_Num.f90

@@ -485,7 +485,7 @@ END SUBROUTINE ConvertUnitsToEngr
 !> This function computes the cross product of two 3-element arrays (resulting in a vector): \n
 !! cross_product = Vector1 \f$\times\f$ Vector2 \n
 !! Use cross_product (nwtc_num::cross_product) instead of directly calling a specific routine in the generic interface.
- FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
+ PURE FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
 
  ! Argument declarations.
 
@@ -505,7 +505,7 @@ FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
  END FUNCTION Cross_ProductR4
 !=======================================================================
 !> \copydoc nwtc_num::cross_productr4
- FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
+ PURE FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
 
  ! Argument declarations.
 
@@ -525,7 +525,7 @@ FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
  END FUNCTION Cross_ProductR4R8
 !=======================================================================
 !> \copydoc nwtc_num::cross_productr4
- FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
+ PURE FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
 
  ! Argument declarations.
 
@@ -545,7 +545,7 @@ FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
  END FUNCTION Cross_ProductR8
 !=======================================================================
 !> \copydoc nwtc_num::cross_productr4
- FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd)
+ PURE FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd)
 
  ! Argument declarations.