Skip to content

Commit

Permalink
AD: Using lin_xIndx for UA linearization
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebranlard
ebranlard committed Jul 18, 2023
1 parent 601abf4 commit 8ce7e8d
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 121 deletions.
173 changes: 53 additions & 120 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6267,34 +6267,14 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg,

if (p%BEMT%UA%lin_nx>0) then
! --- ADENV
!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


! --- ADLEG
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) ! state, State=4 for Oye, otherwise state \in [1,4]
x_op(index) = x%BEMT%UA%element(i,j)%x(k)

index = index + 1
end do

end if

Expand Down Expand Up @@ -6347,33 +6327,14 @@ SUBROUTINE RotGetOP( t, u, p, p_AD, x, xd, z, OtherState, y, m, ErrStat, ErrMsg,
end if

if (p%BEMT%UA%lin_nx>0) then
! --- ADENV
!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) = dxdt%BEMT%UA%element(i,j)%x(k)
!
! index = index + 1
!end do
! --- ADLEG
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) ! state, State=4 for Oye, otherwise state \in [1,4]
x_op(index) = dxdt%BEMT%UA%element(i,j)%x(k)

index = index + 1
end do
end if

if (p%BEMT%lin_nx>0) then
Expand Down Expand Up @@ -6942,24 +6903,24 @@ SUBROUTINE Init_Jacobian_x( p, InitOut, ErrStat, ErrMsg)

k = 1 + p%BEMT%DBEMT%lin_nx
! ---- ADENV
!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)

! 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
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) ! State=4 for Oye, otherwise state \in [1,4]

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

! ---- ADLEG
do j=1,p%NumBlades ! size(x%BEMT%DBEMT%element,2)
Expand Down Expand Up @@ -7181,28 +7142,19 @@ SUBROUTINE Perturb_x( p, n, perturb_sign, x, dx )


! --- ADENV
!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
! StateIndex = n_tmp - p%BEMT%UA%lin_nx
! x%BEMT%V_w(StateIndex) = x%BEMT%V_w(StateIndex) + dx * perturb_sign
!end if

! --- AD LEG
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, State=4 for Oye, otherwise state \in [1,4]

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
Expand Down Expand Up @@ -7305,34 +7257,15 @@ SUBROUTINE Compute_dX(p, x_p, x_m, delta_p, delta_m, dX)
end if

if (p%BEMT%UA%lin_nx>0) then
! --- ADENV
!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)
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) ! state, State=4 for Oye, otherwise state \in [1,4]
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
indx_first = indx_first + 1
end do

! --- ADLEG

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

end if

if (p%BEMT%lin_nx>0) then ! skewWake
Expand Down
7 changes: 6 additions & 1 deletion modules/aerodyn/src/UnsteadyAero.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -821,7 +821,12 @@ subroutine UA_SetParameters( dt, InitInp, p, AFInfo, AFIndx, ErrStat, ErrMsg )
do k=1,UA_NumLinStates
p%lin_xIndx(n,1) = i ! node
p%lin_xIndx(n,2) = j ! blade
p%lin_xIndx(n,3) = k ! state

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
Expand Down

0 comments on commit 8ce7e8d

Please sign in to comment.