View Single Post
edge333
edge333 is offline
#3
Feb24-12, 03:19 PM
P: 12
Thanks, I don't think that's it though. Here's a more simplified version of the code that doesn't involve solving for implicit methods that generate tridiagonal matrices. I think the problems are the same in terms of not getting the final spatial interval and time step. With this explicit scheme it might be easier to see what is wrong.

Also, when I try to run the program at smaller spatial intervals (where dx = 20 m, ndx = 1000), I get an arithmetic operations results in overflow error. Not sure what this means.

module global ! Module is used to define arrays and use in main program blocks

implicit none
real len, dx, k, Dl
real, allocatable, dimension (:, :) :: conc ! I used allocatable arrays so that the size of
! the array can be changed, not necessary for this
! simple of a problem
real, allocatable, dimension (:) :: tim
integer ndx,i,n

end module global

program main
use global

real dt , tottime, u

! This is where you write out what you're going to be doing. You can just link to where a
! source file is located, but make sure you do it.

dt = 100. ! unit is seconds
len = 20000. ! unit is meters
tottime = 43200. ! unit is seconds (43200 seconds in half a day)
ndx = 1000 ! number of grid points - we want 10 cells so we start at 1.
u = 0.25 ! unit is m/s
Dl = 50. ! unit is m^2/sec
k = 0.1 / 86400. ! unit is sec^-1
nts = (tottime) / dt ! number of time steps
dx = len / real(ndx) ! grid spacing

! Allocate the arrays, conc(entration), and tim(e).
allocate (conc ( ndx , nts ), tim ( nts ))

conc = 0.0 ! set initial concentration elements to 0 mg/L
tim = 0.0 ! initial time is zero
conc (1, :) = 100. ! set BC on LHS to 100 mg/l for all time

! Run the internal finite difference equation and populate the array conc ()

do n = 1, nts - 1
tim (n) = dt * (n - 1)
do i = 2, ndx - 1
conc (i, n + 1) = conc (i, n) - u * (dt / dx) * (conc (i, n) - conc (i - 1, n)) &
+ Dl * (dt / (dx**2)) * (conc (i + 1, n) - 2 * conc (i, n) + conc (i - 1, n)) - k * dt * conc (i, n)
end do
conc (ndx, n + 1) = conc (ndx, n) - u * (dt / dx) * (conc (ndx, n) - conc (ndx - 1, n)) &
+ Dl * (dt / dx**2) * (2 * conc (ndx - 1, n) - 2 * conc (ndx, n)) - k * dt * conc (ndx, n)
end do

! Set BC on RHS where c (ndx + 1, nts) = c (ndx - 1, nts)
! Write the array to a file

open (10, file = 'C:\g95\fortran\CEE 573\Problem Set 5\upwind20.csv', status = 'unknown')

! Write the column headings
write (10, 100) (tim (n) / 86400., n = 1, nts, 10)
100 format ('X, ',2x, 10000 ('t=', f8.4, ','))

! Write the columns
do i = 1, ndx
x = dx * (i - 1)
write (10, 110) x, (conc (i, n), n = 1, nts, 10)
110 format (f10.1, ',', 10000 (f8.4, ','))
end do

end program main