Skip to content

Commit e28e1a0

Browse files
andrew-plattandrew-platt
authored
Merge pull request OpenFAST#2399 from deslaughter/f/bd-perf
BeamDyn performance improvements
2 parents cca59dc + e1d3cdb commit e28e1a0

File tree

8 files changed

+525
-444
lines changed

8 files changed

+525
-444
lines changed

‎modules/beamdyn/src/BeamDyn.f90

+396-334
Large diffs are not rendered by default.

‎modules/beamdyn/src/BeamDyn_IO.f90

+10-6
Original file line numberDiff line numberDiff line change
@@ -2472,20 +2472,24 @@ SUBROUTINE Perturb_x( p, fieldIndx, node, dof, perturb_sign, x, dx )
24722472

24732473
REAL(R8Ki) :: orientation(3,3)
24742474
REAL(R8Ki) :: rotation(3,3)
2475+
REAL(R8Ki) :: CrvPerturb(3), CrvBase(3)
24752476

24762477
dx = p%dx(dof)
24772478

24782479
if (fieldIndx==1) then
24792480
if (dof < 4) then ! translational displacement
24802481
x%q( dof, node ) = x%q( dof, node ) + dx * perturb_sign
24812482
else ! w-m parameters
2482-
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
2483-
orientation = transpose(rotation)
24842483

2485-
CALL PerturbOrientationMatrix( orientation, dx * perturb_sign, dof-3 ) ! NOTE: call not using DCM_logmap
2486-
2487-
rotation = transpose(orientation)
2488-
call BD_CrvExtractCrv( rotation, x%q( 4:6, node ), ErrStat2, ErrMsg2 ) ! return the w-m parameters of the new orientation
2484+
! Calculate perturbation in WM parameters
2485+
CrvPerturb = 0.0_R8Ki
2486+
CrvPerturb(dof-3) = 4.0_R8Ki * tan(dx * perturb_sign / 4.0_R8Ki)
2487+
2488+
! Get base rotation in WM parameters
2489+
CrvBase = x%q(4:6, node)
2490+
2491+
! Compose pertubation and base rotation and store in state
2492+
call BD_CrvCompose(x%q(4:6, node), CrvPerturb, CrvBase, FLAG_R1R2)
24892493
end if
24902494
else
24912495
x%dqdt( dof, node ) = x%dqdt( dof, node ) + dx * perturb_sign

‎modules/beamdyn/src/BeamDyn_Subs.f90

+25-30
Original file line numberDiff line numberDiff line change
@@ -242,32 +242,31 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag)
242242
REAL(BDKi), INTENT( OUT):: rr(3) !< Composed rotation
243243
REAL(BDKi), INTENT(IN ):: pp(3) !< Input rotation 1
244244
REAL(BDKi), INTENT(IN ):: qq(3) !< Input rotation 2
245-
INTEGER ,INTENT(IN ):: flag !< Option flag
245+
INTEGER, INTENT(IN ):: flag !< Option flag
246246

247-
REAL(BDKi) :: pp0
248-
REAL(BDKi) :: p(3)
249-
REAL(BDKi) :: qq0
250-
REAL(BDKi) :: q(3)
247+
REAL(BDKi) :: pp0, p(3)
248+
REAL(BDKi) :: qq0, q(3)
251249
REAL(BDKi) :: tr1
252-
REAL(BDKi) :: Delta1
253-
REAL(BDKi) :: Delta2
250+
REAL(BDKi) :: Delta1, Delta2
254251
REAL(BDKi) :: dd1
255252
REAL(BDKi) :: dd2
256253

257254

