diff --git a/.gitignore b/.gitignore index 74be88339b..56647ad7a2 100644 --- a/.gitignore +++ b/.gitignore @@ -69,3 +69,12 @@ openfast_io/openfast_io/_version.py *.code-workspace docker-compose.yml .devcontainer/ + +# AI steering files for specific users +.kiro/* + +# vs-build files +.vs +*.user +*.u2d +gitVersionInfo.h diff --git a/modules/aerodyn/src/AeroAcoustics.f90 b/modules/aerodyn/src/AeroAcoustics.f90 index 708ab8fe41..4b8c52fe91 100644 --- a/modules/aerodyn/src/AeroAcoustics.f90 +++ b/modules/aerodyn/src/AeroAcoustics.f90 @@ -35,7 +35,6 @@ module AeroAcoustics use AeroAcoustics_Types use AeroAcoustics_IO use NWTC_LAPACK - USE NWTC_FFTPACK implicit none private @@ -49,6 +48,9 @@ module AeroAcoustics REAL(ReKi), parameter :: AA_u_min = 0.1_ReKi REAL(ReKi), parameter :: AA_EPSILON = 1.E-16 ! EPSILON(AA_EPSILON) + REAL(ReKi), parameter :: AA_max_exp = log10(huge(1.0_ReKi))-8.0 ! bjj: I picked 8 here somewhat randomly...in debugging, later calculations + ! with the results of this function multiplied by factors of 1E6 + ! so I removed a couple more orders of magnitude (1E8) REAL(ReKi), parameter :: RotorRegionAlph_delta = 60.0_ReKi ! degrees : size of bin, must be a number that evenly divides 360 degrees REAL(ReKi), parameter :: RotorRegionRad_delta = 5.0_ReKi ! meters : size of bin along blade span (rotor radius) @@ -204,8 +206,8 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ! if tno is on or one of the guidati models is on, check if we have airfoil coordinates DO k=1,size(AFInfo) ! if any of the airfoil coordinates are missing change calculation method IF( AFInfo(k)%NumCoords .lt. 5 )then - CALL WrScr( 'Airfoil coordinates are missing: If Full or Simplified Guidati or Bl Calculation is on coordinates are needed ' ) - CALL WrScr( 'Calculation methods enforced as BPM for TBLTE and only Amiet for inflow ' ) + CALL WrScr( 'Airfoil coordinates are missing: If Full or Simplified Guidati or Bl Calculation are used, coordinates are needed.' ) + CALL WrScr( 'Calculation methods have been changed: BPM for TBLTE and only Amiet (BPM) for inflow ' ) p%ITURB = ITURB_BPM p%IInflow = IInflow_BPM exit ! stop checking do loop @@ -272,20 +274,32 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) p%BlSpn = InitInp%BlSpn p%BlChord = InitInp%BlChord - ! Calculate last element size as percentage of blade span + ! Calculate last element size as percentage of blade span, and warn if BldPrcnt is smaller than that. + ! bjj: This check only has meaning for a blade with more than one node. A single node is an explicitly supported + ! configuration (AA_Init requires only NumBlNds >= 1, and the standalone AeroAcoustics driver uses exactly + ! one node): that one element spans the whole blade, so there is no "last element" for BldPrcnt to be smaller + ! than. LastElemPct is 100.0 for a single node, so the warning fired on every single-node + ! run with BldPrcnt < 100 and told the user to refine an AeroDyn blade mesh that has nothing to refine. + LastElemPct = 100.0_ReKi ! a single element spans the entire blade IF (p%NumBlNds > 1) THEN - LastElemPct = 100.0 * (p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1)) / p%BlSpn(p%NumBlNds,1) - ELSE - LastElemPct = 100.0 ! Single node means element spans entire blade - ENDIF + + ! bjj: guard the divide. This would be rare, and an indication that the calling program has an error. + IF (p%BlSpn(p%NumBlNds,1) > 0.0_ReKi) THEN + LastElemPct = 100.0 * (p%BlSpn(p%NumBlNds,1) - p%BlSpn(p%NumBlNds-1,1)) / p%BlSpn(p%NumBlNds,1) + ELSE + CALL SetErrStat(ErrID_Fatal, 'The outermost blade node is at zero span, which is invalid.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + ENDIF - IF (InputFileData%AA_Bl_Prcntge .lt. LastElemPct) THEN - CALL SetErrStat(ErrID_Warn, 'BldPrcnt is smaller than the last blade element size. '// & - 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & - 'Either increase BldPrcnt in your aeroacoustic input file '// & - 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) - if(Failed()) return - endif + if (InputFileData%AA_Bl_Prcntge < LastElemPct) then + CALL SetErrStat(ErrID_Warn, 'BldPrcnt is smaller than the last blade element size. '// & + 'OpenFAST will move on assuming the last blade element, which is the minimum blade span used for noise calculations. '// & + 'Either increase BldPrcnt in your aeroacoustic input file '// & + 'or refine the aerodynamic mesh in your blade AeroDyn input file.', ErrStat2, ErrMsg2, RoutineName ) + if(Failed()) return + endif + + ENDIF p%startnode = min(p%NumBlNds, 2) BladeSpanUsedForNoise = p%BlSpn(p%NumBlNds,1)*(1.0 - InputFileData%AA_Bl_Prcntge/100.0) @@ -301,11 +315,12 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) DO I = 1,p%numBlades DO J = p%startnode,p%NumBlNds ! starts loop from startnode. IF (J < 2) THEN + ! bjj: J == 1 is reachable only when NumBlNds == 1 (otherwise startnode >= 2). p%BlElemSpn(J,I) = p%BlSpn(J,I) !assume this is the innermost node ELSEIF (J .EQ. p%NumBlNds) THEN p%BlElemSpn(J,I) = p%BlSpn(J,I)-p%BlSpn(J-1,I) ELSE - p%BlElemSpn(J,I) = (p%BlSpn(J,I)-p%BlSpn(J-1,I))/2 + (p%BlSpn(J+1,I)-p%BlSpn(J,I))/2 ! this is the average element size around this node, equivalent to (p%BlSpn(J+1,I) - p%BlSpn(J-1,I))/2 + p%BlElemSpn(J,I) = (p%BlSpn(J+1,I) - p%BlSpn(J-1,I)) / 2 ! this is the average element size around this node ENDIF end do end do @@ -368,7 +383,12 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) DO k=1,size(AFInfo) ! for each airfoil interpolation - ! find index where LE is found + ! NOTE: AFInfo(k)%X_Coord(1)/Y_Coord(1) is the aerodynamic center, NOT a point on the airfoil surface. + ! The airfoil surface coordinates are therefore indices 2:NumCoords, and the surface is traversed + ! from the trailing edge, around the leading edge, and back to the trailing edge. + + ! find index where LE is found (i.e., the first index where the x-coordinate starts increasing again) + iLE = -1 ! initialize to an invalid index so we can tell whether the search below succeeded DO i=3,size(AFInfo(k)%X_Coord) IF (AFInfo(k)%X_Coord(i) - AFInfo(k)%X_Coord(i-1) > 0.) THEN iLE = i @@ -376,7 +396,21 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ENDIF ENDDO - ! From LE toward TE + ! bjj: guard against the search failing; without this, iLE (and the indices derived from it) would be + ! used uninitialized, which produces garbage airfoil thicknesses (and thus garbage inflow noise) + ! that changes from run to run. + IF ( iLE < 3 .OR. iLE >= size(AFInfo(k)%X_Coord) ) THEN + call SetErrStat( ErrID_Fatal, 'Unable to locate the leading edge in the coordinates of airfoil '// & + trim(num2lstr(k))//'. The Full Guidati inflow-noise model (InflowMod=2) requires airfoil coordinates '// & + 'that start at the trailing edge, wrap around the leading edge, and return to the trailing edge.', & + ErrStat, ErrMsg, RoutineName ) + return + ENDIF + + ! From LE toward TE: find the points closest to 1% and 10% chord. + ! Seed the search with the LE point itself so that i1_1/i10_1 are always valid indices. + i1_1 = iLE + i10_1 = iLE dist1 = ABS( AFInfo(k)%X_Coord(iLE) - 0.01) dist10 = ABS( AFInfo(k)%X_Coord(iLE) - 0.10) DO i=iLE+1,size(AFInfo(k)%X_Coord) @@ -390,10 +424,15 @@ subroutine SetParameters( InitInp, InputFileData, p, AFInfo, ErrStat, ErrMsg ) ENDIF ENDDO - ! From TE to LE - dist1 = 0.99 - dist10 = 0.90 - DO i=1,iLE-1 + ! From TE to LE: find the points closest to 1% and 10% chord. + ! Start at index 2 (index 1 is the aerodynamic center, not a surface point) and seed the search with + ! that point so that i1_2/i10_2 are always valid indices. (The previous code seeded the search with + ! the arbitrary thresholds 0.99/0.90 and could leave these indices undefined.) + i1_2 = 2 + i10_2 = 2 + dist1 = ABS( AFInfo(k)%X_Coord(2) - 0.01) + dist10 = ABS( AFInfo(k)%X_Coord(2) - 0.10) + DO i=3,iLE-1 IF (ABS(AFInfo(k)%X_Coord(i) - 0.01) < dist1) THEN i1_2 = i dist1 = ABS(AFInfo(k)%X_Coord(i) - 0.01) @@ -490,8 +529,8 @@ subroutine Init_y(y, m, p, errStat, errMsg) end if call AllocAry(y%WriteOutput , p%numOutsAll(1), 'y%WriteOutput' , errStat2 , errMsg2); if(Failed()) return - call AllocAry(y%WriteOutputSep , p%numOutsAll(3), 'y%WriteOutputSep' , errStat2 , errMsg2); if(Failed()) return call AllocAry(y%WriteOutputForPE , p%numOutsAll(2), 'y%WriteOutputForPE' , errStat2 , errMsg2); if(Failed()) return + call AllocAry(y%WriteOutputSep , p%numOutsAll(3), 'y%WriteOutputSep' , errStat2 , errMsg2); if(Failed()) return call AllocAry(y%WriteOutputNodes , p%numOutsAll(4), 'y%WriteOutputSepFreq' , errStat2 , errMsg2); if(Failed()) return y%WriteOutput = 0.0_reki @@ -542,17 +581,35 @@ subroutine Init_MiscVars(m, p, errStat, errMsg) call AllocAry(m%PtotalFreq , size(p%FreqList) , p%NrObsLoc , 'm%PtotalFreq' , errStat2 , errMsg2); if(Failed()) return call AllocAry(m%OASPL , p%NrObsLoc , p%NumBlNds , p%NumBlades , 'm%OASPL' , errStat2 , errMsg2); if(Failed()) return + ! bjj: initialize every allocated array here. Most of these are overwritten on each call to AA_CalcOutput before + ! they are read, but m%ChordAngleLE/m%SpanAngleLE and the SPL* arrays were previously left with whatever happened + ! to be in the freshly allocated memory, which is a latent source of run-to-run differences. m%ChordAngleTE = 0.0_ReKi m%SpanAngleTE = 0.0_ReKi + m%ChordAngleLE = 0.0_ReKi + m%SpanAngleLE = 0.0_ReKi m%rTEtoObserve = 0.0_ReKi m%rLEtoObserve = 0.0_ReKi + m%SPLLBL = 0.0_ReKi + m%SPLP = 0.0_ReKi + m%SPLS = 0.0_ReKi + m%SPLALPH = 0.0_ReKi + m%SPLBLUNT = 0.0_ReKi + m%SPLTIP = 0.0_ReKi + m%SPLTI = 0.0_ReKi m%SPLTIGui = 0.0_ReKi + m%CfVar = 0.0_ReKi m%d99Var = 0.0_ReKi m%dstarVar = 0.0_ReKi m%EdgeVelVar = 0.0_ReKi m%LE_Location = 0.0_ReKi + + m%DirectiviOutput = 0.0_ReKi + m%SumSpecNoiseSep = 0.0_ReKi + m%PtotalFreq = 0.0_ReKi + m%OASPL = 0.0_ReKi contains logical function Failed() call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -644,6 +701,14 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) OtherState%allregcounter(k_minus1,rco_minus1) = OtherState%allregcounter(k_minus1,rco_minus1) + 1 ! increase the sample amount in that specific bin + ! bjj: This counter only ever increases, so on a long enough simulation it would overflow IntKi and go + ! negative, which would make the circular-buffer index below zero or negative (an out-of-bounds access). + ! Once we are past the first full window, subtracting exactly Num_total_sampleTI keeps the counter + ! bounded without changing either the test below or the mod() phase of the buffer index. + if ( OtherState%allregcounter(k_minus1,rco_minus1) > 2*p%Num_total_sampleTI ) then + OtherState%allregcounter(k_minus1,rco_minus1) = OtherState%allregcounter(k_minus1,rco_minus1) - p%Num_total_sampleTI + end if + InflowNorm = TwoNorm( u%Inflow(:,j,i) ) !note: p%Num_total_sampleTI = size(xd%RegVxStor,1) ! with storage region dependent moving average and TI @@ -651,7 +716,10 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) xd%RegVxStor(OtherState%allregcounter(k_minus1,rco_minus1),k_minus1,rco_minus1) = InflowNorm xd%TIVx(j,i) = 0 ELSE - xd%RegVxStor( mod( OtherState%allregcounter(k_minus1,rco_minus1), p%Num_total_sampleTI )+1, k_minus1, rco_minus1)=InflowNorm + ! bjj: use mod(counter-1, N)+1 so that the circular buffer advances in lock step with the counter. + ! The previous expression, mod(counter, N)+1, skipped slot 1 on the first pass, so the very first + ! sample stayed in the window for two full windows instead of one. + xd%RegVxStor( mod( OtherState%allregcounter(k_minus1,rco_minus1)-1, p%Num_total_sampleTI )+1, k_minus1, rco_minus1)=InflowNorm meanInflow = SUM( xd%RegVxStor(:,k_minus1,rco_minus1) ) /p%Num_total_sampleTI if ( EqualRealNos(meanInflow,0.0_ReKi)) then @@ -669,7 +737,12 @@ subroutine AA_UpdateStates( t, n, m, u, p, xd, OtherState, errStat, errMsg ) do j=1,p%NumBlNds ! We scale the incident turbulence intensity by the ratio of average to incident wind speed ! The scaled TI is used by the Amiet model - xd%TIVx(j,i)=p%TI * p%avgV/u%Vrel(J,I) + ! bjj: Vrel can be zero (e.g., before the rotor starts spinning, or at an inboard node), which would + ! make TIVx infinite. InflowNoise() takes 10*log10(TIVx**2), so an Inf here shows up as a huge spike in + ! the inflow-turbulence noise. Apply the same lower bound on relative velocity that AA_CalcOutput uses. + InflowNorm = u%Vrel(J,I) + IF (abs(InflowNorm) < AA_u_min) InflowNorm = SIGN(AA_u_min, InflowNorm) + xd%TIVx(j,i)=p%TI * p%avgV/InflowNorm enddo enddo endif @@ -779,6 +852,26 @@ REAL(ReKi) FUNCTION Log10AA(X) RESULT(F) F = LOG10( MAX(AA_EPSILON, X) ) END FUNCTION Log10AA +!---------------------------------------------------------------------------------------------------------------------------------- +REAL(ReKi) FUNCTION Power10AA(X) RESULT(F) + REAL(ReKi),INTENT(IN) :: X + + !note that this "SAVE" is a deviation from standard OpenFAST code. Because it is used only to throw a warning if we ever + ! trigger numerical overflow, I'm not concerned about it getting reset per instance. If we don't like this behavior, + ! we could store this FirstWarning as a misc var and thread it into every call, or just remove the warning alltogether. + LOGICAL, SAVE :: FirstWarning = .TRUE. + + if (FirstWarning) then + if (X > AA_max_exp) then + call WrScr(' WARNING: noise predictions are extremely large. Use caution in interpreting results.') + FirstWarning = .FALSE. + end if + end if + + + F = 10.0_ReKi ** min(X, AA_max_exp) + +END FUNCTION Power10AA !----------------------------------------------------------------------------------------------------------------------------------! SUBROUTINE Calc_LE_Location_Array(p,m,u) TYPE(AA_ParameterType), intent(in ) :: p !< Parameters @@ -926,11 +1019,7 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,errStat,errMsg) DO I = 1,p%numBlades DO J = p%startnode,p%NumBlNds ! starts loop from startnode. !------------------------------!!------------------------------!!------------------------------!!------------------------------! - - Unoise = u%Vrel(J,I) - IF (abs(Unoise) < AA_u_min) then - Unoise = SIGN(AA_u_min, Unoise) - ENDIF + Unoise = max(abs(u%Vrel(J,I)), AA_u_min) AlphaNoise= u%AoANoise(J,I) call MPi2Pi(AlphaNoise) ! make sure this is in an appropriate range [-pi,pi] @@ -991,7 +1080,7 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,errStat,errMsg) !--------Tip Noise--------------------------------------------------------------! - IF ( (p%ITIP == ITIP_ON) .AND. (J .EQ. p%NumBlNds) ) THEN ! calculate m%SPLTIP(1:nFreq) + IF ( (p%ITIP == ITIP_ON) .AND. (J == p%NumBlNds) ) THEN ! calculate m%SPLTIP(1:nFreq) CALL TIPNOIS(AlphaNoise_Deg,p%ALpRAT,p%BlChord(J,I),UNoise,m%ChordAngleTE(K,J,I),m%SpanAngleTE(K,J,I), & m%rTEtoObserve(K,J,I), p, m%SPLTIP) @@ -1091,7 +1180,7 @@ subroutine TotalContributionFromType(SPL,Ptotal,NoiseMech) do III=1,size(p%FreqList) ! Loops through each 1/3rd octave center frequency - Pt = 10.0_ReKi**(SPL(III)/10.0_ReKi) ! SPL to P Conversion for III Frequency + Pt = Power10AA(SPL(III)/10.0_ReKi) ! SPL to P Conversion for III Frequency P_SumAllFreq = P_SumAllFreq + Pt ! Sum for Running Total m%SumSpecNoiseSep(NoiseMech,III,K) = m%SumSpecNoiseSep(NoiseMech,III,K) + Pt ! Running sum of observer and frequency dependent sound pressure @@ -1137,8 +1226,8 @@ SUBROUTINE LBLVS(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,SPLLAM, integer(intKi) :: I ! I A generic index for DO loops. !compute reynolds number and mach number - M = U / p%SpdSound ! MACH NUMBER - RC = U * C/p%KinVisc ! REYNOLDS NUMBER BASED ON CHORD + M = abs(U / p%SpdSound) ! MACH NUMBER + RC = abs(U * C/p%KinVisc) ! REYNOLDS NUMBER BASED ON CHORD ! compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN @@ -1165,13 +1254,13 @@ SUBROUTINE LBLVS(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,SPLLAM, else ST1PRIM = .28 end if - STPKPRM = 10.**(-.04*ALPSTAR) * ST1PRIM ! Eq 56 from BPM Airfoil Self-noise and Prediction paper + STPKPRM = Power10AA(-.04*ALPSTAR) * ST1PRIM ! Eq 56 from BPM Airfoil Self-noise and Prediction paper ! compute reference reynolds number ! Eq 59 from BPM Airfoil Self-noise and Prediction paper IF (ALPSTAR .LE. 3.0) then - RC0=10.**(.215*ALPSTAR+4.978) + RC0=Power10AA(.215*ALPSTAR+4.978) else - RC0=10.**(.120*ALPSTAR+5.263) + RC0=Power10AA(.120*ALPSTAR+5.263) end if ! compute peak scaled spectrum level @@ -1288,8 +1377,8 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa LOGICAL :: SWITCH !!LOGICAL FOR COMPUTATION OF ANGLE OF ATTACK CONTRIBUTION ! Compute reynolds number and mach number - M = U / p%SpdSound - RC = U * C/p%KinVisc + M = abs(U / p%SpdSound) + RC = abs(U * C/p%KinVisc) ! Compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN @@ -1320,7 +1409,7 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa IF (ALPSTAR .LE. 1.333) then ST2 = ST1 elseif (ALPSTAR .LE. StallVal) then - ST2 = ST1*10.**(.0054*(ALPSTAR-1.333)**2) + ST2 = ST1*Power10AA(.0054*(ALPSTAR-1.333)**2) else ST2 = 4.72 * ST1 end if @@ -1431,9 +1520,9 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,d99Var2,dstarVar1,dstarVar2,StallVa IF (SPLS(I) .LT. -100.) SPLS(I) = -100. ! Similar to Eq 29 of BPM Airfoil Self-noise and Prediction paper IF (SPLALPH(I) .LT. -100.) SPLALPH(I) = -100. ! Eq 30 of BPM Airfoil Self-noise and Prediction paper recommends SPLALPH = 10log(stuff) + A' + K2, where A' is calculated same as A but with x3 Rc - !P1 = 10.**(SPLP(I) / 10.) ! SPL_Pressure - !P2 = 10.**(SPLS(I) / 10.) ! SPL_Suction - !P4 = 10.**(SPLALPH(I) / 10.) ! SPL_AoA + !P1 = Power10AA(SPLP(I) / 10.) ! SPL_Pressure + !P2 = Power10AA(SPLS(I) / 10.) ! SPL_Suction + !P4 = Power10AA(SPLALPH(I) / 10.) ! SPL_AoA !SPLTBL(I) = 10. * LOG10AA(P1 + P2 + P4) ! Eq 24 from BPM Airfoil Self-noise and Prediction paper @@ -1474,7 +1563,7 @@ SUBROUTINE TIPNOIS(ALPHTIP,ALPRAT2,C,U ,THETA,PHI, R,p,SPLTIP) ENDIF !! used to be ALPTIPP = ALPHTIP * ALPRAT2 ALPTIPP = ABS(ALPHTIP) * ALPRAT2 - M = U / p%SpdSound ! MACH NUMBER + M = abs(U / p%SpdSound) ! MACH NUMBER ! Compute directivity function DBARH = DIRECTH_TE(M,THETA,PHI) IF (p%ROUND) THEN @@ -1490,11 +1579,12 @@ SUBROUTINE TIPNOIS(ALPHTIP,ALPRAT2,C,U ,THETA,PHI, R,p,SPLTIP) UM = MM * p%SpdSound ! Eq 65 from BPM Airfoil Self-noise and Prediction paper TERM = M*M*MM**3*L**2*DBARH/R**2 ! TERM = M^2 * M_max^5 *l^2 *D / r^2 according to Semi-Empirical Aeroacoustic Noise Prediction Code for Wind Turbines paper ! Term is correct according to Eq 61 from BPM Airfoil self-noise and Prediction paper - IF (TERM .NE. 0.0) THEN - SCALE = 10.*LOG10(TERM) - ELSE - SCALE = 0.0 - ENDIF + ! bjj: Use Log10AA here instead of a "TERM /= 0" test around LOG10: + ! (1) TERM is negative when U (and therefore M and MM) is negative, and LOG10 of a negative number returns NaN. + ! (2) The old ELSE branch set SCALE = 0, which left SPLTIP near 126 dB when the directivity DBARH was zero -- + ! i.e., it reported a large tip noise level precisely where there should be none. Every other noise mechanism + ! in this module treats DBARH <= 0 as "no contribution", and Log10AA now gives that behavior here too. + SCALE = 10.*Log10AA(TERM) DO I=1,size(p%FreqList) STPP = p%FreqList(I) * L / UM ! Eq 62 from BPM Airfoil Self-noise and Prediction paper SPLTIP(I) = 126.-30.5*(LOG10AA(STPP)+.3)**2 + SCALE ! Eq 61 from BPM Airfoil Self-noise and Prediction paper @@ -1536,7 +1626,7 @@ SUBROUTINE InflowNoise(AlphaNoise,Chord,U,THETA,PHI,d,RObs,TINoise,p,SPLti) INTEGER(intKi) :: I !I A generic index for DO loops. !!!--- NAF NOISE IDENTICAL - Mach = U/p%SpdSound + Mach = abs(U/p%SpdSound) ! This part is recently added for height and surface roughness dependent estimation of turbulence intensity and turbulence scales !%Lturb=300*(Z/300)^(0.46+0.074*log(p%z0_aa)); !% Gives larger length scale @@ -1780,8 +1870,8 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, real(ReKi) :: LogVal ! temp variable to help us not take log10(0) --- ! Reynolds number and mach number - M = U / p%SpdSound - RC = U * C/p%KinVisc + M = abs(U / p%SpdSound) + RC = abs(U * C/p%KinVisc) ! Compute boundary layer thicknesses IF (p%X_BLMethod .eq. X_BLMethod_Tables) THEN DELTAP = d99Var2 @@ -1793,6 +1883,15 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, ! Compute average displacement thickness DSTRAVG = (DSTRS + DSTRP) / 2. + + ! bjj: When the boundary layer properties are read from tables (X_BLMethod_Tables), DSTRS and DSTRP can be zero + ! (e.g., for table entries that were not converged). Without this check, HDSTAR becomes Inf (or NaN if H is also + ! zero), which then propagates through G5COMP into SPLBLUNT and on into the summed noise levels. + IF (DSTRAVG <= 0. .OR. H <= 0.) THEN + SPLBLUNT = -100. ! same floor used by the other noise mechanisms; effectively no contribution + RETURN + ENDIF + HDSTAR = H / DSTRAVG DSTARH = 1. /HDSTAR ! Compute directivity function @@ -1831,11 +1930,11 @@ SUBROUTINE BLUNT(ALPSTAR,C,U ,THETA,PHI,L,R,H,PSI,p,d99Var2,dstarVar1,dstarVar2, G5(I) = G50 + .0714 * PSI * (G514-G50) ! interpolate G5 from G50 and G514 IF (G5(I) .GT. 0.) G5(I) = 0. - G5Sum = 10**(G5(I)/10)+G5Sum ! to be subtracted + G5Sum = Power10AA(G5(I)/10)+G5Sum ! to be subtracted if ( G5Sum .ne. 0) then - LogVal = MAX(AA_EPSILON,1/G5Sum) + LogVal = MAX(AA_EPSILON,1/G5Sum) else - LogVal = 1 + LogVal = 1.0_ReKi ! note that this is reached only when the accumulated bluntness contribution is already ~ 10^-69 end if SPLBLUNT(I) = G4 + G5(I) + SCALE - 10*log10(LogVal) ! equation mentioned there is plus but it is stated subtract, thus ''- 10*log10(1/G5Sum)'' end do @@ -1878,13 +1977,16 @@ REAL(ReKi) FUNCTION G5COMP(HDSTAR,ETA) result(G5) ETA0 = -SQRT((M*M*MU**4)/(6.25+M*M*MU*MU)) ! eq 80 from BPM Airfoil Self-noise and Prediction paper + ! bjj: the arguments of the SQRT calls below are guarded with MAX(0,...): they are all analytically non-negative + ! over the range of each branch, but the third branch in particular has a radicand that is exactly zero at the + ! branch boundary (SQRT(1.5625/1194.99) = 0.03615995), so round-off can make it slightly negative and return NaN. IF (ETA .LE. ETA0) then - K = 2.5*SQRT(1.-(ETA0/MU)**2)-2.5-M*ETA0 ! eq 81 from BPM Airfoil Self-noise and Prediction paper + K = 2.5*SQRT(MAX(0.0_ReKi, 1.-(ETA0/MU)**2))-2.5-M*ETA0 ! eq 81 from BPM Airfoil Self-noise and Prediction paper G5 = M * ETA + K ! begin eq 76 from BPM Airfoil Self-noise and Prediction paper elseif (ETA .LE. 0.) then - G5 = 2.5*SQRT(1.-(ETA/MU)**2)-2.5 + G5 = 2.5*SQRT(MAX(0.0_ReKi, 1.-(ETA/MU)**2))-2.5 elseif (ETA .LE. 0.03615995) then - G5 = SQRT(1.5625-1194.99*ETA**2)-1.25 + G5 = SQRT(MAX(0.0_ReKi, 1.5625-1194.99*ETA**2))-1.25 else G5 = -155.543 * ETA + 4.375 end if @@ -2006,12 +2108,12 @@ SUBROUTINE THICK(C,RC,ALPSTAR,p,DELTAP,DSTRS,DSTRP,StallVal) LogRC = LOG10AA( RC ) ! Boundary layer thickness - DELTA0 = 10.**(1.6569-0.9045*LogRC+0.0596*LogRC**2)*C ! (untripped) Eq. (5) of [1] - IF (p%ITRIP /= ITRIP_None) DELTA0 = 10.**(1.892 -0.9045*LogRC+0.0596*LogRC**2)*C ! (heavily tripped) Eq. (2) of [1] + DELTA0 = Power10AA(1.6569-0.9045*LogRC+0.0596*LogRC**2)*C ! (untripped) Eq. (5) of [1] + IF (p%ITRIP /= ITRIP_None) DELTA0 = Power10AA(1.892 -0.9045*LogRC+0.0596*LogRC**2)*C ! (heavily tripped) Eq. (2) of [1] IF (p%ITRIP .EQ. ITRIP_Light) DELTA0=.6*DELTA0 ! Pressure side boundary layer thickness, Eq (8) of [1] - DELTAP = 10.**(-.04175*ALPSTAR+.00106*ALPSTAR**2)*DELTA0 + DELTAP = Power10AA(-.04175*ALPSTAR+.00106*ALPSTAR**2)*DELTA0 ! Compute zero angle of attack displacement thickness IF (p%ITRIP /= ITRIP_None) THEN @@ -2019,37 +2121,37 @@ SUBROUTINE THICK(C,RC,ALPSTAR,p,DELTAP,DSTRS,DSTRP,StallVal) IF (RC .LE. .3E+06) THEN DSTR0 = .0601 * RC **(-.114)*C ELSE - DSTR0=10.**(3.411-1.5397*LogRC+.1059*LogRC**2)*C + DSTR0=Power10AA(3.411-1.5397*LogRC+.1059*LogRC**2)*C END IF ! Lightly tripped IF (p%ITRIP .EQ. ITRIP_Light) DSTR0 = DSTR0 * .6 ELSE ! Untripped, Eq. (6) of [1] - DSTR0=10.**(3.0187-1.5397*LogRC+.1059*LogRC**2)*C + DSTR0=Power10AA(3.0187-1.5397*LogRC+.1059*LogRC**2)*C ENDIF ! Pressure side displacement thickness, Eq. (9) of [1] - DSTRP = 10.**(-.0432*ALPSTAR+.00113*ALPSTAR**2)*DSTR0 + DSTRP = Power10AA(-.0432*ALPSTAR+.00113*ALPSTAR**2)*DSTR0 ! IF (p%ITRIP .EQ. 3) DSTRP = DSTRP * 1.48 ! commented since itrip is never 3 check if meant 2.(EB_DTU) ! Suction side displacement thickness IF (p%ITRIP .EQ. ITRIP_Heavy) THEN ! Heavily tripped, Eq. (12) of [1] IF (ALPSTAR .LE. 5.) THEN - DSTRS=10.**(.0679*ALPSTAR)*DSTR0 + DSTRS=Power10AA(.0679*ALPSTAR)*DSTR0 ELSEIF (ALPSTAR .LE. StallVal) THEN - DSTRS = 0.381 * 10.**(.1516*ALPSTAR)*DSTR0 + DSTRS = 0.381 * Power10AA(.1516*ALPSTAR)*DSTR0 ELSE - DSTRS = 14.296 * 10.**(.0258*ALPSTAR)*DSTR0 + DSTRS = 14.296 * Power10AA(.0258*ALPSTAR)*DSTR0 ENDIF ELSE ! Untripped or lightly tripped, Eq. (15) of [1] IF (ALPSTAR .LE. 7.5) THEN - DSTRS =10.**(.0679*ALPSTAR)*DSTR0 + DSTRS =Power10AA(.0679*ALPSTAR)*DSTR0 ELSEIF(ALPSTAR .LE. StallVal) THEN - DSTRS = .0162*10.**(.3066*ALPSTAR)*DSTR0 + DSTRS = .0162*Power10AA(.3066*ALPSTAR)*DSTR0 ELSE - DSTRS = 52.42*10.**(.0258*ALPSTAR)*DSTR0 + DSTRS = 52.42*Power10AA(.0258*ALPSTAR)*DSTR0 ENDIF ENDIF @@ -2179,14 +2281,16 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) n_freq = size(p%FreqList) freq = p%FreqList - SPLS = 0.0_ReKi ! initialize in case Cfall(1) <= 0 - SPLP = 0.0_ReKi ! initialize in case Cfall(2) <= 0 + ! bjj: -100 dB (the floor used elsewhere in this module) means "no contribution". The previous value of 0 dB + ! actually adds unit mean-square pressure at every frequency when Cfall <= 0. + SPLS = -100.0_ReKi ! initialize in case Cfall(1) <= 0 + SPLP = -100.0_ReKi ! initialize in case Cfall(2) <= 0 ! Body of TNO band_ratio = 2.**(1./3.) ! Mach number - Mach = U / p%SpdSound + Mach = abs(U / p%SpdSound) ! Directivity function DBARH = DIRECTH_TE(REAL(Mach,ReKi),THETA,PHI) @@ -2194,18 +2298,25 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) do i_omega = 1,n_freq omega = TwoPi*p%FreqList(i_omega) !integration limits + ! bjj: use ABS(Mach) so that the upper limit stays above the lower limit. U (and therefore Mach) carries the + ! sign of Vrel, and a negative upper limit would reverse the integration interval and return a bogus spectrum. int_limits(1) = 0.0e0 int_limits(2) = 10*omega/(Mach*p%SpdSound) ! Convert to third octave band_width = 2. * omega * (sqrt(band_ratio)-1./sqrt(band_ratio)) + ! bjj: The log10 calls below were unguarded. "answer" is the result of a fixed-order Gauss-Kronrod quadrature + ! of an integrand that underflows to zero at high frequency and can also return a small negative value through + ! round-off cancellation. log10(0) is -Inf and log10(negative) is NaN, and the "SPL < -100" floor below does not + ! catch NaN (any comparison with NaN is false), so a NaN would propagate silently into the summed noise levels. + ! Log10AA() floors the argument at AA_EPSILON, which maps both cases to "no contribution". IF (Cfall(1) .GT. 0.) THEN answer = SPL_integrate(omega=omega,limits=int_limits,ISSUCTION=.true., & Mach=Mach,SpdSound=p%SpdSound,AirDens=p%AirDens,KinVisc=p%KinVisc, & Cfall=Cfall,d99all=d99all,EdgeVelAll=EdgeVelAll) Spectrum = D/(4.*pi*R**2)*answer - SPL_suction = 10.*log10(Spectrum*DBARH/2.e-5/2.e-5) - SPLS(i_omega) = SPL_suction + 10.*log10(band_width) + SPL_suction = 10.*Log10AA(Spectrum*DBARH/2.e-5/2.e-5) + SPLS(i_omega) = SPL_suction + 10.*Log10AA(band_width) ENDIF IF (Cfall(2) .GT. 0.) THEN @@ -2213,17 +2324,17 @@ SUBROUTINE TBLTE_TNO(U,THETA,PHI,D,R,Cfall,d99all,EdgeVelAll,p,SPLP,SPLS) Mach=Mach,SpdSound=p%SpdSound,AirDens=p%AirDens,KinVisc=p%KinVisc, & Cfall=Cfall,d99all=d99all,EdgeVelAll=EdgeVelAll) Spectrum = D/(4.*pi*R**2)*answer - SPL_press = 10.*log10(Spectrum*DBARH/2.e-5/2.e-5) - SPLP(i_omega) = SPL_press + 10.*log10(band_width) + SPL_press = 10.*Log10AA(Spectrum*DBARH/2.e-5/2.e-5) + SPLP(i_omega) = SPL_press + 10.*Log10AA(band_width) ENDIF ! Sum the noise sources SPLALPH is BPM value IF (SPLP(i_omega) .LT. -100.) SPLP(i_omega) = -100. IF (SPLS(i_omega) .LT. -100.) SPLS(i_omega) = -100. - !P1 = 10.**(SPLP(i_omega) / 10.) - !P2 = 10.**(SPLS(i_omega) / 10.) - !P4 = 10.**(SPLALPH(i_omega) / 10.) + !P1 = Power10AA(SPLP(i_omega) / 10.) + !P2 = Power10AA(SPLS(i_omega) / 10.) + !P4 = Power10AA(SPLALPH(i_omega) / 10.) ! !SPLTBL(i_omega) = 10. * LOG10(P1 + P2 + P4) enddo @@ -2256,7 +2367,7 @@ SUBROUTINE BL_Param_Interp(p,m,U,AlphaNoise_Deg,C,whichAirfoil) !!!! this if is not used but if necessary two sets of tables can be populated for tripped and untripped cases - RC = U * C/p%KinVisc ! REYNOLDS NUMBER BASED ON CHORD + RC = abs(U * C/p%KinVisc) ! REYNOLDS NUMBER BASED ON CHORD ! find the indices into the arrays representing coordinates of each dimension: diff --git a/modules/aerodyn/src/AeroAcoustics_Driver.f90 b/modules/aerodyn/src/AeroAcoustics_Driver.f90 index ea5bdbf6bd..81040948de 100644 --- a/modules/aerodyn/src/AeroAcoustics_Driver.f90 +++ b/modules/aerodyn/src/AeroAcoustics_Driver.f90 @@ -16,17 +16,6 @@ ! See the License for the specific language governing permissions and ! limitations under the License. ! - -!there will also be various control flags... this may be updated as needed: -!TBLflag = {'BPM','TNO'} -!bluntnessFlag = {'DTU','BPM'} -!BPMBLflag = {'true','false'} -!useOrigModelAtSepOnset = {'true','false'} - - - - - !********************************************************************************************************************************** program AeroAcoustics_Driver use AeroAcoustics_Driver_Subs diff --git a/modules/aerodyn/src/AeroAcoustics_IO.f90 b/modules/aerodyn/src/AeroAcoustics_IO.f90 index 24049bb023..3fa7110495 100644 --- a/modules/aerodyn/src/AeroAcoustics_IO.f90 +++ b/modules/aerodyn/src/AeroAcoustics_IO.f90 @@ -46,7 +46,7 @@ MODULE AeroAcoustics_IO integer(intKi), parameter :: ITURB_None = 0 ! TBLTE noise is not calculated integer(intKi), parameter :: ITURB_BPM = 1 ! TBLTE noise is calculated with BPM - integer(intKi), parameter :: ITURB_TNO = 2 ! TBLTE noise is calculated with TNO + integer(intKi), parameter :: ITURB_TNO = 2 ! TBLTE noise is calculated with TNO (and part of BPM) integer(intKi), parameter :: IInflow_None = 0 ! IInflow noise is not calculated integer(intKi), parameter :: IInflow_BPM = 1 ! IInflow noise is calculated with BPM @@ -452,7 +452,7 @@ SUBROUTINE ValidateInputData( InputFileData, NumBl, ErrStat, ErrMsg ) if (InputFileData%IBLUNT /= IBLUNT_None .and. InputFileData%IBLUNT /= IBLUNT_BPM) then call SetErrStat ( ErrID_Fatal, & - 'IBLUNT must '//trim(num2lstr(IBLUNT_None))//' (none) or '//trim(num2lstr(IBLUNT_BPM))//' (Bluntness noise calculated).', ErrStat, ErrMsg, RoutineName ) + 'IBLUNT must be '//trim(num2lstr(IBLUNT_None))//' (none) or '//trim(num2lstr(IBLUNT_BPM))//' (Bluntness noise calculated).', ErrStat, ErrMsg, RoutineName ) endif if (InputFileData%ILAM /= ILAM_None .and. InputFileData%ilam /= ILAM_BPM) then call SetErrStat ( ErrID_Fatal, 'ILAM must be '//trim(num2lstr(ILAM_None))//' No calculation '//& @@ -467,7 +467,7 @@ SUBROUTINE ValidateInputData( InputFileData, NumBl, ErrStat, ErrMsg ) ' (heavily tripped BL Calculation) or '//trim(num2lstr(ITRIP_Light))//' (lightly tripped BL)' ,ErrStat, ErrMsg, RoutineName ) end if if (InputFileData%ITURB /= ITURB_None .and. InputFileData%ITURB /= ITURB_BPM .and. InputFileData%ITURB /= ITURB_TNO) then - call SetErrStat ( ErrID_Fatal, 'ITURB must be 0 (off) or 1 (BPM) or 2 (TNO) .', ErrStat, ErrMsg, RoutineName ) + call SetErrStat ( ErrID_Fatal, 'ITURB must be 0 (off) or 1 (BPM) or 2 (TNO with BPM alpha) .', ErrStat, ErrMsg, RoutineName ) end if if (InputFileData%IInflow /= IInflow_None .and. InputFileData%IInflow /= IInflow_BPM & .and. InputFileData%IInflow /= IInflow_FullGuidati .and. InputFileData%IInflow /= IInflow_SimpleGuidati ) then diff --git a/modules/aerodyn/src/AeroAcoustics_TNO.f90 b/modules/aerodyn/src/AeroAcoustics_TNO.f90 index 00053096fb..c0d3c46731 100644 --- a/modules/aerodyn/src/AeroAcoustics_TNO.f90 +++ b/modules/aerodyn/src/AeroAcoustics_TNO.f90 @@ -2,7 +2,7 @@ MODULE TNO use NWTC_Library ! ReKi, DBKi, R8Ki - use NWTC_SLATEC ! slatec_qk61 -- which is all that is in that library right now. + !use NWTC_SLATEC ! slatec_qk61: bjj: this has been replaced with a local qk61_ctx that passes a context (TNO_ContextType) to the integrand. See comments below. implicit none PRIVATE @@ -15,27 +15,106 @@ MODULE TNO REAL (TNOKi), PARAMETER :: Cmu = 0.09 ! INTEGER(IntKi),PARAMETER :: limit = 5000 - !TNO variables - REAL (TNOKi) :: Omega_TNO ! NOTE: not a constant and used by function f_int1 and f_int2 + !> Everything the TNO integrands need in order to be evaluated. + !! + !! This type replaces what used to be a set of module-level (and therefore implicitly SAVEd) variables. The reason + !! those existed is that the SLATEC integrator qk61 declares its integrand as "external f" and calls it as f(x), so + !! there was nowhere to thread the flow conditions through. Passing this context as an explicit argument instead + !! makes the module reentrant: every activation of SPL_integrate owns its own context, so concurrent calls (e.g. + !! from an OpenMP region, or two rotors evaluated in parallel) can no longer overwrite each other's state + !! mid-integration. + type :: TNO_ContextType + ! frequency + real(TNOKi) :: Omega = 0.0_TNOKi !< radian frequency + ! atmosphere + real(TNOKi) :: nu = 0.0_TNOKi !< kinematic viscosity + real(TNOKi) :: co = 0.0_TNOKi !< speed of sound + real(TNOKi) :: rho = 0.0_TNOKi !< air density + ! airfoil + real(TNOKi) :: Mach = 0.0_TNOKi !< Mach number of this blade node + logical :: IsSuction = .false. !< .true. for the suction side, .false. for the pressure side + ! blade-node boundary-layer properties; index 1 = suction side, index 2 = pressure side + real(TNOKi) :: d99(2) = 0.0_TNOKi !< boundary layer thickness + real(TNOKi) :: Cf(2) = 0.0_TNOKi !< skin friction coefficient + real(TNOKi) :: edgevel(2) = 0.0_TNOKi !< edge velocity ratio + ! Wavenumbers. Unlike the fields above (which are fixed for the whole integration), these are updated by + ! Pressure() at each point of the OUTER integration over k1 and then read by f_int1 during the INNER + ! integration over x2. That mid-integration coupling is why the context is passed as intent(inout). + real(TNOKi) :: k1 = 0.0_TNOKi + real(TNOKi) :: k3 = 0.0_TNOKi + real(TNOKi) :: k = 0.0_TNOKi + end type TNO_ContextType - !atmosphere variables - REAL (TNOKi) :: nu - REAL (TNOKi) :: co - REAL (TNOKi) :: rho + abstract interface + !> Integrand accepted by qk61_ctx: a function of the integration variable plus a context. + real(TNOKi) function TNO_Integrand(x, ctx) + import :: TNOKi, TNO_ContextType + real(TNOKi), intent(in ) :: x + type(TNO_ContextType), intent(inout) :: ctx + end function TNO_Integrand + end interface - ! Wavenumber variables - REAL (TNOKi) :: k - REAL (TNOKi) :: k1 - REAL (TNOKi) :: k3 + ! ------------------------------------------------------------------------------------------------------------------ + ! Data for the 61-point Gauss-Kronrod rule, transcribed from modules/nwtc-library/src/NetLib/slatec/dqk61.f + ! (SLATEC/QUADPACK; weights and abscissae evaluated with 80-decimal-digit arithmetic by L. W. Fullerton, Bell Labs, + ! Nov. 1981). The abscissae and weights are given for the interval (-1,1); by symmetry only the non-negative + ! abscissae and their corresponding weights are stored. + ! + ! The literals below carry dqk61's full precision and are converted to TNOKi, so a single-precision build rounds to + ! exactly the values hard-coded in qk61.f and a double-precision build matches dqk61.f. That is what lets this + ! routine reproduce the previous results bit-for-bit in both precisions. + ! + ! xgk - abscissae of the 61-point Kronrod rule. + ! xgk(2), xgk(4), ... are the abscissae of the 30-point Gauss rule. + ! xgk(1), xgk(3), ... are the optimally added abscissae. + ! wgk - weights of the 61-point Kronrod rule. + ! wg - weights of the 30-point Gauss rule. + ! ------------------------------------------------------------------------------------------------------------------ + real(TNOKi), parameter :: xgk(31) = (/ & + 0.999484410050490637571325895705811_TNOKi, 0.996893484074649540271630050918695_TNOKi, & + 0.991630996870404594858628366109486_TNOKi, 0.983668123279747209970032581605663_TNOKi, & + 0.973116322501126268374693868423707_TNOKi, 0.960021864968307512216871025581798_TNOKi, & + 0.944374444748559979415831324037439_TNOKi, 0.926200047429274325879324277080474_TNOKi, & + 0.905573307699907798546522558925958_TNOKi, 0.882560535792052681543116462530226_TNOKi, & + 0.857205233546061098958658510658944_TNOKi, 0.829565762382768397442898119732502_TNOKi, & + 0.799727835821839083013668942322683_TNOKi, 0.767777432104826194917977340974503_TNOKi, & + 0.733790062453226804726171131369528_TNOKi, 0.697850494793315796932292388026640_TNOKi, & + 0.660061064126626961370053668149271_TNOKi, 0.620526182989242861140477556431189_TNOKi, & + 0.579345235826361691756024932172540_TNOKi, 0.536624148142019899264169793311073_TNOKi, & + 0.492480467861778574993693061207709_TNOKi, 0.447033769538089176780609900322854_TNOKi, & + 0.400401254830394392535476211542661_TNOKi, 0.352704725530878113471037207089374_TNOKi, & + 0.304073202273625077372677107199257_TNOKi, 0.254636926167889846439805129817805_TNOKi, & + 0.204525116682309891438957671002025_TNOKi, 0.153869913608583546963794672743256_TNOKi, & + 0.102806937966737030147096751318001_TNOKi, 0.051471842555317695833025213166723_TNOKi, & + 0.000000000000000000000000000000000_TNOKi /) - ! Blade params - REAL (TNOKi) :: d99(2) - REAL (TNOKi) :: Cf(2) - REAL (TNOKi) :: edgevel(2) + real(TNOKi), parameter :: wgk(31) = (/ & + 0.001389013698677007624551591226760_TNOKi, 0.003890461127099884051267201844516_TNOKi, & + 0.006630703915931292173319826369750_TNOKi, 0.009273279659517763428441146892024_TNOKi, & + 0.011823015253496341742232898853251_TNOKi, 0.014369729507045804812451432443580_TNOKi, & + 0.016920889189053272627572289420322_TNOKi, 0.019414141193942381173408951050128_TNOKi, & + 0.021828035821609192297167485738339_TNOKi, 0.024191162078080601365686370725232_TNOKi, & + 0.026509954882333101610601709335075_TNOKi, 0.028754048765041292843978785354334_TNOKi, & + 0.030907257562387762472884252943092_TNOKi, 0.032981447057483726031814191016854_TNOKi, & + 0.034979338028060024137499670731468_TNOKi, 0.036882364651821229223911065617136_TNOKi, & + 0.038678945624727592950348651532281_TNOKi, 0.040374538951535959111995279752468_TNOKi, & + 0.041969810215164246147147541285970_TNOKi, 0.043452539701356069316831728117073_TNOKi, & + 0.044814800133162663192355551616723_TNOKi, 0.046059238271006988116271735559374_TNOKi, & + 0.047185546569299153945261478181099_TNOKi, 0.048185861757087129140779492298305_TNOKi, & + 0.049055434555029778887528165367238_TNOKi, 0.049795683427074206357811569379942_TNOKi, & + 0.050405921402782346840893085653585_TNOKi, 0.050881795898749606492297473049805_TNOKi, & + 0.051221547849258772170656282604944_TNOKi, 0.051426128537459025933862879215781_TNOKi, & + 0.051494729429451567558340433647099_TNOKi /) - ! Airfoil - REAL(TNOKi) :: Mach_TNO - LOGICAL :: ISSUCTION_TNO + real(TNOKi), parameter :: wg(15) = (/ & + 0.007968192496166605615465883474674_TNOKi, 0.018466468311090959142302131912047_TNOKi, & + 0.028784707883323369349719179611292_TNOKi, 0.038799192569627049596801936446348_TNOKi, & + 0.048402672830594052902938140422808_TNOKi, 0.057493156217619066481721689402056_TNOKi, & + 0.065974229882180495128128515115962_TNOKi, 0.073755974737705206268243850022191_TNOKi, & + 0.080755895229420215354694938460530_TNOKi, 0.086899787201082979802387530715126_TNOKi, & + 0.092122522237786128717632707087619_TNOKi, 0.096368737174644259639468626351810_TNOKi, & + 0.099593420586795267062780282103569_TNOKi, 0.101762389748405504596428952168554_TNOKi, & + 0.102852652893558840341285636705415_TNOKi /) contains @@ -56,36 +135,144 @@ function SPL_integrate(Omega,limits,ISSUCTION, & real(ReKi), intent(in ) :: EdgeVelAll(2) !< real(ReKi) :: integrand !< integrand result - real(TNOKi) :: answer !< value returned from qk61, NOTE the typing + real(TNOKi) :: answer !< value returned from qk61_ctx, NOTE the typing + + !> All state needed by the integrands. This is a local variable (not module data), which is what makes this + !! routine safe to call from more than one thread at a time. + type(TNO_ContextType) :: ctx ! local variables that are ignored real(TNOKi) :: abserr,resabs,resasc !< accuracy estimates and residuals. Currently ignored - ! Set module values from input - ISSUCTION_TNO = ISSUCTION - Omega_TNO = real(Omega,TNOKi) + ! Set context values from input + ctx%IsSuction = ISSUCTION + ctx%Omega = real(Omega,TNOKi) ! Mach number of segment - Mach_TNO = real(Mach,TNOKi) + ctx%Mach = real(Mach,TNOKi) ! Atmospheric values - co = real(SpdSound, TNOKi) - rho = real(AirDens, TNOKi) - nu = real(KinVisc, TNOKi) + ctx%co = real(SpdSound, TNOKi) + ctx%rho = real(AirDens, TNOKi) + ctx%nu = real(KinVisc, TNOKi) ! Blade node values - Cf = real(Cfall, TNOKi) - d99 = real(d99all, TNOKi) - edgevel = real(ABS(EdgeVelAll),TNOKi) + ctx%Cf = real(Cfall, TNOKi) + ctx%d99 = real(d99all, TNOKi) + ctx%edgevel = real(ABS(EdgeVelAll),TNOKi) - call slatec_qk61(f_int2,limits(1),limits(2),answer,abserr,resabs,resasc) + call qk61_ctx(f_int2,limits(1),limits(2),ctx,answer,abserr,resabs,resasc) integrand = real( answer, ReKi ) end function SPL_integrate +!================================================================================================================================== +!> Integrate FUNC over (a,b) using the 61-point Gauss-Kronrod rule, passing CTX through to the integrand. +!! +!! This is a direct port of the SLATEC routine qk61 (see modules/nwtc-library/src/NetLib/slatec/qk61.f), changed only +!! so that the integrand takes a context argument. The accumulation order is deliberately identical to the original so +!! that the computed RESULT is bit-for-bit the same as before this port. +!! +!! Two things the original could not offer: +!! 1. an explicit interface for the integrand (qk61 used "external f", an implicit interface), and +!! 2. a way to give the integrand its parameters other than global data. +!! +!! Declared RECURSIVE because the TNO model nests two integrations: the outer integral over k1 calls f_int2, which +!! calls Pressure, which integrates over x2. The original relied on per-file compiler flags (-frecursive / +!! -assume recursion, see modules/nwtc-library/CMakeLists.txt) to make qk61's locals automatic; stating RECURSIVE in +!! the source expresses that requirement in the code instead of in the build system. +recursive subroutine qk61_ctx(func, a, b, ctx, result, abserr, resabs, resasc) + procedure(TNO_Integrand) :: func !< integrand, evaluated as func(x, ctx) + real(TNOKi), intent(in ) :: a !< lower limit of integration + real(TNOKi), intent(in ) :: b !< upper limit of integration + type(TNO_ContextType), intent(inout) :: ctx !< context handed to the integrand + real(TNOKi), intent( out) :: result !< approximation to the integral, from the 61-point Kronrod rule + real(TNOKi), intent( out) :: abserr !< estimate of the modulus of the absolute error + real(TNOKi), intent( out) :: resabs !< approximation to the integral of abs(func) + real(TNOKi), intent( out) :: resasc !< approximation to the integral of abs(func-i/(b-a)) + + ! Local variables (names kept from qk61 to keep this reviewable against the original) + real(TNOKi) :: absc ! abscissa + real(TNOKi) :: centr ! mid point of the interval + real(TNOKi) :: dhlgth ! abs(hlgth) + real(TNOKi) :: epmach ! the largest relative spacing + real(TNOKi) :: fc ! function value at the mid point + real(TNOKi) :: fsum + real(TNOKi) :: fval1, fval2 ! function values + real(TNOKi) :: fv1(30), fv2(30) + real(TNOKi) :: hlgth ! half-length of the interval + real(TNOKi) :: resg ! result of the 30-point Gauss rule + real(TNOKi) :: resk ! result of the 61-point Kronrod rule + real(TNOKi) :: reskh ! approximation to the mean value of func over (a,b), i.e. to i/(b-a) + real(TNOKi) :: uflow ! the smallest positive magnitude + integer :: j, jtw, jtwm1 + + ! r1mach(4)/d1mach(4) is the largest relative spacing and r1mach(1)/d1mach(1) the smallest positive magnitude; + ! EPSILON and TINY are the standard intrinsics for exactly those quantities, so the SLATEC *1mach dependency is + ! not needed here. + epmach = epsilon(1.0_TNOKi) + uflow = tiny(1.0_TNOKi) + + centr = 0.5_TNOKi*(b+a) + hlgth = 0.5_TNOKi*(b-a) + dhlgth = abs(hlgth) + + ! Compute the 61-point Kronrod approximation to the integral, and estimate the absolute error. + resg = 0.0_TNOKi + fc = func(centr, ctx) + resk = wgk(31)*fc + resabs = abs(resk) + + do j = 1,15 + jtw = j*2 + absc = hlgth*xgk(jtw) + fval1 = func(centr-absc, ctx) + fval2 = func(centr+absc, ctx) + fv1(jtw) = fval1 + fv2(jtw) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(jtw)*fsum + resabs = resabs+wgk(jtw)*(abs(fval1)+abs(fval2)) + end do + + do j = 1,15 + jtwm1 = j*2-1 + absc = hlgth*xgk(jtwm1) + fval1 = func(centr-absc, ctx) + fval2 = func(centr+absc, ctx) + fv1(jtwm1) = fval1 + fv2(jtwm1) = fval2 + fsum = fval1+fval2 + resk = resk+wgk(jtwm1)*fsum + resabs = resabs+wgk(jtwm1)*(abs(fval1)+abs(fval2)) + end do + + reskh = resk*0.5_TNOKi + resasc = wgk(31)*abs(fc-reskh) + do j = 1,30 + resasc = resasc+wgk(j)*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh)) + end do + + result = resk*hlgth + resabs = resabs*dhlgth + resasc = resasc*dhlgth + abserr = abs((resk-resg)*hlgth) + if (resasc /= 0.0_TNOKi .and. abserr /= 0.0_TNOKi) & + abserr = resasc*min(1.0_TNOKi,(200.0_TNOKi*abserr/resasc)**1.5_TNOKi) + if (resabs > uflow/(50.0_TNOKi*epmach)) & + abserr = max((epmach*50.0_TNOKi)*resabs,abserr) + +end subroutine qk61_ctx +!================================================================================================================================== + + +FUNCTION f_int1(x2, ctx) result(f) + REAL(TNOKi), intent(in) :: x2 + type(TNO_ContextType), intent(inout) :: ctx + REAL(TNOKi):: f -FUNCTION f_int1(x2) REAL(TNOKi):: alpha REAL(TNOKi):: alpha_gauss REAL(TNOKi):: Cfin @@ -105,47 +292,64 @@ FUNCTION f_int1(x2) REAL(TNOKi):: Uc REAL(TNOKi):: Uo REAL(TNOKi):: W - REAL(TNOKi), intent(in) :: x2 - REAL(TNOKi):: f_int1 ! changed and being multiplied with edge velocity taken from xfoil output - ! Uo=Mach_TNO*co ISSUCTION_TNO use edgevel(1) + ! Uo=ctx%Mach*ctx%co ctx%IsSuction use ctx%edgevel(1) !constants from xfoil - if (ISSUCTION_TNO) then + if (ctx%IsSuction) then alpha = 0.45 ! = 0.3 pressure, = 0.45 suction - Cfin = Cf(1) - delta = d99(1) - Uo=Mach_TNO*co*edgevel(1) + Cfin = ctx%Cf(1) + delta = ctx%d99(1) + Uo=ctx%Mach*ctx%co*ctx%edgevel(1) else alpha = 0.30 - Cfin = Cf(2) - delta = d99(2) - Uo=Mach_TNO*co*edgevel(2) + Cfin = ctx%Cf(2) + delta = ctx%d99(2) + Uo=ctx%Mach*ctx%co*ctx%edgevel(2) endif - if (Cfin .le. 0.) then - write(*,*) 'Cf is less than zero, Cf = ',Cfin - stop + ! bjj: Bail out (contributing nothing to the integral) instead of producing NaN/Inf or killing the program. + ! - Cf <= 0 used to execute a bare "stop", which aborts without OpenFAST's error handling and without closing + ! output files. TBLTE_TNO already skips the side whose Cf is non-positive, so this is a backstop. + ! - delta (d99) can be zero for unconverged entries in the pre-tabulated boundary-layer files. That makes L = 0/0 + ! and pi*x2/delta infinite below. + ! - u_star is zero when the edge velocity ratio or the Mach number is zero, and log(u_star*x2/nu) is then -Inf, + ! so U evaluates to 0*(-Inf) = NaN. + ! - L is zero at x2 = 0 (the lower integration limit), which makes ke = sqrt(pi)/L infinite. + if (Cfin .le. 0. .or. delta .le. 0. .or. x2 .le. 0.) then + f = 0. + RETURN endif + u_star = Uo*sqrt(Cfin/2.) + if (u_star .le. 0.) then + f = 0. + RETURN + endif + L = 0.085*delta*tanh(kappa*x2/(0.085*delta)) + if (L .le. 0.) then + f = 0. + RETURN + endif + if (x2 .gt. delta)then U = Uo dudx = 0. - f_int1 = 0. + f = 0. RETURN else W = 1.-cos(pi*x2/delta); - U = u_star*(1./kappa*log(u_star*x2/nu) +Cnuk+ (Uo/u_star-1./kappa*log(u_star*delta/nu)-Cnuk)*0.5*W) - dudx = u_star*(1./(kappa*x2)+(Uo/u_star-1./kappa*log(u_star*delta/nu)-Cnuk)* & + U = u_star*(1./kappa*log(u_star*x2/ctx%nu) +Cnuk+ (Uo/u_star-1./kappa*log(u_star*delta/ctx%nu)-Cnuk)*0.5*W) + dudx = u_star*(1./(kappa*x2)+(Uo/u_star-1./kappa*log(u_star*delta/ctx%nu)-Cnuk)* & 0.5*(pi/delta)*sin(pi*x2/delta)) endif ke=sqrt(pi)/L*0.4213560764 !gamma(5./6.)/gamma(1./3.) - k1_hat = k1/ke - k3_hat = k3/ke + k1_hat = ctx%k1/ke + k3_hat = ctx%k3/ke Nut = (L*kappa)**2.*abs(dudx) kT = sqrt((Nut*dudx)**2./Cmu) @@ -154,45 +358,57 @@ FUNCTION f_int1(x2) Uc = 0.7*U alpha_gauss = 0.05*Uc/L - phim = 1./(alpha_gauss*sqrt(pi))*exp(-((Omega_TNO-Uc*k1)/alpha_gauss)**2.) + phim = 1./(alpha_gauss*sqrt(pi))*exp(-((ctx%Omega-Uc*ctx%k1)/alpha_gauss)**2.) phi22 = 4./9./pi*1/ke**2.*(k1_hat**2.+k3_hat**2.)/(1.+k1_hat**2.+k3_hat**2.)**(7./3.) - f_int1 = L*ums*(dudx)**2*phi22*phim*exp(-2*abs(k)*x2) + f = L*ums*(dudx)**2*phi22*phim*exp(-2*abs(ctx%k)*x2) RETURN END FUNCTION f_int1 -FUNCTION f_int2(k1_in) ! changed name from 'int2' to avoid conflicts with intrinsic of same name +FUNCTION f_int2(k1_in, ctx) result(f) ! changed name from 'int2' to avoid conflicts with intrinsic of same name REAL (TNOKi), intent(in) :: k1_in - REAL (TNOKi) :: f_int2 - f_int2 = Omega_TNO/co/k1_in*Pressure(k1_in) + type(TNO_ContextType), intent(inout) :: ctx + REAL (TNOKi) :: f + + ! bjj: The lower integration limit passed to qk61_ctx is exactly zero. The 61-point Gauss-Kronrod rule does not + ! evaluate the integrand at the interval end points, so k1_in is not zero in practice, but guard the 1/k1_in here + ! (and the k1**2/(k1**2+k3**2) = 0/0 in Pressure) so this does not depend on the internals of the quadrature rule. + if (k1_in .le. 0.0_TNOKi) then + f = 0.0_TNOKi + RETURN + endif + + f = ctx%Omega/ctx%co/k1_in*Pressure(k1_in, ctx) RETURN END FUNCTION f_int2 -FUNCTION Pressure(k1_in) +FUNCTION Pressure(k1_in, ctx) result(P) ! Variables + REAL(TNOKi), intent(in) :: k1_in + type(TNO_ContextType), intent(inout) :: ctx + real(TNOKi) :: P + REAL(TNOKi) :: a,b,answer REAL(TNOKi) :: abserr,resabs,resasc - REAL(TNOKi) :: k1_in - real(TNOKi) :: Pressure - ! Set variables used in f_int1 - k1 = k1_in + ! Set wavenumbers used in f_int1 + ctx%k1 = k1_in - a = 0.0_TNOKi !1e-4*d99(1) - IF (ISSUCTION_TNO)THEN - b = d99(1) + a = 0.0_TNOKi !1e-4*ctx%d99(1) + IF (ctx%IsSuction)THEN + b = ctx%d99(1) ELSE - b = d99(2) + b = ctx%d99(2) ENDIF - k3 = 0. - k= sqrt(k1**2+k3**2) + ctx%k3 = 0. + ctx%k = sqrt(ctx%k1**2+ctx%k3**2) - CALL slatec_qk61(f_int1,a,b,answer,abserr,resabs,resasc) + CALL qk61_ctx(f_int1,a,b,ctx,answer,abserr,resabs,resasc) - Pressure = 4.0_TNOKi*rho**2 * k1**2 / (k1**2 + k3**2)*answer + P = 4.0_TNOKi*ctx%rho**2 * ctx%k1**2 / (ctx%k1**2 + ctx%k3**2)*answer RETURN END FUNCTION Pressure diff --git a/modules/aerodyn/src/AirfoilInfo_Driver.f90 b/modules/aerodyn/src/AirfoilInfo_Driver.f90 new file mode 100644 index 0000000000..3c3662d0b4 --- /dev/null +++ b/modules/aerodyn/src/AirfoilInfo_Driver.f90 @@ -0,0 +1,210 @@ +!********************************************************************************************************************************** +! AFI_Driver: This code tests a stand-alone version of the AFI module +!.................................................................................................................................. +! LICENSING +! Copyright (C) 2018 Envision Energy +! +! This file is part of AirfoilInfo. +! +! Licensed under the Apache License, Version 2.0 (the "License"); +! you may not use this file except in compliance with the License. +! You may obtain a copy of the License at +! +! http://www.apache.org/licenses/LICENSE-2.0 +! +! Unless required by applicable law or agreed to in writing, software +! distributed under the License is distributed on an "AS IS" BASIS, +! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +! See the License for the specific language governing permissions and +! limitations under the License. +! +!********************************************************************************************************************************** + + + +program AFI_Driver + + use NWTC_Library + use VersionInfo + use AirfoilInfo + use AirfoilInfo_Types + + implicit none + + + TYPE(ProgDesc), PARAMETER :: AFI_Ver = ProgDesc( 'AFI_driver', '', '' ) + + + ! Variables + type(AFI_InitInputType) :: AFI_InitInputs ! Input data for initialization + integer, parameter :: NumAFI = 1; + type(AFI_ParameterType) :: AFI_p(NumAFI) ! Parameters + type(AFI_OutputType) :: AFI_interp ! interpolated AFI output values + integer(IntKi) :: ErrStat ! Status of error message + character(1024) :: ErrMsg ! Error message if ErrStat /= ErrID_None + + character(1024) :: afName + character(1024) :: outFileName + integer :: unOutFile(3) = -1 + character(*), parameter :: NumFmt = 'ES16.9E2' + character(*), parameter :: Frmt = '(1x,'//NumFmt//')' + character(*), parameter :: Ext(3) = (/ '.cl','.cd','.cm' /) + character(*), parameter :: delim = ' ' + + integer :: i, j, k, iFile + + real(ReKi), allocatable :: Re(:) + real(ReKi) :: alpha + real(ReKi), parameter :: UserProp = 0.0_ReKi + + + ! Initialize the NWTC library + call NWTC_Init(AFI_Ver%Name) + + ! Initialize error handling variables + ErrMsg = '' + ErrStat = ErrID_None + + + CALL DispNVD(AFI_Ver) + + + ! Check for command line arguments. + afName = '' ! default name for input file + CALL CheckArgs( afName ) + + CALL GetRoot( afName, outFileName ) + outFileName = trim(outFileName)//'.interp.out' + + + ! Setup Airfoil InitInput data structure (should come from an input file): + AFI_InitInputs%AFTabMod = AFITable_1 ! AFITable_2Re ! + AFI_InitInputs%InCol_Alfa = 1 + AFI_InitInputs%InCol_Cl = 2 + AFI_InitInputs%InCol_Cd = 3 + AFI_InitInputs%InCol_Cm = 4 + AFI_InitInputs%InCol_Cpmin = 0 + AFI_InitInputs%FileName = afName + AFI_InitInputs%UAMod = UA_Gonzalez + + + ! Write UA parameters to file: + call AFI_WrHeader(delim, trim(outFileName)//'.AFI.sum', unOutFile(1), ErrStat, ErrMsg) + + do iFile=1,NumAFI + if (NumAFI > 1) then + if (iFile < 10) then + AFI_InitInputs%FileName = 'af00'//trim(num2lstr(iFile))//'.dat' + else + AFI_InitInputs%FileName = 'af0'//trim(num2lstr(iFile))//'.dat' + end if + end if + + ! Initialize the Airfoil Info Params + call AFI_Init ( AFI_InitInputs, AFI_p(iFile), ErrStat, ErrMsg ) + call checkError() + + if (ErrStat < AbortErrLev) then + call AFI_WrData(iFile, unOutFile(1), delim, AFI_p(iFile)) + + call AFI_WrTables(AFI_p(iFile), AFI_InitInputs%UAMod, trim(AFI_InitInputs%FileName) ) + + end if + end do + close(unOutFile(1)) + + iFile = 1 + ! allocate Re array, based on Re in the tables + call allocAry( Re, AFI_p(iFile)%NumTabs*2 + 1, 'Re', ErrStat, ErrMsg ) + call checkError() + + Re(1) = AFI_p(iFile)%Table(1)%Re / 2.0_ReKi + do i = 1, AFI_p(iFile)%NumTabs + Re(2*i) = AFI_p(iFile)%Table(i)%Re + end do + do i = 1, AFI_p(iFile)%NumTabs-1 + Re(2*i+1) = (AFI_p(iFile)%Table(i)%Re + AFI_p(iFile)%Table(i+1)%Re)/2.0_ReKi + end do + Re(size(Re)) = AFI_p(iFile)%Table(AFI_p(iFile)%NumTabs)%Re * 2 + + + ! ------------ + do k=1,size(Ext) + call GetNewUnit( unOutFile(k) ) + call OpenFOutFile ( unOutFile(k), trim(outFileName)//Ext(k), errStat, errMsg ) + call checkError() + + write( unOutFile(k), '('//trim(num2lstr(size(Re)+1))//Frmt//')' ) NaN, Re + end do + + + ! time marching loop + do i = -180, 180 ! alpha + alpha = Real(i, ReKi) ! degrees + + do k=1,size(Ext) + write(unOutFile(k), Frmt, ADVANCE='no') alpha + end do + + alpha = alpha*pi/180.0_ReKi ! radians + + do j=1, size(Re) + + call AFI_ComputeAirfoilCoefs( alpha, Re(j), UserProp, AFI_p(iFile), AFI_interp, ErrStat, ErrMsg) + call checkError() + + write(unOutFile(1), Frmt, ADVANCE='no') AFI_interp%cl + write(unOutFile(2), Frmt, ADVANCE='no') AFI_interp%cd + write(unOutFile(3), Frmt, ADVANCE='no') AFI_interp%cm + + end do + + do k=1,size(Ext) + write (unOutFile(k),'()', IOSTAT=ErrStat) ! write the line return + end do + + end do + + + !------------------------------------------------------------------------------------------------- + ! Close our output files + !------------------------------------------------------------------------------------------------- + + + call Cleanup() + call NormStop() + + contains + + !==================================================================================================== + subroutine Cleanup() + ! The routine closes any open files. + !---------------------------------------------------------------------------------------------------- + integer :: ie + + do ie=1,size(Ext) + if (unOutFile(ie) > 0) close( unOutFile(ie), IOSTAT = ErrStat ) + end do + + end subroutine Cleanup + + !---------------------------------------------------------------------------------------------------- + subroutine checkError() + + if (ErrStat >= AbortErrLev) then + + call Cleanup() + call ProgAbort(ErrMsg) + + elseif ( ErrStat /= ErrID_None ) then + + call WrScr( trim(ErrMsg) ) + + end if + + end subroutine checkError + !---------------------------------------------------------------------------------------------------- + + +end program AFI_Driver + diff --git a/modules/nwtc-library/CMakeLists.txt b/modules/nwtc-library/CMakeLists.txt index c40457fcad..51a43b70f8 100644 --- a/modules/nwtc-library/CMakeLists.txt +++ b/modules/nwtc-library/CMakeLists.txt @@ -106,20 +106,20 @@ set(NWTCLIBS_SOURCES src/NetLib/scalapack/NWTC_ScaLAPACK.f90 # NetLib SLATEC sources - src/NetLib/slatec/NWTC_SLATEC.f90 - src/NetLib/slatec/dqk61.f - src/NetLib/slatec/qk61.f - src/NetLib/slatec/d1mach.f - src/NetLib/slatec/r1mach.f - src/NetLib/slatec/xercnt.f - src/NetLib/slatec/xerhlt.f - src/NetLib/slatec/xerprn.f - src/NetLib/slatec/xersve.f - src/NetLib/slatec/fdump.f - src/NetLib/slatec/i1mach.f - src/NetLib/slatec/j4save.f - src/NetLib/slatec/xgetua.f - src/NetLib/slatec/xermsg.f + #src/NetLib/slatec/NWTC_SLATEC.f90 + #src/NetLib/slatec/dqk61.f + #src/NetLib/slatec/qk61.f + #src/NetLib/slatec/d1mach.f + #src/NetLib/slatec/r1mach.f + #src/NetLib/slatec/xercnt.f + #src/NetLib/slatec/xerhlt.f + #src/NetLib/slatec/xerprn.f + #src/NetLib/slatec/xersve.f + #src/NetLib/slatec/fdump.f + #src/NetLib/slatec/i1mach.f + #src/NetLib/slatec/j4save.f + #src/NetLib/slatec/xgetua.f + #src/NetLib/slatec/xermsg.f ) get_filename_component(FCNAME ${CMAKE_Fortran_COMPILER} NAME) diff --git a/modules/subdyn/src/SubDyn_Types.f90 b/modules/subdyn/src/SubDyn_Types.f90 index 9fee52b5f2..5bdaa7d288 100644 --- a/modules/subdyn/src/SubDyn_Types.f90 +++ b/modules/subdyn/src/SubDyn_Types.f90 @@ -170,7 +170,7 @@ MODULE SubDyn_Types REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberSpin !< Member spin angle about its axis - for rectangular members [rad] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: MemberDivSize !< Optional maximum element length for each member [m] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberNDiv !< Resolved number of finite elements per member [-] - INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] + INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: MemberElemStart !< First element index for each member in p%Elems [-] CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: SSOutList !< List of Output Channels [-] LOGICAL :: OutCOSM = .false. !< Output Cos-matrices Flag [-] LOGICAL :: TabDelim = .false. !< Generate a tab-delimited output file in OutJckF-Flag [-] @@ -196,7 +196,7 @@ MODULE SubDyn_Types REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: K !< System stiffness matrix [-] REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: M !< System mass matrix [-] REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: ElemProps !< Element properties(A, L, Ixx, Iyy, Jzz, Shear, Kappa, E, G, Rho, DirCos(1,1), DirCos(2, 1), ....., DirCos(3, 3) ) [-] - INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and list of nodes making up a member (>2 if subdivided) [-] + INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: MemberNodes !< Member number and endpoint node IDs (interior member nodes reconstructed from connectivity) [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnN !< Nodes that connect to a common node [-] INTEGER(IntKi) , DIMENSION(:,:), ALLOCATABLE :: NodesConnE !< Elements that connect to a common node [-] LOGICAL :: SSSum = .false. !< SubDyn Summary File Flag [-] diff --git a/vs-build/.gitignore b/vs-build/.gitignore deleted file mode 100644 index 5253e137da..0000000000 --- a/vs-build/.gitignore +++ /dev/null @@ -1,4 +0,0 @@ -.vs -*.user -*.u2d -gitVersionInfo.h \ No newline at end of file diff --git a/vs-build/OpenFAST.sln b/vs-build/OpenFAST.sln index 8ec8938fcb..08412001a9 100644 --- a/vs-build/OpenFAST.sln +++ b/vs-build/OpenFAST.sln @@ -589,6 +589,20 @@ Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "SoilDyn_Driver", "drivers\S {F8C26AB7-4B2F-4610-A18B-A8C860926B21} = {F8C26AB7-4B2F-4610-A18B-A8C860926B21} EndProjectSection EndProject +Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "AirfoilInfo_Driver", "drivers\AirfoilInfo_Driver.vfproj", "{3BBE2741-5B28-47BC-9E7F-3E1D172838FB}" + ProjectSection(ProjectDependencies) = postProject + {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} = {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} + {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} = {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} + {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} = {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} + EndProjectSection +EndProject +Project("{6989167D-11E4-40FE-8C1A-2192A86A7E90}") = "AeroAcoustics_Driver", "drivers\AeroAcoustics_Driver.vfproj", "{CFB604E9-0D11-4903-BF4C-97D411400C33}" + ProjectSection(ProjectDependencies) = postProject + {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} = {12DF411B-C7DA-47BA-BB85-7714D5FD2A16} + {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} = {5ADBD025-C654-42C2-BA7C-10F3C3CEEB0E} + {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} = {EAF5E602-E6CD-4194-8CCA-0827AA4CCEC9} + EndProjectSection +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|x64 = Debug|x64 @@ -2390,6 +2404,70 @@ Global {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x64.Build.0 = Release|x64 {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x86.ActiveCfg = Release|x64 {9779535B-8DE4-4484-8B59-6E78344D658A}.Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x64.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x64.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x86.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Debug|x86.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x64.ActiveCfg = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x64.Build.0 = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x86.ActiveCfg = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Debug|x86.Build.0 = Double_Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x64.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x64.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x86.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_OpenMP_Release|x86.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x64.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x64.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x86.ActiveCfg = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Double_Release|x86.Build.0 = Double_Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x64.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x64.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x86.ActiveCfg = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Debug|x86.Build.0 = Debug|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Matlab_Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.OpenMP_Release|x86.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x64.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x64.Build.0 = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x86.ActiveCfg = Release|x64 + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB}.Release|x86.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x64.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x64.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x86.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Debug|x86.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x64.ActiveCfg = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x64.Build.0 = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x86.ActiveCfg = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Debug|x86.Build.0 = Double_Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x64.ActiveCfg = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x64.Build.0 = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x86.ActiveCfg = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_OpenMP_Release|x86.Build.0 = Double_OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x64.ActiveCfg = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x64.Build.0 = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x86.ActiveCfg = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Double_Release|x86.Build.0 = Double_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x64.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x64.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x86.ActiveCfg = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Debug|x86.Build.0 = Debug|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x64.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x64.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x86.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Matlab_Release|x86.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x64.ActiveCfg = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x64.Build.0 = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x86.ActiveCfg = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.OpenMP_Release|x86.Build.0 = OpenMP_Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x64.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x64.Build.0 = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x86.ActiveCfg = Release|x64 + {CFB604E9-0D11-4903-BF4C-97D411400C33}.Release|x86.Build.0 = Release|x64 EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE @@ -2453,6 +2531,8 @@ Global {4AB5A895-5B6B-406D-AE36-85F5F5940974} = {76DE56EB-E48F-4F1A-8961-F97A71E7937B} {F8C26AB7-4B2F-4610-A18B-A8C860926B21} = {272B8080-A022-4F4A-BDD6-835871E44C23} {9779535B-8DE4-4484-8B59-6E78344D658A} = {3517E990-350F-4471-A518-8B0BC77CFDDB} + {3BBE2741-5B28-47BC-9E7F-3E1D172838FB} = {3517E990-350F-4471-A518-8B0BC77CFDDB} + {CFB604E9-0D11-4903-BF4C-97D411400C33} = {3517E990-350F-4471-A518-8B0BC77CFDDB} EndGlobalSection GlobalSection(ExtensibilityGlobals) = postSolution SolutionGuid = {B362252D-3254-4C68-B527-CC85CE3CCF75} diff --git a/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj b/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj index 126c78c66f..b61ec0074c 100644 --- a/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj +++ b/vs-build/c-bindings/AeroDyn_Inflow_C_Binding.vfproj @@ -74,23 +74,23 @@ - - + + + + + - + - + - + - - - - - + + diff --git a/vs-build/drivers/AeroAcoustics_Driver.vfproj b/vs-build/drivers/AeroAcoustics_Driver.vfproj new file mode 100644 index 0000000000..c408345078 --- /dev/null +++ b/vs-build/drivers/AeroAcoustics_Driver.vfproj @@ -0,0 +1,82 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/vs-build/drivers/AirfoilInfo_Driver.vfproj b/vs-build/drivers/AirfoilInfo_Driver.vfproj new file mode 100644 index 0000000000..6319d2a2ee --- /dev/null +++ b/vs-build/drivers/AirfoilInfo_Driver.vfproj @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/vs-build/modules/NWTC-Library.vfproj b/vs-build/modules/NWTC-Library.vfproj index 22bb1fce7c..db76bb4258 100644 --- a/vs-build/modules/NWTC-Library.vfproj +++ b/vs-build/modules/NWTC-Library.vfproj @@ -115,102 +115,39 @@ - - - - - - + + - + - + - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - + - - - - - - - + - + - + - - - - - - + + + + + @@ -218,6 +155,8 @@ + + @@ -257,8 +196,6 @@ - -