Fortran 95 - ADE Implicit Scheme Tridag. Matrix

1. The problem statement, all variables and given/known data

The problem statement is to use Fortran 95 to code a forward time, centered space numerical solution to the 1-D (x-direction) Advection-Dispersion Equation:

dc/dt = u(dc/dx) + D(d2c/dx2) - kc

where c is the concentration of contaminant, u is the advecting velocity, D is the dispersion coefficient, and k is the decay rate.

I am trying to figure out why the following code is writing the solutions out in excel without including the final time step at 0.5 days and the final spatial interval at x = 20,000 meters.

2. Relevant equations

Initial condition:

C(x, t=0) = 0 mg/L
C(x=0, t) = 100 mg/L

dc/dx(X=L, t) = 0

3. The attempt at a solution

I decided not to post the whole code because I did not write the program, so it may not make complete sense. If it doesn't make sense or more coding information is needed, just let me know.

nts is the total number of time steps all the way up to 0.5 days where dt =100 sec.
ndx is the total number of spatial intervals all the way up to L = 20,000 m and dx = 200 m

do n=1,nts-1 ! change from nts-1 to nts to account for last missing time step?
tim(n)=dt*(n-1)
r=0.

do i=1,ndx

do j=1,ndx
r(i)=r(i)+d(i,j)*conc(j,n)

end do

r(i)=r(i)+q(i) ! This is the known RHS matrix elements

end do

call tridag ! This solves the banded tridiagonal matrix using the Thomas algorithm - see my class notes

end do

open(10,file='CNimplicitFTCS.xls',status='unknown')

write(10,100)(tim(n)/86400.,n=1,nts,10)
100 format('X,',2x,10000('t=',f8.4,','))
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

stop

end program main
 PhysOrg.com science news on PhysOrg.com >> King Richard III found in 'untidy lozenge-shaped grave'>> Google Drive sports new view and scan enhancements>> Researcher admits mistakes in stem cell study
 You are writing every 10th result. Are the total number of results exactly divisible by 10? If not, you miss the last one.
 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

Fortran 95 - ADE Implicit Scheme Tridag. Matrix

What kind of integration are you using? It looks like it is explicit in nature. By that I mean that answers at the new time depend solely on the answers from the previous time step. This type of integration can be unstable if time step is too large. A fully implicit scheme is unconditionally stable.