Floating Point Exception (SIGFPE) for iprobes =2

[rickyspaceguy]

SIGFPE generated for input_HOS-NWT.dat in the the bin folder for iprobes=2.
No problems if iprobes=1
Program compiled with flags -fcheck=all -ffpe-trap=invalid,zero,overflow in gfortran 5.4.0 20160609

Console output just before the exception

CPU time for Output time step 2.76 9 0
volume = 2.3105E-04 t= 5.4000E+00 timestep: 1.1926E-02 486 0
H2 operator
KFSBC/DFSBC max. initial = 6.8966151695545354E-004 , reconstructed = 1.1740512391557623E-002
Max relative error = 0.22180840568248397 ... direct method

SIGFPE reads as

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0 0x7FAC399DDE08
#1 0x7FAC399DCF90
#2 0x7FAC3910E4AF
#3 0x4FF3CA in __velocities_MOD_creatematrix_coscos_sinh at velocities.f90:1517 (discriminator 56)
#4 0x5034B3 in __velocities_MOD_velocity_modes_direct at velocities.f90:1299
#5 0x5038EA in __velocities_MOD_hosvel2_direct at velocities.f90:1252
#6 0x45E8B6 in runge_kutta_MOD_runge4_adapt2 at runge_kutta.f90:342
#7 0x464D25 in MAIN at HOS-NWT.f90:520

input_HOS-NWT.txt

[gducrozet]

There was indeed a problem in the current version of the code. The treatment of the case k=0 was done afterwards (without preventing a division over 0). It is actually only detected thanks to the compilation option '-ffpe=invalid' that I was not using.

By the way, to correct it, you should change lines 1510 to 1527 of velocities.f90 to the following:

DO i2=1,1
DO i1=1,n1
index1=i1+n1*(i2-1)
DO ii2=1,1
DO ii1=1,n1
index2=ii1+n1*(ii2-1)
IF(k(ii1,ii2).EQ.0.0_rp) THEN ! Specific treatment of k=0
M(index1,index2) = 1.0_rp
ELSEIF(k(ii1,ii2).LT.50.0_rp) THEN
M(index1,index2) = COS(kx(ii1)*x(i1,i2))*COS(ky(ii2)*y(i1,i2))SINH(k(ii1,ii2)(eta(i1,i2)+1.0_rp))/SINH(k(ii1,ii2))
ELSE
M(index1,index2) = COS(kx(ii1)*x(i1,i2))*COS(ky(ii2)*y(i1,i2))*EXP(k(ii1,ii2)*eta(i1,i2))
ENDIF
ENDDO
ENDDO
ENDDO
ENDDO

[rickyspaceguy]

Shouldn't it be k(ii1,ii2).LE.tiny

[gducrozet]

Yes indeed, would be safer.