Description
I have been experiencing a segmentation fault in the dragcof() subroutine of the Fitch model, which causes a segmentation fault for certain value of the wind speed. The bug depends on the solution fields and most of the times it does not even occur. I fixed the bug by correcting the initialization of nu and nb in the dragcof() function at line 294 of module_wind_fitch.F
To Reproduce
The behavior is very difficult to reproduce. I have been experiencing it only with certain types of turbines for the date of 2017/02/23 at 15:00 UTC in the Belgian Bight of the North Sea. The simulation runs fine if wind turbines are removed.
The bug still appears if the power and thrust coeff on the last row of the turbine table are set to zero.
To Fix
The bug was fixed by changing line 294 of module_wind_fitch.F from
DO k=1,maxvals2
to
DO k=2,maxvals2
The rationale behind this change is that if the loop exits at k=1, then nu=1 and nb=0, which is invalid in FORTRAN. There are guards to prevent the wind speed to be below cut in, but I don't think they are sufficient.
On a different note, the current implementation of the explicit turbine source terms is numerically unstable, because the source term can cause the velocity to be negative if the delta in velocity is larger than udt or vdt. Hence, I propose the following changes, change lines 228 to 230 from
! output u tendency
du(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
! output v tendency
dv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
to
! output u tendency (limited to prevent sign reversal / numerical instability)
du_drag = -.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
if (abs(du_drag)*dt .gt. abs(u(i,k,j))) du_drag = -u(i,k,j)/dt
du(i,k,j) = du(i,k,j) + du_drag
! output v tendency (limited to prevent sign reversal / numerical instability)
dv_drag = -.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec
if (abs(dv_drag)*dt .gt. abs(v(i,k,j))) dv_drag = -v(i,k,j)/dt
dv(i,k,j) = dv(i,k,j) + dv_drag
where du_drag and dv_drag should be declare at the beginning as real :: du_drag,dv_drag
Expected behavior
This would avoid segmentation fault in the dragcof function (the actual bug) and make the introduction of tendencies through the explicit term numerically stable.
Description
I have been experiencing a segmentation fault in the dragcof() subroutine of the Fitch model, which causes a segmentation fault for certain value of the wind speed. The bug depends on the solution fields and most of the times it does not even occur. I fixed the bug by correcting the initialization of nu and nb in the dragcof() function at line 294 of module_wind_fitch.F
To Reproduce
The behavior is very difficult to reproduce. I have been experiencing it only with certain types of turbines for the date of 2017/02/23 at 15:00 UTC in the Belgian Bight of the North Sea. The simulation runs fine if wind turbines are removed.
The bug still appears if the power and thrust coeff on the last row of the turbine table are set to zero.
To Fix
The bug was fixed by changing line 294 of module_wind_fitch.F from
DO k=1,maxvals2to
DO k=2,maxvals2The rationale behind this change is that if the loop exits at k=1, then nu=1 and nb=0, which is invalid in FORTRAN. There are guards to prevent the wind speed to be below cut in, but I don't think they are sufficient.
On a different note, the current implementation of the explicit turbine source terms is numerically unstable, because the source term can cause the velocity to be negative if the delta in velocity is larger than udt or vdt. Hence, I propose the following changes, change lines 228 to 230 from
! output u tendencydu(i,k,j) = du(i,k,j)-.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ec! output v tendencydv(i,k,j) = dv(i,k,j)-.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ecto
! output u tendency (limited to prevent sign reversal / numerical instability)du_drag = -.5*u(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ecif (abs(du_drag)*dt .gt. abs(u(i,k,j))) du_drag = -u(i,k,j)/dtdu(i,k,j) = du(i,k,j) + du_drag! output v tendency (limited to prevent sign reversal / numerical instability)dv_drag = -.5*v(i,k,j)*thrcof*speed*tarea/dz(i,k,j)*ecif (abs(dv_drag)*dt .gt. abs(v(i,k,j))) dv_drag = -v(i,k,j)/dtdv(i,k,j) = dv(i,k,j) + dv_dragwhere du_drag and dv_drag should be declare at the beginning as
real :: du_drag,dv_dragExpected behavior
This would avoid segmentation fault in the dragcof function (the actual bug) and make the introduction of tendencies through the explicit term numerically stable.