Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Expand and improve BeamDyn unit tests #560

Merged
merged 9 commits into from
Oct 14, 2020
2 changes: 1 addition & 1 deletion modules/beamdyn/src/BeamDyn.f90
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I

ELSEIF(p%quadrature .EQ. TRAP_QUADRATURE) THEN

CALL BD_TrapezoidalPointWeight(p, InputFileData) ! computes p%QPtN and p%QPtWeight
CALL BD_TrapezoidalPointWeight(p, InputFileData%InpBl%station_eta, InputFileData%InpBl%station_total) ! computes p%QPtN and p%QPtWeight

ENDIF

Expand Down
19 changes: 9 additions & 10 deletions modules/beamdyn/src/BeamDyn_Subs.f90
Original file line number Diff line number Diff line change
Expand Up @@ -580,10 +580,11 @@ SUBROUTINE BD_GaussPointWeight(n, x, w, ErrStat, ErrMsg)
END SUBROUTINE BD_GaussPointWeight
!-----------------------------------------------------------------------------------------------------------------------------------
! This subroutine computes trapezoidal quadrature points and weights, p%QPtN and p%QPtWeight
SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData)
SUBROUTINE BD_TrapezoidalPointWeight(p, station_eta, station_total)

TYPE(BD_ParameterType),INTENT(INOUT):: p !< BeamDyn parameters
TYPE(BD_InputFile), INTENT(IN ):: InputFileData !< BeamDyn input-file data
Integer(IntKi),INTENT(IN ) :: station_total
REAL(BDKi),INTENT(IN ) :: station_eta(:)

! local variables
REAL(BDKi) :: denom ! denominator for quadrature weight computations
Expand All @@ -593,20 +594,20 @@ SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData)
INTEGER(IntKi) :: id1, j

!bjj: this assumes there is only one member



! compute the trapezoidal quadrature points, p%QPtN, and scale to range [-1,1]:
! If there is refinement, this will add new points between the specified ones. If p%refine == 1, can skip this.
p%QPtN(1) = InputFileData%InpBl%station_eta(1)
p%QPtN(1) = station_eta(1)
DO j = 2,p%nqp
indx = 1+(j-2_IntKi)/p%refine ! note use of integer math here --> (J-2)/p%refine may not be integer.
p%QPtN(j) = InputFileData%InpBl%station_eta(indx) + &
((InputFileData%InpBl%station_eta(indx+1) - InputFileData%InpBl%station_eta(indx))/p%refine) * (MOD(j-2,p%refine) + 1)
p%QPtN(j) = station_eta(indx) + &
((station_eta(indx+1) - station_eta(indx))/p%refine) * (MOD(j-2,p%refine) + 1)
ENDDO
p%QPtN = 2.0_BDKi*p%QPtN - 1.0_BDKi ! rescale range from [0, 1] to [-1,1]


! compute the trapezoidal quadrature weights, p%QPtWeight:
id1 = InputFileData%InpBl%station_total
id1 = station_total
temp_id0 = (id0 - 1)*p%refine + 1 ! Starting index in QPtN --> always going to be 1
temp_id1 = (id1 - 1)*p%refine + 1 ! ending index in QPtN --> will be size(p%QPtN)
denom = p%QPtN(temp_id1) - p%QPtN(temp_id0) ! This is the range of QPtN --> for single member, is always == 2
Expand All @@ -617,8 +618,6 @@ SUBROUTINE BD_TrapezoidalPointWeight(p, InputFileData)
ENDDO
p%QPtWeight(p%nqp) = (p%QPtN(temp_id1 ) - p%QPtN(temp_id1-1 ))/denom



END SUBROUTINE BD_TrapezoidalPointWeight

!-----------------------------------------------------------------------------------------------------------------------------------
Expand Down
4 changes: 3 additions & 1 deletion modules/beamdyn/tests/test_BD_DistrLoadCopy.F90
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ subroutine test_BD_DistrLoadCopy()
testname = "static simple beam under gravity:"

! build the parametertype, inputtype, and miscvartype
parametertype = simpleParameterType()
parametertype = simpleParameterType(1,16,16,0,1)
miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total)
inputtype = simpleInputType(parametertype%nqp, parametertype%elem_total)

Expand All @@ -41,4 +41,6 @@ subroutine test_BD_DistrLoadCopy()
@assertEqual((/ -9*(j-1)-3*(i-1)-1, -9*(j-1)-3*(i-1)-2, -9*(j-1)-3*(i-1)-3 /), miscvartype%DistrLoad_QP(4:6,i,j))
end do
end do

