Fortran 95 - ADE Implicit Scheme Tridag. Matrix

[edge333]

Homework Statement

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.

Homework Equations

Initial condition:

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

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

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

 

[LawrenceC]

You are writing every 10th result. Are the total number of results exactly divisible by 10? If not, you miss the last one.

 

[edge333]

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

 

[LawrenceC]

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.

 

[DeeNos]

Thank you for sharing your code and problem statement. Based on the information provided, it is difficult to determine the exact cause of the issue you are experiencing. However, there are a few things that you can check to try and identify the problem.

Firstly, it is important to check the boundary conditions for your problem. In the code snippet provided, the boundary condition at x = L is set to zero, which may not be correct depending on the problem setup. It is also important to check if the boundary condition at x = 0 is being applied correctly.

Secondly, it is important to check the values of nts and ndx to ensure that they are appropriate for your problem. If these values are not correct, it could lead to incorrect output.

Additionally, it is important to check if the values of dt and dx are appropriate for your problem. These values can greatly affect the accuracy of your solution.

Lastly, it may be helpful to check the tridag subroutine to ensure that it is correctly solving the banded tridiagonal matrix.

I hope these suggestions are helpful in identifying the issue with your code. If you need further assistance, please do not hesitate to ask.