!------------------------------------
! Michael McMahon
! Dept. of Mechanical, Aerospace, and Manufacturing Eng.
! Computational Fluid Dynamics
! Syracuse Univ.     1997
! mmcmahon@npac.syr.edu
!------------------------------------
! This program simulates quasi 1-D inviscid flow through a 
!divering nozzle.  Depending on the back pressure given,
!the code will output mach number data that shows either
! a smooth progression or a shock wave in the middle of the
!nozzle.  The equations governing flow have been cast in
!conservative form so that the shock wave may be captured,
!and the mesh is clustered around x=0.5 to capture the
!sudden discontinuity. To get the data that shows the 
!shock wave, use a back pressure of 0.88; to get the
!unshocked data, use a back pressure of 0.5 or so.
!The program gives the user a choice of second order scheme (Runge 
!kutta) or first order scheme (Lax).  The Runge Kutta takes
!more iterations, especially in the shock wave case.
!This program makes use of the FORALL statement frequently;
!in fact there are only two DO loops in the entire code.
!This construct is more succint, and because it guarantees
!no dependencies make conversion to parallel implementation
!quite easy.
!
!Summary of data files:    mach.dat--> axial location, mach no.
!                          laxiter.dat --> iteration, RMS of res for Lax
!                          jamiter.dat --> "                   " for R-K
!                          mfr.dat --> axial location, mass flow rate
!                          enth.dat --> axial location, enthalpy
!                          tropy.dat --> axial location, entropy
!  This program must be compiled in fixed format

      program nozzle
      Implicit none
      Integer, Parameter :: nx=140
      Integer, Parameter :: dim=3
      Real V(dim,0:nx)
      Real xloc(0:nx)
      Integer scheme,i

      Write(6,*) 'The following schemes are available:'
      Write(6,*) 'Lax                   scheme=1'
      Write(6,*) 'Jameson 4th order RK  scheme=2'
      
      Write(6,*) 'Enter the scheme you wish to run:'
      Read(5,*) scheme

      Call Mesh(xloc,nx)
     
      Call Init_Cons(V,xloc,dim,nx)
      If (scheme .eq. 1 ) then
         Call Solver_Lax(V,xloc,dim,nx,scheme)
		Else if(scheme .eq. 2) then
		Call Solver_James(V,xloc,dim,nx,scheme)
      End if

      End Program

      Subroutine Mesh(xloc,nx)
      Implicit none
      Integer nx
      Real xloc(0:nx)
      Real zeta,xmin,xmax
      Real pi
      Integer i

      zeta=0.0
      xmin=-0.2
      xmax=1.2
      pi=4.0*atan(1.0)

! Create a mesh clustered at x=0.5, near the shock wave

      Do i=0,nx/2
         xloc(i)=(1.0-2.0*xmin)/2.0*sin(pi*zeta) + xmin
         zeta=zeta+ 1.0/nx
      End Do

      Do i=nx/2+1,nx
         xloc(i) = xmax-(2.0*xmax-1)/2.0*sin(pi*zeta)
         zeta=zeta+1.0/nx
      End do

      End Subroutine

      Subroutine Init_Cons(V,xloc,dim,nx)
      Implicit none
      Integer dim,nx
      Real V(dim,0:nx)
      Real xloc(0:nx)
      Integer j
!Set initial conditions

      Forall(j=0:nx) V(1,j)=1.0
      Forall(j=0:nx) V(2,j)=1.0
      Forall(j=0:nx) V(3,j)=1.976
      End Subroutine



      Subroutine Solver_Lax(V,xloc,dim,nx,scheme)
      Implicit none	
      Integer dim,nx,scheme
      Real V(dim,0:nx)
      Real h(dim,0:nx)
      Real Res(dim,0:nx),ef(dim,0:nx),eb(dim,0:nx)
      Real df(dim,0:nx),db(dim,0:nx)
      Real dt(1:nx-1)
      Real s(0:nx)
      Real xloc(0:nx)
      Real p(0:nx),test(0:nx)
      Real rms
      Real cfl,iter,pback
      Integer i,j

      Write(6,*) 'Please enter P_back (0.1 < p  < 0.99)'
      read(5,*) pback
      Write(6,*) 'Please enter CFL number (0.1 < x < 0.9)'
      read(5,*) cfl
      Forall(i=1:nx-1) p(i)=0.0
      p(0)=0.5903
      p(nx)=0.0
	  Forall(i=1:dim,j=0:nx) h(i,j) =0.0
		test=1.0
		res=0.0	
       Call shape_noz(s,xloc,nx)
		iter=0.0
		open(unit=12,file='laxiter.dat',status='unknown')

