Issues with a fortran code in linux debian, solving the flow inside a double convergent-divergent nozzle

Help Homework
1

So I am trying to solve the Euler 1-dimensional equations of fluid flow through a double convergent-divergent nozzle. This code supposedly should work fine, and write the file results.dat with the expected results for density, velocity and pressure. However, all I get is NaN in all cases. I use ifort as a compiler.

Here is the entire code:


module constants

double precision ::GAM, GAMM1 ,rkconst(4),conA,conB,conC,conK,R, T0
double precision, allocatable :: xc(:),xf(:),Sc(:),Sf(:), Scp(: )

end module constants

Program Euler1d

  use constants

  Implicit None

  double precision, allocatable :: U(:,:),U0(:,:),cell_flux(:,:),P(: )

  double precision :: XMIN, XMAX, DX, DT,TIME

  double precision :: UCL(3),UCR(3),UPL(3),UPR(3), NORM,flux(3)

  integer          :: i, j, k, NX,NXC,NTIME, NTIMETOT

R = 287.05;

T0 = 292;

GAM=1.4; GAMM1=GAM - 1

conA= 0.456

conB= 1.09

conC= 0.91

conK= 0.19

!Define runge kutta constants

  rkconst(1) = 0.1084d0

  rkconst(2) = 0.2602d0

  rkconst(3) = 0.5052d0

  rkconst(4) = 1.d0

!Input Parameters

  NTIMETOT = 5000            ! Total timesteps

  XMIN = 0.0                 ! XMIN

  XMAX = 2.0                 ! XMAX

  DX = 1.0e-2                ! DX

  DT = 1.0e-6                ! DT

!Finite Volume Scheme-Nodes represent faces

!Number of points/faces

  NX = int((XMAX-XMIN)/DX)

!Calculate new DX for accurate division by NX

  DX = (XMAX-XMIN)/(NX-1)

  !write (*,*) 'New DX =', DX

!Number of cell centers

  NXC= NX -1

! 0: represents boundary cell at start and NXC+1 boundary cell at end

! U is the solution vector (Dens, Dens *U, Dens *E), U0 the vector after the 1st Runge-Kutta step

! cell_flux is the sum of the fluxes for each cell

  allocate (U(3,0:NXC+1),U0(3,0:NXC+1))

  allocate (cell_flux(1:3,0:NXC+1))
  
  allocate (P(1:NXC))
 
allocate(xf(1:NX))
allocate(xc(1:NXC+1))
allocate(Sf(1:NX))
allocate(Sc(0:NXC+1))
allocate(Scp(1:NXC+1))

do i=1,NX
xf(i)= XMIN + (i-1)*DX

xc(i)= xf(i) + DX/2
Sf(i)= conK+conA*((xf(i)-conC)**2)*((xf(i)-conC)**2-conB**2)

Sc(i)= conK+conA*((xc(i)-conC)**2)*((xc(i)-conC)**2-conB**2)
Scp(i)=4*conA*(xc(i)-conC)**3-(2*conA*conB**2)*(xc(i)-conC)
enddo
Sc(0)= Sc(1)

TIME=0

!Set initial conditions

  call ic

  do ntime = 1, NTIMETOT

     k=1

!k 1:4 four step RK

     do while (k.le.4) 

!Set boundary conditions

        call bc

        cell_flux=0.d0

!For all Faces

        do i=1,NX

           !Left cell

           UCL = U(1:3,i-1)        

           !Right cell

           UCR = U(1:3,i)        

           !Get primitive variables

           call co_to_pri(UCL,UPL,Sf(i))

           call co_to_pri(UCR,UPR,Sf(i))

           NORM=1

           !calculate Roe-flux

           call calc_flux(UPL, UPR, UCL, UCR,NORM,Flux)

           !Normal points from left to right

           !For bc's it's obsolete

           cell_flux(1:3,i-1) = cell_flux(1:3,i-1)  + Flux  

           cell_flux(1:3,i)   = cell_flux(1:3,i)    - Flux

        enddo

!    Time integration

       if(k.eq.1) U0(1:3,1:NXC)=U(1:3,1:NXC)
       do i=1,NXC
        P(i) = GAMM1 * (U(3,i) - 0.5d0 / U(1,i) * U(2,i)**2/Sf(i))
        U(1:3,i) = U0(1:3,i) - DT/DX *rkconst(k)*(cell_flux(1:3,i) - cell_flux(1:3,i-1))
        U(2,i) = U(2,i) + (DT*P(i))*Scp(i)
       enddo
        k=k+1

     enddo

     TIME =TIME + DT

     write(*,*) 'NTIME= ',NTIME,'TIME= ', TIME

!write results 

     if (mod(NTIME,100).eq.0) call writeres(TIME)

  enddo

contains

!Set initial conditions

Subroutine ic

 Implicit None

 double precision :: UC(3),UP(3)

!Sod's shock tube initial conditions

! set pressure, velocity, density and then transform to conservative

 do i=1,NXC

       UP(3)= 3.2e5; UP(2) = 0.d0; UP(1) = UP(3)/(R*T0);

       call pri_to_co(UC,UP,Sf(i))

       U(1:3,i) = UC(1:3)

 enddo

