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')
if (scheme==3) then
end if

close(8)

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