258-
! Set the local values pp and qq allowing for the transpose
259-
260-
IF(flag==FLAG_R1TR2 .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of first rotation parameter
261-
p = -pp
262-
ELSE
263-
p = pp
264-
ENDIF
265-
266-
IF(flag==FLAG_R1R2T .OR. flag==FLAG_R1TR2T) THEN ! "transpose" (negative) of second rotation parameter
267-
q = -qq
268-
ELSE
269-
q = qq
270-
ENDIF
255+
! Set the local values pp (R1) and qq (R2) and apply transpose based on flag value
256+
select case (flag)
257+
case (FLAG_R1R2)
258+
p = pp ! R1
259+
q = qq ! R2
260+
case (FLAG_R1R2T)
261+
p = pp ! R1
262+
q = -qq ! R2^T
263+
case (FLAG_R1TR2)
264+
p = -pp ! R1^T
265+
q = qq ! R2
266+
case (FLAG_R1TR2T)
267+
p = -pp ! R1^T
268+
q = -qq ! R2^T
269+
end select
271270

272271
!> ## Composing the resulting Wiener-Milenkovic parameter
273272
!!
@@ -289,27 +288,23 @@ SUBROUTINE BD_CrvCompose( rr, pp, qq, flag)
289288
!!
290289
!!
291290

292-
293291
! Calculate pp0 and qq0. See Bauchau for the mathematics here (equations 8 to 9 and interviening text)
294292

295293
pp0 = 2.0_BDKi - dot_product(p,p) / 8.0_BDKi ! p_0
296294
qq0 = 2.0_BDKi - dot_product(q,q) / 8.0_BDKi ! q_0
297295

298296
Delta1 = (4.0_BDKi - pp0) * (4.0_BDKi - qq0) ! Delta_1 in Bauchau
299297
Delta2 = pp0 * qq0 - dot_product(p,q) ! Delta_2 in Bauchau
300-
dd1 = Delta1 + Delta2 ! Denomimator term for \Delta_2 >= 0
301-
dd2 = Delta1 - Delta2 ! Denomimator term for \Delta_2 < 0
302298

303-
! Rescaling to remove singularities at +/- 2 \pi
304-
! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation.
305-
IF ( Delta2 >= 0.0_BDKi ) THEN
306-
tr1 = 4.0_BDKi / dd1
299+
! Rescaling to remove singularities at +/- 2 \pi
300+
! Note: changed this to test on \Delta_2 (instead of dd1 > dd2) for better consistency with documentation.
301+
IF (Delta2 >= 0.0_BDKi) THEN
302+
tr1 = 4.0_BDKi / (Delta1 + Delta2)
307303
ELSE
308-
tr1 = -4.0_BDKi / dd2
304+
tr1 = -4.0_BDKi / (Delta1 - Delta2)
309305
ENDIF
310306

311-
rr = tr1 * (qq0*p + pp0*q + cross_product(p,q))
312-
307+
rr = tr1 * (qq0*p + pp0*q + Cross_Product(p,q))
313308

314309
END SUBROUTINE BD_CrvCompose
315310

‎modules/beamdyn/tests/test_tools.F90

+2-2
Original file line numberDiff line numberDiff line change
@@ -173,8 +173,8 @@ type(BD_MiscVarType) function simpleMiscVarType(nqp, dof_node, elem_total, nodes
173173
call AllocAry(m%qp%RR0mEta, 3, nqp, elem_total, 'qp_RR0mEta', ErrStat, ErrMsg)
174174
call AllocAry(m%DistrLoad_QP, 6, nqp, elem_total, 'DistrLoad_QP', ErrStat, ErrMsg)
175175

176-
CALL AllocAry(m%qp%uuu, dof_node ,nqp,elem_total, 'm%qp%uuu displacement at quadrature point',ErrStat,ErrMsg)
177-
CALL AllocAry(m%qp%uup, dof_node/2,nqp,elem_total, 'm%qp%uup displacement prime at quadrature point',ErrStat,ErrMsg)
176+
call AllocAry(m%qp%uuu, dof_node, nqp, elem_total, 'm%qp%uuu displacement at quadrature point', ErrStat, ErrMsg)
177+
call AllocAry(m%qp%uup, dof_node, nqp, elem_total, 'm%qp%uup displacement prime at quadrature point', ErrStat, ErrMsg)
178178

179179
! E1, kappa -- used in force calculations
180180
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

+72-42
Original file line numberDiff line numberDiff line change
@@ -2020,33 +2020,73 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
20202020

20212021
IF (.NOT. SrcMesh%Initialized) RETURN !bjj: maybe we should first CALL MeshDestroy(DestMesh,ErrStat, ErrMess)
20222022

2023-
IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_SIBLING .OR. CtrlCode .EQ. MESH_COUSIN ) THEN
2024-
2025-
IF (CtrlCode .EQ. MESH_NEWCOPY) THEN
2026-
IOS_l = SrcMesh%IOS
2027-
Force_l = SrcMesh%FieldMask(MASKID_FORCE)
2028-
Moment_l = SrcMesh%FieldMask(MASKID_MOMENT)
2029-
Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION)
2030-
TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP)
2031-
TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL)
2032-
RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL)
2033-
TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC)
2034-
RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC)
2035-
nScalars_l = SrcMesh%nScalars
2036-
ELSE ! Sibling or cousin
2037-
IOS_l = SrcMesh%IOS ; IF ( PRESENT(IOS) ) IOS_l = IOS
2038-
Force_l = .FALSE. ; IF ( PRESENT(Force) ) Force_l = Force
2039-
Moment_l = .FALSE. ; IF ( PRESENT(Moment) ) Moment_l = Moment
2040-
Orientation_l = .FALSE. ; IF ( PRESENT(Orientation) ) Orientation_l = Orientation
2041-
TranslationDisp_l = .FALSE. ; IF ( PRESENT(TranslationDisp) ) TranslationDisp_l = TranslationDisp
2042-
TranslationVel_l = .FALSE. ; IF ( PRESENT(TranslationVel) ) TranslationVel_l = TranslationVel
2043-
RotationVel_l = .FALSE. ; IF ( PRESENT(RotationVel) ) RotationVel_l = RotationVel
2044-
TranslationAcc_l = .FALSE. ; IF ( PRESENT(TranslationAcc) ) TranslationAcc_l = TranslationAcc
2045-
RotationAcc_l = .FALSE. ; IF ( PRESENT(RotationAcc) ) RotationAcc_l = RotationAcc
2046-
nScalars_l = 0 ; IF ( PRESENT(nScalars) ) nScalars_l = nScalars
2047-
END IF
2048-
2049-
IF ( CtrlCode .EQ. MESH_NEWCOPY .OR. CtrlCode .EQ. MESH_COUSIN ) THEN
2023+
select case (CtrlCode)
2024+
case (MESH_NEWCOPY)
2025+
IOS_l = SrcMesh%IOS
2026+
Force_l = SrcMesh%FieldMask(MASKID_FORCE)
2027+
Moment_l = SrcMesh%FieldMask(MASKID_MOMENT)
2028+
Orientation_l = SrcMesh%FieldMask(MASKID_ORIENTATION)
2029+
TranslationDisp_l = SrcMesh%FieldMask(MASKID_TRANSLATIONDISP)
2030+
TranslationVel_l = SrcMesh%FieldMask(MASKID_TRANSLATIONVEL)
2031+
RotationVel_l = SrcMesh%FieldMask(MASKID_ROTATIONVEL)
2032+
TranslationAcc_l = SrcMesh%FieldMask(MASKID_TRANSLATIONACC)
2033+
RotationAcc_l = SrcMesh%FieldMask(MASKID_ROTATIONACC)
2034+
nScalars_l = SrcMesh%nScalars
2035+
case (MESH_SIBLING, MESH_COUSIN)
2036+
IF ( PRESENT(IOS) ) then
2037+
IOS_l = IOS
2038+
else
2039+
IOS_l = SrcMesh%IOS
2040+
end if
2041+
IF ( PRESENT(Force) ) then
2042+
Force_l = Force
2043+
else
2044+
Force_l = .FALSE.
2045+
end if
2046+
IF ( PRESENT(Moment) ) then
2047+
Moment_l = Moment
2048+
else
2049+
Moment_l = .FALSE.
2050+
end if
2051+
IF ( PRESENT(Orientation) ) then
2052+
Orientation_l = Orientation
2053+
else
2054+
Orientation_l = .FALSE.
2055+
end if
2056+
IF ( PRESENT(TranslationDisp) ) then
2057+
TranslationDisp_l = TranslationDisp
2058+
else
2059+
TranslationDisp_l = .FALSE.
2060+
end if
2061+
IF ( PRESENT(TranslationVel) ) then
2062+
TranslationVel_l = TranslationVel
2063+
else
2064+
TranslationVel_l = .FALSE.
2065+
end if
2066+
IF ( PRESENT(RotationVel) ) then
2067+
RotationVel_l = RotationVel
2068+
else
2069+
RotationVel_l = .FALSE.
2070+
end if
2071+
IF ( PRESENT(TranslationAcc) ) then
2072+
TranslationAcc_l = TranslationAcc
2073+
else
2074+
TranslationAcc_l = .FALSE.
2075+
end if
2076+
IF ( PRESENT(RotationAcc) ) then
2077+
RotationAcc_l = RotationAcc
2078+
else
2079+
RotationAcc_l = .FALSE.
2080+
end if
2081+
IF ( PRESENT(nScalars) ) then
2082+
nScalars_l = nScalars
2083+
else
2084+
nScalars_l = 0
2085+
end if
2086+
end select
2087+
2088+
select case (CtrlCode)
2089+
case (MESH_NEWCOPY, MESH_COUSIN)
20502090

