Skip to content

Performance improvements in dev-tc #3125

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
deslaughter wants to merge 11 commits into
base: dev-tc
Choose a base branch
from
Open

Performance improvements in dev-tc #3125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
73173d0
Remove unused error handling from DCM_logmap
deslaughter Dec 30, 2025
92441db
Add optional flag to Morison_CalcOutput to skip expensive hydrostatic…
deslaughter Dec 30, 2025
dade1a0
Modified SetErrStat to overwrite contents of ErrMess if ErrStat is Er…
deslaughter Dec 30, 2025
228d47e
Remove some instances of ErrMsg = ""
deslaughter Dec 30, 2025
3721499
Restructure getWaterKin in attempt to be more performant
deslaughter Dec 30, 2025
682224b
Fix bad refactored logic in getWaterKin
deslaughter Dec 31, 2025
dea020b
Minor optimization in FVW ui_seg (~6%)
deslaughter Dec 31, 2025
eddeb74
Remove recursion from Morison RefineElementHstLds_Cyl
deslaughter Dec 31, 2025
9715b1e
Minor performance improvements in OLAF
deslaughter Dec 31, 2025
a50f2cf
Apply suggestions from code review
deslaughter Dec 31, 2025
d14cdf4
Add stack size check in RefineElementHstLds_Cyl
deslaughter Dec 31, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
188 changes: 100 additions & 88 deletions modules/aerodyn/src/FVW_BiotSavart.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,55 +50,61 @@ subroutine ui_seg_11(DeltaPa, DeltaPb, SegGamma, RegFunction, RegParam1, Uind)
real(ReKi) :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb
real(ReKi) :: exp_value !<
!
Uind(1:3)=0.0_ReKi
Uind(1:3) = 0.0_ReKi
xa=DeltaPa(1); ya=DeltaPa(2); za=DeltaPa(3)
xb=DeltaPb(1); yb=DeltaPb(2); zb=DeltaPb(3)
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb) ! |r1|*|r2|*(|r1|*|r2| + r1.r2)
if (denominator>PRECISION_UI) then
crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi
norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
if (norm2_r0>PRECISION_UI) then ! segment of zero length
! --- Far field TODO
! --- Regularization (close field)
norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
select case (RegFunction) !
case ( idRegNone ) ! No vortex core model
Kv=1.0_ReKi
case ( idRegRankine ) ! Rankine
r_bar2 = norm2_orth/ RegParam1**2
if (r_bar2<1) then
Kv=r_bar2
else
Kv=1.0_ReKi
end if
case ( idRegLambOseen ) ! Lamb-Oseen
r_bar2 = norm2_orth/ RegParam1**2
exp_value = -1.25643_ReKi*r_bar2
if(exp_value<MIN_EXP_VALUE) then ! Remove me when Far distance implemented
Kv = 1.0_ReKi
else
Kv = 1.0_ReKi-exp(exp_value)
endif
case ( idRegVatistas ) ! Vatistas n=2
r_bar2 = norm2_orth/ RegParam1**2
Kv = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
case ( idRegOffset ) ! Cut-off radius
Kv = 1.0_ReKi
denominator=denominator+RegParam1**2*norm2_r0
case default
print*,'Unknown SgmtReg', RegFunction
STOP ! Will never happen
Kv=1.0_ReKi !< Should be an error
end select
Kv=SegGamma*fourpi_inv*Kv*(norm_a+norm_b)/denominator
Uind(1:3) = Kv*crossprod(1:3)
endif

if (denominator <= PRECISION_UI) return
crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2

! On the singularity, Uind(1:3)=0.0_ReKi
if (norm2_orth <= PRECISION_UI) return
norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) + (za-zb)*(za-zb)

! segment of zero length
if (norm2_r0 <= PRECISION_UI) return

! --- Far field TODO
! --- Regularization (close field)
norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0

select case (RegFunction) !
case ( idRegNone ) ! No vortex core model
Kv=1.0_ReKi
case ( idRegRankine ) ! Rankine
r_bar2 = norm2_orth/ RegParam1**2
if (r_bar2<1) then
Kv=r_bar2
else
Kv=1.0_ReKi
end if
case ( idRegLambOseen ) ! Lamb-Oseen
r_bar2 = norm2_orth/ RegParam1**2
exp_value = -1.25643_ReKi*r_bar2
if(exp_value<MIN_EXP_VALUE) then ! Remove me when Far distance implemented
Kv = 1.0_ReKi
else
Kv = 1.0_ReKi-exp(exp_value)
endif
endif
case ( idRegVatistas ) ! Vatistas n=2
r_bar2 = norm2_orth/ RegParam1**2
Kv = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
case ( idRegOffset ) ! Cut-off radius
Kv = 1.0_ReKi
denominator=denominator+RegParam1**2*norm2_r0
case default
print*,'Unknown SgmtReg', RegFunction
STOP ! Will never happen
Kv=1.0_ReKi !< Should be an error
end select