call BD_DestroyParam(parametertype, ErrStat, ErrMsg)
end subroutine
6 changes: 4 additions & 2 deletions modules/beamdyn/tests/test_BD_GravityForce.F90
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ subroutine test_BD_GravityForce()
baseline(4:6) = (/ 0.0, 0.0, 0.0 /)

! allocate and build the custom types
parametertype = simpleParameterType()
parametertype = simpleParameterType(1,16,16,0,1)
miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total)

gravity = getGravityInZ()
Expand All @@ -40,5 +40,7 @@ subroutine test_BD_GravityForce()

! test the values
@assertEqual(baseline, miscvartype%qp%Fg(:,1,1), tolerance, testname)


call BD_DestroyParam(parametertype, ErrStat, ErrMsg)

end subroutine
288 changes: 135 additions & 153 deletions modules/beamdyn/tests/test_BD_InitShpDerJaco.F90

Large diffs are not rendered by default.

87 changes: 87 additions & 0 deletions modules/beamdyn/tests/test_BD_MemberEta.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
module test_BD_MemberEta

! tests routine that calculates the length of the beam's reference line
! also finds element boundaries in eta coordinates

use pFUnit_mod
use BeamDyn
use NWTC_Num
use test_tools

implicit none

integer(IntKi) :: nqp
integer(IntKi) :: member_total
real(BDKi) :: total_length
real(BDKi), allocatable :: baseline_jac(:,:)
real(BDKi), allocatable :: baseline_QPtW(:)
real(BDKi), allocatable :: baseline_member_eta(:)
real(BDKi), allocatable :: test_member_eta(:)
real(BDKi) :: baseline_total_length
real(BDKi) :: test_total_length
integer(IntKi) :: ErrStat
character :: ErrMsg
character(1024) :: testname
real(BDKi) :: tolerance

contains

@test
subroutine test_BD_MemberEta_5node()

! initialize NWTC_Num constants
call SetConstants()

tolerance = 1e-14

! --------------------------------------------------------------------------
testname = "test_bd_member_eta_5node"

! this test problem is for a beam reference axis defined by 5 points
! x = 0., 0.16237631096713473, 0.25, -0.30523345382427747, -1.
! y = 0., 0.17578464768961147, 1., 2.4670724951675314, 3.5
! z = 0., 0.1481911137890286, 1.1875, 2.953849702537502, 4.
! where the underlying forth-order polynomial is
! fx[u_] = u - 2 u^3;
! fy[u_] = u/2 + 3 u^2;
! fz[u_] = 5 u^2 - u^4;
!
! exact (to mp) length is 5.627175237247959
! The Jacobian values below are calculate based on a 5-node Legendre Spectral Element
! While we give baseline Jacobian here, we check "test_BD_InitShpDerJaco.F90" elsewhere

nqp = 5 ! number of quadrature points
member_total = 1 ! 1 element

call AllocAry(baseline_jac , nqp, member_total, 'Reference Jacobian', ErrStat, ErrMsg)
call AllocAry(baseline_QPtW, nqp, 'Reference QPtWeight', ErrStat, ErrMsg)
call AllocAry(baseline_member_eta, member_total, 'Reference member_eta', ErrStat, ErrMsg)
call AllocAry(test_member_eta, member_total, 'test member_eta', ErrStat, ErrMsg)

! 5-point Guass-Legendre quadrature; see https://pomax.github.io/bezierinfo/legendre-gauss.html
baseline_QPtW(1:nqp) = (/ 0.2369268850561891, 0.4786286704993665, 0.5688888888888889, 0.4786286704993665, 0.2369268850561891 /)

! assign baseline jacobian based; these values were calculated in separate mathematica script
baseline_jac(1:nqp,1) = (/ 0.6715870058501458, 1.509599209717604, 2.861380785564901, 4.097191592895223, 4.880926263217582 /)

! total length of beam calculated in mathematica; note that for this curved beam, GL quadrature is APPROXIMATE, but
! converged rapidly; TR quadrature is not nearly as good; commented-out values are useful for understanding quadrature performance
! for curved beams.
!baseline_total_length = 5.627175237247959 ! this is actual length based on mathematica
!baseline_total_length = 5.634413547964786 ! this is approximation with 9-point Trapezoidal quadrature; 0.13% error
!baseline_total_length = 5.627202424388781 ! this is approximation with 7-point gauss quadrature; 0.0005% error
baseline_total_length = 5.626918236484061 ! this is approximation with 5-point gauss quadrature; 0.005% error (tested here)
baseline_member_eta(1) = 1. ! just one element; so member_length / total_length = 1

call BD_MemberEta(member_total, baseline_QPtW, baseline_jac, test_member_eta, test_total_length)

@assertEqual(baseline_total_length, test_total_length, tolerance, testname)
@assertEqual(baseline_member_eta, test_member_eta, tolerance, testname)