!______________Start looping___________
!General form of the loop:
!        -Compute time step at each axial location
!        -Compute flux vectors for each "volume"
!        -Compute dissipation vectors
!        -Compute residual
!        -Update solution
!	 -Output iteration & RMS of residual

       Do While(maxval(abs(test)) .gt. 3.8e-2) 
         Call Time_step(V,dt,p,xloc,cfl,nx,dim)
         Call Flux(V,ef,eb,p,xloc,dim,nx)
         Call Dissp(V,df,db,xloc,dt,p,dim,nx,scheme)
         Forall(i=1:dim,j=1:nx-1) res(i,j)= -1.0*((ef(i,j)
     $   -df(i,j))
     $ *0.5*(s(j+1)+s(j)) - (eb(i,j)-db(i,j))
     $ *0.5*(s(j)+s(j-1)))
     $ /(s(j)*(xloc(j+1)-xloc(j-1))*0.5) 
     $    + h(i,j)
        Forall(i=1:dim,j=1:nx-1) V(i,j)=V(i,j)+dt(j)*Res(i,j)
       
        Forall(i=1:nx-1) p(i)=0.4*(V(3,i)-0.5*V(2,i)*V(2,i)/V(1,i))
        Forall(i=1:nx-1) h(2,i)= (p(i)/s(i))*(s(i+1)-s(i-1))
     $   /(xloc(i+1)-xloc(i-1)) 
         Call Boun_cond(V,p,pback,dim,nx)
		iter=iter+1
		ForalL(i=0:nx)   	test(i)=res(1,i)
		write(6,*) maxval(abs(test)) ,iter, p(nx)
		if (mod(int(iter),10) . eq. 0.0) then
		rms=sqrt(sum(test*test)/141)
		write(12,*) iter,rms
		end if
         End do
         Call Output(V,xloc,p,s,dim,nx)
         End subroutine

     	Subroutine Solver_James(V,xloc,dim,nx,scheme)
        Implicit none
        Integer dim,nx,scheme
        Real V(dim,0:nx)
        Real h(dim,0:nx)
        Real Res(dim,0:nx),ef(dim,0:nx),eb(dim,0:nx)
        Real df(dim,0:nx),db(dim,0:nx),Q(dim,0:nx)
        Real dt(1:nx-1)
        Real s(0:nx)
        Real xloc(0:nx)
        Real p(0:nx)
        Real test(0:nx)
        Real cfl,iter,pback,rms
        Integer i,j

        Write(6,*) 'Please enter P_back (0.1 < p < 0.9)'
        read(5,*) pback
        Write(6,*) 'Please enter CFL number (0.1 < x < 0.9)'
        read(5,*) cfl
        p(0)=0.5903

        Forall(i=1:dim,j=0:nx) h(i,j) =0.0
        Forall(i=1:nx) p(i)=0.1
		res=0.0
          test(5)=1.0
		df=0.0
		db=0.0
       Call Shape_noz(s,xloc,nx)
		Q=V
		open(unit=14,file='jamiter.dat',status='unknown')
		
