Skip to content

Initial AeroMap changes for ElastoDyn and BeamDyn #1631

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.

Sign up for GitHub

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

Merged
andrew-platt merged 7 commits into from
Jun 20, 2023
Merged

Initial AeroMap changes for ElastoDyn and BeamDyn #1631

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
4c3185f
ED: update for AeroMaps
bjonkman Jun 16, 2023
6248bc3
BD: update for AeroMaps
bjonkman Jun 16, 2023
a874c64
update logic in ED Nodal output
bjonkman Jun 14, 2023
752fb23
ED: replace code with call to FindValidChannelIndx()
bjonkman Jun 14, 2023
605a0be
Merge remote-tracking branch 'NREL/dev' into f/AeroMap_BD_ED
bjonkman Jun 16, 2023
743e1a0
update comments for FVW
bjonkman Oct 6, 2021
ef8e99b
revert r-test change
bjonkman Jun 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions modules/aerodyn/src/FVW_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -462,6 +462,7 @@ logical function have_nan(p, m, x, z, u, label)
character(len=*), intent(in) :: label !< label for print
integer :: iW
have_nan=.False.
!bjj: If we used Is_NaN (or a version of it for this data type) instead of isnan, I'd get fewer compiler warnings that "Fortran 2003 does not allow this intrinsic procedure."
do iW = 1,size(p%W)
if (any(isnan(x%W(iW)%r_NW))) then
print*,trim(label),'NaN in W(iW)%r_NW'//trim(num2lstr(iW))
Expand All @@ -480,11 +481,11 @@ logical function have_nan(p, m, x, z, u, label)
have_nan=.True.
endif
if (any(isnan(x%W(iW)%Eps_NW))) then
print*,trim(label),'NaN in G_FW'//trim(num2lstr(iW))
print*,trim(label),'NaN in E_NW'
have_nan=.True.
endif
if (any(isnan(x%W(iW)%Eps_FW))) then
print*,trim(label),'NaN in G_FW'//trim(num2lstr(iW))
print*,trim(label),'NaN in E_FW'
have_nan=.True.
endif
if (any(isnan(z%W(iW)%Gamma_LL))) then
Expand Down
175 changes: 101 additions & 74 deletions modules/beamdyn/src/BeamDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -261,6 +261,10 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I
! allocate and initialize misc vars (do this after initializing input and output meshes):
call Init_MiscVars(p, u, y, MiscVar, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
if (ErrStat >= AbortErrLev) then
call cleanup()
return
end if



Expand Down Expand Up @@ -330,7 +334,7 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I
!............................................................................................
! Initialize Jacobian:
!............................................................................................
if (InitInp%Linearize) then
if (InitInp%Linearize .or. p%CompAeroMaps) then
call Init_Jacobian( p, u, y, MiscVar, InitOut, ErrStat2, ErrMsg2)
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
end if
Expand Down Expand Up @@ -922,7 +926,8 @@ subroutine SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg)
ErrStat = ErrID_None
ErrMsg = ""


p%CompAeroMaps = InitInp%CompAeroMaps

! Global position vector
p%GlbPos = InitInp%GlbPos

Expand Down Expand Up @@ -1037,6 +1042,14 @@ subroutine SetParameters(InitInp, InputFileData, p, ErrStat, ErrMsg)
if (ErrStat >= AbortErrLev) return


if (p%CompAeroMaps) then
if (p%BldMotionNodeLoc /= BD_MESH_FE) then
! call SetErrStat(ErrID_Warn, "BeamDyn aero maps must have outputs at FEA nodes; this is different than time-series behavior.", ErrStat, ErrMsg, RoutineName )
p%BldMotionNodeLoc = BD_MESH_FE
call SetErrStat(ErrID_Fatal, "BeamDyn aero maps must have outputs at FEA nodes, which requires Gaussian quadrature. Update the input file.", ErrStat, ErrMsg, RoutineName )
return
end if
end if

!...............................................
! Set start and end node index for each elements
Expand Down Expand Up @@ -2895,21 +2908,21 @@ SUBROUTINE BD_QPDataVelocity( p, x, m )

elem_start = p%node_elem_idx(nelem,1)

DO idx_qp=1,p%nqp
DO idx_qp=1,p%nqp

!> Calculate the values for the
!> Calculate the values for the

! Initialize to zero for summation
m%qp%vvv(:,idx_qp,nelem) = 0.0_BDKi
m%qp%vvp(:,idx_qp,nelem) = 0.0_BDKi
! Initialize to zero for summation
m%qp%vvv(:,idx_qp,nelem) = 0.0_BDKi
m%qp%vvp(:,idx_qp,nelem) = 0.0_BDKi

