* A. Salih, Dept. of Mechanical Engg., NIT - Trichy, India. *
*************************************************************************
* THIS PROGRAM SOLVES ONE-DIMENSIONAL LINEAR WAVE EQUATION WITH *
* PERIODIC AND NONPERIODIC BOUNDARY CONDITIONS USING *
* FINITE DIFFERENCE METHOD *
* -Schemes: 1. Euler FTFS *
* 2. Euler FTCS *
* 3. Upwind Scheme *
* 4. Lax-Friedrichs Scheme *
* 5. Lax-Wendroff Scheme *
* 6. Lax-Wendroff Scheme (multi-step) *
* 7. MacCormack Method *
* 8. Beam-Warming scheme *
* 9. Fromm scheme *
* 10. ENO-2 Scheme *
*************************************************************************
*
c Note: For the initial conditions RAMP FUNCTION and STEP FUNCTION
c non-periodic boundary conditions may be used. For other types of
c initial conditions periodic boundary conditions may be used.
PROGRAM LINEAR_WAVE
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
include 'scalars/parameters.inc'
open(unit=11,file='input.dat',status='unknown')
open(unit=21,file='output.dat',status='unknown')
open(unit=31,file='outputs/x.dat',status='unknown')
open(unit=32,file='outputs/u.dat',status='unknown')
open(unit=33,file='outputs/time.dat',status='unknown')
open(unit=34,file='outputs/time2.dat',status='unknown')
open(unit=35,file='outputs/u_trans.dat',status='unknown')
c
pi = 4*datan(1d0)
c
CALL READ_IN ! read-in data
close (unit=11)
CALL GRID ! setting up the grid points in the domain
CALL INITIAL_COND ! setting up initial condition
CALL SOLVE ! setting up system of equations
CALL PRINTOUT ! printing out the computational results
CALL WRITE_OUT ! write-out the basic data used for computation
c
stop
end
*
***SUBROUTINE: READING IN THE BASIC DATA *******************************
*
subroutine READ_IN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
c
read(11,*) alx ! actual length of domain
read(11,*) bc_type ! type of BC (1-Periodic, 2-Non-Periodic)
read(11,*) ic_type ! initial U profile (1-RAMP FUNCTION,
! 2-STEP DISCONTINUITY, 3-SQUAREWAVE,
! 4-SINEWAVE)
read(11,*) a ! value of advection velocity
read(11,*) Co ! Courant number for calculating time step
read(11,*) tfinal ! time at which the solution is desired
read(11,*) scheme ! 1-FTFS, 2-FTCS, 3-Upwind, 4-Lax-Friedrichs,
! 5-Lax-Wendroff, 6-Lax-Wendroff(multi-step),
! 7-MacCormack, 8-Beam-Warming, 9-Fromm,
! 10-ENO-2
read(11,*) print_freq !
c
c calculation of time-step dt based on the stability condition
dx = alx /(m-1)
dt = Co*dx/a
maxntimestp = tfinal /dt
dt = tfinal /maxntimestp ! revised time-step
Co = a*dt/dx ! revised Courant number
if (tfinal < dt) then
print*, 'Warning: final time is less than time-step!'
stop
endif
c
return
end
*
***SUBROUTINE: WRITING OUT THE BASIC DATA USED FOR COMPUTATION**********
*
subroutine WRITE_OUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
c
write(21,*) 'Linear, 1-dimensional wave equation with periodic bou
$ndary conditions:'
c
if (ic_type == 1) then
write(21,*) 'The initial profile is ramp function'
elseif (ic_type == 2) then
write(21,*) 'The initial profile is step function'
elseif (ic_type == 3) then
write(21,*) 'The initial profile is square wave'
elseif (ic_type == 4) then
write(21,*) 'The initial profile is sine wave'
endif
c
if (scheme == 1) then
write(21,*) 'FTFS Scheme'
elseif (scheme == 2) then
write(21,*) 'FTCS Scheme'
elseif (scheme == 3) then
write(21,*) 'Upwind Scheme'
elseif (scheme == 4) then
write(21,*) 'Lax-Friedrichs Scheme'
elseif (scheme == 5) then
write(21,*) 'Lax-Wendroff Scheme'
elseif (scheme == 6) then
write(21,*) 'Lax-Wendroff Scheme (multi-step)'
elseif (scheme == 7) then
write(21,*) 'MacCormack Method'
elseif (scheme == 8) then
write(21,*) 'Beam-Warming Scheme'
elseif (scheme == 9) then
write(21,*) 'Fromm Scheme'
elseif (scheme == 10) then
write(21,*) 'Essentially NonOscillatory Scheme'
endif
c
write(21,*) 'Advection velocity, a =', a
write(21,*) 'Courant number, Co =', Co
write(21,*) 'Length of the domain', alx
write(21,*) 'Number of grids (including boundaries):', m
write(21,*) 'dx =', dx
write(21,*) 'dt =', dt
write(21,*) 'No. of time step = ', ntimestp
write(21,*) 'Final time = ', time
c
return
end
*
****SUBROUTINE: GRID****************************************************
*
SUBROUTINE GRID
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/grid.inc'
c
c m = no. of grid points (including boundary nodes) in x-direction
dx = alx /(m-1)
x(1) = 0
do i = 2, m
x(i) = x(i-1) + dx
enddo
c
return
end
*
***SUBROUTINE: INITIAL CONDITIONS FOR U ********************************
*
SUBROUTINE INITIAL_COND
implicit doubleprecision (a-h,o-z)
include 'scalars/integers.inc'
c
c defining the initial conditions
if (ic_type == 1) then
CALL RAMPFN
elseif (ic_type == 2) then
CALL STEPFN
elseif (ic_type == 3) then
CALL SQUAREWAVE
elseif (ic_type == 4) then
CALL SINEWAVE
endif
c
return
end
*
****SUBROUTINE: INITIAL CONDITION - RAMP FUNCTION
*
SUBROUTINE RAMPFN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/u.inc'
c
vdisp = 0.1d0*alx
uleft = vdisp + 0.1d0*alx
uright = vdisp
xleft = 0.3d0*alx
xright = 0.4d0*alx
slope = (uleft-uright) /(xleft-xright)
do i = 1, m
if (x(i) < xleft) then
U(i) = uleft
elseif (x(i) < xright) then
U(i) = uleft + slope*(x(i) - xleft)
else
U(i) = uright
endif
enddo
c
return
end
*
****SUBROUTINE: INITIAL CONDITION - STEP FUNCTION
*
SUBROUTINE STEPFN
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/u.inc'
c
vdisp = 0.1d0*alx
uleft = vdisp + 0.1d0*alx
uright = vdisp
xstep = 0.3d0*alx
do i = 1, m
if (x(i) <= xstep) then
U(i) = uleft
else
U(i) = uright
endif
enddo
c
return
end
*
****SUBROUTINE: INITIAL CONDITION - SQUARE WAVE
*
SUBROUTINE SQUAREWAVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/u.inc'
c
amp = 1d0
vdisp = 5d-1
mleft = m/2+1 - m/10
mright = m/2 + m/10
do i = 1, m
if ( (i <= mleft).or.(i >= mright) ) then
U(i) = vdisp
else
U(i) = vdisp + amp
endif
enddo
c
return
end
*
****SUBROUTINE: INITIAL CONDITION - SINE WAVE
*
SUBROUTINE SINEWAVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/u.inc'
c
amp = 1d0
vdisp = 5d-1
mleft = m/2+1 - m/10
mright = m/2 + m/10
denom = mright-mleft
do i = 1, m
if ( (i <= mleft).or.(i >= mright) ) then
U(i) = vdisp
else
U(i) = vdisp + amp*dsin(pi*(i-mleft)/denom)
endif
enddo
c
return
end
*
***SUBROUTINE: NUMERICAL SOLUTION USING DIFFERENT SCHEMES***************
*
subroutine SOLVE
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'scalars/parameters.inc'
include 'arrays/grid.inc'
include 'arrays/u.inc'
doubleprecision UN(-1:m1+1), UPRED(-1:m1+1)
doubleprecision minmod_m, minmod_p
c
c storing the values of u at nth time level
UN = U
c
ntimestp = 0
maxdif = 1d9
time = 0d0
write(33,101) time ! continuous time
write(34,101) time ! intermittent time
monitor = m/4
write(35,*) U(monitor)
c
10 if (ntimestp < maxntimestp) then
if (bc_type == 1) then
c periodic boundary conditions
c west
do i = 0,-1,-1
UN(i) = UN(m+i-1)
enddo
c east
do i = 1,2
UN(m+i) = UN(i+1)
enddo
elseif (bc_type == 2) then
c nonperiodic boundary conditions
c west
do i = 0,-1,-1
UN(i) = UN(1)
enddo
c east
do i = 1,2
UN(m+i) = UN(m)
enddo
endif
c
if (scheme == 1) then ! FTFS
do i = 1,m
U(i) = UN(i) - Co*(UN(i+1) - UN(i))
enddo
c
elseif (scheme == 2) then ! FTCS
do i = 1,m
U(i) = UN(i) - Co*(UN(i+1) - UN(i-1))/2d0
enddo
c
elseif (scheme == 3) then ! Upwind
do i = 1,m
dUdx_m = (UN(i) - UN(i-1)) /dx
dUdx_p = (UN(i+1) - UN(i)) /dx
U(i) = UN(i) - (max(a,0d0)*dUdx_m + min(a,0d0)*dUdx_p)*dt
enddo
c
elseif (scheme == 4) then ! Lax-Friedrichs
do i = 1,m
U(i) = (UN(i+1) + UN(i-1))/2d0 -
$ Co*(UN(i+1) - UN(i-1))/2d0
enddo
c
elseif (scheme == 5) then ! Lax-Wendroff
do i = 1,m
U(i) = UN(i) - Co*(UN(i+1) - UN(i-1))/2d0 +
$ Co**2*(UN(i+1) - 2*UN(i) + UN(i-1))/2d0
enddo
c
elseif (scheme == 6) then ! Lax-Wendroff (multi-step)
c predictor step
do i = 1,m
UPRED(i) = (UN(i+1) + UN(i))/2 - Co*(UN(i+1) - UN(i))/2
enddo
if (bc_type == 1) then
do i = 0,-1,-1
UPRED(i) = UPRED(m+i-1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(i+1)
enddo
elseif (bc_type == 2) then
do i = 0,-1,-1
UPRED(i) = UPRED(1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(m)
enddo
endif
c corrector step
do i = 1,m
U(i) = UN(i) - Co*(UPRED(i) - UPRED(i-1))
enddo
c
elseif (scheme == 7) then ! MacCormack
c predictor step
do i = 1,m
UPRED(i) = UN(i) - Co*(UN(i+1) - UN(i))
enddo
if (bc_type == 1) then
do i = 0,-1,-1
UPRED(i) = UPRED(m+i-1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(i+1)
enddo
elseif (bc_type == 2) then
do i = 0,-1,-1
UPRED(i) = UPRED(1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(m)
enddo
endif
c corrector step
do i = 1,m
U(i) = (UN(i) + UPRED(i))/2d0 -
$ Co*(UPRED(i) - UPRED(i-1))/2d0
enddo
c
elseif (scheme == 8) then ! Beam-Warming (2nd order Upwind)
c predictor step ! for a > 0 only
do i = 1,m
UPRED(i) = UN(i) - Co*(UN(i) - UN(i-1))
enddo
if (bc_type == 1) then
do i = 0,-1,-1
UPRED(i) = UPRED(m+i-1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(i+1)
enddo
elseif (bc_type == 2) then
do i = 0,-1,-1
UPRED(i) = UPRED(1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(m)
enddo
endif
c corrector step
do i = 1,m
U(i) = (UN(i) + UPRED(i))/2d0 -
$ Co*(UPRED(i) - UPRED(i-1))/2d0 -
$ Co*(UN(i) - 2*UN(i-1) + UN(i-2))/2d0
enddo
c
elseif (scheme == 9) then ! Fromm (for a > 0 only)
c predictor step
do i = 1,m
UPRED(i) = UN(i) - Co*(UN(i+1) - UN(i-1))/2d0
enddo
if (bc_type == 1) then
do i = 0,-1,-1
UPRED(i) = UPRED(m+i-1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(i+1)
enddo
elseif (bc_type == 2) then
do i = 0,-1,-1
UPRED(i) = UPRED(1)
enddo
do i = 1,2
UPRED(m+i) = UPRED(m)
enddo
endif
c corrector step
do i = 1,m
U(i) = (UN(i) + UPRED(i))/2d0 -
$ Co*(UPRED(i) - UPRED(i-1))/2d0 -
$ Co*(UN(i) - 2*UN(i-1) + UN(i-2))/4d0
enddo
c
elseif (scheme == 10) then ! ENO-2
c second-order Essentially NonOscillatory scheme
do i = 1,m
dUdxxm = (UN(i) - 2*UN(i-1) + UN(i-2)) /(dx**2)
dUdxx = (UN(i+1) - 2*UN(i) + UN(i-1)) /(dx**2)
dUdxxp = (UN(i+2) - 2*UN(i+1) + UN(i)) /(dx**2)
c
if (dUdxx *dUdxxm > 0) then
minmod_m = dUdxx/(dabs(dUdxx)+1d-30) *
$ min(dabs(dUdxx), dabs(dUdxxm))
else
minmod_m = 0
endif
c
if (dUdxx *dUdxxp > 0) then
minmod_p = dUdxx/(dabs(dUdxx)+1d-30) *
$ min(dabs(dUdxx), dabs(dUdxxp))
else
minmod_p = 0
endif
c
dUdx_m = (UN(i) - UN(i-1))/dx + (dx/2)*minmod_m
dUdx_p = (UN(i+1) - UN(i))/dx - (dx/2)*minmod_p
U(i) = UN(i) - (max(a,0d0)*dUdx_m + min(a,0d0)*dUdx_p)*dt
enddo
endif
c
time = time + dt
ntimestp = ntimestp + 1
c
c updating u values at fictitious nodes
if (bc_type == 1) then
do i = 0,-1,-1
U(i) = U(m+i-1)
enddo
do i = 1,2
U(m+i) = U(i+1)
enddo
elseif (bc_type == 2) then
do i = 0,-1,-1
U(i) = U(1)
enddo
do i = 1,2
U(m+i) = U(m)
enddo
endif
c
UN = U
c
write(33,101) time ! continuous time
write(6,102) ntimestp
c
c printing out transient u at selected point of domain
j = jcounter + 1
jcounter = mod(j,print_freq)
if (jcounter == 0) then
write(34,101) (time + dt) ! intermittent time
write(35,*) U(monitor)
endif
go to 10
endif
c
101 format(e14.7)
102 format(i7,5x,e14.7)
*
return
end
*
***SUBROUTINE: PRINTING-OUT THE OUTPUT DATA*****************************
*
SUBROUTINE PRINTOUT
implicit doubleprecision (a-h,o-z)
include 'array_dimension.inc'
include 'scalars/integers.inc'
include 'scalars/reals.inc'
include 'arrays/u.inc'
include 'arrays/grid.inc'
c
write(31,101) (x(i), i=1,m)
write(32,102) (U(i), i=1,m)
101 format(e13.6)
102 format(e13.6)
c
return
end
|