deallocate(baseline_Jac)
deallocate(baseline_QPtW)
deallocate(baseline_member_eta)
deallocate(test_member_eta)

end subroutine
end module
3 changes: 2 additions & 1 deletion modules/beamdyn/tests/test_BD_QPData_mEta_rho.F90
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ subroutine test_BD_QPData_mEta_rho()
baselineRR0mEta = (/ 0.0, 0.0, 0.0 /)

! allocate and build the custom input types
parametertype = simpleParameterType()
parametertype = simpleParameterType(1,16,16,0,1)
miscvartype = simpleMiscVarType(parametertype%nqp, parametertype%elem_total)

! allocate the results
Expand All @@ -46,4 +46,5 @@ subroutine test_BD_QPData_mEta_rho()
@assertEqual(baselineRR0mEta, miscvartype%qp%RR0mEta(:,i,j), tolerance, testname)
end do
end do
call BD_DestroyParam(parametertype, ErrStat, ErrMsg)
end subroutine
113 changes: 113 additions & 0 deletions modules/beamdyn/tests/test_BD_TrapezoidalPointWeight.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
module test_BD_TrapezoidalPointWeight

use pFUnit_mod
use BeamDyn
use NWTC_Num
use test_tools

implicit none

type(BD_ParameterType) :: p

integer(IntKi) :: nqp
integer(IntKi) :: refine
integer(IntKi) :: station_total
real(BDKi), allocatable :: station_eta(:)
real(BDKi), allocatable :: baseline_QPtN(:)
real(BDKi), allocatable :: baseline_QPtW(:)

integer(IntKi) :: ErrStat
character :: ErrMsg
character(1024) :: testname
real(BDKi) :: tolerance

contains

@test
subroutine test_BD_TrapezoidalPointWeight_2station()

! initialize NWTC_Num constants
call SetConstants()

tolerance = 1e-14

! --------------------------------------------------------------------------
testname = "test_BD_TrapezoidalPointWeight_2station_8refine"

station_total = 2

call AllocAry(station_eta, station_total, "station_eta", ErrStat, ErrMsg)

! simple case where we have enpoints only; typical with constant cross sections
station_eta(1:station_total) = (/0., 1./)

refine = 8

nqp = (station_total - 1)*refine + 1

p=simpleParameterType(1,1,nqp,0,refine)

call AllocAry(baseline_QPtN, nqp, "baseline_QPtN", ErrStat, ErrMsg)
call AllocAry(baseline_QPtW, nqp, "baseline_QPtW", ErrStat, ErrMsg)

baseline_QPtN(1:nqp) = (/ -1., -0.75, -0.5, -0.25, 0., 0.25, 0.5, 0.75, 1. /)
baseline_QPtW(1:nqp) = (/ 0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.125 /)

call BD_TrapezoidalPointWeight(p, station_eta, station_total)

@assertEqual(baseline_QPtN, p%QPtN, tolerance, testname)
@assertEqual(baseline_QPtW, p%QPtWeight, tolerance, testname)

deallocate(station_eta)
deallocate(baseline_QPtN)
deallocate(baseline_QPtW)

call BD_DestroyParam(p, ErrStat, ErrMsg)

end subroutine

@test
subroutine test_BD_TrapezoidalPointWeight_3station()

! initialize NWTC_Num constants
call SetConstants()

tolerance = 1e-14

! --------------------------------------------------------------------------
testname = "test_BD_TrapezoidalPointWeight_3station_2refine"

! provide three stations, unequally distributed
! refine by factor of two

station_total = 3

call AllocAry(station_eta, station_total, "station_eta", ErrStat, ErrMsg)

station_eta(1:station_total) = (/0., 0.25, 1./)

refine = 2

nqp = (station_total - 1)*refine + 1

p=simpleParameterType(1,1,nqp,0,refine)

call AllocAry(baseline_QPtN, nqp, "baseline_QPtN", ErrStat, ErrMsg)
call AllocAry(baseline_QPtW, nqp, "baseline_QPtW", ErrStat, ErrMsg)

baseline_QPtN(1:nqp) = (/ -1., -0.75, -0.5, 0.25, 1. /)
baseline_QPtW(1:nqp) = (/ 0.125, 0.25, 0.125+0.5*0.75, 0.75, 0.5*0.75/)

call BD_TrapezoidalPointWeight(p, station_eta, station_total)

@assertEqual(baseline_QPtN, p%QPtN, tolerance, testname)
@assertEqual(baseline_QPtW, p%QPtWeight, tolerance, testname)

deallocate(station_eta)
deallocate(baseline_QPtN)
deallocate(baseline_QPtW)

call BD_DestroyParam(p, ErrStat, ErrMsg)

end subroutine
end module
Loading