Catch divide-by-zero error in wind speed estimator (#410)

* ability to load rosco yaml on rank 0 alone * Catch divide-by-zero in WSE --------- Co-authored-by: ptrbortolotti <[email protected]>
NREL · Dec 27, 2024 · fcc4394 · fcc4394
1 parent 50e8247
commit fcc4394
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 2 deletions.
diff --git a/Examples/28_tower_resonance.py b/Examples/28_tower_resonance.py
@@ -65,7 +65,7 @@ def main():
     # Ramp: good demo of functionality, short for CI
     r.wind_case_fcn = cl.ramp  
     r.wind_case_opts    = {
-        'U_start': 4,  # from 10 to 15 m/s
+        'U_start': 0,  # from 10 to 15 m/s
         'U_end': 10,
         't_start': 100,
         't_end': 300

diff --git a/rosco/controller/src/ControllerBlocks.f90 b/rosco/controller/src/ControllerBlocks.f90
@@ -282,6 +282,7 @@ SUBROUTINE WindSpeedEstimator(LocalVar, CntrPar, objInst, PerfData, DebugVar, Er
                 LocalVar%WE%v_t = 0.0
                 LocalVar%WE%v_m = max(LocalVar%HorWindV, 3.0_DbKi)   ! avoid divide by 0 below if HorWindV is 0, which some AMRWind setups create
                 LocalVar%WE%v_h = max(LocalVar%HorWindV, 3.0_DbKi)   ! avoid divide by 0 below if HorWindV is 0, which some AMRWind setups create
+                LocalVar%WE_Vw = LocalVar%WE%v_m + LocalVar%WE%v_t   ! Initialize WE_Vw to aviod divide by zero
                 lambda = WE_Inp_Speed * CntrPar%WE_BladeRadius/LocalVar%WE%v_h
                 LocalVar%WE%xh = RESHAPE((/LocalVar%WE%om_r, LocalVar%WE%v_t, LocalVar%WE%v_m/),(/3,1/))
                 LocalVar%WE%P = RESHAPE((/0.01, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 1.0/),(/3,3/))

diff --git a/rosco/controller/src/Functions.f90 b/rosco/controller/src/Functions.f90
@@ -575,12 +575,14 @@ REAL(DbKi) FUNCTION AeroDynTorque(RotSpeed, BldPitch, LocalVar, CntrPar, PerfDat
         REAL(DbKi) :: RotorArea
         REAL(DbKi) :: Cp
         REAL(DbKi) :: Lambda
+        REAL(DbKi) :: WindSpeed
 
         CHARACTER(*), PARAMETER                 :: RoutineName = 'AeroDynTorque'
 
         ! Find Torque
         RotorArea = PI*CntrPar%WE_BladeRadius**2
-        Lambda = RotSpeed*CntrPar%WE_BladeRadius/LocalVar%WE_Vw
+        WindSpeed = MAX(LocalVar%WE_Vw,EPSILON(1.0_DbKi))
+        Lambda = RotSpeed*CntrPar%WE_BladeRadius/WindSpeed
 
         ! Compute Cp
         Cp = interp2d(PerfData%Beta_vec,PerfData%TSR_vec,PerfData%Cp_mat, BldPitch*R2D, Lambda, ErrVar)