20512091
CALL MeshCreate( DestMesh, IOS=IOS_l, Nnodes=SrcMesh%Nnodes, ErrStat=ErrStat, ErrMess=ErrMess &
20522092
,Force=Force_l &
@@ -2123,7 +2163,7 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
21232163

21242164
DestMesh%RemapFlag = SrcMesh%RemapFlag
21252165

2126-
ELSE IF ( CtrlCode .EQ. MESH_SIBLING ) THEN
2166+
case (MESH_SIBLING)
21272167
!bjj: we should make sure the mesh has been committed, otherwise the element lists haven't been created, yet (and thus not shared)
21282168
IF ( ASSOCIATED(SrcMesh%SiblingMesh) ) THEN
21292169
ErrStat = ErrID_Fatal
@@ -2165,25 +2205,15 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
21652205
DestMesh%maxelemlist = SrcMesh%maxelemlist
21662206
DestMesh%nextelem = SrcMesh%nextelem
21672207

2168-
2169-
ENDIF
2170-
2171-
DO i = 1, NELEMKINDS
2172-
IF ( ASSOCIATED(SrcMesh%ElemTable) ) THEN
2173-
ENDIF
2174-
IF ( ASSOCIATED(DestMesh%ElemTable) ) THEN
2175-
ENDIF
2176-
ENDDO
2177-
2178-
ELSE IF ( CtrlCode .EQ. MESH_UPDATECOPY ) THEN
2208+
case (MESH_UPDATECOPY)
21792209

21802210
IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN
21812211
ErrStat = ErrID_Fatal
21822212
ErrMess = "MeshCopy: MESH_UPDATECOPY of meshes with different numbers of nodes."
21832213
RETURN
21842214
ENDIF
21852215

2186-
ELSE IF ( CtrlCode .EQ. MESH_UPDATEREFERENCE ) THEN
2216+
case (MESH_UPDATEREFERENCE)
21872217

21882218
IF ( SrcMesh%nNodes .NE. DestMesh%nNodes ) THEN
21892219
ErrStat = ErrID_Fatal
@@ -2195,11 +2225,11 @@ SUBROUTINE MeshCopy( SrcMesh, DestMesh, CtrlCode, ErrStat , ErrMess &
21952225
DestMesh%RefOrientation = SrcMesh%RefOrientation
21962226
DestMesh%RemapFlag = SrcMesh%RemapFlag
21972227

2198-
ELSE
2228+
case default
21992229
ErrStat = ErrID_Fatal
22002230
ErrMess = 'MeshCopy: Invalid CtrlCode.'
22012231
RETURN
2202-
ENDIF
2232+
end select
22032233

22042234
! These aren't shared between siblings, so they get copied, no matter what the CtrlCode:
22052235

‎modules/nwtc-library/src/ModMesh_Mapping.f90

+15-25
Original file line numberDiff line numberDiff line change
@@ -3262,38 +3262,28 @@ FUNCTION GetLoadsScaleFactor( Src )
32623262

32633263
! LOCAL:
32643264
INTEGER :: I, j
3265-
REAL(ReKi) :: MaxLoad
3265+
REAL(ReKi) :: MaxLoad, MaxForce, MaxMoment
32663266

3267+
IF ( Src%FIELDMASK( MASKID_FORCE ) ) then
3268+
MaxForce = maxval(abs(src%Force))
3269+
else
3270+
MaxForce = 0.0_ReKi
3271+
end if
32673272

3268-
GetLoadsScaleFactor = 1.0
3269-
MaxLoad = 0.0
3270-
3271-
IF ( Src%FIELDMASK( MASKID_FORCE ) ) THEN
3272-
3273-
DO I=1,Src%Nnodes
3274-
DO J=1,3
3275-
MaxLoad = MAX(MaxLoad, ABS(Src%Force(j,I) ) )
3276-
END DO
3277-
END DO
3278-
3279-
END IF
3280-
3273+
IF ( Src%FIELDMASK( MASKID_MOMENT ) ) then
3274+
MaxMoment = maxval(abs(src%Moment))
3275+
else
3276+
MaxMoment = 0.0_ReKi
3277+
end if
3278+
3279+
MaxLoad = max(MaxForce, MaxMoment)
32813280

3282-
IF ( Src%FIELDMASK( MASKID_MOMENT ) ) THEN
3283-
3284-
DO I=1,Src%Nnodes
3285-
DO J=1,3
3286-
MaxLoad = MAX(MaxLoad, ABS(Src%Moment(j,I) ) )
3287-
END DO
3288-
END DO
3289-
3290-
END IF
3291-
32923281
IF ( MaxLoad > 10. ) THEN
32933282
GetLoadsScaleFactor = 10**MIN( NINT(log10(MaxLoad)), 15 ) ! Let's not get carried away and cause overflow; 10E15 is as far as we'll go
3283+
else
3284+
GetLoadsScaleFactor = 1.0_ReKi
32943285
END IF
32953286

3296-
32973287
END FUNCTION GetLoadsScaleFactor
32983288
!----------------------------------------------------------------------------------------------------------------------------------
32993289
SUBROUTINE CreateLoadMap_P_to_P( Src, Dest, MeshMap, ErrStat, ErrMsg )

‎modules/nwtc-library/src/NWTC_Num.f90

+4-4
Original file line numberDiff line numberDiff line change
@@ -485,7 +485,7 @@ END SUBROUTINE ConvertUnitsToEngr
485485
!> This function computes the cross product of two 3-element arrays (resulting in a vector): \n
486486
!! cross_product = Vector1 \f$\times\f$ Vector2 \n
487487
!! Use cross_product (nwtc_num::cross_product) instead of directly calling a specific routine in the generic interface.
488-
FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
488+
PURE FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
489489

490490
! Argument declarations.
491491

@@ -505,7 +505,7 @@ FUNCTION Cross_ProductR4(Vector1, Vector2) result(CProd)
505505
END FUNCTION Cross_ProductR4
506506
!=======================================================================
507507
!> \copydoc nwtc_num::cross_productr4
508-
FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
508+
PURE FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
509509

510510
! Argument declarations.
511511

@@ -525,7 +525,7 @@ FUNCTION Cross_ProductR4R8(Vector1, Vector2) result(CProd)
525525
END FUNCTION Cross_ProductR4R8
526526
!=======================================================================
527527
!> \copydoc nwtc_num::cross_productr4
528-
FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
528+
PURE FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
529529

530530
! Argument declarations.
531531

@@ -545,7 +545,7 @@ FUNCTION Cross_ProductR8(Vector1, Vector2) result(CProd)
545545
END FUNCTION Cross_ProductR8
546546
!=======================================================================
547547
!> \copydoc nwtc_num::cross_productr4
548-
FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd)
548+
PURE FUNCTION Cross_ProductR8R4(Vector1, Vector2) result(CProd)
549549

550550
! Argument declarations.
551551

‎reg_tests/r-test

Submodule r-test updated 20 files

0 commit comments

Comments
 (0)