!-----------Start looping-------------
!General form of the loop:
!        -Compute time step at each axial location
!        -Compute flux vectors for each "volume"
!        -Compute dissipation vectors
!        -Compute residual
!	 -Compute stage of Runge Kutta
!        -Repeat for all 4 stages
!	 -Update solution 
!	 -Output iteration, RMS of residual

	   Do While(maxval(abs(test)) .gt. 2.8e-3)
	   Call Time_step(V,dt,p,xloc,cfl,nx,dim)
        Call Flux(V,ef,eb,p,xloc,dim,nx)
        Call Dissp(V,df,db,xloc,dt,p,dim,nx,scheme)
        Forall(i=1:dim,j=1:nx-1) res(i,j)=-1.0*((ef(i,j)
     $  -df(i,j))
     $ *0.5*(s(j+1)+s(j)) - (eb(i,j)-db(i,j))
     $  *0.5*(s(j)+s(j-1)))
     $ /(s(j)*(xloc(j+1)-xloc(j-1))*0.5)
     $    + h(i,j)
        Forall(i=1:dim,j=1:nx-1) Q(i,j)=V(i,j)+0.25*dt(j)*Res(i,j)
        Forall(i=1:nx-1) p(i)=0.4*(Q(3,i)-0.5*Q(2,i)*Q(2,i)/Q(1,i))
        Forall(i=1:nx-1) h(2,i)= p(i)/s(i)*(s(i+1)-s(i-1))
     $   /(xloc(i+1)-xloc(i-1))
        Call Boun_cond(Q,p,pback,dim,nx)
		Call Time_step(Q,dt,p,xloc,cfl,nx,dim)
         Call Flux(Q,ef,eb,p,xloc,dim,nx)
        Forall(i=1:dim,j=1:nx-1) res(i,j)=-1.0*((ef(i,j)
     $ -df(i,j))
     $ *0.5*(s(j+1)+s(j)) - (eb(i,j)-db(i,j))
     $ *0.5*(s(j)+s(j-1)))
     $ /(s(j)*(xloc(j+1)-xloc(j-1))*0.5)
     $    + h(i,j)
        Forall(i=1:dim,j=1:nx-1) Q(i,j)=V(i,j)+1.0/3.0*dt(j)*Res(i,j)
        Forall(i=1:nx-1) p(i)=0.4*(Q(3,i)-0.5*Q(2,i)*Q(2,i)/Q(1,i))
        Forall(i=1:nx-1) h(2,i)= p(i)/s(i)*(s(i+1)-s(i-1))
     $   /(xloc(i+1)-xloc(i-1))
        Call Boun_cond(Q,p,pback,dim,nx)
		call Time_step(Q,dt,p,xloc,cfl,dim,nx)	
         Call Flux(Q,ef,eb,p,xloc,dim,nx)

         Forall(i=1:dim,j=1:nx-1) res(i,j)=-1.0*((ef(i,j)
     $  -df(i,j))
     $ *0.5*(s(j+1)+s(j)) - (eb(i,j)-db(i,j))
     $  *0.5*(s(j)+s(j-1)))
     $ /(s(j)*(xloc(j+1)-xloc(j-1))*0.5)
     $    + h(i,j)

       Forall(i=1:dim,j=1:nx-1) Q(i,j)=V(i,j) + 0.5*dt(j)*Res(i,j)
       Forall(i=1:nx-1) p(i)=0.4*(Q(3,i)-0.5*Q(2,i)*Q(2,i)/Q(1,i))
       Forall(i=1:nx-1) h(2,i)= p(i)/s(i)*(s(i+1)-s(i-1))
     $   /(xloc(i+1)-xloc(i-1))
       Call Boun_cond(Q,p,pback,dim,nx)
		Call Time_step(Q,dt,p,xloc,cfl,dim,nx)
         Call Flux(Q,ef,eb,p,xloc,dim,nx)
         Forall(i=1:dim,j=1:nx-1) res(i,j)=-1.0*((ef(i,j)-df(i,j))
     $ *0.5*(s(j+1)+s(j)) - (eb(i,j)-db(i,j))*0.5*(s(j)+s(j-1)))
     $ /(s(j)*(xloc(j+1)-xloc(j-1))*0.5)
     $    + h(i,j)

       Forall(i=1:dim,j=1:nx-1) V(i,j)=V(i,j) + (dt(j))*Res(i,j)
       Forall(i=1:nx-1) p(i)=0.4*(V(3,i)-0.5*V(2,i)*V(2,i)/V(1,i))
        Forall(i=1:nx-1) h(2,i)= p(i)/s(i)*(s(i+1)-s(i-1))
     $   /(xloc(i+1)-xloc(i-1))
       Call Boun_cond(V,p,pback,dim,nx)
		iter=iter+1

		 Forall(i=0:nx) test(i)=res(1,i)
             write(6,*) maxval(abs(test)),iter,p(nx)
		if (mod(int(iter),10) .eq. 0) then
		rms=sqrt(sum(test*test)/141)
		write(14,*) iter,rms
		end if
         End do	
		Call Output(V,xloc,p,s,dim,nx)
		End Subroutine


       Subroutine Flux(V,ef,eb,p,xloc,dim,nx)
       Implicit none
       Integer dim,nx
       Real V(dim,0:nx)
       Real ef(dim,0:nx),eb(dim,0:nx)
       Real xloc(0:nx)
       Real p(0:nx)
       Integer i

!Compute mass flux, momentum flux and enthalpy flux into volume
		ef=0.0
		eb=0.0
       Forall(i=1:nx-1) ef(1,i)=0.5*(V(2,i+1)+V(2,i))
       Forall(i=1:nx-1) ef(2,i)=0.5*(V(2,i+1)**2 / V(1,i+1)
     $ +p(i+1) + V(2,i)**2 / V(1,i) + p(i))
       Forall(i=1:nx-1) ef(3,i)=0.5*((V(3,i+1)+p(i+1))*V(2,i+1)
     $ /V(1,i+1)  + (V(3,i)+p(i))*V(2,i)/V(1,i))

  	    Forall(i=1:nx-1) eb(1,i)=0.5*(V(2,i-1)+V(2,i))
       Forall(i=1:nx-1) eb(2,i)=0.5*(V(2,i-1)**2 / V(1,i-1)
     $ +p(i-1) + V(2,i)**2 / V(1,i) + p(i))
       Forall(i=1:nx-1) eb(3,i)=0.5*((V(3,i-1)+p(i-1))*V(2,i-1)
     $  /V(1,i-1) + (V(3,i)+p(i))*V(2,i)/V(1,i)) 
       End Subroutine


       Subroutine Dissp(V,df,db,xloc,dt,p,dim,nx,scheme)
       Implicit none	
       Integer dim,nx,scheme
       Real V(dim,0:nx)
       Real df(dim,0:nx),db(dim,0:nx)
       Real xloc(0:nx),dt(1:nx-1)
       Real eps2(0:nx),eps4(0:nx)	
       Real p(0:nx)
       Real kp2,kp4
       Integer i,j

! Compute artificial dissipation vectors, necessary for stability
!of scheme

        df=0.0
        db=0.0
        eps2=0.0
        eps4=0.0
        kp2=1.0/4.0
        kp4=1.0/256.0
       If (scheme .eq. 1) then
       Forall(i=1:nx-1) df(1,i)=0.5*(xloc(i+1)-xloc(i))
     $  *(V(1,i+1)-V(1,i))/dt(i)
       Forall(i=1:nx-1) df(2,i)=0.5*(xloc(i+1)-xloc(i))
     $  *(V(2,i+1)-V(2,i))/dt(i)
       Forall(i=1:nx-1) df(3,i)=0.5*(xloc(i+1)-xloc(i))
     $  *(V(3,i+1)-V(3,i))/dt(i)
  
        Forall(i=1:nx-1) db(1,i)=0.5*(xloc(i)-xloc(i-1))
     $  *(V(1,i)-V(1,i-1))/dt(i)
       Forall(i=1:nx-1) db(2,i)=0.5*(xloc(i)-xloc(i-1))
     $  *(V(2,i)-V(2,i-1))/dt(i)
       Forall(i=1:nx-1) db(3,i)=0.5*(xloc(i)-xloc(i-1))
     $  *(V(3,i)-V(3,i-1))/dt(i)
		End if

		If (scheme .eq. 2) then
		Forall(i=2:nx-2) eps2(i)=kp2*max(abs((p(i+1)-2.0*p(i)
     $  +p(i-1))/(p(i+1)+2.0*p(i)+p(i-1))),abs((p(i+2)
     $  -2.0*p(i+1) +p(i))/(p(i+2)+2.0*p(i+1)+p(i))))	
		
		Forall(i=2:nx-2) eps4(i)=max(0.0,(kp4-eps2(i)))
        
		Forall(i=1:dim,j=2:nx-2) df(i,j)=eps2(j)*0.5
     $ *(abs(V(2,j+1)/V(1,j+1)+sqrt(1.4*p(j+1)/V(1,j+1)))
     $ +abs(V(2,j)/V(1,j)+sqrt(1.4*p(j)/V(1,j))))*(V(i,j+1)-V(i,j))
     $ -eps4(j)*0.5
     $ * abs(V(2,j+1)/V(1,j+1)+sqrt(1.4*p(j+1)/V(1,j+1))
     $ + V(2,j)/V(1,j)+sqrt(1.4*p(j)/V(1,j)))*(V(i,j+2)-3.0*V(i,j+1)
     $ + 3.0*V(i,j) -V(i,j-1))

	    Forall(i=2:nx-2) eps2(i)=kp2*max(abs((p(i)-2.0*p(i-1)
     $  +p(i-2))/(p(i)+2.0*p(i-1)+p(i-2))),abs((p(i+1)
     $  -2.0*p(i) +p(i-1))/(p(i+1)+2.0*p(i)+p(i-1))))

        Forall(i=2:nx-2) eps4(i)=max(0.0,(kp4-eps2(i)))

        Forall(i=1:dim,j=2:nx-2) db(i,j)=eps2(j)*0.5
     $ * (abs(V(2,j)/V(1,j)+sqrt(1.4*p(j)/V(1,j))
     $ + V(2,j-1)/V(1,j-1)+sqrt(1.4*p(j-1)/V(1,j-1))))
     $ * (V(i,j)-V(i,j-1))
     $ -eps4(j)*0.5
     $ *( abs(V(2,j)/V(1,j)+sqrt(1.4*p(j)/V(1,j))
     $ + V(2,j-1)/V(1,j-1)+sqrt(1.4*p(j-1)/V(1,j-1))))
     $ *(V(i,j+1)-3.0*V(i,j)+3.0*V(i,j-1) -V(i,j-2))
		End if

       End subroutine


       Subroutine Boun_cond(V,p,pback,dim,nx)
       Implicit none
       Integer dim,nx
       Real V(dim,0:nx)
       Real p(0:nx)
       Real pback
       Integer i,j