! Calculate the velocity term, velocity prime (derivative of velocity with respect to X-axis), and acceleration terms
DO idx_node=1,p%nodes_per_elem
m%qp%vvv(:,idx_qp,nelem) = m%qp%vvv(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * x%dqdt(:,elem_start-1+idx_node)
m%qp%vvp(:,idx_qp,nelem) = m%qp%vvp(:,idx_qp,nelem) + p%ShpDer(idx_node,idx_qp)/p%Jacobian(idx_qp,nelem) * x%dqdt(:,elem_start-1+idx_node)
ENDDO
! Calculate the velocity term, velocity prime (derivative of velocity with respect to X-axis), and acceleration terms
DO idx_node=1,p%nodes_per_elem
m%qp%vvv(:,idx_qp,nelem) = m%qp%vvv(:,idx_qp,nelem) + p%Shp(idx_node,idx_qp) * x%dqdt(:,elem_start-1+idx_node)
m%qp%vvp(:,idx_qp,nelem) = m%qp%vvp(:,idx_qp,nelem) + p%ShpDer(idx_node,idx_qp)/p%Jacobian(idx_qp,nelem) * x%dqdt(:,elem_start-1+idx_node)
ENDDO

ENDDO
ENDDO

ENDDO

Expand Down Expand Up @@ -5906,59 +5919,64 @@ SUBROUTINE BD_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrM
end if
end if

if (p%CompAeroMaps) then
dYdu = 0.0_R8Ki
else

! make a copy of outputs because we will need two for the central difference computations (with orientations)
call BD_CopyOutput( y, y_p, MESH_NEWCOPY, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call BD_CopyOutput( y, y_m, MESH_NEWCOPY, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if
! make a copy of outputs because we will need two for the central difference computations (with orientations)
call BD_CopyOutput( y, y_p, MESH_NEWCOPY, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call BD_CopyOutput( y, y_m, MESH_NEWCOPY, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if

do i=1,size(p%Jac_u_indx,1)
do i=1,size(p%Jac_u_indx,1)

! get u_op + delta_p u
call BD_CopyInput( u, u_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
call Perturb_u( p, i, 1, u_perturb, delta_p )
! get u_op + delta_p u
call BD_CopyInput( u, u_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
call Perturb_u( p, i, 1, u_perturb, delta_p )

! compute y at u_op + delta_p u
call BD_CalcOutput( t, u_perturb, p, x, xd, z, OtherState, y_p, m, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
! compute y at u_op + delta_p u
call BD_CalcOutput( t, u_perturb, p, x, xd, z, OtherState, y_p, m, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later

! get u_op - delta_m u
call BD_CopyInput( u, u_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
call Perturb_u( p, i, -1, u_perturb, delta_m )
! get u_op - delta_m u
call BD_CopyInput( u, u_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
call Perturb_u( p, i, -1, u_perturb, delta_m )

! compute y at u_op - delta_m u
call BD_CalcOutput( t, u_perturb, p, x, xd, z, OtherState, y_m, m, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later
! compute y at u_op - delta_m u
call BD_CalcOutput( t, u_perturb, p, x, xd, z, OtherState, y_m, m, ErrStat2, ErrMsg2 )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later

! get central difference:
call Compute_dY( p, y_p, y_m, delta_p, dYdu(:,i) )
! get central difference:
call Compute_dY( p, y_p, y_m, delta_p, dYdu(:,i) )

end do
end do


if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if
call BD_DestroyOutput( y_p, ErrStat2, ErrMsg2 ) ! we don't need this any more
call BD_DestroyOutput( y_m, ErrStat2, ErrMsg2 ) ! we don't need this any more
if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if
call BD_DestroyOutput( y_p, ErrStat2, ErrMsg2 ) ! we don't need this any more
call BD_DestroyOutput( y_m, ErrStat2, ErrMsg2 ) ! we don't need this any more

if (p%RelStates) then
call BD_JacobianPContState_noRotate( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, dYdx=m%lin_C )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if
dYdu = dYdu + matmul(m%lin_C, RelState_x)
end if
if (p%RelStates) then
call BD_JacobianPContState_noRotate( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, dYdx=m%lin_C )
call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
if (ErrStat>=AbortErrLev) then
call cleanup()
return
end if
dYdu = dYdu + matmul(m%lin_C, RelState_x)
end if

end if ! CompAeroMaps

END IF

Expand Down Expand Up @@ -6572,16 +6590,19 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,


index = 1
FieldMask = .false.
FieldMask(MASKID_TranslationDisp) = .true.
FieldMask(MASKID_Orientation) = .true.
FieldMask(MASKID_TranslationVel) = .true.
FieldMask(MASKID_RotationVel) = .true.
FieldMask(MASKID_TranslationAcc) = .true.
FieldMask(MASKID_RotationAcc) = .true.
call PackMotionMesh(u%RootMotion, u_op, index, FieldMask=FieldMask)
if (.not. p%CompAeroMaps) then
FieldMask = .false.
FieldMask(MASKID_TranslationDisp) = .true.
FieldMask(MASKID_Orientation) = .true.
FieldMask(MASKID_TranslationVel) = .true.
FieldMask(MASKID_RotationVel) = .true.
FieldMask(MASKID_TranslationAcc) = .true.
FieldMask(MASKID_RotationAcc) = .true.
call PackMotionMesh(u%RootMotion, u_op, index, FieldMask=FieldMask)

call PackLoadMesh(u%PointLoad, u_op, index)
call PackLoadMesh(u%PointLoad, u_op, index)
end if

call PackLoadMesh(u%DistrLoad, u_op, index)

END IF
Expand All @@ -6606,22 +6627,28 @@ SUBROUTINE BD_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op,
if (ReturnTrimOP) y_op = 0.0_ReKi ! initialize in case we are returning packed orientations and don't fill the entire array

index = 1
call PackLoadMesh(y%ReactionForce, y_op, index)

FieldMask = .false.
FieldMask(MASKID_TranslationDisp) = .true.
FieldMask(MASKID_Orientation) = .true.
FieldMask(MASKID_TranslationVel) = .true.
FieldMask(MASKID_RotationVel) = .true.
FieldMask(MASKID_TranslationAcc) = .true.
FieldMask(MASKID_RotationAcc) = .true.

if (.not. p%CompAeroMaps) then

call PackLoadMesh(y%ReactionForce, y_op, index)

FieldMask(MASKID_RotationVel) = .true.
FieldMask(MASKID_TranslationAcc) = .true.
FieldMask(MASKID_RotationAcc) = .true.
end if
call PackMotionMesh(y%BldMotion, y_op, index, FieldMask=FieldMask, TrimOP=ReturnTrimOP)

index = index - 1
do i=1,p%NumOuts + p%BldNd_TotNumOuts
y_op(i+index) = y%WriteOutput(i)
end do

if (.not. p%CompAeroMaps) then
index = index - 1
do i=1,p%NumOuts + p%BldNd_TotNumOuts
y_op(i+index) = y%WriteOutput(i)
end do
end if


END IF

Expand Down
6 changes: 3 additions & 3 deletions modules/beamdyn/src/BeamDyn_BldNdOuts_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -322,7 +322,7 @@ SUBROUTINE Calc_WriteBldNdOutput( p, m, y, ErrStat, ErrMsg )

! Set the root rotation DCM relative to the reference.
! NOTE: the orientations used in this routine are DCM's. These are directly from the mesh.
call LAPACK_DGEMM('T', 'N', 1.0_BDKi, m%u2%RootMotion%Orientation(:,:,1), m%u2%RootMotion%RefOrientation(:,:,1), 0.0_BDKi, RootRelOrient, ErrStat2, ErrMsg2 )
call LAPACK_GEMM('T', 'N', 1.0_BDKi, m%u2%RootMotion%Orientation(:,:,1), m%u2%RootMotion%RefOrientation(:,:,1), 0.0_BDKi, RootRelOrient, ErrStat2, ErrMsg2 )


! Loop over the channel sets
Expand Down Expand Up @@ -406,8 +406,8 @@ SUBROUTINE Calc_WriteBldNdOutput( p, m, y, ErrStat, ErrMsg )
!-------------------------
!FIXME: we are not trapping errors here. Do we need to?
! Sectional angular/rotational deflection Wiener-Milenkovic parameter (relative to the undeflected orientation) expressed in r
call LAPACK_DGEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,idx_node), RootRelOrient, 0.0_BDKi, Tmp33b, ErrStat2, ErrMsg2 )
call LAPACK_DGEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,idx_node), Tmp33b, 0.0_BDKi, Tmp33a, ErrStat2, ErrMsg2 )
call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,idx_node), RootRelOrient, 0.0_BDKi, Tmp33b, ErrStat2, ErrMsg2 )
call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,idx_node), Tmp33b, 0.0_BDKi, Tmp33a, ErrStat2, ErrMsg2 )
call BD_CrvExtractCrv(Tmp33a,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the node's angular/rotational defelctions
WM_ParamRD = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2) ! Rotate the parameters to the correct coordinate system for output

Expand Down
Loading