End Subroutine ic

!Set boundary conditions at ends

Subroutine bc

double precision :: UP(3),UC(3)

!REFLECTIVE BOUNDARY CONDITIONS (as defined by the problem)

!give bc in a primitive form 

!start

 UP(1) =   4e5/(R*T0) !Density

 UP(2) =  U(2,1)/U(1,1) !Velocity

 UP(3) =   4e5   !Pressure

 call pri_to_co(UC,UP,Sc(0))

 U(1,0) =  UC(1) !Density

 U(2,0) =  UC(2) !Velocity

 U(3,0) =  UC(3) !Pressure

!end

!OUTLET

 UP(1) =  U(1,NXC)/(conK+conA*((2-conC)**2)*((2-conC)**2-conB**2)) !Density

 UP(2) =  U(2,NXC)/U(1,NXC) !Velocity

 UP(3) =  3.2e5     !Pressure

 call pri_to_co(UC,UP,Sc(NXC+1))

!end

 U(1,NXC+1) =  UC(1) !Density

 U(2,NXC+1) =  UC(2) !Velocity

 U(3,NXC+1) =  UC(3) !Pressure

end Subroutine bc

!write results in primitive variable form

Subroutine writeres(TIME)

double precision :: TIME,UP(3),UC(3)

  open (1,file='results.dat')

  write(1,'(a)') '#Time, X, Density, Velocity, Pressure'

  do i=1,NXC

     UC(1:3) = U(1:3,i)

     call co_to_pri(UC,UP,Sc(i))

     write(1,'(5(e28.17,1x))') TIME,XMIN + (i-1) * DX + 0.5 * DX, UP(1:3)

  enddo

  write(1,*) 

End Subroutine writeres

! Flux:Calculated Roe flux

Subroutine calc_flux(UPL, UPR, UCL, UCR,NORM,Flux)
use constants

  Implicit None

  double precision, intent(in)  :: UPL(3), UPR(3),UCL(3),UCR(3),NORM

  double precision, intent(out) :: Flux(3)

  double precision              :: Flux_left(3), Flux_right(3), Aroe(3,3), DF1(3)

!1D fluxes

  call flux1d(UPL,UCL,NORM,Flux_left,Sf(i))

  call flux1d(UPR,UCR,NORM,Flux_right,Sf(i))

!Roe matric

  call roematrix(UPL,UCL,UPR,UCR,NORM,Aroe,Sf(i))

  DF1(1:3) = UCR(1:3)-UCL(1:3)

!Final flux 

 Flux(1 : 3) = 0.5d0 * (Flux_left(1:3) + Flux_right(1:3) - matmul(Aroe,DF1))

End subroutine calc_flux

End Program Euler1d

!==============================================================================

!Subroutine to  convert conservative (UC) to primitive (UP)

Subroutine co_to_pri(UC,UP,S)

  use constants

  Implicit None

  double precision, intent(in)  :: UC(3),S

  double precision, intent(out) :: UP(3)

  UP(1) = UC(1)/S

  UP(2) = UC(2)/UC(1)

  UP(3) = GAMM1 * (UC(3) - 0.5d0 / UC(1) * UC(2)**2)/S

End Subroutine co_to_pri

!Subroutine to  convert primitive(UP) to conservative (UC)

Subroutine pri_to_co(UC,UP,S)

  use constants

  Implicit None

  double precision, intent(out)  :: UC(3)

  double precision, intent(in)   :: UP(3),S

  UC(1) = UP(1)*S

  UC(2) = UP(2)*UP(1)*S

  UC(3) = (UP(3) / GAMM1 + 0.5d0 * UP(1) * UP(2)**2)*S

End Subroutine pri_to_co

!Subroutine to calculate Roe flux

! UPL :Primitive Variable left 

! UPR :Primitive Variable right

! UCL :Conservative Variable left

! UCR :Conservative Variable Right

! NORM:Normal

!1D euler equaations flux functions

! UP :Primitive Variable

! UC :Conservative Variable 

! NORM:Normal

Subroutine flux1d(UP,UC,NORM,Flux,S)

  Implicit None

  double precision, intent(in)  :: UP(3), UC(3),NORM, S

  double precision, intent(out) :: Flux(3)

  double precision :: VN

  VN =  UP(2) * NORM

  Flux(1) = UP(1) * VN * S

  Flux(2) = (UP(1) * UP(2)* VN + UP(3) * NORM) * S

  Flux(3) = (UC(3) + UP(3)) * VN * S

End subroutine flux1d

!Subroutine to calculate Roe matrix

! UPL :Primitive Variable left 

! UPR :Primitive Variable right

! UCL :Conservative Variable left

! UCR :Conservative Variable Right

! NORM:Normal

! Aroe: ROe matrix in conservative variables

Subroutine roematrix(UPL,UCL,UPR,UCR,NORM,Aroe,S)

  use constants

  Implicit None

  double precision, intent(in)  :: UPL(3),UCL(3),UPR(3),UCR(3),NORM, S

  double precision, intent(out) :: Aroe(3,3)

  double precision              :: DensRoe, Uroe, Hl,Hr, Hroe,croe,L1,L2,L3,delta

  double precision              :: left(3,3),right(3,3),eigen(3,3)

