diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index 0efcf60232..b1a6eb931e 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -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)) @@ -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 diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index d2354776ef..21a9d52c0a 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/modules/beamdyn/src/BeamDyn_BldNdOuts_IO.f90 b/modules/beamdyn/src/BeamDyn_BldNdOuts_IO.f90 index 7cc98b7a03..1c420d1fc7 100644 --- a/modules/beamdyn/src/BeamDyn_BldNdOuts_IO.f90 +++ b/modules/beamdyn/src/BeamDyn_BldNdOuts_IO.f90 @@ -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 @@ -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 diff --git a/modules/beamdyn/src/BeamDyn_IO.f90 b/modules/beamdyn/src/BeamDyn_IO.f90 index af3548c49e..43c284d8f3 100644 --- a/modules/beamdyn/src/BeamDyn_IO.f90 +++ b/modules/beamdyn/src/BeamDyn_IO.f90 @@ -1006,11 +1006,13 @@ SUBROUTINE BD_ReadPrimaryFile(InputFile,InputFileData,OutFileRoot,UnEc,ErrStat,E ! OutList - List of user-requested output channels at each node(-): CALL ReadOutputList ( UnIn, InputFile, InputFileData%BldNd_OutList, InputFileData%BldNd_NumOuts, 'BldNd_OutList', "List of user-requested output channels", ErrStat2, ErrMsg2, UnEc ) ! Routine in NWTC Subroutine Library - IF ( ErrStat2 >= AbortErrLev ) THEN + IF ( ErrStat2 >= AbortErrLev .and. InputFileData%BldNd_NumOuts < 1) THEN InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoBldNdOuts) ) CALL Cleanup() RETURN + ELSE + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ENDIF !---------------------- END OF FILE ----------------------------------------- @@ -1684,7 +1686,7 @@ SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg, CalcWriteOutput ! compute the root relative orientation, RootRelOrient, which is used in several calculations below ! RootRelOrient = matmul( transpose(m%u2%RootMotion%Orientation(:,:,1)), m%u2%RootMotion%RefOrientation(:,:,1)) - 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 ) !------------------------------------ ! Tip translational deflection (relative to the undeflected position) expressed in r @@ -1699,8 +1701,8 @@ SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg, CalcWriteOutput !------------------------- ! Tip 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(:,:,y%BldMotion%NNodes), RootRelOrient, 0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 ) - call LAPACK_DGEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,y%BldMotion%NNodes), temp33_2, 0.0_BDKi, temp33, ErrStat2, ErrMsg2 ) + call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,y%BldMotion%NNodes), RootRelOrient, 0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 ) + call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,y%BldMotion%NNodes), temp33_2, 0.0_BDKi, temp33, ErrStat2, ErrMsg2 ) call BD_CrvExtractCrv(temp33,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the tip angular/rotational defelctions CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -1785,8 +1787,8 @@ SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg, CalcWriteOutput !------------------------- ! 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(:,:,j_BldMotion), RootRelOrient, 0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 ) - call LAPACK_DGEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,j_BldMotion), temp33_2, 0.0_BDKi, temp33, ErrStat2, ErrMsg2 ) + call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,j_BldMotion), RootRelOrient, 0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 ) + call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation( :,:,j_BldMotion), temp33_2, 0.0_BDKi, temp33, ErrStat2, ErrMsg2 ) call BD_CrvExtractCrv(temp33,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the node's angular/rotational defelctions CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return @@ -2088,10 +2090,14 @@ SUBROUTINE Init_Jacobian( p, u, y, m, InitOut, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return ! determine how many inputs there are in the Jacobians - nu = u%RootMotion%NNodes * 18 & ! 3 Translation Displacements + 3 orientations + 6 velocities (rotation+translation) + 6 accelerations at each node - + u%PointLoad%NNodes * 6 & ! 3 forces + 3 moments at each node - + u%DistrLoad%NNodes * 6 ! 3 forces + 3 moments at each node - + if (p%CompAeroMaps) then + nu = u%DistrLoad%NNodes * 6 ! 3 forces + 3 moments at each node + else + nu = u%RootMotion%NNodes * 18 & ! 3 Translation Displacements + 3 orientations + 6 velocities (rotation+translation) + 6 accelerations at each node + + u%PointLoad%NNodes * 6 & ! 3 forces + 3 moments at each node + + u%DistrLoad%NNodes * 6 ! 3 forces + 3 moments at each node + end if + ! all other inputs (e.g., hub motion) ignored !............................ @@ -2116,29 +2122,31 @@ SUBROUTINE Init_Jacobian( p, u, y, m, InitOut, ErrStat, ErrMsg) !Module/Mesh/Field: u%RootMotion%RotationVel = 4; !Module/Mesh/Field: u%RootMotion%TranslationAcc = 5; !Module/Mesh/Field: u%RootMotion%RotationAcc = 6; - do i_meshField = 1,6 - do i=1,u%RootMotion%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField - p%Jac_u_indx(index,2) = j !component index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do + if (.not. p%CompAeroMaps) then + do i_meshField = 1,6 + do i=1,u%RootMotion%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField + p%Jac_u_indx(index,2) = j !component index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do - !Module/Mesh/Field: u%PointLoad%Force = 7; - !Module/Mesh/Field: u%PointLoad%Moment = 8; - do i_meshField = 7,8 - do i=1,u%PointLoad%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField - p%Jac_u_indx(index,2) = j !component index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do + !Module/Mesh/Field: u%PointLoad%Force = 7; + !Module/Mesh/Field: u%PointLoad%Moment = 8; + do i_meshField = 7,8 + do i=1,u%PointLoad%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField + p%Jac_u_indx(index,2) = j !component index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do + end if !Module/Mesh/Field: u%DistrLoad%Force = 9; !Module/Mesh/Field: u%DistrLoad%Moment = 10; @@ -2194,10 +2202,12 @@ SUBROUTINE Init_Jacobian( p, u, y, m, InitOut, ErrStat, ErrMsg) InitOut%RotFrame_u = .false. ! every input is on a mesh, which stores values in the global (not rotating) frame index = 1 - call PackMotionMesh_Names(u%RootMotion, 'RootMotion', InitOut%LinNames_u, index) ! all 6 motion fields - InitOut%IsLoad_u(1:index-1) = .false. ! the RootMotion inputs are not loads - InitOut%IsLoad_u(index:) = .true. ! the remaining inputs are loads - call PackLoadMesh_Names( u%PointLoad, 'PointLoad', InitOut%LinNames_u, index) + InitOut%IsLoad_u = .true. ! initialize all inputs as loads, and overwrite for the RootMotion mesh, below: + if (.not. p%CompAeroMaps) then + call PackMotionMesh_Names(u%RootMotion, 'RootMotion', InitOut%LinNames_u, index) ! all 6 motion fields + InitOut%IsLoad_u(1:index-1) = .false. ! the RootMotion inputs are not loads + call PackLoadMesh_Names( u%PointLoad, 'PointLoad', InitOut%LinNames_u, index) + end if call PackLoadMesh_Names( u%DistrLoad, 'DistrLoad', InitOut%LinNames_u, index) @@ -2224,16 +2234,20 @@ SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) CHARACTER(ChanLen) :: ChannelName LOGICAL :: isRotating + LOGICAL :: BladeMask(FIELDMASK_SIZE) ! flags to determine if this field is part of the packing ErrStat = ErrID_None ErrMsg = "" + if (p%CompAeroMaps) then + p%Jac_ny = y%BldMotion%NNodes * 12 ! 6 displacements (translation, rotation) + 6 velocities + else - ! determine how many outputs there are in the Jacobians - p%Jac_ny = y%ReactionForce%NNodes * 6 & ! 3 forces + 3 moments at each node - + y%BldMotion%NNodes * 18 & ! 6 displacements (translation, rotation) + 6 velocities + 6 accelerations at each node - + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values - + ! determine how many outputs there are in the Jacobians + p%Jac_ny = y%ReactionForce%NNodes * 6 & ! 3 forces + 3 moments at each node + + y%BldMotion%NNodes * 18 & ! 6 displacements (translation, rotation) + 6 velocities + 6 accelerations at each node + + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values + end if ! get the names of the linearized outputs: call AllocAry(InitOut%LinNames_y, p%Jac_ny,'LinNames_y',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) @@ -2244,50 +2258,58 @@ SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) InitOut%RotFrame_y = .false. ! need to set all the values in the global system to .false index_next = 1 - call PackLoadMesh_Names( y%ReactionForce, 'Reaction force', InitOut%LinNames_y, index_next) - call PackMotionMesh_Names(y%BldMotion, 'Blade motion', InitOut%LinNames_y, index_next) + if (p%CompAeroMaps) then + BladeMask = .true. ! default is all the fields + BladeMask(MASKID_TRANSLATIONACC) = .false. + BladeMask(MASKID_ROTATIONACC) = .false. + + call PackMotionMesh_Names(y%BldMotion, 'Blade motion', InitOut%LinNames_y, index_next, FieldMask=BladeMask) + else + call PackLoadMesh_Names( y%ReactionForce, 'Reaction force', InitOut%LinNames_y, index_next) + call PackMotionMesh_Names(y%BldMotion, 'Blade motion', InitOut%LinNames_y, index_next) - do i=1,p%NumOuts + p%BldNd_TotNumOuts - InitOut%LinNames_y(i+index_next-1) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) - end do + do i=1,p%NumOuts + p%BldNd_TotNumOuts + InitOut%LinNames_y(i+index_next-1) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) + end do - AllOut = .true. ! all output values except those specifically in the global system are in the rotating system - AllOut(TipTVXg) = .false. - AllOut(TipTVYg) = .false. - AllOut(TipTVZg) = .false. - AllOut(TipRVXg) = .false. - AllOut(TipRVYg) = .false. - AllOut(TipRVZg) = .false. + AllOut = .true. ! all output values except those specifically in the global system are in the rotating system + AllOut(TipTVXg) = .false. + AllOut(TipTVYg) = .false. + AllOut(TipTVZg) = .false. + AllOut(TipRVXg) = .false. + AllOut(TipRVYg) = .false. + AllOut(TipRVZg) = .false. - do j=1,9 - do i=1,3 !x,y,z - AllOut(NTVg(j,i)) = .false. - AllOut(NRVg(j,i)) = .false. + do j=1,9 + do i=1,3 !x,y,z + AllOut(NTVg(j,i)) = .false. + AllOut(NRVg(j,i)) = .false. + end do end do - end do - do i=1,p%NumOuts - if (p%OutParam(i)%Indx == 0 ) then - InitOut%RotFrame_y(i+index_next-1) = .false. - else - InitOut%RotFrame_y(i+index_next-1) = AllOut( p%OutParam(i)%Indx ) - end if - end do + do i=1,p%NumOuts + if (p%OutParam(i)%Indx == 0 ) then + InitOut%RotFrame_y(i+index_next-1) = .false. + else + InitOut%RotFrame_y(i+index_next-1) = AllOut( p%OutParam(i)%Indx ) + end if + end do - ! set outputs for all nodes out: - index_next = index_next + p%NumOuts - DO i=1,p%BldNd_NumOuts - ChannelName = p%BldNd_OutParam(i)%Name - call Conv2UC(ChannelName) - if ( ChannelName( LEN_TRIM(ChannelName):LEN_TRIM(ChannelName) ) == 'G') then ! channel is in global coordinate system - isRotating = .false. - else - isRotating = .true. - end if - InitOut%RotFrame_y(index_next : index_next+size(p%BldNd_BlOutNd)-1 ) = isRotating - index_next = index_next + size(p%BldNd_BlOutNd) - ENDDO + ! set outputs for all nodes out: + index_next = index_next + p%NumOuts + DO i=1,p%BldNd_NumOuts + ChannelName = p%BldNd_OutParam(i)%Name + call Conv2UC(ChannelName) + if ( ChannelName( LEN_TRIM(ChannelName):LEN_TRIM(ChannelName) ) == 'G') then ! channel is in global coordinate system + isRotating = .false. + else + isRotating = .true. + end if + InitOut%RotFrame_y(index_next : index_next+size(p%BldNd_BlOutNd)-1 ) = isRotating + index_next = index_next + size(p%BldNd_BlOutNd) + ENDDO + end if END SUBROUTINE Init_Jacobian_y @@ -2438,14 +2460,22 @@ SUBROUTINE Compute_dY(p, y_p, y_m, delta, dY) ! local variables: INTEGER(IntKi) :: i ! loop over outputs INTEGER(IntKi) :: indx_first ! index indicating next value of dY to be filled - + LOGICAL :: Mask(FIELDMASK_SIZE) ! flags to determine if this field is part of the packing + indx_first = 1 - call PackLoadMesh_dY( y_p%ReactionForce, y_m%ReactionForce, dY, indx_first) - call PackMotionMesh_dY(y_p%BldMotion, y_m%BldMotion, dY, indx_first) ! all 6 motion fields + if (p%CompAeroMaps) then + Mask = .true. + Mask(MASKID_TRANSLATIONACC) = .false. + Mask(MASKID_ROTATIONACC) = .false. + call PackMotionMesh_dY(y_p%BldMotion, y_m%BldMotion, dY, indx_first, FieldMask=Mask) ! 4 motion fields + else + call PackLoadMesh_dY( y_p%ReactionForce, y_m%ReactionForce, dY, indx_first) + call PackMotionMesh_dY(y_p%BldMotion, y_m%BldMotion, dY, indx_first) ! all 6 motion fields - do i=1,p%NumOuts + p%BldNd_TotNumOuts - dY(i+indx_first-1) = y_p%WriteOutput(i) - y_m%WriteOutput(i) - end do + do i=1,p%NumOuts + p%BldNd_TotNumOuts + dY(i+indx_first-1) = y_p%WriteOutput(i) - y_m%WriteOutput(i) + end do + end if dY = dY / (2.0_R8Ki*delta) diff --git a/modules/beamdyn/src/BeamDyn_Types.f90 b/modules/beamdyn/src/BeamDyn_Types.f90 index 2e38e89c90..196f71caae 100644 --- a/modules/beamdyn/src/BeamDyn_Types.f90 +++ b/modules/beamdyn/src/BeamDyn_Types.f90 @@ -53,6 +53,7 @@ MODULE BeamDyn_Types REAL(R8Ki) , DIMENSION(1:3,1:3) :: HubRot !< Initial Hub direction cosine matrix [-] LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-] LOGICAL :: DynamicSolve = .TRUE. !< Use dynamic solve option. Set to False for static solving (handled by glue code or driver code). [-] + LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false) [-] END TYPE BD_InitInputType ! ======================= ! ========= BD_InitOutputType ======= @@ -237,6 +238,7 @@ MODULE BeamDyn_Types INTEGER(IntKi) :: Jac_nx !< half the number of continuous states in jacobian matrix [-] LOGICAL :: RotStates !< Orient states in rotating frame during linearization? (flag) [-] LOGICAL :: RelStates !< Define states relative to root motion during linearization? (flag) [-] + LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false) [-] END TYPE BD_ParameterType ! ======================= ! ========= BD_InputType ======= @@ -364,6 +366,7 @@ SUBROUTINE BD_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrSt DstInitInputData%HubRot = SrcInitInputData%HubRot DstInitInputData%Linearize = SrcInitInputData%Linearize DstInitInputData%DynamicSolve = SrcInitInputData%DynamicSolve + DstInitInputData%CompAeroMaps = SrcInitInputData%CompAeroMaps END SUBROUTINE BD_CopyInitInput SUBROUTINE BD_DestroyInitInput( InitInputData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -436,6 +439,7 @@ SUBROUTINE BD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Db_BufSz = Db_BufSz + SIZE(InData%HubRot) ! HubRot Int_BufSz = Int_BufSz + 1 ! Linearize Int_BufSz = Int_BufSz + 1 ! DynamicSolve + Int_BufSz = Int_BufSz + 1 ! CompAeroMaps IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -513,6 +517,8 @@ SUBROUTINE BD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Int_Xferred = Int_Xferred + 1 IntKiBuf(Int_Xferred) = TRANSFER(InData%DynamicSolve, IntKiBuf(1)) Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1)) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE BD_PackInitInput SUBROUTINE BD_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -617,6 +623,8 @@ SUBROUTINE BD_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err Int_Xferred = Int_Xferred + 1 OutData%DynamicSolve = TRANSFER(IntKiBuf(Int_Xferred), OutData%DynamicSolve) Int_Xferred = Int_Xferred + 1 + OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE BD_UnPackInitInput SUBROUTINE BD_CopyInitOutput( SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg ) @@ -4132,6 +4140,7 @@ SUBROUTINE BD_CopyParam( SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg ) DstParamData%Jac_nx = SrcParamData%Jac_nx DstParamData%RotStates = SrcParamData%RotStates DstParamData%RelStates = SrcParamData%RelStates + DstParamData%CompAeroMaps = SrcParamData%CompAeroMaps END SUBROUTINE BD_CopyParam SUBROUTINE BD_DestroyParam( ParamData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -4531,6 +4540,7 @@ SUBROUTINE BD_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si Int_BufSz = Int_BufSz + 1 ! Jac_nx Int_BufSz = Int_BufSz + 1 ! RotStates Int_BufSz = Int_BufSz + 1 ! RelStates + Int_BufSz = Int_BufSz + 1 ! CompAeroMaps IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -5340,6 +5350,8 @@ SUBROUTINE BD_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si Int_Xferred = Int_Xferred + 1 IntKiBuf(Int_Xferred) = TRANSFER(InData%RelStates, IntKiBuf(1)) Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1)) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE BD_PackParam SUBROUTINE BD_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -6302,6 +6314,8 @@ SUBROUTINE BD_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg Int_Xferred = Int_Xferred + 1 OutData%RelStates = TRANSFER(IntKiBuf(Int_Xferred), OutData%RelStates) Int_Xferred = Int_Xferred + 1 + OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE BD_UnPackParam SUBROUTINE BD_CopyInput( SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg ) diff --git a/modules/beamdyn/src/Driver_Beam.f90 b/modules/beamdyn/src/Driver_Beam.f90 index eed92052de..d437714c9f 100644 --- a/modules/beamdyn/src/Driver_Beam.f90 +++ b/modules/beamdyn/src/Driver_Beam.f90 @@ -102,7 +102,8 @@ PROGRAM BeamDyn_Driver_Program BD_InitInput%RootName = TRIM(RootName)//'.BD' BD_InitInput%RootDisp = MATMUL(BD_InitInput%GlbPos(:),DvrData%RootRelInit) - BD_InitInput%GlbPos(:) BD_InitInput%DynamicSolve = DvrData%DynamicSolve ! QuasiStatic options handled within the BD code. - + BD_InitInput%CompAeroMaps = .false. + t_global = DvrData%t_initial n_t_final = ((DvrData%t_final - DvrData%t_initial) / dt_global ) diff --git a/modules/beamdyn/src/Registry_BeamDyn.txt b/modules/beamdyn/src/Registry_BeamDyn.txt index 768160a515..975965e002 100644 --- a/modules/beamdyn/src/Registry_BeamDyn.txt +++ b/modules/beamdyn/src/Registry_BeamDyn.txt @@ -34,6 +34,7 @@ typedef ^ InitInputType ReKi HubPos {3} - - "Initial typedef ^ InitInputType R8Ki HubRot {3}{3} - - "Initial Hub direction cosine matrix" typedef ^ InitInputType Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." - typedef ^ InitInputType Logical DynamicSolve - .TRUE. - "Use dynamic solve option. Set to False for static solving (handled by glue code or driver code)." - +typedef ^ InitInputType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false)" - # Define outputs that the initialization routine may need here: @@ -249,6 +250,7 @@ typedef ^ ParameterType Integer Jac_ny - typedef ^ ParameterType Integer Jac_nx - - - "half the number of continuous states in jacobian matrix" - typedef ^ ParameterType logical RotStates - - - "Orient states in rotating frame during linearization? (flag)" - typedef ^ ParameterType Logical RelStates - - - "Define states relative to root motion during linearization? (flag)" - +typedef ^ ParameterType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false)" - # ..... Inputs diff --git a/modules/elastodyn/src/ElastoDyn.f90 b/modules/elastodyn/src/ElastoDyn.f90 index 257f40b673..36d7222ce6 100644 --- a/modules/elastodyn/src/ElastoDyn.f90 +++ b/modules/elastodyn/src/ElastoDyn.f90 @@ -123,7 +123,7 @@ SUBROUTINE ED_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut p%UseAD14 = LEN_TRIM(InitInp%ADInputFile) > 0 ! if we're using AD14, we need to use the AD14 input files p%RootName = InitInp%RootName ! FAST already adds '.ED' to the root name - + p%CompAeroMaps = InitInp%CompAeroMaps p%Gravity = InitInp%Gravity CALL ED_ReadInput( InitInp%InputFile, InitInp%ADInputFile, InputFileData, p%BD4Blades, Interval, p%RootName, ErrStat2, ErrMsg2 ) @@ -148,6 +148,54 @@ SUBROUTINE ED_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut END IF + IF (p%CompAeroMaps) THEN + InputFileData%DT = Interval + p%Gravity = 0.0_ReKi + + ! DEGREES OF FREEDOM + InputFileData%TeetDOF = .false. + InputFileData%DrTrDOF = .false. + InputFileData%GenDOF = .false. + InputFileData%YawDOF = .false. + InputFileData%TwFADOF1 = .false. + InputFileData%TwFADOF2 = .false. + InputFileData%TwSSDOF1 = .false. + InputFileData%TwSSDOF2 = .false. + InputFileData%PtfmSgDOF= .false. + InputFileData%PtfmSwDOF= .false. + InputFileData%PtfmHvDOF= .false. + InputFileData%PtfmRDOF = .false. + InputFileData%PtfmPDOF = .false. + InputFileData%PtfmYDOF = .false. + + ! INITIAL CONDITIONS + InputFileData%RotSpeed = InitInp%RotSpeed + InputFileData%OoPDefl = 0.0_ReKi + InputFileData%IPDefl = 0.0_ReKi + InputFileData%BlPitch(1) = 0.0_ReKi + InputFileData%BlPitch(2) = 0.0_ReKi + InputFileData%BlPitch(3) = 0.0_ReKi + InputFileData%TeetDefl = 0.0_ReKi + InputFileData%Azimuth = 0.0_ReKi + InputFileData%NacYaw = 0.0_ReKi + InputFileData%TTDspFA = 0.0_ReKi + InputFileData%TTDspSS = 0.0_ReKi + InputFileData%PtfmSurge = 0.0_ReKi + InputFileData%PtfmSway = 0.0_ReKi + InputFileData%PtfmHeave = 0.0_ReKi + InputFileData%PtfmRoll = 0.0_ReKi + InputFileData%PtfmPitch = 0.0_ReKi + InputFileData%PtfmYaw = 0.0_ReKi + + ! TURBINE CONFIGURATION + ! CHECK THAT precone is same for all blades??? + InputFileData%ShftTilt = 0.0_ReKi + + ! CHECK THAT BldFile is same for all blades??? + + InputFileData%TeetMod = 0 + + END IF CALL ED_ValidateInput( InputFileData, p%BD4Blades, InitInp%Linearize, InitInp%MHK, ErrStat2, ErrMsg2 ) CALL CheckError( ErrStat2, ErrMsg2 ) @@ -257,6 +305,7 @@ SUBROUTINE ED_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut InitOut%HubRad = p%HubRad InitOut%RotSpeed = p%RotSpeed InitOut%isFixed_GenDOF = .not. InputFileData%GenDOF + InitOut%GearBox_index = DOF_GeAz ! for steady-state solver changing rotor speed if (.not. p%BD4Blades) then @@ -286,7 +335,7 @@ SUBROUTINE ED_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! set up data needed for linearization analysis !............................................................................................ - if (InitInp%Linearize) then + if (InitInp%Linearize .or. p%CompAeroMaps) then call ED_Init_Jacobian(p, u, y, InitOut, ErrStat2, ErrMsg2) call CheckError( ErrStat2, ErrMsg2 ) if (ErrStat >= AbortErrLev) return @@ -3869,7 +3918,7 @@ SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg ) INTEGER :: I ! Generic loop-counting index INTEGER :: J ! Generic loop-counting index INTEGER :: INDX ! Index for valid arrays - INTEGER :: startIndx ! Index for BeamDyn + INTEGER :: startIndx ! Index for using BeamDyn for Blades LOGICAL :: CheckOutListAgain ! Flag used to determine if output parameter starting with "M" is valid (or the negative of another parameter) LOGICAL :: InvalidOutput(0:MaxOutPts) ! This array determines if the output channel is valid for this configuration @@ -4490,50 +4539,20 @@ SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg ) DO I = 1,p%NumOuts p%OutParam(I)%Name = OutList(I) - OutListTmp = OutList(I) - - ! Reverse the sign (+/-) of the output channel if the user prefixed the - ! channel name with a "-", "_", "m", or "M" character indicating "minus". - - - CheckOutListAgain = .FALSE. - - IF ( INDEX( "-_", OutListTmp(1:1) ) > 0 ) THEN - p%OutParam(I)%SignM = -1 ! ex, "-TipDxc1" causes the sign of TipDxc1 to be switched. - OutListTmp = OutListTmp(2:) - ELSE IF ( INDEX( "mM", OutListTmp(1:1) ) > 0 ) THEN ! We'll assume this is a variable name for now, (if not, we will check later if OutListTmp(2:) is also a variable name) - CheckOutListAgain = .TRUE. - p%OutParam(I)%SignM = 1 - ELSE - p%OutParam(I)%SignM = 1 - END IF - - CALL Conv2UC( OutListTmp ) ! Convert OutListTmp to upper case - - - Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry ) - - - ! If it started with an "M" (CheckOutListAgain) we didn't find the value in our list (Indx < 1) - - IF ( CheckOutListAgain .AND. Indx < 1 ) THEN ! Let's assume that "M" really meant "minus" and then test again - p%OutParam(I)%SignM = -1 ! ex, "MTipDxc1" causes the sign of TipDxc1 to be switched. - OutListTmp = OutListTmp(2:) - - Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry ) - END IF + Indx = FindValidChannelIndx(OutList(I), ValidParamAry, p%OutParam(I)%SignM) IF ( Indx > 0 ) THEN ! we found the channel name - p%OutParam(I)%Indx = ParamIndxAry(Indx) IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN ! but, it isn't valid for these settings + p%OutParam(I)%Indx = 0 ! pick any valid channel (I just picked "Time=0" here because it's universal) p%OutParam(I)%Units = "INVALID" p%OutParam(I)%SignM = 0 ELSE + p%OutParam(I)%Indx = ParamIndxAry(Indx) p%OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output END IF ELSE ! this channel isn't valid - p%OutParam(I)%Indx = Time ! pick any valid channel (I just picked "Time" here because it's universal) + p%OutParam(I)%Indx = 0 ! pick any valid channel (I just picked "Time=0" here because it's universal) p%OutParam(I)%Units = "INVALID" p%OutParam(I)%SignM = 0 ! multiply all results by zero @@ -10233,7 +10252,7 @@ SUBROUTINE ED_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrM ! allocate dYdu if necessary if (.not. allocated(dYdu)) then - call AllocAry(dYdu, p%Jac_ny, size(p%Jac_u_indx,1)+1, 'dYdu', ErrStat2, ErrMsg2) + call AllocAry(dYdu, p%Jac_ny, size(p%Jac_u_indx,1)+p%NumExtendedInputs, 'dYdu', ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) then call cleanup() @@ -10241,57 +10260,62 @@ SUBROUTINE ED_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrM end if end if - ! make a copy of outputs because we will need two for the central difference computations (with orientations) - call ED_CopyOutput( y, y_p, MESH_NEWCOPY, ErrStat2, ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - call ED_CopyOutput( y, y_m, MESH_NEWCOPY, ErrStat2, ErrMsg2) - call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - if (ErrStat>=AbortErrLev) then - call cleanup() - return - 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 ED_CopyOutput( y, y_p, MESH_NEWCOPY, ErrStat2, ErrMsg2) + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + call ED_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 u - call ED_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 ED_Perturb_u( p, i, 1, u_perturb, delta ) - - ! compute y at u_op + delta u - call ED_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 u + call ED_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 ED_Perturb_u( p, i, 1, u_perturb, delta ) + + ! compute y at u_op + delta u + call ED_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 u - call ED_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 ED_Perturb_u( p, i, -1, u_perturb, delta ) + ! get u_op - delta u + call ED_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 ED_Perturb_u( p, i, -1, u_perturb, delta ) - ! compute y at u_op - delta u - call ED_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 u + call ED_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, dYdu(:,i) ) + ! get central difference: + call Compute_dY( p, y_p, y_m, delta, dYdu(:,i) ) - end do - - ! now do the extended input: sum the p%NumBl blade pitch columns - dYdu(:,size(p%Jac_u_indx,1)+1) = dYdu(:,size(p%Jac_u_indx,1)-p%NumBl-1) ! last NumBl+2 columns are: GenTrq, YawMom, and BlPitchCom - do i=2,p%NumBl - dYdu(:,size(p%Jac_u_indx,1)+1) = dYdu(:,size(p%Jac_u_indx,1)+1) + dYdu(:,size(p%Jac_u_indx,1)-p%NumBl-2+i) - end do + end do + ! now do the extended input: sum the p%NumBl blade pitch columns + if (p%NumExtendedInputs > 0) then + dYdu(:,size(p%Jac_u_indx,1)+1) = dYdu(:,size(p%Jac_u_indx,1)-p%NumBl-1) ! last NumBl+2 columns are: GenTrq, YawMom, and BlPitchCom + do i=2,p%NumBl + dYdu(:,size(p%Jac_u_indx,1)+1) = dYdu(:,size(p%Jac_u_indx,1)+1) + dYdu(:,size(p%Jac_u_indx,1)-p%NumBl-2+i) + end do + end if - if (ErrStat>=AbortErrLev) then - call cleanup() - return - end if - call ED_DestroyOutput( y_p, ErrStat2, ErrMsg2 ) ! we don't need this any more - call ED_DestroyOutput( y_m, ErrStat2, ErrMsg2 ) ! we don't need this any more + if (ErrStat>=AbortErrLev) then + call cleanup() + return + end if + call ED_DestroyOutput( y_p, ErrStat2, ErrMsg2 ) ! we don't need this any more + call ED_DestroyOutput( y_m, ErrStat2, ErrMsg2 ) ! we don't need this any more + end if !CompAeroMaps END IF @@ -10302,7 +10326,7 @@ SUBROUTINE ED_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrM ! allocate dXdu if necessary if (.not. allocated(dXdu)) then - call AllocAry(dXdu, p%DOFs%NActvDOF * 2, size(p%Jac_u_indx,1)+1, 'dXdu', ErrStat2, ErrMsg2) + call AllocAry(dXdu, p%NActvDOF_Lin + p%NActvVelDOF_Lin, size(p%Jac_u_indx,1)+p%NumExtendedInputs, 'dXdu', ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) then call cleanup() @@ -10334,31 +10358,25 @@ SUBROUTINE ED_JacobianPInput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrM call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - ! get central difference: - ! we may have had an error allocating memory, so we'll check if (ErrStat>=AbortErrLev) then call cleanup() return - end if - - do j=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dXdu(j, i) = x_p%QT( p%DOFs%PS(j) ) - x_m%QT( p%DOFs%PS(j) ) - end do - do j=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dXdu(j+p%DOFs%NActvDOF, i) = x_p%QDT( p%DOFs%PS(j) ) - x_m%QDT( p%DOFs%PS(j) ) - end do - dXdu(:,i) = dXdu(:,i) / (2*delta) + end if + ! get central difference: + call Compute_dX( p, x_p, x_m, delta, dXdu(:,i) ) + end do ! now do the extended input: sum the p%NumBl blade pitch columns - dXdu(:,size(p%Jac_u_indx,1)+1) = dXdu(:,size(p%Jac_u_indx,1)-p%NumBl-1) ! last NumBl+2 columns are: GenTrq, YawMom, and BlPitchCom - do i=2,p%NumBl - dXdu(:,size(p%Jac_u_indx,1)+1) = dXdu(:,size(p%Jac_u_indx,1)+1) + dXdu(:,size(p%Jac_u_indx,1)-p%NumBl-2+i) - end do - + if (p%NumExtendedInputs > 0) then + dXdu(:,size(p%Jac_u_indx,1)+1) = dXdu(:,size(p%Jac_u_indx,1)-p%NumBl-1) ! last NumBl+2 columns are: GenTrq, YawMom, and BlPitchCom + do i=2,p%NumBl + dXdu(:,size(p%Jac_u_indx,1)+1) = dXdu(:,size(p%Jac_u_indx,1)+1) + dXdu(:,size(p%Jac_u_indx,1)-p%NumBl-2+i) + end do + end if call ED_DestroyContState( x_p, ErrStat2, ErrMsg2 ) ! we don't need this any more call ED_DestroyContState( x_m, ErrStat2, ErrMsg2 ) ! we don't need this any more @@ -10453,7 +10471,7 @@ SUBROUTINE ED_JacobianPContState( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ! allocate dYdx if necessary if (.not. allocated(dYdx)) then - call AllocAry(dYdx, p%Jac_ny, p%DOFs%NActvDOF*2, 'dYdx', ErrStat2, ErrMsg2) + call AllocAry(dYdx, p%Jac_ny, p%NActvDOF_Lin + p%NActvVelDOF_Lin, 'dYdx', ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) then call cleanup() @@ -10472,7 +10490,7 @@ SUBROUTINE ED_JacobianPContState( t, u, p, x, xd, z, OtherState, y, m, ErrStat, end if - do i=1,p%DOFs%NActvDOF*2 + do i=1,p%NActvDOF_Lin + p%NActvVelDOF_Lin ! get x_op + delta x call ED_CopyContState( x, x_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 ) @@ -10514,7 +10532,7 @@ SUBROUTINE ED_JacobianPContState( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ! allocate dXdx if necessary if (.not. allocated(dXdx)) then - call AllocAry(dXdx, p%DOFs%NActvDOF * 2, p%DOFs%NActvDOF * 2, 'dXdx', ErrStat2, ErrMsg2) + call AllocAry(dXdx, p%NActvDOF_Lin + p%NActvVelDOF_Lin, p%NActvDOF_Lin + p%NActvVelDOF_Lin, 'dXdx', ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) then call cleanup() @@ -10522,7 +10540,7 @@ SUBROUTINE ED_JacobianPContState( t, u, p, x, xd, z, OtherState, y, m, ErrStat, end if end if - do i=1,p%DOFs%NActvDOF * 2 + do i=1,p%NActvDOF_Lin + p%NActvVelDOF_Lin ! get x_op + delta x call ED_CopyContState( x, x_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 ) @@ -10544,22 +10562,17 @@ SUBROUTINE ED_JacobianPContState( t, u, p, x, xd, z, OtherState, y, m, ErrStat, call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - ! get central difference: - ! we may have had an error allocating memory, so we'll check if (ErrStat>=AbortErrLev) then call cleanup() return - end if + end if - do j=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dXdx(j, i) = x_p%QT( p%DOFs%PS(j) ) - x_m%QT( p%DOFs%PS(j) ) - end do - do j=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dXdx(j+p%DOFs%NActvDOF, i) = x_p%QDT( p%DOFs%PS(j) ) - x_m%QDT( p%DOFs%PS(j) ) - end do - dXdx(:,i) = dXdx(:,i) / (2*delta) + ! get central difference: + + call Compute_dX( p, x_p, x_m, delta, dXdx(:,i) ) + end do call ED_DestroyContState( x_p, ErrStat2, ErrMsg2 ) ! we don't need this any more @@ -10749,6 +10762,7 @@ SUBROUTINE ED_Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) CHARACTER(ErrMsgLen) :: ErrMsg2 CHARACTER(*), PARAMETER :: RoutineName = 'ED_Init_Jacobian_y' LOGICAL :: Mask(FIELDMASK_SIZE) ! flags to determine if this field is part of the packing + LOGICAL :: BladeMask(FIELDMASK_SIZE) ! flags to determine if this field is part of the packing logical, allocatable :: AllOut(:) @@ -10758,25 +10772,37 @@ SUBROUTINE ED_Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) ! determine how many outputs there are in the Jacobians - p%Jac_ny = 0 - if (allocated(y%BladeLn2Mesh)) then - do i=1,p%NumBl - p%Jac_ny = p%Jac_ny + y%BladeLn2Mesh(i)%NNodes * 18 ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node on each blade - end do - end if + p%Jac_ny = 0 + BladeMask = .true. ! default is all the fields + if (p%CompAeroMaps) then + if (allocated(y%BladeLn2Mesh)) then + do i=1,p%NumBl_Lin + p%Jac_ny = p%Jac_ny + y%BladeLn2Mesh(i)%NNodes * 12 ! 3 TranslationDisp, Orientation, TranslationVel, and RotationVel at each node on each blade (skip accelerations) + end do + end if + BladeMask(MASKID_TRANSLATIONACC) = .false. + BladeMask(MASKID_ROTATIONACC) = .false. + else - p%Jac_ny = p%Jac_ny & - + y%PlatformPtMesh%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node - + y%TowerLn2Mesh%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node - + y%HubPtMotion%NNodes * 9 & ! 3 TranslationDisp, Orientation, and RotationVel at each node - + y%NacelleMotion%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node - + 3 & ! Yaw, YawRate, and HSS_Spd - + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values + if (allocated(y%BladeLn2Mesh)) then + do i=1,p%NumBl_Lin + p%Jac_ny = p%Jac_ny + y%BladeLn2Mesh(i)%NNodes * 18 ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node on each blade + end do + end if + + p%Jac_ny = p%Jac_ny & + + y%PlatformPtMesh%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node + + y%TowerLn2Mesh%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node + + y%HubPtMotion%NNodes * 9 & ! 3 TranslationDisp, Orientation, and RotationVel at each node + + y%NacelleMotion%NNodes * 18 & ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each node + + 3 & ! Yaw, YawRate, and HSS_Spd + + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values - do i=1,p%NumBl - p%Jac_ny = p%Jac_ny + y%BladeRootMotion(i)%NNodes * 18 ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each (1) node on each blade - end do + do i=1,p%NumBl_Lin + p%Jac_ny = p%Jac_ny + y%BladeRootMotion(i)%NNodes * 18 ! 3 TranslationDisp, Orientation, TranslationVel, RotationVel, TranslationAcc, and RotationAcc at each (1) node on each blade + end do + end if !................. ! set linearization output names: @@ -10787,101 +10813,103 @@ SUBROUTINE ED_Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) InitOut%RotFrame_y = .false. ! note that meshes are in the global, not rotating frame - ! note that this Mask is for the y%HubPtMotion mesh ONLY. The others pack *all* of the motion fields - Mask = .false. - Mask(MASKID_TRANSLATIONDISP) = .true. - Mask(MASKID_ORIENTATION) = .true. - Mask(MASKID_ROTATIONVEL) = .true. index_next = 1 if (allocated(y%BladeLn2Mesh)) then index_last = index_next - do i=1,p%NumBl - call PackMotionMesh_Names(y%BladeLn2Mesh(i), 'Blade '//trim(num2lstr(i)), InitOut%LinNames_y, index_next) + do i=1,p%NumBl_Lin + call PackMotionMesh_Names(y%BladeLn2Mesh(i), 'Blade '//trim(num2lstr(i)), InitOut%LinNames_y, index_next, FieldMask=BladeMask) end do - !InitOut%RotFrame_y(index_last:index_next-1) = .true. ! values on the mesh are in global, not rotating frame end if - call PackMotionMesh_Names(y%PlatformPtMesh, 'Platform', InitOut%LinNames_y, index_next) - call PackMotionMesh_Names(y%TowerLn2Mesh, 'Tower', InitOut%LinNames_y, index_next) - call PackMotionMesh_Names(y%HubPtMotion, 'Hub', InitOut%LinNames_y, index_next, FieldMask=Mask) - index_last = index_next - do i=1,p%NumBl - call PackMotionMesh_Names(y%BladeRootMotion(i), 'Blade root '//trim(num2lstr(i)), InitOut%LinNames_y, index_next) - end do - !InitOut%RotFrame_y(index_last:index_next-1) = .true. ! values on the mesh are in global, not rotating frame - - call PackMotionMesh_Names(y%NacelleMotion, 'Nacelle', InitOut%LinNames_y, index_next) - InitOut%LinNames_y(index_next) = 'Yaw, rad'; index_next = index_next+1 - InitOut%LinNames_y(index_next) = 'YawRate, rad/s'; index_next = index_next+1 - InitOut%LinNames_y(index_next) = 'HSS_Spd, rad/s' + + if (.not. p%CompAeroMaps) then + call PackMotionMesh_Names(y%PlatformPtMesh, 'Platform', InitOut%LinNames_y, index_next) + call PackMotionMesh_Names(y%TowerLn2Mesh, 'Tower', InitOut%LinNames_y, index_next) + + ! note that this Mask is for the y%HubPtMotion mesh ONLY. The others pack *all* of the motion fields + Mask = .false. + Mask(MASKID_TRANSLATIONDISP) = .true. + Mask(MASKID_ORIENTATION) = .true. + Mask(MASKID_ROTATIONVEL) = .true. + + call PackMotionMesh_Names(y%HubPtMotion, 'Hub', InitOut%LinNames_y, index_next, FieldMask=Mask) + index_last = index_next + do i=1,p%NumBl_Lin + call PackMotionMesh_Names(y%BladeRootMotion(i), 'Blade root '//trim(num2lstr(i)), InitOut%LinNames_y, index_next) + end do + + call PackMotionMesh_Names(y%NacelleMotion, 'Nacelle', InitOut%LinNames_y, index_next) + InitOut%LinNames_y(index_next) = 'Yaw, rad'; index_next = index_next+1 + InitOut%LinNames_y(index_next) = 'YawRate, rad/s'; index_next = index_next+1 + InitOut%LinNames_y(index_next) = 'HSS_Spd, rad/s' - do i=1,p%NumOuts + p%BldNd_TotNumOuts - InitOut%LinNames_y(i+index_next) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) !trim(p%OutParam(i)%Name)//', '//p%OutParam(i)%Units - end do + do i=1,p%NumOuts + p%BldNd_TotNumOuts + InitOut%LinNames_y(i+index_next) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) !trim(p%OutParam(i)%Name)//', '//p%OutParam(i)%Units + end do - !! check for AllOuts in rotating frame - allocate( AllOut(0:MaxOutPts), STAT=ErrStat2 ) ! allocate starting at zero to account for invalid output channels - if (ErrStat2 /=0 ) then - call SetErrStat(ErrID_Info, 'error allocating temporary space for AllOut',ErrStat,ErrMsg,RoutineName) - return; - end if + !! check for AllOuts in rotating frame + allocate( AllOut(0:MaxOutPts), STAT=ErrStat2 ) ! allocate starting at zero to account for invalid output channels + if (ErrStat2 /=0 ) then + call SetErrStat(ErrID_Info, 'error allocating temporary space for AllOut',ErrStat,ErrMsg,RoutineName) + return; + end if - AllOut = .false. - do k=1,3 - AllOut(TipDxc( k)) = .true. - AllOut(TipDyc( k)) = .true. - AllOut(TipDzc( k)) = .true. - AllOut(TipDxb( k)) = .true. - AllOut(TipDyb( k)) = .true. - AllOut(TipALxb( k)) = .true. - AllOut(TipALyb( k)) = .true. - AllOut(TipALzb( k)) = .true. - AllOut(TipRDxb( k)) = .true. - AllOut(TipRDyb( k)) = .true. - AllOut(TipRDzc( k)) = .true. - AllOut(TipClrnc(k)) = .true. - AllOut(PtchPMzc(k)) = .true. - AllOut(RootFxc( k)) = .true. - AllOut(RootFyc( k)) = .true. - AllOut(RootFzc( k)) = .true. - AllOut(RootFxb( k)) = .true. - AllOut(RootFyb( k)) = .true. - AllOut(RootMxc( k)) = .true. - AllOut(RootMyc( k)) = .true. - AllOut(RootMzc( k)) = .true. - AllOut(RootMxb( k)) = .true. - AllOut(RootMyb( k)) = .true. - - do j=1,9 - AllOut(SpnALxb( j,k)) = .true. - AllOut(SpnALyb( j,k)) = .true. - AllOut(SpnALzb( j,k)) = .true. - AllOut(SpnFLxb( j,k)) = .true. - AllOut(SpnFLyb( j,k)) = .true. - AllOut(SpnFLzb( j,k)) = .true. - AllOut(SpnMLxb( j,k)) = .true. - AllOut(SpnMLyb( j,k)) = .true. - AllOut(SpnMLzb( j,k)) = .true. - AllOut(SpnTDxb( j,k)) = .true. - AllOut(SpnTDyb( j,k)) = .true. - AllOut(SpnTDzb( j,k)) = .true. - AllOut(SpnRDxb( j,k)) = .true. - AllOut(SpnRDyb( j,k)) = .true. - AllOut(SpnRDzb( j,k)) = .true. + AllOut = .false. + do k=1,3 + AllOut(TipDxc( k)) = .true. + AllOut(TipDyc( k)) = .true. + AllOut(TipDzc( k)) = .true. + AllOut(TipDxb( k)) = .true. + AllOut(TipDyb( k)) = .true. + AllOut(TipALxb( k)) = .true. + AllOut(TipALyb( k)) = .true. + AllOut(TipALzb( k)) = .true. + AllOut(TipRDxb( k)) = .true. + AllOut(TipRDyb( k)) = .true. + AllOut(TipRDzc( k)) = .true. + AllOut(TipClrnc(k)) = .true. + AllOut(PtchPMzc(k)) = .true. + AllOut(RootFxc( k)) = .true. + AllOut(RootFyc( k)) = .true. + AllOut(RootFzc( k)) = .true. + AllOut(RootFxb( k)) = .true. + AllOut(RootFyb( k)) = .true. + AllOut(RootMxc( k)) = .true. + AllOut(RootMyc( k)) = .true. + AllOut(RootMzc( k)) = .true. + AllOut(RootMxb( k)) = .true. + AllOut(RootMyb( k)) = .true. + + do j=1,9 + AllOut(SpnALxb( j,k)) = .true. + AllOut(SpnALyb( j,k)) = .true. + AllOut(SpnALzb( j,k)) = .true. + AllOut(SpnFLxb( j,k)) = .true. + AllOut(SpnFLyb( j,k)) = .true. + AllOut(SpnFLzb( j,k)) = .true. + AllOut(SpnMLxb( j,k)) = .true. + AllOut(SpnMLyb( j,k)) = .true. + AllOut(SpnMLzb( j,k)) = .true. + AllOut(SpnTDxb( j,k)) = .true. + AllOut(SpnTDyb( j,k)) = .true. + AllOut(SpnTDzb( j,k)) = .true. + AllOut(SpnRDxb( j,k)) = .true. + AllOut(SpnRDyb( j,k)) = .true. + AllOut(SpnRDzb( j,k)) = .true. + end do end do - end do - do i=1,p%NumOuts - InitOut%RotFrame_y(i+index_next) = AllOut( p%OutParam(i)%Indx ) - end do + do i=1,p%NumOuts + InitOut%RotFrame_y(i+index_next) = AllOut( p%OutParam(i)%Indx ) + end do - do i=1, p%BldNd_TotNumOuts - InitOut%RotFrame_y(i+p%NumOuts+index_next) = .true. - end do - - deallocate(AllOut) + do i=1, p%BldNd_TotNumOuts + InitOut%RotFrame_y(i+p%NumOuts+index_next) = .true. + end do + deallocate(AllOut) + end if !.not. p%CompAeroMaps END SUBROUTINE ED_Init_Jacobian_y !---------------------------------------------------------------------------------------------------------------------------------- @@ -10899,17 +10927,26 @@ SUBROUTINE ED_Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) CHARACTER(*), PARAMETER :: RoutineName = 'ED_Init_Jacobian_x' ! local variables: - INTEGER(IntKi) :: i + INTEGER(IntKi) :: i, indx ErrStat = ErrID_None ErrMsg = "" + if (p%CompAeroMaps) then + p%NActvDOF_Lin = p%DOFs%NActvDOF / p%NumBl ! we have only blade DOFs, and we are going to use only 1 of the blades + p%NActvDOF_Stride = p%NumBl + p%NActvVelDOF_Lin = 0 ! we do NOT have velocity states + else + p%NActvDOF_Lin = p%DOFs%NActvDOF + p%NActvDOF_Stride = 1 + p%NActvVelDOF_Lin = p%NActvDOF_Lin ! we have velocity states + end if ! allocate space for the row/column names and for perturbation sizes - call allocAry(p%dx, p%NDof, 'p%dx', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - CALL AllocAry(InitOut%LinNames_x, p%DOFs%NActvDOF*2, 'LinNames_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - CALL AllocAry(InitOut%RotFrame_x, p%DOFs%NActvDOF*2, 'RotFrame_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - CALL AllocAry(InitOut%DerivOrder_x, p%DOFs%NActvDOF*2, 'DerivOrder_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call allocAry(p%dx, p%NDof, 'p%dx', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + CALL AllocAry(InitOut%LinNames_x, p%NActvDOF_Lin + p%NActvVelDOF_Lin,'LinNames_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + CALL AllocAry(InitOut%RotFrame_x, p%NActvDOF_Lin + p%NActvVelDOF_Lin,'RotFrame_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + CALL AllocAry(InitOut%DerivOrder_x, p%NActvDOF_Lin + p%NActvVelDOF_Lin,'DerivOrder_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return ! All Elastodyn continuous states are max order = 2 @@ -10939,26 +10976,33 @@ SUBROUTINE ED_Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) p%dx(i) = max(p%dx(i), MinPerturb) end do - InitOut%RotFrame_x = .false. - do i=1,p%DOFs%NActvDOF - if ( p%DOFs%PS(i) >= DOF_BF(1,1) ) then - if ( p%NumBl == 2 ) then - InitOut%RotFrame_x(i) = p%DOFs%PS(i) < DOF_Teet - else - InitOut%RotFrame_x(i) = .true. ! = p%DOFs%PS(i) <= DOF_BF (MaxBl,NumBF) + if (p%CompAeroMaps) then + InitOut%RotFrame_x = .true. + else + InitOut%RotFrame_x = .false. + do i=1,p%DOFs%NActvDOF + if ( p%DOFs%PS(i) >= DOF_BF(1,1) ) then + if ( p%NumBl == 2 ) then + InitOut%RotFrame_x(i) = p%DOFs%PS(i) < DOF_Teet + else + InitOut%RotFrame_x(i) = .true. ! = p%DOFs%PS(i) <= DOF_BF (MaxBl,NumBF) + end if end if - end if - end do + end do + end if ! set linearization output names: - do i=1,p%DOFs%NActvDOF - InitOut%LinNames_x(i) = p%DOF_Desc( p%DOFs%PS(i) ) + indx = 0 + do i=1,p%DOFs%NActvDOF,p%NActvDOF_Stride + indx = indx + 1 + InitOut%LinNames_x(indx) = p%DOF_Desc( p%DOFs%PS(i) ) end do - do i=1,p%DOFs%NActvDOF - InitOut%LinNames_x(i+p%DOFs%NActvDOF) = 'First time derivative of '//trim(InitOut%LinNames_x(i))//'/s' - InitOut%RotFrame_x(i+p%DOFs%NActvDOF) = InitOut%RotFrame_x(i) - end do + + do i=1,p%NActvVelDOF_Lin + InitOut%LinNames_x(i+p%NActvDOF_Lin) = 'First time derivative of '//trim(InitOut%LinNames_x(i))//'/s' + InitOut%RotFrame_x(i+p%NActvDOF_Lin) = InitOut%RotFrame_x(i) + end do END SUBROUTINE ED_Init_Jacobian_x !---------------------------------------------------------------------------------------------------------------------------------- @@ -10983,10 +11027,15 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) REAL(R8Ki) :: ScaleLength - ErrStat = ErrID_None ErrMsg = "" - + + if (p%CompAeroMaps) then + p%NumBl_Lin = 1 + else + p%NumBl_Lin = p%NumBl + end if + call ED_Init_Jacobian_y( p, y, InitOut, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -10999,18 +11048,23 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) ! determine how many inputs there are in the Jacobians nu = 0; if (allocated(u%BladePtLoads)) then - do i=1,p%NumBl + do i=1,p%NumBl_Lin nu = nu + u%BladePtLoads(i)%NNodes * 6 ! 3 forces + 3 moments at each node on each blade end do end if - nu = nu & - + u%PlatformPtMesh%NNodes * 6 & ! 3 forces + 3 moments at each node - + u%TowerPtLoads%NNodes * 6 & ! 3 forces + 3 moments at each node - + u%HubPtLoad%NNodes * 6 & ! 3 forces + 3 moments at each node - + u%NacelleLoads%NNodes * 6 & ! 3 forces + 3 moments at each node - + p%NumBl & ! blade pitch command (BlPitchCom) - + 2 ! YawMom and GenTrq - + + if (p%CompAeroMaps) then + p%NumExtendedInputs = 0 + else + nu = nu & + + u%PlatformPtMesh%NNodes * 6 & ! 3 forces + 3 moments at each node + + u%TowerPtLoads%NNodes * 6 & ! 3 forces + 3 moments at each node + + u%HubPtLoad%NNodes * 6 & ! 3 forces + 3 moments at each node + + u%NacelleLoads%NNodes * 6 & ! 3 forces + 3 moments at each node + + p%NumBl & ! blade pitch command (BlPitchCom) + + 2 ! YawMom and GenTrq + p%NumExtendedInputs = 1 + end if ! note: all other inputs are ignored !.................... @@ -11037,7 +11091,7 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) !Module/Mesh/Field: u%BladePtLoads(2)%Moment = 4; !Module/Mesh/Field: u%BladePtLoads(3)%Force = 5; !Module/Mesh/Field: u%BladePtLoads(3)%Moment = 6; - do k=1,p%NumBl + do k=1,p%NumBl_Lin do i_meshField = 1,2 do i=1,u%BladePtLoads(k)%NNodes @@ -11054,64 +11108,66 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) end if - !if MaxBl ever changes (i.e., MaxBl /=3), we need to modify this accordingly: - do i_meshField = 7,8 - do i=1,u%PlatformPtMesh%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%PlatformPtMesh%Force = 7; u%PlatformPtMesh%Moment = 8; - p%Jac_u_indx(index,2) = j !index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do + if (.not. p%CompAeroMaps) then + !if MaxBl ever changes (i.e., MaxBl /=3), we need to modify this accordingly: + do i_meshField = 7,8 + do i=1,u%PlatformPtMesh%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%PlatformPtMesh%Force = 7; u%PlatformPtMesh%Moment = 8; + p%Jac_u_indx(index,2) = j !index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do - do i_meshField = 9,10 - do i=1,u%TowerPtLoads%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%TowerPtLoads%Force = 9; u%TowerPtLoads%Moment = 10; - p%Jac_u_indx(index,2) = j !index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do + do i_meshField = 9,10 + do i=1,u%TowerPtLoads%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%TowerPtLoads%Force = 9; u%TowerPtLoads%Moment = 10; + p%Jac_u_indx(index,2) = j !index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do - do i_meshField = 11,12 - do i=1,u%HubPtLoad%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%HubPtLoad%Force = 11; u%HubPtLoad%Moment = 12; - p%Jac_u_indx(index,2) = j !index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do - - do i_meshField = 13,14 - do i=1,u%NacelleLoads%NNodes - do j=1,3 - p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%NacelleLoads%Force = 13; u%NacelleLoads%Moment = 14; - p%Jac_u_indx(index,2) = j !index: j - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !j - end do !i - end do + do i_meshField = 11,12 + do i=1,u%HubPtLoad%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%HubPtLoad%Force = 11; u%HubPtLoad%Moment = 12; + p%Jac_u_indx(index,2) = j !index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do + + do i_meshField = 13,14 + do i=1,u%NacelleLoads%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%NacelleLoads%Force = 13; u%NacelleLoads%Moment = 14; + p%Jac_u_indx(index,2) = j !index: j + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !j + end do !i + end do - do i_meshField = 1,p%NumBl ! scalars - p%Jac_u_indx(index,1) = 15 !Module/Mesh/Field: u%BlPitchCom = 15; - p%Jac_u_indx(index,2) = 1 !index: n/a - p%Jac_u_indx(index,3) = i_meshField !Node: blade - index = index + 1 - end do + do i_meshField = 1,p%NumBl ! scalars + p%Jac_u_indx(index,1) = 15 !Module/Mesh/Field: u%BlPitchCom = 15; + p%Jac_u_indx(index,2) = 1 !index: n/a + p%Jac_u_indx(index,3) = i_meshField !Node: blade + index = index + 1 + end do - do i_meshField = 16,17 ! scalars - p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%YawMom = 16; u%GenTrq = 17; - p%Jac_u_indx(index,2) = 1 !index: j - p%Jac_u_indx(index,3) = 1 !Node: i - index = index + 1 - end do + do i_meshField = 16,17 ! scalars + p%Jac_u_indx(index,1) = i_meshField !Module/Mesh/Field: u%YawMom = 16; u%GenTrq = 17; + p%Jac_u_indx(index,2) = 1 !index: j + p%Jac_u_indx(index,3) = 1 !Node: i + index = index + 1 + end do + end if ! .not. p%CompAeroMaps !................ ! input perturbations, du: @@ -11154,9 +11210,9 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) !................ ! names of the columns, InitOut%LinNames_u: !................ - call AllocAry(InitOut%LinNames_u, nu+1, 'LinNames_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - call AllocAry(InitOut%RotFrame_u, nu+1, 'RotFrame_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - call AllocAry(InitOut%IsLoad_u, nu+1, 'IsLoad_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call AllocAry(InitOut%LinNames_u, nu+p%NumExtendedInputs, 'LinNames_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call AllocAry(InitOut%RotFrame_u, nu+p%NumExtendedInputs, 'RotFrame_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call AllocAry(InitOut%IsLoad_u, nu+p%NumExtendedInputs, 'IsLoad_u', ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return InitOut%IsLoad_u = .true. ! most of ED's inputs are loads; we will override the non-load inputs below. @@ -11164,27 +11220,29 @@ SUBROUTINE ED_Init_Jacobian( p, u, y, InitOut, ErrStat, ErrMsg) index = 1 if (allocated(u%BladePtLoads)) then index_last = index - do k=1,p%NumBl + do k=1,p%NumBl_Lin call PackLoadMesh_Names(u%BladePtLoads(k), 'Blade '//trim(num2lstr(k)), InitOut%LinNames_u, index) end do !InitOut%RotFrame_u(index_last:index-1) = .true. ! values on the mesh are in global, not rotating frame end if - call PackLoadMesh_Names(u%PlatformPtMesh, 'Platform', InitOut%LinNames_u, index) - call PackLoadMesh_Names(u%TowerPtLoads, 'Tower', InitOut%LinNames_u, index) - call PackLoadMesh_Names(u%HubPtLoad, 'Hub', InitOut%LinNames_u, index) - call PackLoadMesh_Names(u%NacelleLoads, 'Nacelle', InitOut%LinNames_u, index) + if (.not. p%CompAeroMaps) then + call PackLoadMesh_Names(u%PlatformPtMesh, 'Platform', InitOut%LinNames_u, index) + call PackLoadMesh_Names(u%TowerPtLoads, 'Tower', InitOut%LinNames_u, index) + call PackLoadMesh_Names(u%HubPtLoad, 'Hub', InitOut%LinNames_u, index) + call PackLoadMesh_Names(u%NacelleLoads, 'Nacelle', InitOut%LinNames_u, index) - do k = 1,p%NumBl ! scalars - InitOut%LinNames_u(index) = 'Blade '//trim(num2lstr(k))//' pitch command, rad' - InitOut%IsLoad_u( index) = .false. - InitOut%RotFrame_u(index) = .true. - index = index + 1 - end do + do k = 1,p%NumBl ! scalars + InitOut%LinNames_u(index) = 'Blade '//trim(num2lstr(k))//' pitch command, rad' + InitOut%IsLoad_u( index) = .false. + InitOut%RotFrame_u(index) = .true. + index = index + 1 + end do - InitOut%LinNames_u(index) = 'Yaw moment, Nm' ; index = index + 1 - InitOut%LinNames_u(index) = 'Generator torque, Nm' ; index = index + 1 - InitOut%LinNames_u(index) = 'Extended input: collective blade-pitch command, rad' - InitOut%IsLoad_u( index) = .false. + InitOut%LinNames_u(index) = 'Yaw moment, Nm' ; index = index + 1 + InitOut%LinNames_u(index) = 'Generator torque, Nm' ; index = index + 1 + InitOut%LinNames_u(index) = 'Extended input: collective blade-pitch command, rad' + InitOut%IsLoad_u( index) = .false. + end if END SUBROUTINE ED_Init_Jacobian !---------------------------------------------------------------------------------------------------------------------------------- @@ -11258,10 +11316,10 @@ END SUBROUTINE ED_Perturb_u !---------------------------------------------------------------------------------------------------------------------------------- !> This routine perturbs the nth element of the continuous state array. !! Do not change this without making sure subroutine elastodyn::ed_init_jacobian is consistant with this routine! -SUBROUTINE ED_Perturb_x( p, n, perturb_sign, x, dx ) +SUBROUTINE ED_Perturb_x( p, n_in, perturb_sign, x, dx ) TYPE(ED_ParameterType) , INTENT(IN ) :: p !< parameters - INTEGER( IntKi ) , INTENT(IN ) :: n !< number of array element to use + INTEGER( IntKi ) , INTENT(IN ) :: n_in !< number of array element to use INTEGER( IntKi ) , INTENT(IN ) :: perturb_sign !< +1 or -1 (value to multiply perturbation by; positive or negative difference) TYPE(ED_ContinuousStateType) , INTENT(INOUT) :: x !< perturbed ED states REAL( R8Ki ) , INTENT( OUT) :: dx !< amount that specific state was perturbed @@ -11269,14 +11327,19 @@ SUBROUTINE ED_Perturb_x( p, n, perturb_sign, x, dx ) ! local variables integer(intKi) :: indx + integer(intKi) :: n + n = (n_in - 1) * p%NActvDOF_Stride + 1 if (n > p%DOFs%NActvDOF) then + indx = p%DOFs%PS(n-p%DOFs%NActvDOF) dx = p%dx( indx ) - x%QDT( indx ) = x%QDT( indx ) + dx * perturb_sign + x%QDT( indx ) = x%QDT( indx ) + dx * perturb_sign + else + indx = p%DOFs%PS(n) dx = p%dx( indx ) @@ -11301,40 +11364,81 @@ SUBROUTINE Compute_dY(p, y_p, y_m, delta, dY) LOGICAL :: Mask(FIELDMASK_SIZE) ! flags to determine if this field is part of the packing - Mask = .false. - Mask(MASKID_TRANSLATIONDISP) = .true. - Mask(MASKID_ORIENTATION) = .true. - Mask(MASKID_ROTATIONVEL) = .true. - - - indx_first = 1 + indx_first = 1 if (allocated(y_p%BladeLn2Mesh)) then - do k=1,p%NumBl - call PackMotionMesh_dY(y_p%BladeLn2Mesh(k), y_m%BladeLn2Mesh(k), dY, indx_first) - end do + Mask = .true. + if (p%CompAeroMaps) then + Mask(MASKID_TRANSLATIONACC) = .false. + Mask(MASKID_ROTATIONACC) = .false. + end if + + do k=1,p%NumBl_Lin + call PackMotionMesh_dY(y_p%BladeLn2Mesh(k), y_m%BladeLn2Mesh(k), dY, indx_first, FieldMask=Mask) + end do end if - call PackMotionMesh_dY(y_p%PlatformPtMesh, y_m%PlatformPtMesh, dY, indx_first, UseSmlAngle=.true.) - call PackMotionMesh_dY(y_p%TowerLn2Mesh, y_m%TowerLn2Mesh, dY, indx_first, UseSmlAngle=.true.) - call PackMotionMesh_dY(y_p%HubPtMotion, y_m%HubPtMotion, dY, indx_first, FieldMask=Mask) - do k=1,p%NumBl - call PackMotionMesh_dY(y_p%BladeRootMotion(k), y_m%BladeRootMotion(k), dY, indx_first) - end do - call PackMotionMesh_dY(y_p%NacelleMotion, y_m%NacelleMotion, dY, indx_first) + if (.not. p%CompAeroMaps) then + call PackMotionMesh_dY(y_p%PlatformPtMesh, y_m%PlatformPtMesh, dY, indx_first, UseSmlAngle=.true.) + call PackMotionMesh_dY(y_p%TowerLn2Mesh, y_m%TowerLn2Mesh, dY, indx_first, UseSmlAngle=.true.) + + Mask = .false. + Mask(MASKID_TRANSLATIONDISP) = .true. + Mask(MASKID_ORIENTATION) = .true. + Mask(MASKID_ROTATIONVEL) = .true. + call PackMotionMesh_dY(y_p%HubPtMotion, y_m%HubPtMotion, dY, indx_first, FieldMask=Mask) + + do k=1,p%NumBl_Lin + call PackMotionMesh_dY(y_p%BladeRootMotion(k), y_m%BladeRootMotion(k), dY, indx_first) + end do + call PackMotionMesh_dY(y_p%NacelleMotion, y_m%NacelleMotion, dY, indx_first) - dY(indx_first) = y_p%Yaw - y_m%Yaw; indx_first = indx_first + 1 - dY(indx_first) = y_p%YawRate - y_m%YawRate; indx_first = indx_first + 1 - dY(indx_first) = y_p%HSS_Spd - y_m%HSS_Spd; indx_first = indx_first + 1 + dY(indx_first) = y_p%Yaw - y_m%Yaw; indx_first = indx_first + 1 + dY(indx_first) = y_p%YawRate - y_m%YawRate; indx_first = indx_first + 1 + dY(indx_first) = y_p%HSS_Spd - y_m%HSS_Spd; indx_first = indx_first + 1 - !indx_last = indx_first + p%NumOuts - 1 - do k=1,p%NumOuts + p%BldNd_TotNumOuts - dY(k+indx_first-1) = y_p%WriteOutput(k) - y_m%WriteOutput(k) - end do + !indx_last = indx_first + p%NumOuts - 1 + do k=1,p%NumOuts + p%BldNd_TotNumOuts + dY(k+indx_first-1) = y_p%WriteOutput(k) - y_m%WriteOutput(k) + end do + end if dY = dY / (2.0_R8Ki*delta) END SUBROUTINE Compute_dY !---------------------------------------------------------------------------------------------------------------------------------- +!> This routine uses values of two continuous state types to compute an array of differences. +!! Do not change this packing without making sure subroutine elastodyn::init_jacobian is consistant with this routine! +SUBROUTINE Compute_dX(p, x_p, x_m, delta, dX) + + TYPE(ED_ParameterType) , INTENT(IN ) :: p !< parameters + TYPE(ED_ContinuousStateType) , INTENT(IN ) :: x_p !< ED continuous states at \f$ u + \Delta_p u \f$ or \f$ x + \Delta_p x \f$ (p=plus) + TYPE(ED_ContinuousStateType) , INTENT(IN ) :: x_m !< ED continuous states at \f$ u - \Delta_m u \f$ or \f$ x - \Delta_m x \f$ (m=minus) + REAL(R8Ki) , INTENT(IN ) :: delta !< difference in inputs or states \f$ delta = \Delta u \f$ or \f$ delta = \Delta_p x \f$ + REAL(R8Ki) , INTENT(INOUT) :: dX(:) !< column of dXdu or dXdx: \f$ \frac{\partial Y}{\partial u_i} = \frac{y_p - y_m}{2 \, \Delta u}\f$ or \f$ \frac{\partial Y}{\partial x_i} = \frac{y_p - y_m}{2 \, \Delta x}\f$ + + ! local variables: + INTEGER(IntKi) :: i ! loop over blade nodes + INTEGER(IntKi) :: j ! loop over blades + INTEGER(IntKi) :: indx_first ! index indicating next value of dY to be filled + + indx_first = 0 + + if (p%NActvVelDOF_Lin > 0) then + do j=1,p%DOFs%NActvDOF, p%NActvDOF_Stride ! Loop through all active (enabled) DOFs for linearization + indx_first = indx_first + 1 + dX(indx_first) = x_p%QT( p%DOFs%PS(j) ) - x_m%QT( p%DOFs%PS(j) ) + end do + end if + + do j=1,p%DOFs%NActvDOF, p%NActvDOF_Stride ! Loop through all active (enabled) DOFs for linearization + indx_first = indx_first + 1 + dX(indx_first) = x_p%QDT( p%DOFs%PS(j) ) - x_m%QDT( p%DOFs%PS(j) ) + end do + + dX = dX / (2*delta) ! whole array operation + +END SUBROUTINE Compute_dX +!---------------------------------------------------------------------------------------------------------------------------------- !> Routine to pack the data structures representing the operating points into arrays for linearization. SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, y_op, x_op, dx_op, xd_op, z_op, NeedTrimOP ) @@ -11373,43 +11477,45 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, ErrStat = ErrID_None ErrMsg = '' - + !.................................. IF ( PRESENT( u_op ) ) THEN if (.not. allocated(u_op)) then - call AllocAry(u_op, size(p%Jac_u_indx,1)+1,'u_op',ErrStat2,ErrMsg2) ! +1 for extended input here + call AllocAry(u_op, size(p%Jac_u_indx,1)+p%NumExtendedInputs,'u_op',ErrStat2,ErrMsg2) ! +1 for extended input here call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return end if index = 1 if (allocated(u%BladePtLoads)) then - do k=1,p%NumBl + do k=1,p%NumBl_Lin call PackLoadMesh(u%BladePtLoads(k), u_op, index) end do end if - call PackLoadMesh(u%PlatformPtMesh, u_op, index) - call PackLoadMesh(u%TowerPtLoads, u_op, index) - call PackLoadMesh(u%HubPtLoad, u_op, index) - call PackLoadMesh(u%NacelleLoads, u_op, index) - - do k = 1,p%NumBl ! scalars - u_op(index) = u%BlPitchCom(k) - index = index + 1 - end do - u_op(index) = u%YawMom ; index = index + 1 - u_op(index) = u%GenTrq ; index = index + 1 - - ! extended input: - u_op(index) = u%BlPitchCom(1) - - do k = 2,p%NumBl - if (.not. EqualRealNos( u%BlPitchCom(1), u%BlPitchCom(k) ) ) then - call SetErrStat(ErrID_Info,"Operating point of collective pitch extended input is invalid because "// & - "the commanded blade pitch angles are not the same for each blade.", ErrStat, ErrMsg, RoutineName) - exit - end if - end do + if (.not. p%CompAeroMaps) then + call PackLoadMesh(u%PlatformPtMesh, u_op, index) + call PackLoadMesh(u%TowerPtLoads, u_op, index) + call PackLoadMesh(u%HubPtLoad, u_op, index) + call PackLoadMesh(u%NacelleLoads, u_op, index) + + do k = 1,p%NumBl_Lin ! scalars + u_op(index) = u%BlPitchCom(k) + index = index + 1 + end do + u_op(index) = u%YawMom ; index = index + 1 + u_op(index) = u%GenTrq ; index = index + 1 + + ! extended input: ! note this happens only if .not. p%CompAeroMaps, so p%NumExtendedInputs > 0 + u_op(index) = u%BlPitchCom(1) + + do k = 2,p%NumBl_Lin + if (.not. EqualRealNos( u%BlPitchCom(1), u%BlPitchCom(k) ) ) then + call SetErrStat(ErrID_Info,"Operating point of collective pitch extended input is invalid because "// & + "the commanded blade pitch angles are not the same for each blade.", ErrStat, ErrMsg, RoutineName) + exit + end if + end do + end if END IF @@ -11423,20 +11529,26 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, if (.not. allocated(y_op)) then ! our operating point includes DCM (orientation) matrices, not just small angles like the perturbation matrices do - ny = p%Jac_ny + y%PlatformPtMesh%NNodes * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node + if (p%CompAeroMaps) then + ny = p%Jac_ny + else + ny = p%Jac_ny + y%PlatformPtMesh%NNodes * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node + y%TowerLn2Mesh%NNodes * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node + y%HubPtMotion%NNodes * 6 & ! Jac_ny has 3 for Orientation, but we need 9 at each node + y%NacelleMotion%NNodes * 6 ! Jac_ny has 3 for Orientation, but we need 9 at each node - + + do k=1,p%NumBl_Lin + ny = ny + y%BladeRootMotion(k)%NNodes * 6 ! Jac_ny has 3 for Orientation, but we need 9 at each node on each blade + end do + + end if + if (allocated(y%BladeLn2Mesh)) then - do k=1,p%NumBl + do k=1,p%NumBl_Lin ny = ny + y%BladeLn2Mesh(k)%NNodes * 6 ! Jac_ny has 3 for Orientation, but we need 9 (at each node on each blade) end do end if - do k=1,p%NumBl - ny = ny + y%BladeRootMotion(k)%NNodes * 6 ! Jac_ny has 3 for Orientation, but we need 9 at each node on each blade - end do - + call AllocAry(y_op, ny,'y_op',ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return @@ -11445,33 +11557,44 @@ SUBROUTINE ED_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 - Mask = .false. - Mask(MASKID_TRANSLATIONDISP) = .true. - Mask(MASKID_ORIENTATION) = .true. - Mask(MASKID_ROTATIONVEL) = .true. - + if ( p%CompAeroMaps ) then + Mask = .false. + Mask(MASKID_TRANSLATIONDISP) = .true. + Mask(MASKID_ORIENTATION) = .true. + Mask(MASKID_TRANSLATIONVEL) = .true. + Mask(MASKID_ROTATIONVEL) = .true. + else + Mask = .true. + end if + index = 1 if (allocated(y%BladeLn2Mesh)) then - do k=1,p%NumBl - call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index, TrimOP=ReturnTrimOP) + do k=1,p%NumBl_Lin + call PackMotionMesh(y%BladeLn2Mesh(k), y_op, index, FieldMask=Mask, TrimOP=ReturnTrimOP) end do end if - call PackMotionMesh(y%PlatformPtMesh, y_op, index, TrimOP=ReturnTrimOP) - call PackMotionMesh(y%TowerLn2Mesh, y_op, index, TrimOP=ReturnTrimOP) - call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask, TrimOP=ReturnTrimOP) - - do k=1,p%NumBl - call PackMotionMesh(y%BladeRootMotion(k), y_op, index, TrimOP=ReturnTrimOP) - end do - call PackMotionMesh(y%NacelleMotion, y_op, index, TrimOP=ReturnTrimOP) - - y_op(index) = y%Yaw ; index = index + 1 - y_op(index) = y%YawRate ; index = index + 1 - y_op(index) = y%HSS_Spd - - do i=1,p%NumOuts + p%BldNd_TotNumOuts - y_op(i+index) = y%WriteOutput(i) - end do + if (.not. p%CompAeroMaps) then + Mask = .false. + Mask(MASKID_TRANSLATIONDISP) = .true. + Mask(MASKID_ORIENTATION) = .true. + Mask(MASKID_ROTATIONVEL) = .true. + + call PackMotionMesh(y%PlatformPtMesh, y_op, index, TrimOP=ReturnTrimOP) + call PackMotionMesh(y%TowerLn2Mesh, y_op, index, TrimOP=ReturnTrimOP) + call PackMotionMesh(y%HubPtMotion, y_op, index, FieldMask=Mask, TrimOP=ReturnTrimOP) + do k=1,p%NumBl_Lin + call PackMotionMesh(y%BladeRootMotion(k), y_op, index, TrimOP=ReturnTrimOP) + end do + call PackMotionMesh(y%NacelleMotion, y_op, index, TrimOP=ReturnTrimOP) + + y_op(index) = y%Yaw ; index = index + 1 + y_op(index) = y%YawRate ; index = index + 1 + y_op(index) = y%HSS_Spd + + do i=1,p%NumOuts + p%BldNd_TotNumOuts + y_op(i+index) = y%WriteOutput(i) + end do + end if END IF @@ -11479,17 +11602,23 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, IF ( PRESENT( x_op ) ) THEN if (.not. allocated(x_op)) then - call AllocAry(x_op, p%DOFs%NActvDOF * 2,'x_op',ErrStat2,ErrMsg2) + call AllocAry(x_op, p%NActvDOF_Lin + p%NActvVelDOF_Lin,'x_op',ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return end if - do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - x_op(i) = x%QT( p%DOFs%PS(i) ) + index = 0 + do i=1,p%DOFs%NActvDOF,p%NActvDOF_Stride ! Loop through all active (enabled) DOFs in the Jacobian + index = index + 1 + x_op(index) = x%QT( p%DOFs%PS(i) ) end do - do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - x_op(i+p%DOFs%NActvDOF) = x%QDT( p%DOFs%PS(i) ) - end do + + if (p%NActvVelDOF_Lin > 0) then ! .not. p%CompAeroMaps + do i=1,p%DOFs%NActvDOF,p%NActvDOF_Stride ! Loop through all active (enabled) DOFs in the Jacobian + index = index + 1 + x_op(index) = x%QDT( p%DOFs%PS(i) ) + end do + end if END IF @@ -11497,7 +11626,7 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, IF ( PRESENT( dx_op ) ) THEN if (.not. allocated(dx_op)) then - call AllocAry(dx_op, p%DOFs%NActvDOF * 2,'dx_op',ErrStat2,ErrMsg2) + call AllocAry(dx_op, p%NActvDOF_Lin + p%NActvVelDOF_Lin,'dx_op',ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return end if @@ -11509,12 +11638,18 @@ SUBROUTINE ED_GetOP( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, u_op, return end if - do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dx_op(i) = dx%QT( p%DOFs%PS(i) ) + index = 0 + if (p%NActvVelDOF_Lin > 0) then ! p%CompAeroMaps + do i=1,p%DOFs%NActvDOF,p%NActvDOF_Stride ! Loop through all active (enabled) DOFs in the Jacobian + index = index + 1 + dx_op(index) = dx%QT( p%DOFs%PS(i) ) + end do + end if + + do i=1,p%DOFs%NActvDOF,p%NActvDOF_Stride ! Loop through all active (enabled) DOFs in the Jacobian + index = index + 1 + dx_op(index) = dx%QDT( p%DOFs%PS(i) ) end do - do i=1,p%DOFs%NActvDOF ! Loop through all active (enabled) DOFs - dx_op(i+p%DOFs%NActvDOF) = dx%QDT( p%DOFs%PS(i) ) - end do call ED_DestroyContState( dx, ErrStat2, ErrMsg2) diff --git a/modules/elastodyn/src/ElastoDyn_AllBldNdOuts_IO.f90 b/modules/elastodyn/src/ElastoDyn_AllBldNdOuts_IO.f90 index adbbd13201..2bc7764364 100644 --- a/modules/elastodyn/src/ElastoDyn_AllBldNdOuts_IO.f90 +++ b/modules/elastodyn/src/ElastoDyn_AllBldNdOuts_IO.f90 @@ -482,9 +482,11 @@ SUBROUTINE AllBldNdOuts_SetParameters( p, InputFileData, ErrStat, ErrMsg ) ErrMsg = "" ! Check if the requested blades exist - IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) .OR. (InputFileData%BldNd_BladesOut > p%NumBl) ) THEN - CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all bladed nodes (BldNd_BladesOut) must be between 0 and "//TRIM(Num2LStr(p%NumBl))//".", ErrStat, ErrMsg, RoutineName) + IF ( (InputFileData%BldNd_BladesOut < 0_IntKi) ) THEN p%BldNd_BladesOut = 0_IntKi + ELSE IF ((InputFileData%BldNd_BladesOut > p%NumBl) ) THEN + CALL SetErrStat( ErrID_Warn, " Number of blades to output data at all blade nodes (BldNd_BladesOut) must be less than "//TRIM(Num2LStr(p%NumBl))//".", ErrStat, ErrMsg, RoutineName) + p%BldNd_BladesOut = p%NumBl ! NOTE: we are forgiving and plateau to numBlades ELSE p%BldNd_BladesOut = InputFileData%BldNd_BladesOut ENDIF @@ -503,6 +505,7 @@ SUBROUTINE AllBldNdOuts_SetParameters( p, InputFileData, ErrStat, ErrMsg ) ELSE p%BldNd_NumOuts = InputFileData%BldNd_NumOuts ENDIF + if (p%BldNd_BladesOut==0) p%BldNd_NumOuts = 0 ! Set the total number of outputs ( requested channel groups * number requested nodes * number requested blades ) p%BldNd_TotNumOuts = p%BldNodes*p%BldNd_BladesOut*p%BldNd_NumOuts !p%BldNd_NumOuts * size(p%BldNd_BlOutNd) * size(p%BldNd_BladesOut) diff --git a/modules/elastodyn/src/ElastoDyn_IO.f90 b/modules/elastodyn/src/ElastoDyn_IO.f90 index 2a27e52742..45348664e6 100644 --- a/modules/elastodyn/src/ElastoDyn_IO.f90 +++ b/modules/elastodyn/src/ElastoDyn_IO.f90 @@ -3574,12 +3574,12 @@ SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, BldFile, FurlFile, TwrFile !----------- OUTLIST ----------------------------------------------------------- ! In case there is something ill-formed in the additional nodal outputs section, we will simply ignore it and assume that this section does not exist. ErrMsg_NoAllBldNdOuts='Nodal outputs section of ElastoDyn input file not found or improperly formatted.' + InputFileData%BldNd_NumOuts = 0 ! initialize in case of error + InputFileData%BldNd_BladesOut = 0 ! initialize in case of error !----------- OUTLIST for BldNd ----------------------------------------------------------- CALL ReadCom( UnIn, InputFile, 'Section Header: OutList for Blade node channels', ErrStat2, ErrMsg2, UnEc ) IF ( ErrStat2 >= AbortErrLev ) THEN - InputFileData%BldNd_BladesOut = 0 - InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) CALL Cleanup() RETURN @@ -3592,7 +3592,6 @@ SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, BldFile, FurlFile, TwrFile CALL ReadVar( UnIn, InputFile, InputFileData%BldNd_BladesOut, 'BldNd_BladesOut', 'Which blades to output node data on.'//TRIM(Num2Lstr(I)), ErrStat2, ErrMsg2, UnEc ) IF ( ErrStat2 >= AbortErrLev ) THEN InputFileData%BldNd_BladesOut = 0 - InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) CALL Cleanup() RETURN @@ -3603,8 +3602,6 @@ SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, BldFile, FurlFile, TwrFile ! TODO: Parse this string into an array of nodes to output at (one idea is to set an array of boolean to T/F for which nodes to output). At present, we ignore it entirely. CALL ReadVar( UnIn, InputFile, InputFileData%BldNd_BlOutNd_Str, 'BldNd_BlOutNd_Str', 'Which nodes to output node data on.'//TRIM(Num2Lstr(I)), ErrStat2, ErrMsg2, UnEc ) IF ( ErrStat2 >= AbortErrLev ) THEN - InputFileData%BldNd_BladesOut = 0 - InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) CALL Cleanup() RETURN @@ -3614,8 +3611,6 @@ SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, BldFile, FurlFile, TwrFile ! Section header for outlist CALL ReadCom( UnIn, InputFile, 'Section Header: OutList', ErrStat2, ErrMsg2, UnEc ) IF ( ErrStat2 >= AbortErrLev ) THEN - InputFileData%BldNd_BladesOut = 0 - InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) CALL Cleanup() RETURN @@ -3624,12 +3619,13 @@ SUBROUTINE ReadPrimaryFile( InputFile, InputFileData, BldFile, FurlFile, TwrFile ! OutList - List of user-requested output channels at each node(-): CALL ReadOutputList ( UnIn, InputFile, InputFileData%BldNd_OutList, InputFileData%BldNd_NumOuts, 'BldNd_OutList', "List of user-requested output channels", ErrStat2, ErrMsg2, UnEc ) ! Routine in NWTC Subroutine Library - IF ( ErrStat2 >= AbortErrLev ) THEN - InputFileData%BldNd_BladesOut = 0 + IF ( ErrStat2 >= AbortErrLev .and. InputFileData%BldNd_NumOuts < 1) THEN InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) CALL Cleanup() RETURN + ELSE + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ENDIF !---------------------- END OF FILE ----------------------------------------- @@ -4131,9 +4127,7 @@ SUBROUTINE ValidatePrimaryData( InputFileData, BD4Blades, Linearize, MHK, ErrSta REAL(ReKi) :: SmallAngleLimit_Rad ! Largest input angle considered "small" (check in input file), radians INTEGER(IntKi) :: I ! loop counter INTEGER(IntKi) :: K ! blade number - INTEGER(IntKi) :: FmtWidth ! width of the field returned by the specified OutFmt - INTEGER(IntKi) :: ErrStat2 ! Temporary error status - CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary rror message + !!INTEGER(IntKi) :: FmtWidth ! width of the field returned by the specified OutFmt CHARACTER(*), PARAMETER :: RoutineName = 'ValidatePrimaryData' ! Initialize error status and angle limit defined locally (in correct units) @@ -4212,7 +4206,7 @@ SUBROUTINE ValidatePrimaryData( InputFileData, BD4Blades, Linearize, MHK, ErrSta ! Check that the gearbox efficiency is valid: IF ( ( InputFileData%GBoxEff <= 0.0_ReKi ) .OR. ( InputFileData%GBoxEff > 1.0_ReKi ) ) THEN CALL SetErrStat( ErrID_Fatal, 'GBoxEff must be in the range (0,1] (i.e., (0,100] percent).',ErrStat,ErrMsg,RoutineName ) - ENDIF + ENDIF ! warn if 2nd modes are enabled without their corresponding 1st modes diff --git a/modules/elastodyn/src/ElastoDyn_Registry.txt b/modules/elastodyn/src/ElastoDyn_Registry.txt index d18ab670cd..4d1038aebf 100644 --- a/modules/elastodyn/src/ElastoDyn_Registry.txt +++ b/modules/elastodyn/src/ElastoDyn_Registry.txt @@ -26,6 +26,8 @@ typedef ^ InitInputType CHARACTER(1024) RootName - - - "RootName for writing out typedef ^ InitInputType ReKi Gravity - - - "Gravitational acceleration" m/s^2 typedef ^ InitInputType IntKi MHK - - - "MHK turbine type switch" - typedef ^ InitInputType ReKi WtrDpth - - - "Water depth" m +typedef ^ InitInputType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if ElastoDyn is computing aero maps (true) or running a normal simulation (false)" - +typedef ^ InitInputType ReKi RotSpeed - - - "Rotor speed used when ElastoDyn is computing aero maps" "rad/s" # Define outputs from the initialization routine here: typedef ^ InitOutputType CHARACTER(ChanLen) WriteOutputHdr {:} - - "Names of the output-to-file channels" - @@ -55,6 +57,7 @@ typedef ^ InitOutputType LOGICAL RotFrame_x {:} - - "Flag that tells FAST/MBC3 i typedef ^ InitOutputType IntKi DerivOrder_x {:} - - "Integer that tells FAST/MBC3 the maximum derivative order of continuous states used in linearization" - typedef ^ InitOutputType LOGICAL RotFrame_u {:} - - "Flag that tells FAST/MBC3 if the inputs used in linearization are in the rotating frame" - typedef ^ InitOutputType LOGICAL IsLoad_u {:} - - "Flag that tells FAST if the inputs used in linearization are loads (for preconditioning matrix)" - +typedef ^ InitOutputType IntKi GearBox_index - - - "Index to gearbox rotation in state array (for steady-state calculations)" - # ..... Blade Input file data ........................................................................................................... typedef ElastoDyn/ED BladeInputData IntKi NBlInpSt - - - "Number of blade input stations" - @@ -748,6 +751,12 @@ typedef ^ ParameterType Integer Jac_u_indx {:}{:} - - "matrix to help fill/pack typedef ^ ParameterType R8Ki du {:} - - "vector that determines size of perturbation for u (inputs)" typedef ^ ParameterType R8Ki dx {:} - - "vector that determines size of perturbation for x (continuous states)" typedef ^ ParameterType Integer Jac_ny - - - "number of outputs in jacobian matrix" - +typedef ^ ParameterType Logical CompAeroMaps - - - "number of outputs in jacobian matrix" - +typedef ^ ParameterType Integer NumExtendedInputs - - - "number of extended inputs for linearization" - +typedef ^ ParameterType Integer NumBl_Lin - - - "number of blades in the jacobian" - +typedef ^ ParameterType Integer NActvVelDOF_Lin - - - "number of velocity states in the jacobian" - +typedef ^ ParameterType Integer NActvDOF_Lin - - - "number of active DOFs to use in the jacobian" - +typedef ^ ParameterType Integer NActvDOF_Stride - - - "stride for active DOFs to use in the jacobian" - # ..... Inputs .................................................................................................................... # Define inputs that are contained on the mesh here: diff --git a/modules/elastodyn/src/ElastoDyn_Types.f90 b/modules/elastodyn/src/ElastoDyn_Types.f90 index 6f19f87bdb..e4831db34a 100644 --- a/modules/elastodyn/src/ElastoDyn_Types.f90 +++ b/modules/elastodyn/src/ElastoDyn_Types.f90 @@ -44,6 +44,8 @@ MODULE ElastoDyn_Types REAL(ReKi) :: Gravity !< Gravitational acceleration [m/s^2] INTEGER(IntKi) :: MHK !< MHK turbine type switch [-] REAL(ReKi) :: WtrDpth !< Water depth [m] + LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if ElastoDyn is computing aero maps (true) or running a normal simulation (false) [-] + REAL(ReKi) :: RotSpeed !< Rotor speed used when ElastoDyn is computing aero maps [rad/s] END TYPE ED_InitInputType ! ======================= ! ========= ED_InitOutputType ======= @@ -75,6 +77,7 @@ MODULE ElastoDyn_Types INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: DerivOrder_x !< Integer that tells FAST/MBC3 the maximum derivative order of continuous states used in linearization [-] LOGICAL , DIMENSION(:), ALLOCATABLE :: RotFrame_u !< Flag that tells FAST/MBC3 if the inputs used in linearization are in the rotating frame [-] LOGICAL , DIMENSION(:), ALLOCATABLE :: IsLoad_u !< Flag that tells FAST if the inputs used in linearization are loads (for preconditioning matrix) [-] + INTEGER(IntKi) :: GearBox_index !< Index to gearbox rotation in state array (for steady-state calculations) [-] END TYPE ED_InitOutputType ! ======================= ! ========= BladeInputData ======= @@ -752,6 +755,12 @@ MODULE ElastoDyn_Types REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: du !< vector that determines size of perturbation for u (inputs) [-] REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: dx !< vector that determines size of perturbation for x (continuous states) [-] INTEGER(IntKi) :: Jac_ny !< number of outputs in jacobian matrix [-] + LOGICAL :: CompAeroMaps !< number of outputs in jacobian matrix [-] + INTEGER(IntKi) :: NumExtendedInputs !< number of extended inputs for linearization [-] + INTEGER(IntKi) :: NumBl_Lin !< number of blades in the jacobian [-] + INTEGER(IntKi) :: NActvVelDOF_Lin !< number of velocity states in the jacobian [-] + INTEGER(IntKi) :: NActvDOF_Lin !< number of active DOFs to use in the jacobian [-] + INTEGER(IntKi) :: NActvDOF_Stride !< stride for active DOFs to use in the jacobian [-] END TYPE ED_ParameterType ! ======================= ! ========= ED_InputType ======= @@ -841,6 +850,8 @@ SUBROUTINE ED_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrSt DstInitInputData%Gravity = SrcInitInputData%Gravity DstInitInputData%MHK = SrcInitInputData%MHK DstInitInputData%WtrDpth = SrcInitInputData%WtrDpth + DstInitInputData%CompAeroMaps = SrcInitInputData%CompAeroMaps + DstInitInputData%RotSpeed = SrcInitInputData%RotSpeed END SUBROUTINE ED_CopyInitInput SUBROUTINE ED_DestroyInitInput( InitInputData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -909,6 +920,8 @@ SUBROUTINE ED_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Re_BufSz = Re_BufSz + 1 ! Gravity Int_BufSz = Int_BufSz + 1 ! MHK Re_BufSz = Re_BufSz + 1 ! WtrDpth + Int_BufSz = Int_BufSz + 1 ! CompAeroMaps + Re_BufSz = Re_BufSz + 1 ! RotSpeed IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -958,6 +971,10 @@ SUBROUTINE ED_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg Int_Xferred = Int_Xferred + 1 ReKiBuf(Re_Xferred) = InData%WtrDpth Re_Xferred = Re_Xferred + 1 + IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1)) + Int_Xferred = Int_Xferred + 1 + ReKiBuf(Re_Xferred) = InData%RotSpeed + Re_Xferred = Re_Xferred + 1 END SUBROUTINE ED_PackInitInput SUBROUTINE ED_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -1013,6 +1030,10 @@ SUBROUTINE ED_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err Int_Xferred = Int_Xferred + 1 OutData%WtrDpth = ReKiBuf(Re_Xferred) Re_Xferred = Re_Xferred + 1 + OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps) + Int_Xferred = Int_Xferred + 1 + OutData%RotSpeed = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 END SUBROUTINE ED_UnPackInitInput SUBROUTINE ED_CopyInitOutput( SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg ) @@ -1203,6 +1224,7 @@ SUBROUTINE ED_CopyInitOutput( SrcInitOutputData, DstInitOutputData, CtrlCode, Er END IF DstInitOutputData%IsLoad_u = SrcInitOutputData%IsLoad_u ENDIF + DstInitOutputData%GearBox_index = SrcInitOutputData%GearBox_index END SUBROUTINE ED_CopyInitOutput SUBROUTINE ED_DestroyInitOutput( InitOutputData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -1400,6 +1422,7 @@ SUBROUTINE ED_PackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs Int_BufSz = Int_BufSz + 2*1 ! IsLoad_u upper/lower bounds for each dimension Int_BufSz = Int_BufSz + SIZE(InData%IsLoad_u) ! IsLoad_u END IF + Int_BufSz = Int_BufSz + 1 ! GearBox_index IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -1700,6 +1723,8 @@ SUBROUTINE ED_PackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs Int_Xferred = Int_Xferred + 1 END DO END IF + IntKiBuf(Int_Xferred) = InData%GearBox_index + Int_Xferred = Int_Xferred + 1 END SUBROUTINE ED_PackInitOutput SUBROUTINE ED_UnPackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -2068,6 +2093,8 @@ SUBROUTINE ED_UnPackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Er Int_Xferred = Int_Xferred + 1 END DO END IF + OutData%GearBox_index = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE ED_UnPackInitOutput SUBROUTINE ED_CopyBladeInputData( SrcBladeInputDataData, DstBladeInputDataData, CtrlCode, ErrStat, ErrMsg ) @@ -16436,6 +16463,12 @@ SUBROUTINE ED_CopyParam( SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg ) DstParamData%dx = SrcParamData%dx ENDIF DstParamData%Jac_ny = SrcParamData%Jac_ny + DstParamData%CompAeroMaps = SrcParamData%CompAeroMaps + DstParamData%NumExtendedInputs = SrcParamData%NumExtendedInputs + DstParamData%NumBl_Lin = SrcParamData%NumBl_Lin + DstParamData%NActvVelDOF_Lin = SrcParamData%NActvVelDOF_Lin + DstParamData%NActvDOF_Lin = SrcParamData%NActvDOF_Lin + DstParamData%NActvDOF_Stride = SrcParamData%NActvDOF_Stride END SUBROUTINE ED_CopyParam SUBROUTINE ED_DestroyParam( ParamData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -17183,6 +17216,12 @@ SUBROUTINE ED_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si Db_BufSz = Db_BufSz + SIZE(InData%dx) ! dx END IF Int_BufSz = Int_BufSz + 1 ! Jac_ny + Int_BufSz = Int_BufSz + 1 ! CompAeroMaps + Int_BufSz = Int_BufSz + 1 ! NumExtendedInputs + Int_BufSz = Int_BufSz + 1 ! NumBl_Lin + Int_BufSz = Int_BufSz + 1 ! NActvVelDOF_Lin + Int_BufSz = Int_BufSz + 1 ! NActvDOF_Lin + Int_BufSz = Int_BufSz + 1 ! NActvDOF_Stride IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -18743,6 +18782,18 @@ SUBROUTINE ED_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si END IF IntKiBuf(Int_Xferred) = InData%Jac_ny Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1)) + Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = InData%NumExtendedInputs + Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = InData%NumBl_Lin + Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = InData%NActvVelDOF_Lin + Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = InData%NActvDOF_Lin + Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = InData%NActvDOF_Stride + Int_Xferred = Int_Xferred + 1 END SUBROUTINE ED_PackParam SUBROUTINE ED_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -20551,6 +20602,18 @@ SUBROUTINE ED_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg END IF OutData%Jac_ny = IntKiBuf(Int_Xferred) Int_Xferred = Int_Xferred + 1 + OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps) + Int_Xferred = Int_Xferred + 1 + OutData%NumExtendedInputs = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 + OutData%NumBl_Lin = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 + OutData%NActvVelDOF_Lin = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 + OutData%NActvDOF_Lin = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 + OutData%NActvDOF_Stride = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 END SUBROUTINE ED_UnPackParam SUBROUTINE ED_CopyInput( SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg ) diff --git a/modules/openfast-library/src/FAST_Subs.f90 b/modules/openfast-library/src/FAST_Subs.f90 index 80094ca781..79c95888f8 100644 --- a/modules/openfast-library/src/FAST_Subs.f90 +++ b/modules/openfast-library/src/FAST_Subs.f90 @@ -237,6 +237,8 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, BD, SrvD, END IF Init%InData_ED%Linearize = p_FAST%Linearize + Init%InData_ED%CompAeroMaps = .FALSE. + Init%InData_ED%RotSpeed = 0.0 ! will set this in a future commit that includes the OpenFAST AeroMap/Steady-State solver Init%InData_ED%InputFile = p_FAST%EDFile IF ( p_FAST%CompAero == Module_AD14 ) THEN Init%InData_ED%ADInputFile = p_FAST%AeroFile @@ -333,6 +335,7 @@ SUBROUTINE FAST_InitializeAll( t_initial, p_FAST, y_FAST, m_FAST, ED, BD, SrvD, Init%InData_BD%DynamicSolve = .TRUE. ! FAST can only couple to BeamDyn when dynamic solve is used. Init%InData_BD%Linearize = p_FAST%Linearize + Init%InData_BD%CompAeroMaps = .FALSE. Init%InData_BD%gravity = (/ 0.0_ReKi, 0.0_ReKi, -p_FAST%Gravity /) ! "Gravitational acceleration" m/s^2 ! now initialize BeamDyn for all beams diff --git a/reg_tests/r-test b/reg_tests/r-test index 95e2af3fac..77b73ebca9 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 95e2af3fac290a5fa2c0078566950faf88cd7792 +Subproject commit 77b73ebca9c8923d9e3cc532e1eb0bfb32a1c37d