diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index c867913b1a..ad37cbbb96 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -341,7 +341,26 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut !FIXME: add handling for passing of blade files and other types of files. call ReadInputFiles( InitInp%InputFile, InputFileData, interval, p%RootName, NumBlades, AeroProjMod, UnEcho, ErrStat2, ErrMsg2 ) if (Failed()) return; - + + ! override some parameters to simplify for aero maps + ! bjj: do we put a warning here if any of these values aren't currently set this way? + if (InitInp%CompAeroMaps) then + InputFileData%DTAero = interval ! we're not using this, so set it to something "safe" + do iR = 1, nRotors + InputFileData%AFAeroMod = AFAeroMod_Steady + InputFileData%TwrPotent = TwrPotent_none + InputFileData%TwrShadow = TwrShadow_none + InputFileData%TwrAero = .false. + InputFileData%FrozenWake = .false. + !InputFileData%CavitCheck = .false. + end do + + if (InputFileData%WakeMod == WakeMod_DBEMT) then + ! these models (DBEMT and BEMT) should be the same at the first time step, so we'll simplify here + InputFileData%WakeMod = WakeMod_BEMT + end if + end if + ! Validate the inputs call ValidateInputData( InitInp, InputFileData, NumBlades, ErrStat2, ErrMsg2 ) if (Failed()) return; @@ -487,7 +506,7 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut !............................................................................................ ! Initialize Jacobian: !............................................................................................ - if (InitInp%Linearize) then + if (InitInp%Linearize .or. InitInp%CompAeroMaps) then do iR = 1, nRotors call Init_Jacobian(InputFileData%rotors(iR), p%rotors(iR), p, u%rotors(iR), y%rotors(iR), m%rotors(iR), InitOut%rotors(iR), errStat2, errMsg2) if (Failed()) return; @@ -1061,6 +1080,7 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er u%InflowOnHub = 0.0_ReKi u%InflowOnNacelle = 0.0_ReKi u%InflowOnTailFin = 0.0_ReKi + u%AvgDiskVel = 0.0_ReKi ! Meshes for motion inputs (ElastoDyn and/or BeamDyn) !................ @@ -1224,6 +1244,11 @@ subroutine Init_u( u, p, p_AD, InputFileData, MHK, WtrDpth, InitInp, errStat, er u%BladeMotion(k)%RotationVel = 0.0_ReKi u%BladeMotion(k)%TranslationAcc = 0.0_ReKi + if (p_AD%CompAeroMaps) then + do j=1,InputFileData%BladeProps(k)%NumBlNds + u%BladeMotion(k)%TranslationVel(:,j) = cross_product(u%HubMotion%RefOrientation(1,:,1)*InitInp%RotSpeed, u%BladeMotion(k)%Position(:,j)-u%HubMotion%Position(:,1)) + end do + end if end do !k=numBlades @@ -4137,6 +4162,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x call SetErrStat( ErrID_Fatal, "DTAero was changed in Init_BEMTmodule(); this is not allowed.", ErrStat2, ErrMsg2, RoutineName) !m%UseFrozenWake = .FALSE. !BJJ: set this in BEMT + if (p_AD%CompAeroMaps) p%BEMT%lin_nx = 0 ! we are going to ignore this call Cleanup() return @@ -6069,7 +6095,7 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: xd_op(:) !< values of linearized discrete states REAL(ReKi), ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: z_op(:) !< values of linearized constraint states - INTEGER(IntKi) :: index, i, j, k + INTEGER(IntKi) :: index, i, j, k, n INTEGER(IntKi) :: nu INTEGER(IntKi) :: ErrStat2 CHARACTER(ErrMsgLen) :: ErrMsg2 @@ -6084,14 +6110,19 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, ErrMsg = '' IF ( PRESENT( u_op ) ) THEN - - nu = size(p%Jac_u_indx,1) + u%TowerMotion%NNodes * 6 & ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM - + u%hubMotion%NNodes * 6 ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM - do i=1,p%NumBlades - nu = nu + u%BladeMotion(i)%NNodes * 6 & ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM - + u%BladeRootMotion(i)%NNodes * 6 ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM + nu = size(p%Jac_u_indx,1) + do i=1,p%NumBl_Lin + nu = nu + u%BladeMotion(i)%NNodes * 6 ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM end do - + + if (.not. p_AD%CompAeroMaps) then + nu = nu + u%TowerMotion%NNodes * 6 & ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM + + u%hubMotion%NNodes * 6 ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM + do i=1,p%NumBlades + nu = nu + u%BladeRootMotion(i)%NNodes * 6 ! Jac_u_indx has 3 orientation angles, but the OP needs the full 9 elements of the DCM + end do + end if + if (.not. allocated(u_op)) then call AllocAry(u_op, nu, 'u_op', ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -6100,55 +6131,70 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, index = 1 - FieldMask = .false. - FieldMask(MASKID_TRANSLATIONDISP) = .true. - FieldMask(MASKID_Orientation) = .true. - FieldMask(MASKID_TRANSLATIONVel) = .true. - call PackMotionMesh(u%TowerMotion, u_op, index, FieldMask=FieldMask) - - FieldMask(MASKID_TRANSLATIONVel) = .false. - FieldMask(MASKID_RotationVel) = .true. - call PackMotionMesh(u%HubMotion, u_op, index, FieldMask=FieldMask) - - FieldMask = .false. - FieldMask(MASKID_Orientation) = .true. - do k = 1,p%NumBlades - call PackMotionMesh(u%BladeRootMotion(k), u_op, index, FieldMask=FieldMask) - end do + if (.not. p_AD%CompAeroMaps) then + FieldMask = .false. + FieldMask(MASKID_TRANSLATIONDISP) = .true. + FieldMask(MASKID_Orientation) = .true. + FieldMask(MASKID_TRANSLATIONVel) = .true. + call PackMotionMesh(u%TowerMotion, u_op, index, FieldMask=FieldMask) + + FieldMask(MASKID_TRANSLATIONVel) = .false. + FieldMask(MASKID_RotationVel) = .true. + call PackMotionMesh(u%HubMotion, u_op, index, FieldMask=FieldMask) + + FieldMask = .false. + FieldMask(MASKID_Orientation) = .true. + do k = 1,p%NumBlades + call PackMotionMesh(u%BladeRootMotion(k), u_op, index, FieldMask=FieldMask) + end do - FieldMask(MASKID_TRANSLATIONDISP) = .true. - FieldMask(MASKID_Orientation) = .true. - FieldMask(MASKID_TRANSLATIONVel) = .true. - FieldMask(MASKID_RotationVel) = .true. - FieldMask(MASKID_TRANSLATIONAcc) = .true. - do k=1,p%NumBlades + FieldMask(MASKID_TRANSLATIONDISP) = .true. + FieldMask(MASKID_Orientation) = .true. + FieldMask(MASKID_TRANSLATIONVel) = .true. + FieldMask(MASKID_RotationVel) = .true. + FieldMask(MASKID_TRANSLATIONAcc) = .true. + else + FieldMask = .false. + FieldMask(MASKID_TRANSLATIONDISP) = .true. + FieldMask(MASKID_Orientation) = .true. + FieldMask(MASKID_TRANSLATIONVel) = .true. + end if + + do k=1,p%NumBl_Lin call PackMotionMesh(u%BladeMotion(k), u_op, index, FieldMask=FieldMask) end do - do k=1,p%NumBlades - do i=1,p%NumBlNds + if (.not. p_AD%CompAeroMaps) then + do k=1,p%NumBlades + do i=1,p%NumBlNds + do j=1,3 + u_op(index) = u%InflowOnBlade(j,i,k) + index = index + 1 + end do + end do + end do + + do i=1,p%NumTwrNds do j=1,3 - u_op(index) = u%InflowOnBlade(j,i,k) + u_op(index) = u%InflowOnTower(j,i) index = index + 1 end do end do - end do - - do i=1,p%NumTwrNds - do j=1,3 - u_op(index) = u%InflowOnTower(j,i) - index = index + 1 - end do - end do - - do k=1,p%NumBlades - do j = 1, size(u%UserProp,1) ! Number of nodes for a blade - u_op(index) = u%UserProp(j,k) - index = index + 1 + ! UserProp + do k=1,p%NumBlades + do j = 1, size(u%UserProp,1) ! Number of nodes for a blade + u_op(index) = u%UserProp(j,k) + index = index + 1 + end do end do - end do - - ! I'm not including this in the linearization yet + + ! AvgDiskVel + !do i=1,3 + ! u_op(index) = u%AvgDiskVel(i) + ! index = index + 1 + !end do + + ! I'm not including this in the linearization yet !do i=1,u%NacelleMotion%NNodes ! 1 or 0 ! do j=1,3 ! u_op(index) = u%InflowOnNacelle(j) @@ -6163,6 +6209,7 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, ! end do !end do + end if END IF IF ( PRESENT( y_op ) ) THEN @@ -6176,23 +6223,24 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, index = 1 - call PackLoadMesh(y%TowerLoad, y_op, index) - do k=1,p%NumBlades + if (.not. p_AD%CompAeroMaps) call PackLoadMesh(y%TowerLoad, y_op, index) + do k=1,p%NumBl_Lin call PackLoadMesh(y%BladeLoad(k), y_op, index) end do - index = index - 1 - do i=1,p%NumOuts + p%BldNd_TotNumOuts - y_op(i+index) = y%WriteOutput(i) - end do - + if (.not. p_AD%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 IF ( PRESENT( x_op ) ) THEN if (.not. allocated(x_op)) then - call AllocAry(x_op, p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx,'x_op',ErrStat2,ErrMsg2) + call AllocAry(x_op, p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx + p%BEMT%lin_nx,'x_op',ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return end if @@ -6219,34 +6267,32 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, end do end if - + ! UA states if (p%BEMT%UA%lin_nx>0) then - if (p%BEMT%UA%UAMod==UA_OYE) then - do j=1,p%NumBlades ! size(x%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(x%BEMT%UA%element,1) - x_op(index) = x%BEMT%UA%element(i,j)%x(4) - index = index + 1 - end do - end do - else - do j=1,p%NumBlades ! size(x%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(x%BEMT%UA%element,1) - do k=1,4 !size(x%BEMT%UA%element(i,j)%x) !linearize only first 4 states (5th is vortex) - x_op(index) = x%BEMT%UA%element(i,j)%x(k) - index = index + 1 - end do - end do - end do - endif + do n=1,p%BEMT%UA%lin_nx + i = p%BEMT%UA%lin_xIndx(n,1) + j = p%BEMT%UA%lin_xIndx(n,2) + k = p%BEMT%UA%lin_xIndx(n,3) + x_op(index) = x%BEMT%UA%element(i,j)%x(k) + + index = index + 1 + end do end if - + ! BEMT states + if (p%BEMT%lin_nx>0) then + !do k = 1,size(x%BEMT%V_w) + ! x_op(index) = x%BEMT%v_w(k) + ! index = index + 1 + !end do + end if + END IF IF ( PRESENT( dx_op ) ) THEN if (.not. allocated(dx_op)) then - call AllocAry(dx_op, p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx,'dx_op',ErrStat2,ErrMsg2) + call AllocAry(dx_op, p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx + p%BEMT%lin_nx,'dx_op',ErrStat2,ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) if (ErrStat>=AbortErrLev) return end if @@ -6281,25 +6327,25 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, end do end if - + ! UA states derivatives if (p%BEMT%UA%lin_nx>0) then - if (p%BEMT%UA%UAMod==UA_OYE) then - do j=1,p%NumBlades ! size(dxdt%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(dxdt%BEMT%UA%element,1) - dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(4) - index = index + 1 - end do - end do - else - do j=1,p%NumBlades ! size(dxdt%BEMT%UA%element,2) - do i=1,p%NumBlNds ! size(dxdt%BEMT%UA%element,1) - do k=1,4 !size(dxdt%BEMT%UA%element(i,j)%x) don't linearize 5th state - dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(k) - index = index + 1 - end do - end do - end do - endif + do n=1,p%BEMT%UA%lin_nx + i = p%BEMT%UA%lin_xIndx(n,1) + j = p%BEMT%UA%lin_xIndx(n,2) + k = p%BEMT%UA%lin_xIndx(n,3) + dx_op(index) = dxdt%BEMT%UA%element(i,j)%x(k) + + index = index + 1 + end do + end if + ! BEMT states derivatives + if (p%BEMT%lin_nx>0) then + call SetErrStat(ErrID_Fatal,'Number of lin states for bem should be zero for now.', ErrStat, ErrMsg, RoutineName) + return + !do k = 1,size(x%BEMT%V_w) + ! dx_op(index) = dxdt%BEMT%v_w(k) + ! index = index + 1 + !end do end if call AD_DestroyRotContinuousStateType( dxdt, ErrStat2, ErrMsg2) @@ -6331,9 +6377,10 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, END SUBROUTINE RotGetOP !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) +SUBROUTINE Init_Jacobian_y( p, p_AD, y, InitOut, ErrStat, ErrMsg) TYPE(RotParameterType) , INTENT(INOUT) :: p !< parameters + TYPE(AD_ParameterType) , INTENT(INOUT) :: p_AD !< parameters TYPE(RotOutputType) , INTENT(IN ) :: y !< outputs TYPE(RotInitOutputType) , INTENT(INOUT) :: InitOut !< Initialization output data (for Jacobian row/column names) @@ -6352,11 +6399,15 @@ SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) ErrMsg = "" - ! determine how many outputs there are in the Jacobians - p%Jac_ny = y%TowerLoad%NNodes * 6 & ! 3 forces + 3 moments at each node - + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values - - do k=1,p%NumBlades + ! determine how many outputs there are in the Jacobians + if (p_AD%CompAeroMaps) then + p%Jac_ny = 0 ! we skip tower and writeOutput values in the solve (note: y%TowerLoad%NNodes=0) + else + p%Jac_ny = y%TowerLoad%NNodes * 6 & ! 3 forces + 3 moments at each node + + p%NumOuts + p%BldNd_TotNumOuts ! WriteOutput values + end if + + do k=1,p%NumBl_Lin p%Jac_ny = p%Jac_ny + y%BladeLoad(k)%NNodes * 6 ! 3 forces + 3 moments at each node end do @@ -6369,97 +6420,102 @@ SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) InitOut%RotFrame_y = .false. ! default all to false, then set the true ones below indx_next = 1 - call PackLoadMesh_Names(y%TowerLoad, 'Tower', InitOut%LinNames_y, indx_next) + if (.not. p_AD%CompAeroMaps) call PackLoadMesh_Names(y%TowerLoad, 'Tower', InitOut%LinNames_y, indx_next) ! note: y%TowerLoad%NNodes=0 for aeroMaps indx_last = indx_next - do k=1,p%NumBlades + do k=1,p%NumBl_Lin call PackLoadMesh_Names(y%BladeLoad(k), 'Blade '//trim(num2lstr(k)), InitOut%LinNames_y, indx_next) end do ! InitOut%RotFrame_y(indx_last:indx_next-1) = .true. ! The mesh fields are in the global frame, so are not in the rotating frame - do i=1,p%NumOuts + p%BldNd_TotNumOuts - InitOut%LinNames_y(i+indx_next-1) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) !trim(p%OutParam(i)%Name)//', '//p%OutParam(i)%Units - end do + if (.not. p_AD%CompAeroMaps) then - - ! check for all the WriteOutput values that are functions of blade number: - 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 + do i=1,p%NumOuts + p%BldNd_TotNumOuts + InitOut%LinNames_y(i+indx_next-1) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i)) !trim(p%OutParam(i)%Name)//', '//p%OutParam(i)%Units + end do - AllOut = .false. - do k=1,3 - AllOut( BAzimuth(k)) = .true. - AllOut( BPitch (k)) = .true. - - ! AllOut( BFldFx( k)) = .true. - ! AllOut( BFldFy( k)) = .true. - ! AllOut( BFldFz( k)) = .true. - ! AllOut( BFldMx( k)) = .true. - ! AllOut( BFldMy( k)) = .true. - ! AllOut( BFldMz( k)) = .true. - - do j=1,9 - AllOut(BNVUndx(j,k)) = .true. - AllOut(BNVUndy(j,k)) = .true. - AllOut(BNVUndz(j,k)) = .true. - AllOut(BNVDisx(j,k)) = .true. - AllOut(BNVDisy(j,k)) = .true. - AllOut(BNVDisz(j,k)) = .true. - AllOut(BNSTVx (j,k)) = .true. - AllOut(BNSTVy (j,k)) = .true. - AllOut(BNSTVz (j,k)) = .true. - AllOut(BNVRel (j,k)) = .true. - AllOut(BNDynP (j,k)) = .true. - AllOut(BNRe (j,k)) = .true. - AllOut(BNM (j,k)) = .true. - AllOut(BNVIndx(j,k)) = .true. - AllOut(BNVIndy(j,k)) = .true. - AllOut(BNAxInd(j,k)) = .true. - AllOut(BNTnInd(j,k)) = .true. - AllOut(BNAlpha(j,k)) = .true. - AllOut(BNTheta(j,k)) = .true. - AllOut(BNPhi (j,k)) = .true. - AllOut(BNCurve(j,k)) = .true. - AllOut(BNCl (j,k)) = .true. - AllOut(BNCd (j,k)) = .true. - AllOut(BNCm (j,k)) = .true. - AllOut(BNCx (j,k)) = .true. - AllOut(BNCy (j,k)) = .true. - AllOut(BNCn (j,k)) = .true. - AllOut(BNCt (j,k)) = .true. - AllOut(BNFl (j,k)) = .true. - AllOut(BNFd (j,k)) = .true. - AllOut(BNMm (j,k)) = .true. - AllOut(BNFx (j,k)) = .true. - AllOut(BNFy (j,k)) = .true. - AllOut(BNFn (j,k)) = .true. - AllOut(BNFt (j,k)) = .true. - AllOut(BNClrnc(j,k)) = .true. + ! check for all the WriteOutput values that are functions of blade number: + 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( BAzimuth(k)) = .true. + AllOut( BPitch (k)) = .true. + + AllOut( BAeroFx( k)) = .true. + AllOut( BAeroFy( k)) = .true. + AllOut( BAeroFz( k)) = .true. + AllOut( BAeroMx( k)) = .true. + AllOut( BAeroMy( k)) = .true. + AllOut( BAeroMz( k)) = .true. + !AllOut( TipClrnc(k)) = .true. + + do j=1,9 + AllOut(BNVUndx(j,k)) = .true. + AllOut(BNVUndy(j,k)) = .true. + AllOut(BNVUndz(j,k)) = .true. + AllOut(BNVDisx(j,k)) = .true. + AllOut(BNVDisy(j,k)) = .true. + AllOut(BNVDisz(j,k)) = .true. + AllOut(BNSTVx (j,k)) = .true. + AllOut(BNSTVy (j,k)) = .true. + AllOut(BNSTVz (j,k)) = .true. + AllOut(BNVRel (j,k)) = .true. + AllOut(BNDynP (j,k)) = .true. + AllOut(BNRe (j,k)) = .true. + AllOut(BNM (j,k)) = .true. + AllOut(BNVIndx(j,k)) = .true. + AllOut(BNVIndy(j,k)) = .true. + AllOut(BNAxInd(j,k)) = .true. + AllOut(BNTnInd(j,k)) = .true. + AllOut(BNAlpha(j,k)) = .true. + AllOut(BNTheta(j,k)) = .true. + AllOut(BNPhi (j,k)) = .true. + AllOut(BNCurve(j,k)) = .true. + AllOut(BNCl (j,k)) = .true. + AllOut(BNCd (j,k)) = .true. + AllOut(BNCm (j,k)) = .true. + AllOut(BNCx (j,k)) = .true. + AllOut(BNCy (j,k)) = .true. + AllOut(BNCn (j,k)) = .true. + AllOut(BNCt (j,k)) = .true. + AllOut(BNFl (j,k)) = .true. + AllOut(BNFd (j,k)) = .true. + AllOut(BNMm (j,k)) = .true. + AllOut(BNFx (j,k)) = .true. + AllOut(BNFy (j,k)) = .true. + AllOut(BNFn (j,k)) = .true. + AllOut(BNFt (j,k)) = .true. + AllOut(BNClrnc(j,k)) = .true. + end do end do - end do - do i=1,p%NumOuts - InitOut%RotFrame_y(i+indx_next-1) = AllOut( p%OutParam(i)%Indx ) - end do + do i=1,p%NumOuts + InitOut%RotFrame_y(i+indx_next-1) = AllOut( p%OutParam(i)%Indx ) + end do - do i=1,p%BldNd_TotNumOuts - InitOut%RotFrame_y(i+p%NumOuts+indx_next-1) = .true. - !AbsCant, AbsToe, AbsTwist should probably be set to .false. - end do + do i=1,p%BldNd_TotNumOuts + InitOut%RotFrame_y(i+p%NumOuts+indx_next-1) = .true. + !AbsCant, AbsToe, AbsTwist should probably be set to .false. + end do - deallocate(AllOut) + deallocate(AllOut) + + end if END SUBROUTINE Init_Jacobian_y !---------------------------------------------------------------------------------------------------------------------------------- -SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) +SUBROUTINE Init_Jacobian_u( InputFileData, p, p_AD, u, InitOut, ErrStat, ErrMsg) TYPE(RotInputFile) , INTENT(IN ) :: InputFileData !< input file data (for default blade perturbation) TYPE(RotParameterType) , INTENT(INOUT) :: p !< parameters + TYPE(AD_ParameterType) , INTENT(INOUT) :: p_AD !< parameters TYPE(RotInputType) , INTENT(IN ) :: u !< inputs TYPE(RotInitOutputType) , INTENT(INOUT) :: InitOut !< Initialization output data (for Jacobian row/column names) @@ -6468,6 +6524,7 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) ! local variables: INTEGER(IntKi) :: i, j, k, index, index_last, nu, i_meshField + INTEGER(IntKi) :: NumFieldsForLinearization REAL(ReKi) :: perturb, perturb_t, perturb_b(MaxBl) LOGICAL :: FieldMask(FIELDMASK_SIZE) CHARACTER(1), PARAMETER :: UVW(3) = (/'U','V','W'/) @@ -6480,17 +6537,28 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) ! determine how many inputs there are in the Jacobians - nu = u%TowerMotion%NNodes * 9 & ! 3 Translation Displacements + 3 orientations + 3 Translation velocities at each node - + u%hubMotion%NNodes * 9 & ! 3 Translation Displacements + 3 orientations + 3 Rotation velocities at each node - + size( u%InflowOnBlade) & - + size( u%InflowOnTower) & !note that we are not passing the inflow on nacelle or hub here - + size( u%UserProp) - - do i=1,p%NumBlades - nu = nu + u%BladeMotion(i)%NNodes * 15 & ! 3 Translation Displacements + 3 orientations + 3 Translation velocities + 3 Rotation velocities + 3 TranslationAcc at each node - + u%BladeRootMotion(i)%NNodes * 3 ! 3 orientations at each node - end do + if (p_AD%CompAeroMaps) then + nu = 0 + NumFieldsForLinearization = 3 ! Translation Displacements + orientations + Translation velocities at each node on the blade mesh + else + nu = u%TowerMotion%NNodes * 9 & ! 3 Translation Displacements + 3 orientations + 3 Translation velocities at each node + + u%hubMotion%NNodes * 9 & ! 3 Translation Displacements + 3 orientations + 3 Rotation velocities at each node + + size( u%InflowOnBlade) & + + size( u%InflowOnTower) & !note that we are not passing the inflow on nacelle or hub here + + size( u%UserProp) + !+ 3 ! 3 velocity components in AvgDiskVel; note that we are not passing the inflow on nacelle or hub here + + NumFieldsForLinearization = 5 ! Translation Displacements + orientations + Translation velocities + Rotation velocities + TranslationAcc at each node on the blade mesh + do i=1,p%NumBlades + nu = nu + u%BladeRootMotion(i)%NNodes * 3 ! 3 orientations at each node + end do + end if + + do i=1,p%NumBl_Lin + nu = nu + u%BladeMotion(i)%NNodes * 3*NumFieldsForLinearization ! 3 components per field + end do + ! all other inputs ignored @@ -6508,54 +6576,59 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) !............... ! AD input mappings stored in p%Jac_u_indx: - !............... + !............... index = 1 - !Module/Mesh/Field: u%TowerMotion%TranslationDisp = 1; - !Module/Mesh/Field: u%TowerMotion%Orientation = 2; - !Module/Mesh/Field: u%TowerMotion%TranslationVel = 3; - do i_meshField = 1,3 - do i=1,u%TowerMotion%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%HubMotion%TranslationDisp = 4; - !Module/Mesh/Field: u%HubMotion%Orientation = 5; - !Module/Mesh/Field: u%HubMotion%RotationVel = 6; - do i_meshField = 4,6 - do i=1,u%HubMotion%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_AD%CompAeroMaps) then - !bjj: if MaxBl (max blades) changes, we need to modify this - !Module/Mesh/Field: u%BladeRootMotion(1)%Orientation = 7; - !Module/Mesh/Field: u%BladeRootMotion(2)%Orientation = 8; - !Module/Mesh/Field: u%BladeRootMotion(3)%Orientation = 9; - do k=1,p%NumBlades - do i_meshField = 6,6 - do i=1,u%BladeRootMotion(k)%NNodes + !Module/Mesh/Field: u%TowerMotion%TranslationDisp = 1; + !Module/Mesh/Field: u%TowerMotion%Orientation = 2; + !Module/Mesh/Field: u%TowerMotion%TranslationVel = 3; + do i_meshField = 1,3 + do i=1,u%TowerMotion%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%HubMotion%TranslationDisp = 4; + !Module/Mesh/Field: u%HubMotion%Orientation = 5; + !Module/Mesh/Field: u%HubMotion%RotationVel = 6; + do i_meshField = 4,6 + do i=1,u%HubMotion%NNodes do j=1,3 - p%Jac_u_indx(index,1) = i_meshField + k + 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 + + !bjj: if MaxBl (max blades) changes, we need to modify this + !Module/Mesh/Field: u%BladeRootMotion(1)%Orientation = 7; + !Module/Mesh/Field: u%BladeRootMotion(2)%Orientation = 8; + !Module/Mesh/Field: u%BladeRootMotion(3)%Orientation = 9; + do k=1,p%NumBlades + do i_meshField = 6,6 + do i=1,u%BladeRootMotion(k)%NNodes + do j=1,3 + p%Jac_u_indx(index,1) = i_meshField + k + 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 !i_meshField - end do !k + end do !i_meshField + end do !k + end if ! .not. compAeroMaps + !bjj: if MaxBl (max blades) changes, we need to modify this !Module/Mesh/Field: u%BladeMotion(1)%TranslationDisp = 10; !Module/Mesh/Field: u%BladeMotion(1)%Orientation = 11; @@ -6574,54 +6647,67 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) !Module/Mesh/Field: u%BladeMotion(3)%TranslationVel = 22; !Module/Mesh/Field: u%BladeMotion(3)%RotationVel = 23; !Module/Mesh/Field: u%BladeMotion(3)%TranslationAcc = 24; - do k=1,p%NumBlades - do i_meshField = 1,5 + do k=1,p%NumBl_Lin + do i_meshField = 1,NumFieldsForLinearization do i=1,u%BladeMotion(k)%NNodes do j=1,3 - p%Jac_u_indx(index,1) = 9 + i_meshField + (k-1)*5 + p%Jac_u_indx(index,1) = 9 + i_meshField + (k-1)*5 ! this should use the MAX possible NumFieldsForLinearization = 5 (so that it's consistent for all cases) 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 !i_meshField + end do !i_meshField end do !k - !Module/Mesh/Field: u%InflowOnBlade(:,:,1) = 25; - !Module/Mesh/Field: u%InflowOnBlade(:,:,2) = 26; - !Module/Mesh/Field: u%InflowOnBlade(:,:,3) = 27; - do k=1,size(u%InflowOnBlade,3) ! p%NumBlades - do i=1,size(u%InflowOnBlade,2) ! numNodes + if (.not. p_AD%CompAeroMaps) then + + !Module/Mesh/Field: u%InflowOnBlade(:,:,1) = 25; + !Module/Mesh/Field: u%InflowOnBlade(:,:,2) = 26; + !Module/Mesh/Field: u%InflowOnBlade(:,:,3) = 27; + do k=1,size(u%InflowOnBlade,3) ! p%NumBlades + do i=1,size(u%InflowOnBlade,2) ! numNodes + do j=1,3 + p%Jac_u_indx(index,1) = 24 + k + 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 !k + + !Module/Mesh/Field: u%InflowOnTower(:,:) = 28; + do i=1,size(u%InflowOnTower,2) ! numNodes do j=1,3 - p%Jac_u_indx(index,1) = 24 + k + p%Jac_u_indx(index,1) = 28 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 !j end do !i - end do !k + + !Module/Mesh/Field: u%UserProp(:,:) = 29,30,31; + do k=1,size(u%UserProp,2) ! p%NumBlades + do i=1,size(u%UserProp,1) ! numNodes + p%Jac_u_indx(index,1) = 28 + k + p%Jac_u_indx(index,2) = 1 !component index: this is a scalar, so 1, but is never used + p%Jac_u_indx(index,3) = i !Node: i + index = index + 1 + end do !i + end do !k + + !Module/Mesh/Field: u%AvgDiskVel(:,:) = 32; + !do j=1,3 + ! p%Jac_u_indx(index,1) = 32 + ! p%Jac_u_indx(index,2) = j !component index: j + ! p%Jac_u_indx(index,3) = 1 !Node: 1 (not really necessary here, since we have only a 1 dimensional array) + ! index = index + 1 + !end do !j + + + end if ! .not. compAeroMaps - !Module/Mesh/Field: u%InflowOnTower(:,:) = 28; - do i=1,size(u%InflowOnTower,2) ! numNodes - do j=1,3 - p%Jac_u_indx(index,1) = 28 - 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 - - !Module/Mesh/Field: u%UserProp(:,:) = 29,30,31; - - do k=1,size(u%UserProp,2) ! p%NumBlades - do i=1,size(u%UserProp,1) ! numNodes - p%Jac_u_indx(index,1) = 28 + k - p%Jac_u_indx(index,2) = 1 !component index: this is a scalar, so 1, but is never used - p%Jac_u_indx(index,3) = i !Node: i - index = index + 1 - end do !i - end do !k !...................................... ! default perturbations, p%du: !...................................... @@ -6660,9 +6746,12 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) p%du(24 + k) = perturb_b(k) ! u%InflowOnBlade(:,:,k) = 24 + k end do p%du(28) = perturb_t ! u%InflowOnTower(:,:) = 28 - do k=1,p%NumBlades + do k=1,p%NumBl_Lin p%du(28+k) = perturb ! u%UserProp(:,:) = 29,30,31 end do + !p%du(32) = minval(perturb_b(1:p%numBlades)) ! u%AvgDiskVel(:) = 32 + + !..................... ! get names of linearized inputs !..................... @@ -6676,59 +6765,76 @@ SUBROUTINE Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat, ErrMsg) InitOut%IsLoad_u = .false. ! None of AeroDyn's inputs are loads InitOut%RotFrame_u = .false. - do k=0,p%NumBlades*p%NumBlNds-1 - InitOut%RotFrame_u(nu - k ) = .true. ! UserProp(:,:) - end do + if (.not. p_AD%CompAeroMaps) then + do k=0,p%NumBl_Lin*p%NumBlNds-1 + InitOut%RotFrame_u(nu - k ) = .true. ! UserProp(:,:) ! TODO TODO TODO add -3 due to DiskAvgVel + end do + endif + index = 1 FieldMask = .false. FieldMask(MASKID_TRANSLATIONDISP) = .true. FieldMask(MASKID_Orientation) = .true. FieldMask(MASKID_TRANSLATIONVel) = .true. - call PackMotionMesh_Names(u%TowerMotion, 'Tower', InitOut%LinNames_u, index, FieldMask=FieldMask) + if (.not. p_AD%CompAeroMaps) call PackMotionMesh_Names(u%TowerMotion, 'Tower', InitOut%LinNames_u, index, FieldMask=FieldMask) FieldMask(MASKID_TRANSLATIONVel) = .false. FieldMask(MASKID_RotationVel) = .true. - call PackMotionMesh_Names(u%HubMotion, 'Hub', InitOut%LinNames_u, index, FieldMask=FieldMask) + if (.not. p_AD%CompAeroMaps) call PackMotionMesh_Names(u%HubMotion, 'Hub', InitOut%LinNames_u, index, FieldMask=FieldMask) index_last = index FieldMask = .false. FieldMask(MASKID_Orientation) = .true. - do k = 1,p%NumBlades - call PackMotionMesh_Names(u%BladeRootMotion(k), 'Blade root '//trim(num2lstr(k)), InitOut%LinNames_u, index, FieldMask=FieldMask) - end do + if (.not. p_AD%CompAeroMaps) then + do k = 1,p%NumBlades + call PackMotionMesh_Names(u%BladeRootMotion(k), 'Blade root '//trim(num2lstr(k)), InitOut%LinNames_u, index, FieldMask=FieldMask) + end do + + + FieldMask(MASKID_RotationVel) = .true. + FieldMask(MASKID_TRANSLATIONAcc) = .true. + end if FieldMask(MASKID_TRANSLATIONDISP) = .true. FieldMask(MASKID_TRANSLATIONVel) = .true. - FieldMask(MASKID_RotationVel) = .true. - FieldMask(MASKID_TRANSLATIONAcc) = .true. - do k=1,p%NumBlades + do k=1,p%NumBl_Lin call PackMotionMesh_Names(u%BladeMotion(k), 'Blade '//trim(num2lstr(k)), InitOut%LinNames_u, index, FieldMask=FieldMask) end do - do k=1,p%NumBlades - do i=1,p%NumBlNds - do j=1,3 - InitOut%LinNames_u(index) = UVW(j)//'-component inflow on blade '//trim(num2lstr(k))//', node '//trim(num2lstr(i))//', m/s' - index = index + 1 + if (.not. p_AD%CompAeroMaps) then + do k=1,p%NumBlades + do i=1,p%NumBlNds + do j=1,3 + InitOut%LinNames_u(index) = UVW(j)//'-component inflow on blade '//trim(num2lstr(k))//', node '//trim(num2lstr(i))//', m/s' + index = index + 1 + end do end do end do - end do - !InitOut%RotFrame_u(index_last:index-1) = .true. ! values on the mesh (and from IfW) are in global coordinates, thus not in the rotating frame + !InitOut%RotFrame_u(index_last:index-1) = .true. ! values on the mesh (and from IfW) are in global coordinates, thus not in the rotating frame - do i=1,p%NumTwrNds - do j=1,3 - InitOut%LinNames_u(index) = UVW(j)//'-component inflow on tower node '//trim(num2lstr(i))//', m/s' - index = index + 1 + do i=1,p%NumTwrNds + do j=1,3 + InitOut%LinNames_u(index) = UVW(j)//'-component inflow on tower node '//trim(num2lstr(i))//', m/s' + index = index + 1 + end do end do - end do - - do k=1,p%NumBlades - do i=1,p%NumBlNds - InitOut%LinNames_u(index) = 'User property on blade '//trim(num2lstr(k))//', node '//trim(num2lstr(i))//', -' - index = index + 1 + + ! UserProp + do k=1,p%NumBl_Lin + do i=1,p%NumBlNds + InitOut%LinNames_u(index) = 'User property on blade '//trim(num2lstr(k))//', node '//trim(num2lstr(i))//', -' + index = index + 1 + end do end do - end do - + + ! AvgDiskVel + !do j=1,3 + ! InitOut%LinNames_u(index) = UVW(j)//'-component inflow of average disk velocity, m/s' + ! index = index + 1 + !end do + + end if + END SUBROUTINE Init_Jacobian_u !---------------------------------------------------------------------------------------------------------------------------------- SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) @@ -6744,7 +6850,7 @@ SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) CHARACTER(*), PARAMETER :: RoutineName = 'Init_Jacobian_x' ! local variables: - INTEGER(IntKi) :: i, j, k + INTEGER(IntKi) :: i, j, k, n, state INTEGER(IntKi) :: nx INTEGER(IntKi) :: nx1 CHARACTER(25) :: NodeTxt @@ -6753,7 +6859,7 @@ SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) ErrMsg = "" - nx = p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx + nx = p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx + p%BEMT%lin_nx ! allocate space for the row/column names and for perturbation sizes ! always allocate this in case it is size zero ... (we use size(p%dx) for many calculations) @@ -6793,34 +6899,46 @@ SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg) InitOut%RotFrame_x(i+nx1) = InitOut%RotFrame_x(i) end do end if - + ! UA states if (p%BEMT%UA%lin_nx>0) then InitOut%DerivOrder_x(1+p%BEMT%DBEMT%lin_nx:nx) = 1 InitOut%RotFrame_x( 1+p%BEMT%DBEMT%lin_nx:nx) = .true. k = 1 + p%BEMT%DBEMT%lin_nx - do j=1,p%NumBlades ! size(x%BEMT%DBEMT%element,2) - do i=1,p%NumBlNds ! size(x%BEMT%DBEMT%element,1) - NodeTxt = 'blade '//trim(num2lstr(j))//', node '//trim(num2lstr(i)) - if (p%BEMT%UA%UAMod/=UA_OYE) then - - InitOut%LinNames_x(k) = 'x1 '//trim(NodeTxt)//', rad' - k = k + 1 + do n=1,p%BEMT%UA%lin_nx + i = p%BEMT%UA%lin_xIndx(n,1) + j = p%BEMT%UA%lin_xIndx(n,2) + state = p%BEMT%UA%lin_xIndx(n,3) - InitOut%LinNames_x(k) = 'x2 '//trim(NodeTxt)//', rad' - k = k + 1 - - InitOut%LinNames_x(k) = 'x3 '//trim(NodeTxt)//', -' - k = k + 1 - endif - - InitOut%LinNames_x(k) = 'x4 '//trim(NodeTxt)//', -' - p%dx(k) = 0.001 ! x4 is a number between 0 and 1, so we need this to be small - k = k + 1 - end do + p%dx(k) = p%BEMT%UA%dx(state) + + NodeTxt = 'x'//trim(num2lstr(state))//' blade '//trim(num2lstr(j))//', node '//trim(num2lstr(i)) + if (state<3) then + InitOut%LinNames_x(k) = trim(NodeTxt)//', rad' ! x1 and x2 are radians + else + InitOut%LinNames_x(k) = trim(NodeTxt)//', -' ! x3, x4 (and x5) are units of cl or cn + end if + InitOut%DerivOrder_x(k) = 1 + InitOut%RotFrame_x(k) = .true. + + k = k + 1 end do end if + ! BEMT states + if (p%BEMT%lin_nx>0) then + call SetErrStat(ErrID_Fatal,'Number of lin states for bem should be zero for now.', ErrStat, ErrMsg, RoutineName) + return + !k = 1 + p%BEMT%DBEMT%lin_nx + p%BEMT%UA%lin_nx + + !InitOut%DerivOrder_x(k:nx) = 1 + !InitOut%RotFrame_x( k:nx) = .false. + ! + !InitOut%LinNames_x(k ) = 'X-component of wake velocity, m/s' + !InitOut%LinNames_x(k+1) = 'Y-component of wake velocity, m/s' + !InitOut%LinNames_x(k+2) = 'Z-component of wake velocity, m/s' + end if + END SUBROUTINE Init_Jacobian_x !---------------------------------------------------------------------------------------------------------------------------------- @@ -6847,8 +6965,13 @@ SUBROUTINE Init_Jacobian( InputFileData, p, p_AD, u, y, m, InitOut, ErrStat, Err ErrStat = ErrID_None ErrMsg = "" -!FIXME: add logic to check that p%NumBlades is not greater than MaxBl. Cannot linearize if that is true. - call Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg) + if (p_AD%CompAeroMaps) then + p%NumBl_Lin = 1 + else + p%NumBl_Lin = p%NumBlades + end if + + call Init_Jacobian_y( p, p_AD, y, InitOut, ErrStat, ErrMsg) ! these matrices will be needed for linearization with frozen wake feature if (p%FrozenWake) then @@ -6856,7 +6979,7 @@ SUBROUTINE Init_Jacobian( InputFileData, p, p_AD, u, y, m, InitOut, ErrStat, Err call AllocAry(m%BEMT%TnInd_op,p%NumBlNds,p%numBlades,'m%BEMT%TnInd_op', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) end if - call Init_Jacobian_u( InputFileData, p, u, InitOut, ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call Init_Jacobian_u( InputFileData, p, p_AD, u, InitOut, ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call Init_Jacobian_x( p, InitOut, ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -6950,12 +7073,17 @@ SUBROUTINE Perturb_u( p, n, perturb_sign, u, du ) CASE (28) !Module/Mesh/Field: u%InflowOnTower(:,:) = 28; u%InflowOnTower(fieldIndx,node) = u%InflowOnTower(fieldIndx,node) + du * perturb_sign + CASE (29) !Module/Mesh/Field: u%UserProp(:,1) = 29; u%UserProp(node,1) = u%UserProp(node,1) + du * perturb_sign CASE (30) !Module/Mesh/Field: u%UserProp(:,2) = 30; u%UserProp(node,2) = u%UserProp(node,2) + du * perturb_sign CASE (31) !Module/Mesh/Field: u%UserProp(:,3) = 31; u%UserProp(node,3) = u%UserProp(node,3) + du * perturb_sign + + !CASE (32) !Module/Mesh/Field: u%AvgDiskVel(:) = 32; + ! u%AvgDiskVel(fieldIndx) = u%AvgDiskVel(fieldIndx) + du * perturb_sign + END SELECT END SUBROUTINE Perturb_u @@ -6974,7 +7102,8 @@ SUBROUTINE Perturb_x( p, n, perturb_sign, x, dx ) ! local variables INTEGER(IntKi) :: Blade ! loop over blade nodes INTEGER(IntKi) :: BladeNode ! loop over blades - INTEGER(IntKi) :: StateIndex ! loop over blades + INTEGER(IntKi) :: StateIndex ! which state we are perturbing + INTEGER(IntKi) :: n_tmp ! dx = p%dx( n ) @@ -6990,16 +7119,19 @@ SUBROUTINE Perturb_x( p, n, perturb_sign, x, dx ) endif else - !call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), size(x%BEMT%UA%element(1,1)%x), Blade, BladeNode, StateIndex ) - if (p%BEMT%UA%UAMod==UA_OYE) then - call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), 1, Blade, BladeNode, StateIndex ) - StateIndex=4 ! Always the 4th one + n_tmp = n - p%BEMT%DBEMT%lin_nx + + if (n_tmp <= p%BEMT%UA%lin_nx) then + BladeNode = p%BEMT%UA%lin_xIndx(n_tmp,1) ! node + Blade = p%BEMT%UA%lin_xIndx(n_tmp,2) ! blade + StateIndex = p%BEMT%UA%lin_xIndx(n_tmp,3) ! state + + x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) = x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) + dx * perturb_sign else - call GetStateIndices( n - p%BEMT%DBEMT%lin_nx, size(x%BEMT%UA%element,2), size(x%BEMT%UA%element,1), 4, Blade, BladeNode, StateIndex ) - endif - x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) = x%BEMT%UA%element(BladeNode,Blade)%x(StateIndex) + dx * perturb_sign - + StateIndex = n_tmp - p%BEMT%UA%lin_nx + x%BEMT%V_w(StateIndex) = x%BEMT%V_w(StateIndex) + dx * perturb_sign + end if end if contains @@ -7046,17 +7178,17 @@ SUBROUTINE Compute_dY(p, p_AD, y_p, y_m, delta_p, delta_m, dY) indx_first = 1 - call PackLoadMesh_dY(y_p%TowerLoad, y_m%TowerLoad, dY, indx_first) + if (.not. p_AD%CompAeroMaps) call PackLoadMesh_dY(y_p%TowerLoad, y_m%TowerLoad, dY, indx_first) - do k=1,p%NumBlades + do k=1,p%NumBl_Lin call PackLoadMesh_dY(y_p%BladeLoad(k), y_m%BladeLoad(k), dY, indx_first) end do - - do k=1,p%NumOuts + p%BldNd_TotNumOuts - dY(k+indx_first-1) = y_p%WriteOutput(k) - y_m%WriteOutput(k) - end do - + if (.not. p_AD%CompAeroMaps) then + 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 / (delta_p + delta_m) @@ -7076,6 +7208,8 @@ SUBROUTINE Compute_dX(p, x_p, x_m, delta_p, delta_m, dX) ! local variables: INTEGER(IntKi) :: i ! loop over blade nodes INTEGER(IntKi) :: j ! loop over blades + INTEGER(IntKi) :: k ! loop over states + INTEGER(IntKi) :: n ! loop over active UA states INTEGER(IntKi) :: indx_first ! index indicating next value of dY to be filled @@ -7100,25 +7234,23 @@ SUBROUTINE Compute_dX(p, x_p, x_m, delta_p, delta_m, dX) end if if (p%BEMT%UA%lin_nx>0) then - - if (p%BEMT%UA%UAMod==UA_OYE) then - do j=1,size(x_p%BEMT%UA%element,2) ! number of blades - do i=1,size(x_p%BEMT%UA%element,1) ! number of nodes per blade - dX(indx_first) = x_p%BEMT%UA%element(i,j)%x(4) - x_m%BEMT%UA%element(i,j)%x(4) - indx_first = indx_first + 1 ! = index_first += 4 - end do - end do - else - do j=1,size(x_p%BEMT%UA%element,2) ! number of blades - do i=1,size(x_p%BEMT%UA%element,1) ! number of nodes per blade - dX(indx_first:indx_first+3) = x_p%BEMT%UA%element(i,j)%x(1:4) - x_m%BEMT%UA%element(i,j)%x(1:4) - indx_first = indx_first + 4 ! = index_first += 4 - end do - end do - endif + do n=1,p%BEMT%UA%lin_nx + i = p%BEMT%UA%lin_xIndx(n,1) + j = p%BEMT%UA%lin_xIndx(n,2) + k = p%BEMT%UA%lin_xIndx(n,3) + dX(indx_first) = x_p%BEMT%UA%element(i,j)%x(k) - x_m%BEMT%UA%element(i,j)%x(k) + indx_first = indx_first + 1 + end do + end if + if (p%BEMT%lin_nx>0) then ! skewWake + !do j=1,size(x_p%BEMT%v_w) + ! dX(indx_first) = x_p%BEMT%v_w(j) - x_m%BEMT%v_w(j) + ! indx_first = indx_first + 1 + !end do + end if dX = dX / (delta_p + delta_m) END SUBROUTINE Compute_dX diff --git a/modules/aerodyn/src/AeroDyn_Registry.txt b/modules/aerodyn/src/AeroDyn_Registry.txt index 5f763c8a6c..87d7233eee 100644 --- a/modules/aerodyn/src/AeroDyn_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Registry.txt @@ -84,6 +84,7 @@ typedef ^ RotInitInputType R8Ki NacellePosition {3} - - "X-Y-Z reference positio typedef ^ RotInitInputType R8Ki NacelleOrientation {3}{3} - - "DCM reference orientation of nacelle" - typedef ^ RotInitInputType IntKi AeroProjMod - 1 - "Flag to switch between different projection models" - typedef ^ RotInitInputType IntKi AeroBEM_Mod - -1 - "Flag to switch between different BEM Model" - +typedef ^ RotInitInputType ReKi RotSpeed - - 0 "Rotor speed used when AeroDyn is computing aero maps" "rad/s" typedef ^ InitInputType RotInitInputType rotors {:} - - "Init Input Types for rotors" - typedef ^ InitInputType CHARACTER(1024) InputFile - - - "Name of the input file" - @@ -91,6 +92,7 @@ typedef ^ InitInputType CHARACTER(1024) RootName - - - "RootName for writing out typedef ^ InitInputType LOGICAL UsePrimaryInputFile - .TRUE. - "Read input file instead of passed data" - typedef ^ InitInputType FileInfoType PassedPrimaryInputData - - - "Primary input file as FileInfoType (set by driver/glue code)" - typedef ^ InitInputType Logical Linearize - .FALSE. - "Flag that tells this module if the glue code wants to linearize." - +typedef ^ InitInputType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false)" - typedef ^ InitInputType ReKi Gravity - - - "Gravity force" Nm/s^2 typedef ^ InitInputType IntKi MHK - - - "MHK turbine type switch" - typedef ^ InitInputType ReKi defFldDens - - - "Default fluid density from the driver; may be overwritten" kg/m^3 @@ -421,6 +423,7 @@ typedef ^ RotInputType ReKi InflowOnTower {:}{:} - - "U,V,W at nodes on the towe typedef ^ RotInputType ReKi InflowOnHub {3} - - "U,V,W at hub" m/s typedef ^ RotInputType ReKi InflowOnNacelle {3} - - "U,V,W at nacelle" m/s typedef ^ RotInputType ReKi InflowOnTailFin {3} - - "U,V,W at tailfin" m/s +typedef ^ RotInputType ReKi AvgDiskVel {3} - 0.0 "disk-averaged U,V,W" m/s typedef ^ RotInputType ReKi UserProp {:}{:} - - "Optional user property for interpolating airfoils (per element per blade)" - diff --git a/modules/aerodyn/src/AeroDyn_Types.f90 b/modules/aerodyn/src/AeroDyn_Types.f90 index c26eef9940..b0957a03a7 100644 --- a/modules/aerodyn/src/AeroDyn_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Types.f90 @@ -101,6 +101,7 @@ MODULE AeroDyn_Types REAL(R8Ki) , DIMENSION(1:3,1:3) :: NacelleOrientation !< DCM reference orientation of nacelle [-] INTEGER(IntKi) :: AeroProjMod = 1 !< Flag to switch between different projection models [-] INTEGER(IntKi) :: AeroBEM_Mod = -1 !< Flag to switch between different BEM Model [-] + REAL(ReKi) :: RotSpeed !< Rotor speed used when AeroDyn is computing aero maps [rad/s] END TYPE RotInitInputType ! ======================= ! ========= AD_InitInputType ======= @@ -111,6 +112,7 @@ MODULE AeroDyn_Types LOGICAL :: UsePrimaryInputFile = .TRUE. !< Read input file instead of passed data [-] TYPE(FileInfoType) :: PassedPrimaryInputData !< Primary input file as FileInfoType (set by driver/glue code) [-] LOGICAL :: Linearize = .FALSE. !< Flag that tells this module if the glue code wants to linearize. [-] + LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false) [-] REAL(ReKi) :: Gravity !< Gravity force [Nm/s^2] INTEGER(IntKi) :: MHK !< MHK turbine type switch [-] REAL(ReKi) :: defFldDens !< Default fluid density from the driver; may be overwritten [kg/m^3] @@ -457,6 +459,7 @@ MODULE AeroDyn_Types REAL(ReKi) , DIMENSION(1:3) :: InflowOnHub !< U,V,W at hub [m/s] REAL(ReKi) , DIMENSION(1:3) :: InflowOnNacelle !< U,V,W at nacelle [m/s] REAL(ReKi) , DIMENSION(1:3) :: InflowOnTailFin !< U,V,W at tailfin [m/s] + REAL(ReKi) , DIMENSION(1:3) :: AvgDiskVel !< disk-averaged U,V,W [m/s] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: UserProp !< Optional user property for interpolating airfoils (per element per blade) [-] END TYPE RotInputType ! ======================= @@ -1436,6 +1439,7 @@ SUBROUTINE AD_CopyRotInitInputType( SrcRotInitInputTypeData, DstRotInitInputType DstRotInitInputTypeData%NacelleOrientation = SrcRotInitInputTypeData%NacelleOrientation DstRotInitInputTypeData%AeroProjMod = SrcRotInitInputTypeData%AeroProjMod DstRotInitInputTypeData%AeroBEM_Mod = SrcRotInitInputTypeData%AeroBEM_Mod + DstRotInitInputTypeData%RotSpeed = SrcRotInitInputTypeData%RotSpeed END SUBROUTINE AD_CopyRotInitInputType SUBROUTINE AD_DestroyRotInitInputType( RotInitInputTypeData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -1519,6 +1523,7 @@ SUBROUTINE AD_PackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Db_BufSz = Db_BufSz + SIZE(InData%NacelleOrientation) ! NacelleOrientation Int_BufSz = Int_BufSz + 1 ! AeroProjMod Int_BufSz = Int_BufSz + 1 ! AeroBEM_Mod + Re_BufSz = Re_BufSz + 1 ! RotSpeed IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -1617,6 +1622,8 @@ SUBROUTINE AD_PackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Int_Xferred = Int_Xferred + 1 IntKiBuf(Int_Xferred) = InData%AeroBEM_Mod Int_Xferred = Int_Xferred + 1 + ReKiBuf(Re_Xferred) = InData%RotSpeed + Re_Xferred = Re_Xferred + 1 END SUBROUTINE AD_PackRotInitInputType SUBROUTINE AD_UnPackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -1737,6 +1744,8 @@ SUBROUTINE AD_UnPackRotInitInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrSt Int_Xferred = Int_Xferred + 1 OutData%AeroBEM_Mod = IntKiBuf(Int_Xferred) Int_Xferred = Int_Xferred + 1 + OutData%RotSpeed = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 END SUBROUTINE AD_UnPackRotInitInputType SUBROUTINE AD_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrStat, ErrMsg ) @@ -1777,6 +1786,7 @@ SUBROUTINE AD_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrSt CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) IF (ErrStat>=AbortErrLev) RETURN DstInitInputData%Linearize = SrcInitInputData%Linearize + DstInitInputData%CompAeroMaps = SrcInitInputData%CompAeroMaps DstInitInputData%Gravity = SrcInitInputData%Gravity DstInitInputData%MHK = SrcInitInputData%MHK DstInitInputData%defFldDens = SrcInitInputData%defFldDens @@ -1900,6 +1910,7 @@ SUBROUTINE AD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg DEALLOCATE(Int_Buf) END IF Int_BufSz = Int_BufSz + 1 ! Linearize + Int_BufSz = Int_BufSz + 1 ! CompAeroMaps Re_BufSz = Re_BufSz + 1 ! Gravity Int_BufSz = Int_BufSz + 1 ! MHK Re_BufSz = Re_BufSz + 1 ! defFldDens @@ -2017,6 +2028,8 @@ SUBROUTINE AD_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg ENDIF IntKiBuf(Int_Xferred) = TRANSFER(InData%Linearize, IntKiBuf(1)) Int_Xferred = Int_Xferred + 1 + IntKiBuf(Int_Xferred) = TRANSFER(InData%CompAeroMaps, IntKiBuf(1)) + Int_Xferred = Int_Xferred + 1 ReKiBuf(Re_Xferred) = InData%Gravity Re_Xferred = Re_Xferred + 1 IntKiBuf(Int_Xferred) = InData%MHK @@ -2172,6 +2185,8 @@ SUBROUTINE AD_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Err IF(ALLOCATED(Int_Buf)) DEALLOCATE(Int_Buf) OutData%Linearize = TRANSFER(IntKiBuf(Int_Xferred), OutData%Linearize) Int_Xferred = Int_Xferred + 1 + OutData%CompAeroMaps = TRANSFER(IntKiBuf(Int_Xferred), OutData%CompAeroMaps) + Int_Xferred = Int_Xferred + 1 OutData%Gravity = ReKiBuf(Re_Xferred) Re_Xferred = Re_Xferred + 1 OutData%MHK = IntKiBuf(Int_Xferred) @@ -15894,6 +15909,7 @@ SUBROUTINE AD_CopyRotInputType( SrcRotInputTypeData, DstRotInputTypeData, CtrlCo DstRotInputTypeData%InflowOnHub = SrcRotInputTypeData%InflowOnHub DstRotInputTypeData%InflowOnNacelle = SrcRotInputTypeData%InflowOnNacelle DstRotInputTypeData%InflowOnTailFin = SrcRotInputTypeData%InflowOnTailFin + DstRotInputTypeData%AvgDiskVel = SrcRotInputTypeData%AvgDiskVel IF (ALLOCATED(SrcRotInputTypeData%UserProp)) THEN i1_l = LBOUND(SrcRotInputTypeData%UserProp,1) i1_u = UBOUND(SrcRotInputTypeData%UserProp,1) @@ -16127,6 +16143,7 @@ SUBROUTINE AD_PackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Err Re_BufSz = Re_BufSz + SIZE(InData%InflowOnHub) ! InflowOnHub Re_BufSz = Re_BufSz + SIZE(InData%InflowOnNacelle) ! InflowOnNacelle Re_BufSz = Re_BufSz + SIZE(InData%InflowOnTailFin) ! InflowOnTailFin + Re_BufSz = Re_BufSz + SIZE(InData%AvgDiskVel) ! AvgDiskVel Int_BufSz = Int_BufSz + 1 ! UserProp allocated yes/no IF ( ALLOCATED(InData%UserProp) ) THEN Int_BufSz = Int_BufSz + 2*2 ! UserProp upper/lower bounds for each dimension @@ -16410,6 +16427,10 @@ SUBROUTINE AD_PackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, Err ReKiBuf(Re_Xferred) = InData%InflowOnTailFin(i1) Re_Xferred = Re_Xferred + 1 END DO + DO i1 = LBOUND(InData%AvgDiskVel,1), UBOUND(InData%AvgDiskVel,1) + ReKiBuf(Re_Xferred) = InData%AvgDiskVel(i1) + Re_Xferred = Re_Xferred + 1 + END DO IF ( .NOT. ALLOCATED(InData%UserProp) ) THEN IntKiBuf( Int_Xferred ) = 0 Int_Xferred = Int_Xferred + 1 @@ -16802,6 +16823,12 @@ SUBROUTINE AD_UnPackRotInputType( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, OutData%InflowOnTailFin(i1) = ReKiBuf(Re_Xferred) Re_Xferred = Re_Xferred + 1 END DO + i1_l = LBOUND(OutData%AvgDiskVel,1) + i1_u = UBOUND(OutData%AvgDiskVel,1) + DO i1 = LBOUND(OutData%AvgDiskVel,1), UBOUND(OutData%AvgDiskVel,1) + OutData%AvgDiskVel(i1) = ReKiBuf(Re_Xferred) + Re_Xferred = Re_Xferred + 1 + END DO IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! UserProp not allocated Int_Xferred = Int_Xferred + 1 ELSE @@ -18311,6 +18338,12 @@ SUBROUTINE AD_Input_ExtrapInterp1(u1, u2, tin, u_out, tin_out, ErrStat, ErrMsg ) END DO ENDDO DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1) + DO i1 = LBOUND(u_out%rotors(i01)%AvgDiskVel,1),UBOUND(u_out%rotors(i01)%AvgDiskVel,1) + b = -(u1%rotors(i01)%AvgDiskVel(i1) - u2%rotors(i01)%AvgDiskVel(i1)) + u_out%rotors(i01)%AvgDiskVel(i1) = u1%rotors(i01)%AvgDiskVel(i1) + b * ScaleFactor + END DO + ENDDO + DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1) IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN DO i2 = LBOUND(u_out%rotors(i01)%UserProp,2),UBOUND(u_out%rotors(i01)%UserProp,2) DO i1 = LBOUND(u_out%rotors(i01)%UserProp,1),UBOUND(u_out%rotors(i01)%UserProp,1) @@ -18469,6 +18502,13 @@ SUBROUTINE AD_Input_ExtrapInterp2(u1, u2, u3, tin, u_out, tin_out, ErrStat, ErrM END DO ENDDO DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1) + DO i1 = LBOUND(u_out%rotors(i01)%AvgDiskVel,1),UBOUND(u_out%rotors(i01)%AvgDiskVel,1) + b = (t(3)**2*(u1%rotors(i01)%AvgDiskVel(i1) - u2%rotors(i01)%AvgDiskVel(i1)) + t(2)**2*(-u1%rotors(i01)%AvgDiskVel(i1) + u3%rotors(i01)%AvgDiskVel(i1)))* scaleFactor + c = ( (t(2)-t(3))*u1%rotors(i01)%AvgDiskVel(i1) + t(3)*u2%rotors(i01)%AvgDiskVel(i1) - t(2)*u3%rotors(i01)%AvgDiskVel(i1) ) * scaleFactor + u_out%rotors(i01)%AvgDiskVel(i1) = u1%rotors(i01)%AvgDiskVel(i1) + b + c * t_out + END DO + ENDDO + DO i01 = LBOUND(u_out%rotors,1),UBOUND(u_out%rotors,1) IF (ALLOCATED(u_out%rotors(i01)%UserProp) .AND. ALLOCATED(u1%rotors(i01)%UserProp)) THEN DO i2 = LBOUND(u_out%rotors(i01)%UserProp,2),UBOUND(u_out%rotors(i01)%UserProp,2) DO i1 = LBOUND(u_out%rotors(i01)%UserProp,1),UBOUND(u_out%rotors(i01)%UserProp,1) diff --git a/modules/aerodyn/src/BEMT.f90 b/modules/aerodyn/src/BEMT.f90 index 33fe5d2821..c71362ca7f 100644 --- a/modules/aerodyn/src/BEMT.f90 +++ b/modules/aerodyn/src/BEMT.f90 @@ -882,7 +882,7 @@ subroutine BEMT_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, AFInfo, !........................ do j = 1,p%numBlades do i = 1,p%numBladeNodes - call DBEMT_UpdateStates(i, j, t, n, m%u_DBEMT, p%DBEMT, x%DBEMT, OtherState%DBEMT, m%DBEMT, errStat2, errMsg2) + call DBEMT_UpdateStates(i, j, t, n, m%u_DBEMT, uTimes, p%DBEMT, x%DBEMT, OtherState%DBEMT, m%DBEMT, errStat2, errMsg2) if (ErrStat2 /= ErrID_None) then call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName//trim(NodeText(i,j))) if (errStat >= AbortErrLev) return diff --git a/modules/aerodyn/src/DBEMT.f90 b/modules/aerodyn/src/DBEMT.f90 index 3c12e5c3b4..a80c2c0a09 100644 --- a/modules/aerodyn/src/DBEMT.f90 +++ b/modules/aerodyn/src/DBEMT.f90 @@ -25,6 +25,7 @@ ! [2] R. Damiani, J.Jonkman ! DBEMT Theory Rev. 3 ! Unpublished +! module DBEMT use NWTC_Library @@ -328,13 +329,14 @@ end subroutine DBEMT_InitStates !!---------------------------------------------------------------------------------------------------------------------------------- !> Loose coupling routine for solving for constraint states, integrating continuous states, and updating discrete and other states. !! Continuous, constraint, discrete, and other states are updated for t + Interval -subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errMsg ) +subroutine DBEMT_UpdateStates( i, j, t, n, u, uTimes, p, x, OtherState, m, errStat, errMsg ) !.................................................................................................................................. integer(IntKi), intent(in ) :: i !< blade node counter integer(IntKi), intent(in ) :: j !< blade counter real(DbKi), intent(in ) :: t !< Current simulation time in seconds integer(IntKi), intent(in ) :: n !< Current simulation time step n = 0,1,... type(DBEMT_InputType), intent(in ) :: u(2) !< Inputs at t and t+dt + real(DbKi), intent(in ) :: uTimes(2) ! Times associated with u(:), in seconds type(DBEMT_ParameterType), intent(in ) :: p !< Parameters type(DBEMT_ContinuousStateType), intent(inout) :: x !< Input: Continuous states at t; !! Output: Continuous states at t + Interval @@ -346,7 +348,6 @@ subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errM ! local variables real(ReKi) :: A, B, C0, k_tau, C0_2 ! tau1_plus1, C_tau1, C, K1 integer(IntKi) :: indx - real(DbKi) :: utimes(2) TYPE(DBEMT_ElementInputType) :: u_elem(2) !< Inputs at utimes @@ -364,9 +365,6 @@ subroutine DBEMT_UpdateStates( i, j, t, n, u, p, x, OtherState, m, errStat, errM call DBEMT_InitStates( i, j, u(1), p, x, OtherState ) if (p%DBEMT_Mod == DBEMT_cont_tauConst) then ! continuous formulation: - utimes(1) = t - utimes(2) = t + p%dt - u_elem(1) = u(1)%element(i,j) u_elem(2) = u(2)%element(i,j) call DBEMT_ABM4( i, j, t, n, u_elem, utimes, p, x, OtherState, m, ErrStat, ErrMsg ) @@ -924,4 +922,4 @@ subroutine DBEMT_End( u, p, x, OtherState, m, ErrStat, ErrMsg ) END SUBROUTINE DBEMT_End -end module DBEMT +end module DBEMT \ No newline at end of file diff --git a/modules/aerodyn/src/FVW_Subs.f90 b/modules/aerodyn/src/FVW_Subs.f90 index b1a6eb931e..bf53755c2b 100644 --- a/modules/aerodyn/src/FVW_Subs.f90 +++ b/modules/aerodyn/src/FVW_Subs.f90 @@ -462,7 +462,6 @@ 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)) @@ -481,11 +480,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 E_NW' + print*,trim(label),'NaN in E_NW'//trim(num2lstr(iW)) have_nan=.True. endif if (any(isnan(x%W(iW)%Eps_FW))) then - print*,trim(label),'NaN in E_FW' + print*,trim(label),'NaN in E_FW'//trim(num2lstr(iW)) have_nan=.True. endif if (any(isnan(z%W(iW)%Gamma_LL))) then diff --git a/modules/aerodyn/src/UnsteadyAero.f90 b/modules/aerodyn/src/UnsteadyAero.f90 index dce8be0f54..f40d8f5d90 100644 --- a/modules/aerodyn/src/UnsteadyAero.f90 +++ b/modules/aerodyn/src/UnsteadyAero.f90 @@ -42,6 +42,7 @@ module UnsteadyAero use NWTC_Library use UnsteadyAero_Types use AirfoilInfo + use NWTC_LAPACK implicit none @@ -723,8 +724,9 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) integer(IntKi) :: ErrStat2 character(*), parameter :: RoutineName = 'UA_SetParameters' logical :: IsUsed(size(AFInfo)) + INTEGER :: UA_NumLinStates ! potentially put in p, number of states per blade node that are linearized - INTEGER(IntKi) :: i, j + INTEGER(IntKi) :: i, j, k, n @@ -746,19 +748,35 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) p%ShedEffect = InitInp%ShedEffect if (p%UAMod==UA_HGM .or. p%UAMod==UA_HGMV) then - p%lin_nx = p%numBlades*p%nNodesPerBlade*4 ! 4 continuous states per node per blade (5th state isn't currently linearizable) + UA_NumLinStates = 4 + ! set the maximum number of states + ! note: we will subtract states for nodes where UA is off for good, below + p%lin_nx = p%numBlades*p%nNodesPerBlade*UA_NumLinStates ! 4 continuous states per node per blade (5th state isn't currently linearizable) else if (p%UAMod==UA_OYE) then - p%lin_nx = p%numBlades*p%nNodesPerBlade*1 ! continuous state per node per blade, but stored at position 4 + UA_NumLinStates = 1 + p%lin_nx = p%numBlades*p%nNodesPerBlade*UA_NumLinStates ! continuous state per node per blade, but stored at position 4 else p%lin_nx = 0 + UA_NumLinStates = 0 end if + ! Compute derivative step size + ! --- ADENV + !p%dx = 0.5_R8Ki * D2R_D + !p%dx(4) = 0.0001_R8Ki + ! --- ADLEG + p%dx = 2.0_R8Ki * D2R_D + p%dx(4) = 0.001 ! x4 is a number between 0 and 1, so we need this to be small + + + p%UA_off_forGood = .false. ! flag that determines if UA should be turned off for the whole simulation if (allocated(InitInp%UAOff_innerNode)) then do j=1,min(size(p%UA_off_forGood,2), size(InitInp%UAOff_innerNode)) !blade do i=1,min(InitInp%UAOff_innerNode(j),size(p%UA_off_forGood,1)) !node ! call WrScr( 'Warning: Turning off Unsteady Aerodynamics on inner node (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' ) p%UA_off_forGood(i,j) = .true. + !p%lin_nx = p%lin_nx - UA_NumLinStates end do end do end if @@ -767,7 +785,10 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) do j=1,min(size(p%UA_off_forGood,2), size(InitInp%UAOff_outerNode)) !blade do i=InitInp%UAOff_outerNode(j), size(p%UA_off_forGood,1) !node ! call WrScr( 'Warning: Turning off Unsteady Aerodynamics on outer node (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' ) - p%UA_off_forGood(i,j) = .true. + if (.not. p%UA_off_forGood(i,j)) then + p%UA_off_forGood(i,j) = .true. + !p%lin_nx = p%lin_nx - UA_NumLinStates + end if end do end do end if @@ -780,12 +801,40 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) if (ErrStat2 > ErrID_None) then call WrScr( 'Warning: Turning off Unsteady Aerodynamics because '//trim(ErrMsg2)//' (node '//trim(num2lstr(i))//', blade '//trim(num2lstr(j))//')' ) p%UA_off_forGood(i,j) = .true. + !p%lin_nx = p%lin_nx - UA_NumLinStates end if end if end do end do + + ! set up index array for linearization + p%lin_nx = max(0, p%lin_nx) + call AllocAry(p%lin_xIndx,p%lin_nx,3,'p%lin_xIndx',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + if (ErrStat >= AbortErrLev) return + + if (p%lin_nx > 0) then + n = 1 + do j=1,size(p%UA_off_forGood,2) !blade + do i=1,size(p%UA_off_forGood,1) !node + !if (.not. p%UA_off_forGood(i,j)) then + do k=1,UA_NumLinStates + p%lin_xIndx(n,1) = i ! node + p%lin_xIndx(n,2) = j ! blade + + if (p%UAMod==UA_OYE) then + p%lin_xIndx(n,3) = 4 ! Hack for UA_Oye, state is 4 instead of 1 for now + else + p%lin_xIndx(n,3) = k ! state + endif + n = n + 1 + end do + !end if + end do + end do + end if + ! check that the airfoils have appropriate data for UA IsUsed = .false. do j=1,size(p%UA_off_forGood,2) !blade @@ -798,7 +847,7 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg ) do i=1,size(AFInfo,1) if (IsUsed(i)) then - call UA_ValidateAFI(InitInp%UAMod, AFInfo(i), ErrStat2, ErrMsg2) + call UA_ValidateAFI(p%UAMod, p%Flookup, AFInfo(i), ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) end if end do @@ -971,6 +1020,8 @@ subroutine UA_ReInit( p, x, xd, OtherState, m, ErrStat, ErrMsg ) do i = 1, size(OtherState%xdot) call UA_CopyContState( x, OtherState%xdot(i), MESH_UPDATECOPY, ErrStat2, ErrMsg2) ! there are no meshes, so the control code is irrelevant call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + call UA_CopyContState( x, OtherState%xHistory(i), MESH_UPDATECOPY, ErrStat2, ErrMsg2) ! there are no meshes, so the control code is irrelevant + call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) end do if (p%UAMod == UA_HGMV) then @@ -1315,8 +1366,8 @@ subroutine UA_Init( InitInp, u, p, x, xd, OtherState, y, m, Interval, & InitOut%WriteOutputUnt(iOffset+41) ='(-)' InitOut%WriteOutputUnt(iOffset+42) ='(-)' InitOut%WriteOutputUnt(iOffset+43) ='(-)' - InitOut%WriteOutputUnt(iOffset+44) ='(deg)' - InitOut%WriteOutputUnt(iOffset+45) ='(-)' + InitOut%WriteOutputUnt(iOffset+44) ='(-)' + InitOut%WriteOutputUnt(iOffset+45) ='(deg)' end if @@ -1388,8 +1439,9 @@ subroutine UA_ValidateInput(InitInp, ErrStat, ErrMsg) end subroutine UA_ValidateInput !============================================================================== -subroutine UA_ValidateAFI(UAMod, AFInfo, ErrStat, ErrMsg) +subroutine UA_ValidateAFI(UAMod, FLookup, AFInfo, ErrStat, ErrMsg) integer(IntKi), intent(in ) :: UAMod ! which UA model we are using + logical, intent(in ) :: FLookup ! lookup table type(AFI_ParameterType), target, intent(in ) :: AFInfo ! The airfoil parameter data integer(IntKi), intent( out) :: ErrStat ! Error status of the operation character(*), intent( out) :: ErrMsg ! Error message if ErrStat /= ErrID_None @@ -1419,6 +1471,8 @@ subroutine UA_ValidateAFI(UAMod, AFInfo, ErrStat, ErrMsg) call SetErrStat(ErrID_Fatal, 'UA St_sh parameter must not be 0 in "'//trim(AFInfo%FileName)//'".', ErrStat, ErrMsg, RoutineName ) end if + ! we won't check alpha1 or alph2 validity if we aren't using them for the lookup (curve fit) + if (.not. Flookup) then if ( AFInfo%Table(j)%UA_BL%alpha1 > pi .or. AFInfo%Table(j)%UA_BL%alpha1 < -pi ) then call SetErrStat(ErrID_Fatal, 'UA alpha1 parameter must be between -180 and 180 degrees in "'//trim(AFInfo%FileName)//'".', ErrStat, ErrMsg, RoutineName ) end if @@ -1438,7 +1492,7 @@ subroutine UA_ValidateAFI(UAMod, AFInfo, ErrStat, ErrMsg) if ( AFInfo%Table(j)%UA_BL%alpha2 > AFInfo%Table(j)%UA_BL%alpha0 ) then call SetErrStat(ErrID_Fatal, 'UA alpha0 parameter must be greater than alpha2 in "'//trim(AFInfo%FileName)//'".', ErrStat, ErrMsg, RoutineName ) end if - + end if ! don't check alpha1 and alpha2 unless they are going to be used if ( AFInfo%Table(j)%UA_BL%filtCutOff < 0.0_ReKi ) then call SetErrStat(ErrID_Fatal, 'UA filtCutOff parameter must be greater than 0 in "'//trim(AFInfo%FileName)//'".', ErrStat, ErrMsg, RoutineName ) end if @@ -1634,7 +1688,6 @@ subroutine UA_UpdateDiscOtherState_BV( i, j, u, p, xd, OtherState, AFInfo, m, Er real(ReKi) :: alpha_minus1 !< 3/4 chord angle of attack at real(ReKi) :: alpha_filt_cur !< real(ReKi) :: alpha_filt_minus1 !< - real(ReKi) :: Tu !< Time constant based on u=Vrel and chord real(ReKi) :: dynamicFilterCutoffHz !< find frequency based on reduced frequency of k = BL_p%filtCutOff real(ReKi) :: LowPassConst !< ! @@ -1758,7 +1811,7 @@ subroutine BV_getAlphas(i, j, u, p, xd, BL_p, tc, alpha_34, alphaE_L, alphaLag_D alphaLag_D = alpha_34 - dalphaD*isgn ! NOTE: not effective alpha yet for drag end subroutine BV_getAlphas !============================================================================== -!> Calculate gamma for lift and drag based rel thickness. See CACTUS BV_DynStall.f95 +!> Calculate gamma for lift and drag based rel thickness. See CACTUS BV_DynStall.f95 subroutine BV_getGammas(tc, umach, gammaL, gammaD) real(ReKi), intent(in) :: tc !< Relative thickness of airfoil real(ReKi), intent(in) :: umach !< Mach number of Urel, = Urel*MinfMinf (freestrem Mach), 0 for incompressible @@ -2214,7 +2267,7 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, type(UA_InputType) :: u_interp_raw ! Input at current timestep, t and t+dt type(UA_InputType) :: u_interp ! Input at current timestep, t and t+dt type(AFI_UA_BL_Type) :: BL_p ! airfoil UA parameters retrieved in Kelvin Chain - real(ReKi) :: Tu + real(R8Ki) :: Tu ! Initialize variables @@ -2235,6 +2288,7 @@ subroutine UA_UpdateStates( i, j, t, n, u, uTimes, p, x, xd, OtherState, AFInfo, call UA_fixInputs(u_interp_raw, u_interp, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + if (p%UAMod == UA_HGM .or. p%UAMod == UA_HGMV) then @@ -2421,7 +2475,7 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) integer(IntKi) :: errStat2 character(*), parameter :: RoutineName = 'HGM_Steady' - real(ReKi) :: Tu + real(R8Ki) :: Tu real(ReKi) :: alphaE real(ReKi) :: alphaF real(ReKi) :: alpha_34 @@ -2486,10 +2540,12 @@ SUBROUTINE HGM_Steady( i, j, u, p, x, AFInfo, ErrStat, ErrMsg ) !call AFI_ComputeAirfoilCoefs( alphaF, u%Re, u%UserProp, AFInfo, AFI_interp, ErrStat, ErrMsg) !x%x(4) = AFI_interp%f_st else - print*,'HGM_steady, should never happen' - STOP + call WrScr('>>> HGM_steady logic error: should never happen.') + call SetErrStat(ErrID_FATAL,"Programming error.",ErrStat,ErrMsg,RoutineName) + return end if + ! calculate x%x(4) = fs_aF = f_st(alphaF): x%x(4) = AFI_interp%f_st x%x(5) = 0.0_R8Ki @@ -2521,14 +2577,14 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx integer(IntKi) :: errStat2 character(*), parameter :: RoutineName = 'UA_CalcContStateDeriv' - real(ReKi) :: Tu + real(R8Ki) :: Tu real(ReKi) :: alphaE real(ReKi) :: alphaF - real(ReKi) :: Clp + real(R8Ki) :: Clp real(ReKi) :: cRate ! slope of the piecewise linear region of fully attached polar real(R8Ki) :: x4 real(ReKi) :: alpha_34 - real(ReKi), parameter :: U_dot = 0.0_ReKi ! at some point we may add this term + real(R8Ki), parameter :: U_dot = 0.0_R8Ki ! at some point we may add this term TYPE(UA_InputType) :: u ! Inputs at t real(R8Ki) :: CnC_dot, One_Plus_Sqrt_x4, cv_dot, CnC @@ -2654,8 +2710,9 @@ subroutine UA_CalcContStateDeriv( i, j, t, u_in, p, x, OtherState, AFInfo, m, dx dxdt%x(5) = cv_dot - x%x(5)/(BL_p%T_V0 * Tu) else - print*,'>>> UA_CalcContStateDeriv, should never happen.' - STOP ! should never happen + call WrScr('>>> UA_CalcContStateDeriv logic error: should never happen.') + call SetErrStat(ErrID_FATAL,"Programming error.",ErrStat,ErrMsg,RoutineName) + return end if END SUBROUTINE UA_CalcContStateDeriv @@ -2668,14 +2725,10 @@ SUBROUTINE Get_HGM_constants(i, j, p, u, x, BL_p, Tu, alpha_34, alphaE) TYPE(UA_ElementContinuousStateType), INTENT(IN ) :: x ! Continuous states at t TYPE(AFI_UA_BL_Type), INTENT(IN ) :: BL_p ! potentially interpolated UA parameters + REAL(R8Ki), INTENT( OUT) :: Tu REAL(ReKi), optional, INTENT( OUT) :: alpha_34 - REAL(ReKi), INTENT( OUT) :: Tu REAL(ReKi), optional, INTENT( OUT) :: alphaE - ! Local variables - real(ReKi) :: vx_34 - - ! Variables derived from inputs !u%u = U_ac = TwoNorm(u%v_ac) ! page 4 definitions Tu = Get_Tu(u%u, p%c(i,j)) @@ -2795,7 +2848,7 @@ SUBROUTINE UA_RK4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat, k1%x = p%dt * k1%x - x_tmp%x = x%element(i,j)%x + 0.5 * k1%x + x_tmp%x = x%element(i,j)%x + 0.5_R8Ki * k1%x ! interpolate u to find u_interp = u(t + dt/2) TPlusHalfDt = t + 0.5_DbKi*p%dt @@ -2808,7 +2861,7 @@ SUBROUTINE UA_RK4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat, k2%x = p%dt * k2%x - x_tmp%x = x%element(i,j)%x + 0.5 * k2%x + x_tmp%x = x%element(i,j)%x + 0.5_R8Ki * k2%x ! find xdot at t + dt/2 (note x_tmp has changed) CALL UA_CalcContStateDeriv( i, j, TPlusHalfDt, u_interp, p, x_tmp, OtherState, AFInfo, m, k3, ErrStat2, ErrMsg2 ) @@ -2828,7 +2881,7 @@ SUBROUTINE UA_RK4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat, k4%x = p%dt * k4%x - x%element(i,j)%x = x%element(i,j)%x + ( k1%x + 2. * k2%x + 2. * k3%x + k4%x ) / 6. + x%element(i,j)%x = x%element(i,j)%x + ( k1%x + 2.0_R8Ki * k2%x + 2.0_R8Ki * k3%x + k4%x ) / 6.0_R8Ki END SUBROUTINE UA_RK4 !---------------------------------------------------------------------------------------------------------------------------------- @@ -2915,9 +2968,8 @@ SUBROUTINE UA_AB4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat, IF ( ErrStat >= AbortErrLev ) RETURN else - x%element(i,j)%x = x%element(i,j)%x + p%DT/24. * ( 55.*OtherState%xdot(1)%element(i,j)%x - 59.*OtherState%xdot(2)%element(i,j)%x & - + 37.*OtherState%xdot(3)%element(i,j)%x - 9.*OtherState%xdot(4)%element(i,j)%x ) - + x%element(i,j)%x = x%element(i,j)%x + p%DT/24.0_R8Ki * ( 55.0_R8Ki*OtherState%xdot(1)%element(i,j)%x - 59.0_R8Ki*OtherState%xdot(2)%element(i,j)%x & + + 37.0_R8Ki*OtherState%xdot(3)%element(i,j)%x - 9.0_R8Ki*OtherState%xdot(4)%element(i,j)%x ) endif @@ -2995,16 +3047,208 @@ SUBROUTINE UA_ABM4( i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat, IF ( ErrStat >= AbortErrLev ) RETURN - x%element(i,j)%x = x_in%x + p%DT/24. * ( 9. * xdot_pred%x + 19. * OtherState%xdot(1)%element(i,j)%x & - - 5. * OtherState%xdot(2)%element(i,j)%x & - + 1. * OtherState%xdot(3)%element(i,j)%x ) + x%element(i,j)%x = x_in%x + p%DT/24.0_R8Ki * ( 9.0_R8Ki * xdot_pred%x + 19.0_R8Ki * OtherState%xdot(1)%element(i,j)%x & + - 5.0_R8Ki * OtherState%xdot(2)%element(i,j)%x & + + 1.0_R8Ki * OtherState%xdot(3)%element(i,j)%x ) endif END SUBROUTINE UA_ABM4 !---------------------------------------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------------------------------------- +!> This subroutine implements a Newton solve of the 2nd-order BDF system for numerically integrating ordinary differential equations: +SUBROUTINE UA_BDF2( i, j, t, n, u_interp, p, x, OtherState, AFInfo, m, ErrStat, ErrMsg ) +!.................................................................................................................................. + integer(IntKi), INTENT(IN ) :: i !< blade node counter + integer(IntKi), INTENT(IN ) :: j !< blade counter + REAL(DbKi), INTENT(IN ) :: t !< Current simulation time in seconds + INTEGER(IntKi), INTENT(IN ) :: n !< time step number + TYPE(UA_InputType), INTENT(IN ) :: u_interp !< Inputs at t + dt + TYPE(UA_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(UA_ContinuousStateType), INTENT(INOUT) :: x !< Continuous states at t on input at t + dt on output + TYPE(UA_OtherStateType), INTENT(INOUT) :: OtherState !< Other states + TYPE(AFI_ParameterType), INTENT(IN ) :: AFInfo ! The airfoil parameter data + TYPE(UA_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + ! local variables + + INTEGER(IntKi) :: k + INTEGER(IntKi), parameter :: KMax = 10 + REAL(R8Ki), parameter :: TolerSquared = (1D-6)**2 + + INTEGER(IntKi) :: ErrStat2 ! local error status + CHARACTER(ErrMsgLen) :: ErrMsg2 ! local error message (ErrMsg) + CHARACTER(*), PARAMETER :: RoutineName = 'UA_BDF2' + REAL(R8Ki) :: JMat(size(x%element(i,j)%x), size(x%element(i,j)%x)) + INTEGER :: iPivot(size(JMat,1)) + REAL(R8Ki) :: x_delta(size(JMat,1)) + REAL(R8Ki) :: x_constant(size(JMat,1)) + REAL(R8Ki) :: err + TYPE(UA_ElementContinuousStateType) :: xdot_pred ! Derivative of continuous states at t + + REAL(R8Ki), parameter :: Beta = 2.0_R8Ki/3.0_R8Ki + REAL(R8Ki), parameter :: Alpha0 = 4.0_R8Ki/3.0_R8Ki + REAL(R8Ki), parameter :: Alpha1 = -1.0_R8Ki/3.0_R8Ki + + !!!! for p=0, we get backward Euler: + !!!REAL(R8Ki), parameter :: Beta = 1.0_R8Ki + !!!REAL(R8Ki), parameter :: Alpha0 = 1.0_R8Ki + !!!REAL(R8Ki), parameter :: Alpha1 = 0.0_R8Ki + + + !NOTE: the error handling here assumes that we do not have any allocatable data in the inputs (u_interp) to be concerned with. + ! Also, We assume that if there is going to be an error in UA_CalcContStateDeriv, it will happen only on the first call + ! to the routine. + + ! Initialize ErrStat + + ErrStat = ErrID_None + ErrMsg = "" + + if (OtherState%n(i,j) < n) then + + OtherState%n(i,j) = n + + OtherState%xHistory(4)%element(i,j) = OtherState%xHistory(3)%element(i,j) + OtherState%xHistory(3)%element(i,j) = OtherState%xHistory(2)%element(i,j) + OtherState%xHistory(2)%element(i,j) = OtherState%xHistory(1)%element(i,j) + + if (n <= 1) then + OtherState%xHistory(2)%element(i,j) = x%element(i,j) + end if + + elseif (OtherState%n(i,j) > n) then + + CALL SetErrStat(ErrID_Fatal,'Backing up in time is not supported with a multistep method.',ErrStat,ErrMsg,RoutineName) + RETURN + + endif + + OtherState%xHistory(1)%element(i,j) = x%element(i,j) + + !!if (n<=1) then ! initialize because we don't have values for x + !! CALL UA_RK4(i, j, t, n, u, utimes, p, x, OtherState, AFInfo, m, ErrStat2, ErrMsg2 ) + !! CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) + !! RETURN + !!end if + + + x_constant = Alpha0 * x%element(i,j)%x + Alpha1 * OtherState%xHistory(2)%element(i,j)%x + + err = HUGE(err) + k = 0 + DO + + IF (K >= KMax) EXIT + + IF (K==0) THEN + ! This Jacobian will change when x changes, only if the values of x1, x2, or x3 are near boundaries of slope changes in + ! the FullyAttached function of the airfoil. At that point, it should be okay if the derivative is computed + ! on a slightly different region anyway. + CALL UA_Jacobian( i, j, t, n, u_interp, p, x, OtherState, AFInfo, m, Beta, JMat, ErrStat2, ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + + ! Get the LU decomposition of this matrix using a LAPACK routine: + ! The result is of the form JMat = P * L * U + + CALL LAPACK_getrf( M=size(JMat,1), N=size(JMat,2), A=JMat, IPIV=iPivot, ErrStat=ErrStat2, ErrMsg=ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF ( ErrStat >= AbortErrLev ) RETURN + END IF + + !------------------------------------------------------------------------------------------------- + ! Solve for delta x: JMat * x_delta = - F = - ( x(t+dt) - x(t) - dt * X(t+dt) + ! using the LAPACK routine + !------------------------------------------------------------------------------------------------- + CALL UA_CalcContStateDeriv( i, j, t, u_interp, p, x%element(i,j), OtherState, AFInfo, m, xdot_pred, ErrStat=ErrStat2, ErrMsg=ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF ( ErrStat >= AbortErrLev ) RETURN + + x_delta = - x%element(i,j)%x + x_constant + p%dt * Beta * xdot_pred%x + CALL LAPACK_getrs( TRANS="N", N=SIZE(JMat,1), A=JMat, & + IPIV=iPivot, B=x_delta, ErrStat=ErrStat2, ErrMsg=ErrMsg2 ) + CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) + IF ( ErrStat >= AbortErrLev ) RETURN + + !------------------------------------------------------------------------------------------------- + ! check for error, update inputs if necessary, and iterate again + !------------------------------------------------------------------------------------------------- + !err_prev = err + err = DOT_PRODUCT(x_delta, x_delta) + IF ( err <= TolerSquared) EXIT + IF (K >= KMax ) EXIT + + !!------------------------------------------------------------------------------------------------- + !! modify states for next iteration + !!------------------------------------------------------------------------------------------------- + !if (err > err_prev ) then + ! u_delta = u_delta * reduction_factor ! don't take a full step if we're getting farther from the solution! + ! err_prev = err_prev * reduction_factor + !end if + + x%element(i,j)%x = x%element(i,j)%x + x_delta + + + K = K + 1 + + END DO ! K + +END SUBROUTINE UA_BDF2 +!---------------------------------------------------------------------------------------------------------------------------------- +SUBROUTINE UA_Jacobian( i, j, t, n, u_interp, p, x, OtherState, AFInfo, m, Beta, JMat, ErrStat, ErrMsg ) + integer(IntKi), INTENT(IN ) :: i !< blade node counter + integer(IntKi), INTENT(IN ) :: j !< blade counter + REAL(DbKi), INTENT(IN ) :: t !< Current simulation time in seconds + INTEGER(IntKi), INTENT(IN ) :: n !< time step number + TYPE(UA_InputType), INTENT(IN ) :: u_interp !< Inputs at utimes + TYPE(UA_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(UA_ContinuousStateType), INTENT(INOUT) :: x !< Continuous states at t on input at t + dt on output + TYPE(UA_OtherStateType), INTENT(INOUT) :: OtherState !< Other states + TYPE(AFI_ParameterType), INTENT(IN ) :: AFInfo ! The airfoil parameter data + TYPE(UA_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + REAL(R8Ki), INTENT(IN ) :: Beta !< Value of Beta for p-th order BDF method + REAL(R8Ki), INTENT( OUT) :: JMat(:,:) !< Jacobian matrix + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + + TYPE(UA_ElementContinuousStateType) :: x_tmp ! Holds temporary modification to x + TYPE(UA_ElementContinuousStateType) :: X_p ! Holds derivative of X + TYPE(UA_ElementContinuousStateType) :: X_m ! Holds derivative of X + + INTEGER :: k + + ! Initialize ErrStat + + ErrStat = ErrID_None + ErrMsg = "" + + + ! compute JMat = I - dt*dXdx + + call eye(JMat, ErrStat, ErrMsg) + + x_tmp%x = x%element(i,j)%x + do k=1,size(p%dx) + x_tmp%x(k) = x%element(i,j)%x(k) + p%dx(k) + CALL UA_CalcContStateDeriv( i, j, t, u_interp, p, x_tmp, OtherState, AFInfo, m, X_p, ErrStat, ErrMsg ) + if (ErrStat >= AbortErrLev) return + + x_tmp%x(k) = x%element(i,j)%x(k) - p%dx(k) + CALL UA_CalcContStateDeriv( i, j, t, u_interp, p, x_tmp, OtherState, AFInfo, m, X_m, ErrStat, ErrMsg ) + if (ErrStat >= AbortErrLev) return + + ! reset + x_tmp%x(k) = x%element(i,j)%x(k) + + ! compute I(:,k) - dt * dXdx(:,k) + JMat(:,k) = JMat(:,k) - p%dt*Beta*(X_p%x - X_m%x) / (2.0_R8Ki * p%dx(k)) + end do + +END SUBROUTINE UA_Jacobian +!---------------------------------------------------------------------------------------------------------------------------------- !============================================================================== @@ -3050,7 +3294,7 @@ subroutine UA_CalcOutput( i, j, t, u_in, p, x, xd, OtherState, AFInfo, y, misc, ! for UA_HGM real(ReKi) :: alphaE - real(ReKi) :: Tu + real(R8Ki) :: Tu real(ReKi) :: alpha_34 real(ReKi) :: fs_aE real(ReKi) :: cl_fs diff --git a/modules/aerodyn/src/UnsteadyAero_Registry.txt b/modules/aerodyn/src/UnsteadyAero_Registry.txt index b71aef051d..d3eb0ef03a 100644 --- a/modules/aerodyn/src/UnsteadyAero_Registry.txt +++ b/modules/aerodyn/src/UnsteadyAero_Registry.txt @@ -161,6 +161,7 @@ typedef ^ OtherStateType ReKi typedef ^ OtherStateType ReKi sigma3 {:}{:} - - "multiplier for T_V" - typedef ^ OtherStateType IntKi n {:}{:} - - "counter for continuous state integration" - typedef ^ OtherStateType UA_ContinuousStateType xdot 4 - - "history states for continuous state integration" - +typedef ^ OtherStateType UA_ContinuousStateType xHistory 4 - - "history states for continuous state integration" - typedef ^ OtherStateType ReKi t_vortexBegin {:}{:} - - "HGMV model: simulation time when vortex lift term became active" s typedef ^ OtherStateType ReKi SignOfOmega {:}{:} - - "HGMV model: sign of omega when vortex lift term became active " s typedef ^ OtherStateType LOGICAL PositivePressure {:}{:} - - "HGMV model: logical flag indicating if the vortex lift became active because of positive pressure (or negative)" - @@ -202,6 +203,8 @@ typedef ^ ^ INTEGER typedef ^ ^ Logical ShedEffect - - - "Include the effect of shed vorticity. If False, the input alpha is assumed to already contain this effect (e.g. vortex methods)" - typedef ^ ParameterType IntKi lin_nx - 0 - "Number of continuous states for linearization" - typedef ^ ^ LOGICAL UA_off_forGood {:}{:} - - "logical flag indicating if UA is off for good" - +typedef ^ ^ INTEGER lin_xIndx {:}{:} - - "array to indicate which state to perturb for UA" - +typedef ^ ^ R8Ki dx {5} - - "array to indicate size of state perturbations" - # ..... Inputs .................................................................................................................... # Define inputs that are contained on the mesh here: diff --git a/modules/aerodyn/src/UnsteadyAero_Types.f90 b/modules/aerodyn/src/UnsteadyAero_Types.f90 index 3a42342068..eb5efeb6e6 100644 --- a/modules/aerodyn/src/UnsteadyAero_Types.f90 +++ b/modules/aerodyn/src/UnsteadyAero_Types.f90 @@ -181,6 +181,7 @@ MODULE UnsteadyAero_Types REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: sigma3 !< multiplier for T_V [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: n !< counter for continuous state integration [-] TYPE(UA_ContinuousStateType) , DIMENSION(1:4) :: xdot !< history states for continuous state integration [-] + TYPE(UA_ContinuousStateType) , DIMENSION(1:4) :: xHistory !< history states for continuous state integration [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: t_vortexBegin !< HGMV model: simulation time when vortex lift term became active [s] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: SignOfOmega !< HGMV model: sign of omega when vortex lift term became active [s] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: PositivePressure !< HGMV model: logical flag indicating if the vortex lift became active because of positive pressure (or negative) [-] @@ -221,6 +222,8 @@ MODULE UnsteadyAero_Types LOGICAL :: ShedEffect !< Include the effect of shed vorticity. If False, the input alpha is assumed to already contain this effect (e.g. vortex methods) [-] INTEGER(IntKi) :: lin_nx = 0 !< Number of continuous states for linearization [-] LOGICAL , DIMENSION(:,:), ALLOCATABLE :: UA_off_forGood !< logical flag indicating if UA is off for good [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: lin_xIndx !< array to indicate which state to perturb for UA [-] + REAL(R8Ki) , DIMENSION(1:5) :: dx !< array to indicate size of state perturbations [-] END TYPE UA_ParameterType ! ======================= ! ========= UA_InputType ======= @@ -4412,6 +4415,11 @@ SUBROUTINE UA_CopyOtherState( SrcOtherStateData, DstOtherStateData, CtrlCode, Er CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) IF (ErrStat>=AbortErrLev) RETURN ENDDO + DO i1 = LBOUND(SrcOtherStateData%xHistory,1), UBOUND(SrcOtherStateData%xHistory,1) + CALL UA_CopyContState( SrcOtherStateData%xHistory(i1), DstOtherStateData%xHistory(i1), CtrlCode, ErrStat2, ErrMsg2 ) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName) + IF (ErrStat>=AbortErrLev) RETURN + ENDDO IF (ALLOCATED(SrcOtherStateData%t_vortexBegin)) THEN i1_l = LBOUND(SrcOtherStateData%t_vortexBegin,1) i1_u = UBOUND(SrcOtherStateData%t_vortexBegin,1) @@ -4555,6 +4563,10 @@ SUBROUTINE UA_DestroyOtherState( OtherStateData, ErrStat, ErrMsg, DEALLOCATEpoin CALL UA_DestroyContState( OtherStateData%xdot(i1), ErrStat2, ErrMsg2, DEALLOCATEpointers_local ) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) ENDDO +DO i1 = LBOUND(OtherStateData%xHistory,1), UBOUND(OtherStateData%xHistory,1) + CALL UA_DestroyContState( OtherStateData%xHistory(i1), ErrStat2, ErrMsg2, DEALLOCATEpointers_local ) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) +ENDDO IF (ALLOCATED(OtherStateData%t_vortexBegin)) THEN DEALLOCATE(OtherStateData%t_vortexBegin) ENDIF @@ -4663,6 +4675,25 @@ SUBROUTINE UA_PackOtherState( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs DEALLOCATE(Int_Buf) END IF END DO + DO i1 = LBOUND(InData%xHistory,1), UBOUND(InData%xHistory,1) + Int_BufSz = Int_BufSz + 3 ! xHistory: size of buffers for each call to pack subtype + CALL UA_PackContState( Re_Buf, Db_Buf, Int_Buf, InData%xHistory(i1), ErrStat2, ErrMsg2, .TRUE. ) ! xHistory + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + IF (ErrStat >= AbortErrLev) RETURN + + IF(ALLOCATED(Re_Buf)) THEN ! xHistory + Re_BufSz = Re_BufSz + SIZE( Re_Buf ) + DEALLOCATE(Re_Buf) + END IF + IF(ALLOCATED(Db_Buf)) THEN ! xHistory + Db_BufSz = Db_BufSz + SIZE( Db_Buf ) + DEALLOCATE(Db_Buf) + END IF + IF(ALLOCATED(Int_Buf)) THEN ! xHistory + Int_BufSz = Int_BufSz + SIZE( Int_Buf ) + DEALLOCATE(Int_Buf) + END IF + END DO Int_BufSz = Int_BufSz + 1 ! t_vortexBegin allocated yes/no IF ( ALLOCATED(InData%t_vortexBegin) ) THEN Int_BufSz = Int_BufSz + 2*2 ! t_vortexBegin upper/lower bounds for each dimension @@ -4875,6 +4906,36 @@ SUBROUTINE UA_PackOtherState( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMs IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1 ENDIF END DO + DO i1 = LBOUND(InData%xHistory,1), UBOUND(InData%xHistory,1) + CALL UA_PackContState( Re_Buf, Db_Buf, Int_Buf, InData%xHistory(i1), ErrStat2, ErrMsg2, OnlySize ) ! xHistory + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + IF (ErrStat >= AbortErrLev) RETURN + + IF(ALLOCATED(Re_Buf)) THEN + IntKiBuf( Int_Xferred ) = SIZE(Re_Buf); Int_Xferred = Int_Xferred + 1 + IF (SIZE(Re_Buf) > 0) ReKiBuf( Re_Xferred:Re_Xferred+SIZE(Re_Buf)-1 ) = Re_Buf + Re_Xferred = Re_Xferred + SIZE(Re_Buf) + DEALLOCATE(Re_Buf) + ELSE + IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1 + ENDIF + IF(ALLOCATED(Db_Buf)) THEN + IntKiBuf( Int_Xferred ) = SIZE(Db_Buf); Int_Xferred = Int_Xferred + 1 + IF (SIZE(Db_Buf) > 0) DbKiBuf( Db_Xferred:Db_Xferred+SIZE(Db_Buf)-1 ) = Db_Buf + Db_Xferred = Db_Xferred + SIZE(Db_Buf) + DEALLOCATE(Db_Buf) + ELSE + IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1 + ENDIF + IF(ALLOCATED(Int_Buf)) THEN + IntKiBuf( Int_Xferred ) = SIZE(Int_Buf); Int_Xferred = Int_Xferred + 1 + IF (SIZE(Int_Buf) > 0) IntKiBuf( Int_Xferred:Int_Xferred+SIZE(Int_Buf)-1 ) = Int_Buf + Int_Xferred = Int_Xferred + SIZE(Int_Buf) + DEALLOCATE(Int_Buf) + ELSE + IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1 + ENDIF + END DO IF ( .NOT. ALLOCATED(InData%t_vortexBegin) ) THEN IntKiBuf( Int_Xferred ) = 0 Int_Xferred = Int_Xferred + 1 @@ -5227,6 +5288,50 @@ SUBROUTINE UA_UnPackOtherState( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, Er IF(ALLOCATED(Db_Buf )) DEALLOCATE(Db_Buf ) IF(ALLOCATED(Int_Buf)) DEALLOCATE(Int_Buf) END DO + i1_l = LBOUND(OutData%xHistory,1) + i1_u = UBOUND(OutData%xHistory,1) + DO i1 = LBOUND(OutData%xHistory,1), UBOUND(OutData%xHistory,1) + Buf_size=IntKiBuf( Int_Xferred ) + Int_Xferred = Int_Xferred + 1 + IF(Buf_size > 0) THEN + ALLOCATE(Re_Buf(Buf_size),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating Re_Buf.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + Re_Buf = ReKiBuf( Re_Xferred:Re_Xferred+Buf_size-1 ) + Re_Xferred = Re_Xferred + Buf_size + END IF + Buf_size=IntKiBuf( Int_Xferred ) + Int_Xferred = Int_Xferred + 1 + IF(Buf_size > 0) THEN + ALLOCATE(Db_Buf(Buf_size),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating Db_Buf.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + Db_Buf = DbKiBuf( Db_Xferred:Db_Xferred+Buf_size-1 ) + Db_Xferred = Db_Xferred + Buf_size + END IF + Buf_size=IntKiBuf( Int_Xferred ) + Int_Xferred = Int_Xferred + 1 + IF(Buf_size > 0) THEN + ALLOCATE(Int_Buf(Buf_size),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating Int_Buf.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + Int_Buf = IntKiBuf( Int_Xferred:Int_Xferred+Buf_size-1 ) + Int_Xferred = Int_Xferred + Buf_size + END IF + CALL UA_UnpackContState( Re_Buf, Db_Buf, Int_Buf, OutData%xHistory(i1), ErrStat2, ErrMsg2 ) ! xHistory + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + IF (ErrStat >= AbortErrLev) RETURN + + IF(ALLOCATED(Re_Buf )) DEALLOCATE(Re_Buf ) + IF(ALLOCATED(Db_Buf )) DEALLOCATE(Db_Buf ) + IF(ALLOCATED(Int_Buf)) DEALLOCATE(Int_Buf) + END DO IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! t_vortexBegin not allocated Int_Xferred = Int_Xferred + 1 ELSE @@ -5991,6 +6096,21 @@ SUBROUTINE UA_CopyParam( SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg ) END IF DstParamData%UA_off_forGood = SrcParamData%UA_off_forGood ENDIF +IF (ALLOCATED(SrcParamData%lin_xIndx)) THEN + i1_l = LBOUND(SrcParamData%lin_xIndx,1) + i1_u = UBOUND(SrcParamData%lin_xIndx,1) + i2_l = LBOUND(SrcParamData%lin_xIndx,2) + i2_u = UBOUND(SrcParamData%lin_xIndx,2) + IF (.NOT. ALLOCATED(DstParamData%lin_xIndx)) THEN + ALLOCATE(DstParamData%lin_xIndx(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%lin_xIndx.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + END IF + DstParamData%lin_xIndx = SrcParamData%lin_xIndx +ENDIF + DstParamData%dx = SrcParamData%dx END SUBROUTINE UA_CopyParam SUBROUTINE UA_DestroyParam( ParamData, ErrStat, ErrMsg, DEALLOCATEpointers ) @@ -6019,6 +6139,9 @@ SUBROUTINE UA_DestroyParam( ParamData, ErrStat, ErrMsg, DEALLOCATEpointers ) ENDIF IF (ALLOCATED(ParamData%UA_off_forGood)) THEN DEALLOCATE(ParamData%UA_off_forGood) +ENDIF +IF (ALLOCATED(ParamData%lin_xIndx)) THEN + DEALLOCATE(ParamData%lin_xIndx) ENDIF END SUBROUTINE UA_DestroyParam @@ -6081,6 +6204,12 @@ SUBROUTINE UA_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si Int_BufSz = Int_BufSz + 2*2 ! UA_off_forGood upper/lower bounds for each dimension Int_BufSz = Int_BufSz + SIZE(InData%UA_off_forGood) ! UA_off_forGood END IF + Int_BufSz = Int_BufSz + 1 ! lin_xIndx allocated yes/no + IF ( ALLOCATED(InData%lin_xIndx) ) THEN + Int_BufSz = Int_BufSz + 2*2 ! lin_xIndx upper/lower bounds for each dimension + Int_BufSz = Int_BufSz + SIZE(InData%lin_xIndx) ! lin_xIndx + END IF + Db_BufSz = Db_BufSz + SIZE(InData%dx) ! dx IF ( Re_BufSz .GT. 0 ) THEN ALLOCATE( ReKiBuf( Re_BufSz ), STAT=ErrStat2 ) IF (ErrStat2 /= 0) THEN @@ -6182,6 +6311,30 @@ SUBROUTINE UA_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, Si END DO END DO END IF + IF ( .NOT. ALLOCATED(InData%lin_xIndx) ) THEN + IntKiBuf( Int_Xferred ) = 0 + Int_Xferred = Int_Xferred + 1 + ELSE + IntKiBuf( Int_Xferred ) = 1 + Int_Xferred = Int_Xferred + 1 + IntKiBuf( Int_Xferred ) = LBOUND(InData%lin_xIndx,1) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%lin_xIndx,1) + Int_Xferred = Int_Xferred + 2 + IntKiBuf( Int_Xferred ) = LBOUND(InData%lin_xIndx,2) + IntKiBuf( Int_Xferred + 1) = UBOUND(InData%lin_xIndx,2) + Int_Xferred = Int_Xferred + 2 + + DO i2 = LBOUND(InData%lin_xIndx,2), UBOUND(InData%lin_xIndx,2) + DO i1 = LBOUND(InData%lin_xIndx,1), UBOUND(InData%lin_xIndx,1) + IntKiBuf(Int_Xferred) = InData%lin_xIndx(i1,i2) + Int_Xferred = Int_Xferred + 1 + END DO + END DO + END IF + DO i1 = LBOUND(InData%dx,1), UBOUND(InData%dx,1) + DbKiBuf(Db_Xferred) = InData%dx(i1) + Db_Xferred = Db_Xferred + 1 + END DO END SUBROUTINE UA_PackParam SUBROUTINE UA_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg ) @@ -6292,6 +6445,35 @@ SUBROUTINE UA_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg END DO END DO END IF + IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN ! lin_xIndx not allocated + Int_Xferred = Int_Xferred + 1 + ELSE + Int_Xferred = Int_Xferred + 1 + i1_l = IntKiBuf( Int_Xferred ) + i1_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + i2_l = IntKiBuf( Int_Xferred ) + i2_u = IntKiBuf( Int_Xferred + 1) + Int_Xferred = Int_Xferred + 2 + IF (ALLOCATED(OutData%lin_xIndx)) DEALLOCATE(OutData%lin_xIndx) + ALLOCATE(OutData%lin_xIndx(i1_l:i1_u,i2_l:i2_u),STAT=ErrStat2) + IF (ErrStat2 /= 0) THEN + CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%lin_xIndx.', ErrStat, ErrMsg,RoutineName) + RETURN + END IF + DO i2 = LBOUND(OutData%lin_xIndx,2), UBOUND(OutData%lin_xIndx,2) + DO i1 = LBOUND(OutData%lin_xIndx,1), UBOUND(OutData%lin_xIndx,1) + OutData%lin_xIndx(i1,i2) = IntKiBuf(Int_Xferred) + Int_Xferred = Int_Xferred + 1 + END DO + END DO + END IF + i1_l = LBOUND(OutData%dx,1) + i1_u = UBOUND(OutData%dx,1) + DO i1 = LBOUND(OutData%dx,1), UBOUND(OutData%dx,1) + OutData%dx(i1) = REAL(DbKiBuf(Db_Xferred), R8Ki) + Db_Xferred = Db_Xferred + 1 + END DO END SUBROUTINE UA_UnPackParam SUBROUTINE UA_CopyInput( SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg )