From f593904cb7dfb3d787e58f5b810a3423bc0cd200 Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Thu, 11 May 2023 11:55:38 -0600 Subject: [PATCH 1/2] AD: Fix divide-by-zero errors --- modules/aerodyn/src/AeroDyn.f90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index d78fff7c1b..0349c719b1 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -2559,13 +2559,23 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) ! Calculate Yaw and Tilt for use in xVelCorr ! Define a vector wrt which the yaw is defined - x_hat_wind = m%V_diskAvg/twonorm(m%V_diskAvg) - ! Yaw + denom = twonorm(m%V_diskAvg) + if (EqualRealNos(denom, 0.0_ReKi)) then + x_hat_wind = 0.0_ReKi + else + x_hat_wind = m%V_diskAvg/denom + end if + ! Yaw tmpD = x_hat_disk tmpD(3) = 0.0 tmpW = x_hat_wind tmpW(3) = 0.0 - yaw = acos(max(-1.0_ReKi,min(1.0_ReKi,dot_product(tmpD,tmpW)/(twonorm(tmpD)*TwoNorm(tmpW))))) + denom = TwoNorm(tmpD)*TwoNorm(tmpW) + if (EqualRealNos(denom, 0.0_ReKi)) then + yaw = 0.0_ReKi + else + yaw = acos(max(-1.0_ReKi,min(1.0_ReKi,dot_product(tmpD,tmpW)/denom))) + end if tmp_skewVec = cross_product(tmpW,tmpD); yaw = sign(yaw,tmp_skewVec(3)) m%Yaw = yaw From c283aa3a355c95848f3e282eb9f7499f06ce3c9b Mon Sep 17 00:00:00 2001 From: Bonnie Jonkman Date: Thu, 11 May 2023 12:18:26 -0600 Subject: [PATCH 2/2] more divide-by-zero checks --- modules/aerodyn/src/AeroDyn.f90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 0349c719b1..2ef0eed64f 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -2585,7 +2585,13 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) tmpD(2) = 0.0 tmpW = x_hat_wind tmpW(2) = 0.0 - tilt = acos(max(-1.0_ReKi,min(1.0_ReKi,dot_product(tmpD,tmpW)/(twonorm(tmpD)*TwoNorm(tmpW))))) + denom = TwoNorm(tmpD)*TwoNorm(tmpW) + if (EqualRealNos(denom, 0.0_Reki)) then + tilt = 0.0_ReKi + else + tilt = acos(max(-1.0_ReKi,min(1.0_ReKi,dot_product(tmpD,tmpW)/denom))) + end if + tmp_skewVec = cross_product(tmpW,tmpD) tilt = sign(tilt,tmp_skewVec(2)) m%tilt = tilt @@ -2709,7 +2715,12 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) dr(:) = elemPosRotorProj(:,j) - elemPosRotorProj(:,j-1) dz(:) = elemPosRelToHub(:,j) - elemPosRelToHub(:,j-1) - m%BEMT_u(indx)%drdz(j,k) = TwoNorm(dr(:)) / TwoNorm(dz(:)) + denom = TwoNorm(dz(:)) + if (EqualRealNos(denom,0.0_ReKi)) then ! this should not happen, but we'll check anyway + m%BEMT_u(indx)%drdz(j,k) = 0.0_ReKi + else + m%BEMT_u(indx)%drdz(j,k) = TwoNorm(dr(:)) / denom + end if end do ! j m%BEMT_u(indx)%drdz(1,k) = m%BEMT_u(indx)%drdz(2,k) end do !k=blades