Kv = SegGamma*fourpi_inv*Kv*(norm_a+norm_b)/denominator
Uind(1:3) = Kv*crossprod(1:3)

end subroutine ui_seg_11


Expand Down Expand Up @@ -134,21 +140,23 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
real(ReKi) :: norm2_r0 !< Squared length of the segment d = (r1xr2)/r0
real(ReKi) :: xa, ya, za, xb, yb, zb !< Coordinates of X-Xa and X-Xb
real(ReKi) :: exp_value !<
real(ReKi) :: CPs_icp(3) !<

! Branching based on regularization model
! NOTE: copy paste of code is done for optimization!
! The only thing changing is the part labelled "regularization"
select case (RegFunction)
case ( idRegNone ) ! No vortex core
!$OMP PARALLEL default(shared)
!$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
!$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
do icp=iCPStart,iCPEnd ! loop on CPs
Uind = 0.0_ReKi
CPs_icp = CPs(:,icp)
do is=iSegStart,iSegEnd ! loop on selected segments
Uind = 0.0_ReKi
P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
P2 = SegPoints(1:3, SegConnct(2,is))
xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
Expand All @@ -162,26 +170,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
! --- Far field TODO
! --- NO Regularization (close field)
Kv = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Kv*crossprod(1:3)
Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
end if
end if
end if
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
end do ! Loop on segments
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
enddo ! Loop on control points
!$OMP END DO
!$OMP END PARALLEL

case ( idRegRankine ) ! Rankine
!$OMP PARALLEL default(shared)
!$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
!$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
do icp=iCPStart,iCPEnd ! loop on CPs
Uind = 0.0_ReKi
CPs_icp = CPs(:,icp)
do is=iSegStart,iSegEnd ! loop on selected segments
Uind = 0.0_ReKi
P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
P2 = SegPoints(1:3, SegConnct(2,is))
xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
Expand All @@ -201,26 +210,27 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
Kv=1.0_ReKi
end if
Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Kv*crossprod(1:3)
Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
end if
end if ! denominator size or distances too small
end if !
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
end do ! Loop on segments
Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
enddo ! Loop on control points
!$OMP END DO
!$OMP END PARALLEL

case ( idRegLambOseen ) ! LambOseen
!$OMP PARALLEL default(shared)
!$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime)
!$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb,exp_value) schedule(runtime)
do icp=iCPStart,iCPEnd ! loop on CPs
Uind = 0.0_ReKi
CPs_icp = CPs(:,icp)
do is=iSegStart,iSegEnd ! loop on selected segments
Uind = 0.0_ReKi
P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
P2 = SegPoints(1:3, SegConnct(2,is))
xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
Expand All @@ -241,61 +251,63 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
Kv = 1.0_ReKi-exp(exp_value)
endif
Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Kv*crossprod(1:3)
Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
endif
end if ! denominator size or distances too small
end if !
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
end do ! Loop on segments
Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
enddo ! Loop on control points
!$OMP END DO
!$OMP END PARALLEL

case ( idRegVatistas ) ! Vatistas n=2
!$OMP PARALLEL default(shared)
!$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
!$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
do icp=iCPStart,iCPEnd ! loop on CPs
Uind = 0.0_ReKi
CPs_icp = CPs(:,icp)
do is=iSegStart,iSegEnd ! loop on selected segments
Uind = 0.0_ReKi
P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
P2 = SegPoints(1:3, SegConnct(2,is))
xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
if (denominator>PRECISION_UI) then
crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
if (norm2_orth>PRECISION_UI) then ! On the singularity, Uind(1:3)=0.0_ReKi
norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
if (norm2_r0>PRECISION_UI) then
! --- Far field TODO
! --- Regularization (close field) --- Vatistas
norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
r_bar2 = norm2_orth/RegParam(is)**2
Kv = r_bar2/sqrt(1+r_bar2**2)
Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Kv*crossprod(1:3)
end if
end if ! denominator size or distances too small
end if !
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
! denominator size or distances too small
if (denominator <= PRECISION_UI) cycle
crossprod(1) = ya*zb-za*yb; crossprod(2) = za*xb-xa*zb; crossprod(3) = xa*yb-ya*xb
norm2_orth = crossprod(1)**2 + crossprod(2)**2 + crossprod(3)**2
! On the singularity, cycle
if (norm2_orth <= PRECISION_UI) cycle
norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
! segment of zero length
if (norm2_r0 <= PRECISION_UI) cycle
! --- Far field TODO
! --- Regularization (close field) --- Vatistas
norm2_orth = norm2_orth/norm2_r0 ! d = (r1xr2)/r0
r_bar2 = norm2_orth/RegParam(is)**2
Kv = r_bar2/sqrt(1.0_ReKi+r_bar2**2)
Kv = SegGamma(is)*fourpi_inv*Kv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
end do ! Loop on segments
Uind_out(1:3,icp) = Uind_out(1:3,icp) + Uind(1:3)
enddo ! Loop on control points
!$OMP END DO
!$OMP END PARALLEL

case ( idRegOffset ) ! Denominator offset
!$OMP PARALLEL default(shared)
!$OMP do private(icp,is,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
!$OMP do private(icp,is,CPs_icp,Uind,P1,P2,crossprod,denominator,r_bar2,Kv,norm_a,norm_b,norm2_r0,norm2_orth,xa,ya,za,xb,yb,zb) schedule(runtime)
do icp=iCPStart,iCPEnd ! loop on CPs
Uind = 0.0_ReKi
CPs_icp = CPs(:,icp)
do is=iSegStart,iSegEnd ! loop on selected segments
Uind = 0.0_ReKi
P1 = SegPoints(1:3, SegConnct(1,is)) ! Segment extremity points
P2 = SegPoints(1:3, SegConnct(2,is))
xa=CPs(1,icp)-P1(1); ya=CPs(2,icp)-P1(2); za=CPs(3,icp)-P1(3);
xb=CPs(1,icp)-P2(1); yb=CPs(2,icp)-P2(2); zb=CPs(3,icp)-P2(3);
xa=CPs_icp(1)-P1(1); ya=CPs_icp(2)-P1(2); za=CPs_icp(3)-P1(3);
xb=CPs_icp(1)-P2(1); yb=CPs_icp(2)-P2(2); zb=CPs_icp(3)-P2(3);
norm_a = sqrt(xa*xa + ya*ya + za*za)
norm_b = sqrt(xb*xb + yb*yb + zb*zb)
denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
Expand All @@ -309,12 +321,12 @@ subroutine ui_seg(iCPStart, iCPEnd, CPs, &
! --- Regularization (close field) -- Offset
denominator = denominator+RegParam(is)**2*norm2_r0
Kv = SegGamma(is)*fourpi_inv*(norm_a+norm_b)/(denominator + MINDENOM)
Uind(1:3) = Kv*crossprod(1:3)
Uind(1:3) = Uind(1:3) + Kv*crossprod(1:3)
end if
end if ! denominator size or distances too small
end if !
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
end do ! Loop on segments
Uind_out(1:3,icp) = Uind_out(1:3,icp)+Uind(1:3)
enddo ! Loop on control points
!$OMP END DO
!$OMP END PARALLEL
Expand Down
16 changes: 9 additions & 7 deletions modules/aerodyn/src/FVW_Wings.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -281,23 +281,25 @@ subroutine Wings_ComputeCirculation(t, z, z_prev, p, x, m, AFInfo, ErrStat, ErrM
GammaScale=1.0_ReKi
endif

if (p%CircSolvMethod==idCircPrescribed) then
select case (p%CircSolvMethod)
! Prescribed circulation
case (idCircPrescribed)
do iW = 1, p%nWings !Loop over lifting lines
z%W(iW)%Gamma_LL(1:p%W(iW)%nSpan) = p%W(iW)%PrescribedCirculation(1:p%W(iW)%nSpan)
m%W(iW)%Vind_CP=-9999._ReKi !< Safety
m%W(iW)%Vtot_CP=-9999._ReKi !< Safety
enddo

else if (p%CircSolvMethod==idCircPolarData) then
! --- Solve for circulation using polar data
! Solve for circulation using polar data
case (idCircPolarData)
CALL Wings_ComputeCirculationPolarData(z, z_prev, p, x, m, AFInfo, GammaScale, ErrStat, ErrMsg, iLabel)

else if (p%CircSolvMethod==idCircNoFlowThrough) then
case (idCircNoFlowThrough)
! --- Solve for circulation using the no-flow through condition
ErrMsg='Circulation method nor implemented'; ErrStat=ErrID_Fatal; return ! should never happen
else
case default
ErrMsg='Circulation method nor implemented'; ErrStat=ErrID_Fatal; return ! should never happen
endif
end select

! Scale circulation (for initial transient)
do iW = 1, p%nWings !Loop over lifting lines
Expand Down Expand Up @@ -428,9 +430,9 @@ subroutine Wings_ComputeCirculationPolarData(z, z_prev, p, x, m, AFInfo, GammaSc
do iDepth=1,p%iNWStart ! Two first panels
! --- Defining a ring
P1=x%W(iW)%r_NW(1:3,iSpan ,iDepth )
P4=x%W(iW)%r_NW(1:3,iSpan ,iDepth+1)
P2=x%W(iW)%r_NW(1:3,iSpan+1,iDepth )
P3=x%W(iW)%r_NW(1:3,iSpan+1,iDepth+1)
P4=x%W(iW)%r_NW(1:3,iSpan ,iDepth+1)
! --- Induced velocity from ring, on all other control points (have to loop on rotors and wings and span again)
kCP2=1
do iWCP=1,p%nWings
Expand Down
Loading