The Algorithm

• Interfacing

The program is designed to read the data from a file called « don » and to write the results into « res.txt ».
 

• The data file DON

The Don files permits the setting of all user defined parameters:

1                domain length
300           number of elements x
0.001        time step dt
1                final time tf
5                number of time step between saving passauv
0.01          viscosity nu
3                scheme
0                initialization
1                boundary type
0.               lambda
0.               mu
1.               flux or u value at 0
-1.              flux or u value at l
0.               q
0.               delta

• The length of the domain (usually 1m)
• The number of elements (number of points-1)
• The time step
• The final time
• The saving frequency (100 saves U to the res.txt file every 100 time steps)
• The scheme type: 1=explicit: this scheme is present for development purposes as it is easy to code. It is nit suitable for most computation. 2=implicit Crank Nicholson (CN) finite difference a 3 point scheme that gives good results in most cases, 3=generalized CN a 4 points scheme, the possibility of setting the weight of each point using the parameters d and q gives the opportunity to test a large array of schemes including the already coded CN finite difference one.
• The type of initial conditions: 3 are available: 0=sinus; 1=shock, 2=K/(1+x). The initial condition setting is related to the value of flxg and flxd.
• The type of boundary conditions: 0=Dirichlet (fixed u value at both extremities, value defines by flxg ar x=0 and flxd at x=1); 1=Neumann (fixed fluxes defined by flxd and flxg).
• The value of l and m for the stability study.
• The boundaries values flxd and flxg.
• The value of d and q for the generalized CN scheme.

• The FORTRAN 90 routine

 
• The use of case structures enables quick adding of new schemes or initialization.
• For implicit schemes' Matrix inversion the Thomas algorithm is used

C     Last change:  SPI  23 Nov 100    2:34 pm

! Start declaration ***********************************************************************
      implicit none
        Real  u1     !velocity vector
        Real  u2     !velocity vector
        Real  dt                        !time step
        Real  dx                        !spatial step
        REAL  x,l                       !position vector, total length
        Real  tf                        !end time
        Real  t                         !current time
        Real  nu                        !viscosity
        Real  r,k                       !r=nu*dt/dx**2, k=0.5*dt/dx
        Real  a,b,c,d,c2,d2             !matrix coefficents for implicit scheme
 Real  e,f,f3,c3,d3,a2,b2        !
 Real  q,delta                   !coefficients for implicit scheme
        REAL  lambda, mu                !coefficients for u and u**3 terms
        REAL flxd,flxg                  !boundaries fluxes or velocity value

        Integer n                       !number of spatial discretization points
        Integer boun,init               !boundary type index, initialization type index
        Integer scheme                  !scheme index
        Integer i,j                     !loop indexes
        INTEGER it                      !current iteration
        Integer passauv                 !saving to file every passauv iterations

        allocatable u1(:),u2(:),x(:),a(:),b(:),c(:),d(:),c2(:),d2(:),
     & e(:),a2(:),b2(:),d3(:),c3(:),f(:),f3(:)

!End declaration *************************************************************************
 

! Reading Data from don file **************************************************************

        open(8,file='don',status='old')
                read (8,*) l
         read(8,*) n
                read(8,*) dt
                read(8,*) tf
                read(8,*) passauv
                read(8,*) nu
                read(8,*) scheme
                READ(8,*) init
                read(8,*) boun
                read (8,*) lambda
                read (8,*) mu
                read (8,*) flxg
                read (8,*) flxd
               if (scheme==3) then
                read(8,*) q
                read(8,*) delta
               end if

                close(8)
! End reading data ***********************************************************************
 

                allocate (u1(n+1),u2(n+1),x(n+1),a(n+1),b(n+1),c(n+1),
     &             d(n+1),c2(n+1),d2(n+1), e(n+1),c3(n+1),d3(n+1),
     &             a2(n+1),b2(n+1),f(n+1),f3(n+1))

! Initialization **************************************************************************

                dx=1./real(n)
         r=nu*dt/(dx*dx)
                k=0.5*dt/dx
 

               open(7,file='res.txt')

        do i=1,n+1
                x(i)=(i-1)*dx*l
                a(i)=0.
                b(i)=0.
                c(i)=0.
                d(i)=0.
                e(i)=0.
                f(i)=0.
                end do
 
 

        t=0.
        it=0

! Setting the initialization type ******************************************************

! init=0 sinusoidal
! init=1 shock
! init=2 1/(1+x)

                select case (init)

                case (0)  ! init=0 sinusoidal

                  do i=1,n+1
               u1(i)=sin(2*3.1415927*(i-1)/n)
                  end do

                case (1) ! init=1 shock

                  do i=1,n+1
                  IF (X(i)>l/2.) THEN
                  u1(i)=flxd
                  else
                  u1(i)=flxg
              END IF
                  end do

                case (2)  ! init=2 1/(1+x)

                  do i=1,n+1
               u1(i)=-2*nu/(x(i)+1)
                  end do

                end select
! End Initialization *******************************************************************
 
 

!****************************************************************************************
! TIME LOOP ******************************************************************************

                do

         if (t>tf) exit !exit when end time is reached

                t=t+dt
                it=it+1
 
 
 

! Start scheme computation **************************************************************

!scheme=1 explicit
!scheme=2 Crank Nicholson
!scheme=3 Generalized Crank Nicholson
 

  select case(scheme)
 

!**************************************************************************************
  case (1) !scheme=1 explicit

                do i=2,n
                u2(i)=r*u1(i+1)+(1-2*r)*u1(i)+r*u1(i-1)
     &            -k*u1(i)*(u1(i+1)-u1(i-1))
                        end do

                do i=1,n
                        u1(i)=u2(i)
                        end do

!***************************************************************************************

                case(2)   !scheme=2 Crank Nicholson

! Matrix coefficients

                do j=2,n+1
                a(j)=-0.25*dt*u1(j-1)/dx-0.5*r
                end do

                do j=1,n+1
                b(j)=1+r
                end do

                do j=1,n
                c(j)=0.25*dt*u1(j+1)/dx-0.5*r
                end do

                do j=2,n
      d(j)=0.5*r*u1(j-1)+(1-r)*u1(j)+0.5*r*u1(j+1)
     &+lambda*u1(j)+mu*u1(j)**3
                 end do

!  Boundaries

               if (boun==0) then !Dirichlet boundaries
        b(1)=1.
        c(1)=0.
        d(1)=flxg
        b(n+1)=1.
        a(n+1)=0.
        d(n+1)=flxd

                else             !Neumann boundaries

        d(1)=(1-r)*u1(1)+2*r*flxg+r*u1(2)+lambda*u1(1)+mu*u1(1)**3
      d(n+1)=(1-r)*u1(n+1)+r*u1(n)+2*r*flxd+lambda*u1(n+1)+mu*u1(n+1)**3
                endif

! Thomas Algorithm

                c2(1)=c(1)/b(1)
                d2(1)=d(1)/b(1)

                do i=2,n+1
                c2(i)=c(i)/(b(i)-a(i)*c2(i-1))
                d2(i)=(d(i)-a(i)*d2(i-1))/(b(i)-a(i)*c2(i-1))
                end do
                u1(n+1)=d2(n+1)
                do i=1,n
                u1(n+1-i)=d2(n+1-i)-u1(n-i+2)*c2(n+1-i)
                end do
 
 

!scheme=3 Generalized Crank Nicholson ****************************************************

  case(3)

! Matrix coefficients

                do j=3,n-1
                if (u1(j)>=0) then
                e(j)=q/6.*dt/dx*u1(j-2)
                else
                e(j)=0.
                END if
  end do

                e(n)=0.
                e(n+1)=0.

                if (u1(2)>=0) then
                a(2)=-(0.25)*dt*u1(1)/dx-0.5*r+delta
                else
                a(2)=-(0.25)*dt*u1(1)/dx-0.5*r+delta
                END if

  do j=3,n-1
                if (u1(j)>=0) then
                a(j)=-(0.25+0.5*q)*dt*u1(j-1)/dx-0.5*r+delta
                else
                a(j)=-(0.25-q/6.)*dt*u1(j-1)/dx-0.5*r+delta
                END if
                end do

                if (u1(n)>=0) then
                a(n)=-(0.25)*dt*u1(n-1)/dx-0.5*r+delta
                else
                a(n)=-(0.25)*dt*u1(n-1)/dx-0.5*r+delta
                END if

                if (u1(n+1)>=0) then
                a(n+1)=-(0.25)*dt*u1(n)/dx-r+2*delta
                else
                a(n+1)=-(0.25)*dt*u1(n)/dx-r+2*delta
                END if
 
 
 

                DO j=1,2
                if (u1(2)>=0) then
                b(j)=1+r-2*delta
                else
                b(j)=1+r-2*delta
                END if
                END do

                do j=3,n+1
                if (u1(j)>=0) then
                b(j)=1+r-2*delta+0.5*q*dt*u1(j)/dx
                else
                b(j)=1+r-2*delta-0.5*q*dt*u1(j)/dx
                END if
                end do

                DO j=n,n+1
                if (u1(2)>=0) then
                b(j)=1+r-2*delta
                else
                b(j)=1+r-2*delta
                END if
                END do

                if (u1(1)>=0) then
                c(1)=(0.25)*dt*u1(2)/dx-r+2*delta
                else
                c(1)=(0.25)*dt*u1(2)/dx-r+2*delta
                END if

                if (u1(2)>=0) then
                c(2)=(0.25)*dt*u1(3)/dx-0.5*r+delta
                else
                c(2)=(0.25)*dt*u1(3)/dx-0.5*r+delta
                END if

                do j=3,n-1
                if (u1(j)>=0) then
                c(j)=(0.25-q/6.)*dt*u1(j+1)/dx-0.5*r+delta
                else
                c(j)=(0.25+0.5*q)*dt*u1(j+1)/dx-0.5*r+delta
                END if
                end do

                if (u1(n)>=0) then
                c(n)=(0.25)*dt*u1(n+1)/dx-0.5*r+delta
                else
                c(n)=(0.25)*dt*u1(n+1)/dx-0.5*r+delta
                END if

                f(1)=0.
                f(2)=0.

                do j=3,n-1
                if (u1(j)>=0) then
                f(j)=0.
                else
                f(j)=-q/6.*dt/dx*u1(j+2)
                END if
                end do
 

                do j=2,n
               d(j)=(0.5*r+delta)*u1(j-1)+(1-r-2*delta)*u1(j)
     & +(delta+0.5*r)*u1(j+1)+lambda*u1(j)+mu*u1(j)**3
                end do

! Boundaries

         if (boun==0) then  !Dirichlet boundaries
         b(1)=1.
         c(1)=0.
         d(1)=flxg
         f(1)=0.

         b(n+1)=1.
         a(n+1)=0.
         d(n+1)=flxd
         e(n+1)=0.

                else          !Neumann boundaries
       d(1)=(1-r-2*delta)*u1(1)
     & +(delta+0.5*r)*(2*u1(2)+2*flxg)+lambda*u1(1)+mu*u1(1)**3
       d(n+1)=(1-r-2*delta)*u1(n+1)
     & +(delta+0.5*r)*(2*u1(n)+2*flxd)+lambda*u1(n+1)+mu*u1(n+1)**3
                endif
 

! Generalized Thomas Algoritm
 

  b2(1)=b(1)
  c2(1)=c(1)
                d2(1)=d(1)

                a2(2)=a(2)
  b2(2)=b(2)
  c2(2)=c(2)
                d2(2)=d(2)
 

                do i=3,n
                a2(i)=a(i)-e(i)*b2(i-1)/a2(i-1)
  b2(i)=b(i)-e(i)*c2(i-1)/a2(i-1)
  c2(i)=c(i)-e(i)*f(i-1)/a2(i-1)
                d2(i)=d(i)-e(i)*d2(i-1)/a2(i-1)
                end do

                a2(n+1)=a(n+1)-e(n+1)*b2(n)/a2(n)
  b2(n+1)=b(n+1)-e(n+1)*c2(n)/a2(n)
  d2(n+1)=d(n+1)-e(n+1)*d2(n)/a2(n)

          c3(1)=c2(1)/b2(1)
  d3(1)=d2(1)/b2(1)
                f3(1)=f(1)/b2(1)

  do i=2,n-1
  c3(i)=(c2(i)-a2(i)*f3(i-1))/(b2(i)-a2(i)*c3(i-1))
  d3(i)=(d2(i)-a2(i)*d3(i-1))/(b2(i)-a2(i)*c3(i-1))
                f3(i)=f(i)/(b2(i)-a2(i)*c3(i-1))
  end do

  c3(n)=(c2(n)-a2(n)*f3(n-1))/(b2(n)-a2(n)*c3(n-1))
  d3(n)=(d2(n)-a2(n)*d3(n-1))/(b2(n)-a2(n)*c3(n-1))

              d3(n+1)=(d2(n+1)-a2(n+1)*d3(n))/(b2(n+1)-a2(n+1)*c3(n))

  u1(n+1)=d3(n+1)
                u1(n)=d3(n)-u1(n+1)*c3(n)

                do i=2,n
            u1(n+1-i)=d3(n+1-i)-u1(n-i+2)*c3(n+1-i)-u1(n+3-i)*f3(n-i+1)
                end do

                end select
!End of scheme computation ****************************************************************

!Saving Data to file res.txt **************************************************************
                if (MOD(it,passauv)==0.) then

              write(*,*) 't=',t,'s',it
              write(7,*) 't=',t,'s',it

              do i=1,n+1
              write(7,*) x(i),u1(i)
              end do

               END if
!End Saving Data***************************************************************************

        end do
! END TIME LOOP ***************************************************************************
*******************************************************************************************

!Saving last computation to file *********************************************************
             write(7,*) 't=',t,'s'

             do i=1,n+1
             write(7,*) x(i),u1(i)
             END do
 
 

!End saving *******************************************************************************
 

        close (7)

        deallocate (u1,u2,x,a,b,c,d,c2,d2)

       end program