The Algorithm


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


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


 

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