!If the mach number at the exit is supersonic, dont specify back pressure
!If the mach number at the exit is subsonic, specify back pressure as given

       V(1,nx)=V(1,nx-1)
       V(2,nx)=V(2,nx-1)
       V(3,nx)=V(3,nx-1)
		If ((V(2,nx-1)/V(1,nx-1))/(sqrt(1.4*p(nx-1)
     $   /V(1,nx-1))) .gt. 1.0 ) then
			p(nx)=p(nx-1)
		End if
		If ((V(2,nx-1)/V(1,nx-1))
     $   /(sqrt(1.4*p(nx-1)/V(1,nx-1))) .le. 1.0 ) then
			p(nx) = pback
		End if


       end subroutine

      Subroutine Time_step(V,dt,p,xloc,cfl,nx,dim)
      Implicit none	
      Integer nx,dim
      Real V(dim,nx)	
      Real dt(1:nx-1),p(0:nx),xloc(0:nx)
      Real cfl
      Integer i

      Do i=1,nx-1
         dt(i)=0.5*cfl/(V(2,i)/V(1,i)+sqrt(1.4*p(i)/V(1,i)))
     $   *(xloc(i+1)-xloc(i-1))
      End do 

      End subroutine

      Subroutine Shape_noz(s,xloc,nx)
	  Implicit none
      Integer nx
      Real s(0:nx),xloc(0:nx)
      Real sexit,pi
      Integer i

!Get shape of the divering nozzle (cos wave between 0< x < 1

      sexit=2.0
      pi=4.0*atan(1.0)
		open(unit=8,file='noz.dat',status='unknown')
      Do i=0,nx
         If (xloc(i) .le. 0.0) then
            s(i)=1.0
         Else if (xloc(i) .gt. 0.0 .and. xloc(i) .le. 1.0)then
            s(i) = 0.5*(1.0-sexit)*cos(pi*xloc(i))  
     $ + 0.5*(1.0+sexit)

         Else if (xloc(i) .gt. 1.0) then
            s(i)=sexit
         End if
		write(8,*) xloc(i),s(i)
       End do

       End subroutine

       Subroutine Output(V,xloc,p,s,dim,nx)
	   implicit none
       Integer dim,nx
       Real V(dim,0:nx)
       Real xloc(0:nx),p(0:nx),s(0:nx)
       Integer i

       Open(unit=8,file='mach.dat',status='unknown')
		Open(unit=9,file='mfr.dat',status='unknown')
	   Open(unit=10,file='enth.dat',status='unknown')
	   Open(unit=11,file='tropy.dat',status='unknown')
       DO i=0,nx
          Write(8,*) xloc(i),V(2,i)/V(1,i)/sqrt(1.4*p(i)/V(1,i))
		  Write(9,*) xloc(i),V(2,i)*s(i)
	      Write(10,*) xloc(i), v(3,i)/V(1,i) + p(i)/V(1,i)
          Write(11,*) xloc(i),(1.4/0.4)*log(p(i)/V(1,i) )
     $   -log(p(i)) - (1.4/0.4)*log(p(0)/V(1,0)) + log(p(0))
       End do

       End subroutine