!Calculate Roe averaged variables (Density,Velocity,Enthalpy)

  DensRoe = sqrt(UPL(1))*sqrt(UPR(1))
  
  Uroe    = (sqrt(UPL(1)) * UPL(2) + sqrt(UPR(1)) * UPR(2)) / ( sqrt(UPL(1)) + sqrt(UPR(1)) )

  Hl     = (UCL(3)/S + UPL(3))/UPL(1)

  Hr     = (UCR(3)/S + UPR(3))/UPR(1)

  Hroe    = (sqrt(UPL(1)) * Hl + sqrt(UPR(1)) * Hr) / ( sqrt(UPL(1)) + sqrt(UPR(1)) )

  croe    =  sqrt( GAMM1 * (Hroe - 0.5d0 * Uroe**2) )                           

!Get left and right eigenvectors in conservative form

  call leftrightcons(DensRoe,croe,Uroe,left,right)

  eigen     = 0.d0

   L1=abs(Uroe)

   L2=abs(Uroe+Croe)

   L3=abs(Uroe-Croe)

!Harten's entropy fix  and fill eigenvalue matrix

   delta = 0.1d0 * 0.5d0*Croe

   if (L1.le.delta) L1 = (L1**2 + delta**2) / (2 * delta)

   if (L2.le.delta) L2 = (L2**2 + delta**2) / (2 * delta)

   if (L3.le.delta) L3 = (L3**2 + delta**2) / (2 * delta)

   eigen=0

   eigen(1,1)= L1

   eigen(2,2)= L2

   eigen(3,3)= L3

!get Roe matrix

   Aroe = matmul(right,eigen)

   Aroe = matmul(Aroe,left)

End Subroutine roematrix

!Subroutine to calculate left and right eigenvectors in conservative form

Subroutine leftrightcons(Dens,C,U,left,right)

   use constants

   Implicit None

   double precision, intent(in) :: Dens,U,C

   double precision, intent(out):: left(3,3),right(3,3)

   left(1,1) =  1.d0 - 0.5d0 * GAMM1 * U**2/C**2

   left(1,2) =  GAMM1 * U /C**2

   left(1,3) = -GAMM1/C**2

   left(2,1) = (0.5d0*GAMM1*U**2 - U*C) * (1.d0/(Dens*C))

   left(2,2) = (C - GAMM1*U) *(1.d0/(Dens*C))

   left(2,3) = GAMM1/(Dens*C)

   left(3,1) = -(0.5d0*GAMM1*U**2 + U*C) * (1.d0/(Dens*C))

   left(3,2) = (C + GAMM1*U) *(1.d0/(Dens*C))

   left(3,3) = -GAMM1/(Dens*C)

   right(1,1) =  1.d0

   right(1,2) =  0.5d0 * Dens/C

   right(1,3) = -0.5d0 * Dens/C
   
   right(2,1) = U

   right(2,2) = 0.5d0 * (U+C) * Dens/C

   right(2,3) =-0.5d0 * (U-C) * Dens/C

   right(3,1) = 0.5d0 * U**2

   right(3,2) = (0.5d0 * U**2 + U*C + C**2/GAMM1) * 0.5d0 * Dens/C

   right(3,3) =-(0.5d0 * U**2 - U*C + C**2/GAMM1) * 0.5d0 * Dens/C

End Subroutine leftrightcons```
2

Welcome to the Discourse @Oresibius !

Can you edit your message and put your whole code between triple quotes (with the name of the language besides the first triplet for syntax colouring):

```fortran

and

    ```

It will help people to read and analyze it.

1 Like
3

Same problem here in GFortran (the NANs appear in the results.dat file):

     0.49999999873762135E-02      0.19849246231155779E+01   NaN    NaN   NaN
     0.49999999873762135E-02      0.19949748743718592E+01   NaN    NaN   NaN
4

with a bit of editing, I found out that the NaN values start when TIME = 0.53999E-4 approximately, currently trying to find where the issue is (maybe division by zero?)

5

Or a square root of a negative number? There are many sqrt() in the Subroutine roematrix.

6

Run-time diagnostics are your friend. Enable the compiler option to trap NaNs. For ifort this is /fpe:0, or if using the Visual Studio GUI it is “Floating Point Exception Handling” => Underflow gives 0.0; Abort on other IEEE exceptions".

Rebuild and run (with a debug build). I get an error and break to the debugger on the line

croe    =  sqrt( GAMM1 * (Hroe - 0.5d0 * Uroe**2) )

with
GAMM1 = 0.399999976158142
Hroe = 172222.047861949
Uroe = -1269.36765662690

This will be in NTIME=54, as the last output to console is

NTIME=           53 TIME=   5.299999986618786E-005

Note: If you move the line “End Program Euler1d” to the end of the file then the remaining subroutines will be in contained in the main program. This will eliminate warnings about missing interfaces.

Edit: on linux the option is -fpe0 fpe

1 Like