forums.silverfrost.com

[curoi]

I have code trying to run multiple equations with 3-D array representing a well-mixed tank that is has some inflow = outflow where the 2-D mass transport of some contaminant mass ( c - concentration ) is defined by the advection-diffusion equation (ADE).

This tank has four boundary conditions that define the sides of the tank where

dc/dx = 0 at x = 0
dc/dx = 0 at x = L
dc/dy = 0 at y = 0
dc/dy = 0 at y = B

These boundary conditions yield 9 different finite difference approximation equations for the ADE ( 4 for the sides, 4 for the corners, and 1 for the internal cells)

The problem I am having is that I need to run all 9 equations simultaneously in an iterative way for all time ( using integer n ) so that the do loop for n = 1 to nts encompasses all 9 equations, I get this error from Silverfrost/Plato IDE at line 53:


"System. Overflow Exception

Arithmetic operation resulted in an overflow"

I had originally applied do = 1, nts loop to the internal cell equation and than another do = 1, nts loop for the other 8 equations and it ran the code even though it wasn't calculating the FDA's in the right order.

I'm wondering how to fix this error so that I can calculate all 9 equations within 1 do n = 1, nts do loop.

Fortran 95 code:

program main

! Declaration of variables

implicit none

real :: tottime, length, width, depth, dx, dy, dt, Q, u, D, mass
real, allocatable, dimension (Smile :: time, x, y
real, allocatable, dimension (:,:,Smile :: conc
integer :: i, j, n, ndx, ndy, nts

! Variable values

length = 21.0 ! Length of the well-mixed tank (in meters)
width = 15.5 ! Width of the well-mixed tank (in meters)
depth = 2.0 ! Depth of the well-mixed tank (in meters)
ndx = 21 ! Number of spatial intervals in x-direction
ndy = 31 ! Number of spatial intervals in y-direction
dx = length / real( ndx ) ! Grid spacing in x-direction (in meters)
dy = width / real ( ndy ) ! Grid spacing in the y-direction (in meters)
Q = -3.0 ! Flow rate of water into and out of the tank (m^3/s)
u = Q / ( 15.0 * depth ) ! Velocity of water into and out of the tank (m/s)
D = 100.0 ! Diffusion coefficient (m^2/s) [DL = DT] - Longitudinal (x-direction) is the same as transverse (y-direction)
tottime = 25.0 ! Total time of the system
dt = 0.0005 ! Time step
nts = (tottime) / dt ! Number of time steps
mass = 100. ! Amount of mass injected into the system at x = 15, y = 7.5, and t = 0

! Allocation of time and concentration

allocate ( time( nts ), x( ndx ), y( ndy ), conc( ndy, ndx, nts ) )

conc = 0.0
time = 0.0

! Initial concentration is injected at x = 15 m (i = 16), y = 7.5 m (j = 16), and at t = 0 (n = 1)

conc( 16, 16, 1 ) = mass / ( dx * dy * depth ) * ( 1000000. / ( 100.**3 ) * 1000. )

! INTERNAL CELLS

do n = 2, nts

! INTERNAL CELLS
do i = 2, ndx - 1

do j = 2, ndy - 1

conc( j, i, n ) = conc( j, i, n - 1 ) - u * dt / dx * ( conc( j, i + 1, n - 1 ) - conc( j, i, n - 1 ) ) &
+ D * dt * ( ( 1 / dx**2 ) * ( conc( j, i + 1, n - 1 ) - 2 * conc( j, i, n - 1 ) + conc( j, i - 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( conc( j + 1, i, n - 1 ) - 2 * conc( j, i, n - 1 ) + conc( j - 1, i, n - 1 ) ) ) <------------------- line 53

end do

end do


! BOUNDARY CONDITIONS

! CORNERS

! Bottom-left corner ! is conc supposed to be i + 1 near boundary?
conc( 1, 1, n ) = conc( 1, 1, n - 1 ) - u * dt / dx * ( conc( 1, 2, n - 1 ) - conc( 1, 1, n - 1 ) ) &
+ D * dt * ( ( 1 / dx**2 ) * ( 2 * conc( 1, 2, n - 1 ) - 2 * conc( 1, 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( 2, 1, n - 1 ) - 2 * conc( 1, 1, n - 1 ) ) )

[curoi]

! Top-left corner ! is conc supposed to be i + 1 near boundary?
conc( ndy, 1, n ) = conc( ndy, 1, n - 1 ) - u * dt / dx * ( conc( ndy, 2, n - 1 ) - conc( ndy, 1, n - 1 ) ) &
+ D * dt * ( ( 1 / dx**2 ) * ( 2 * conc( ndy, 2, n - 1 ) - 2 * conc( ndy, 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( ndy - 1, 1, n - 1 ) - 2 * conc( ndy, 1, n - 1 ) ) )

! Bottom-right corner
conc( 1, ndx, n ) = conc( 1, ndx, n - 1 ) - u * dt / dx * ( conc( 1, ndx, n - 1 ) - conc( 1, ndx - 1, n - 1 ) ) &
+ D * dt * ( ( 1 / dx**2 ) * ( 2 * conc( 1, ndx - 1, n - 1 ) - 2 * conc( 1, 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( 2, ndx, n - 1 ) - 2 * conc( 1, ndx, n - 1 ) ) )

! Top-right corner
conc( ndy, ndx, n ) = conc( ndy, ndx, n - 1 ) - u * dt / dx * ( conc( ndy, ndx, n - 1 ) - conc( ndy, ndx - 1, n - 1 ) ) &
+ D * dt * ( ( 1 / dx**2 ) * ( 2 * conc( ndy, ndx - 1, n - 1 ) - 2 * conc( ndy, ndx, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( ndy - 1, ndx, n - 1 ) - 2 * conc( ndy, ndx, n - 1 ) ) )

! SIDES

do i = 2, ndx - 1

! Bottom side
conc( 1, i, n ) = conc( 1, i, n - 1 ) - u * dt / dx * ( conc( 1, i, n - 1 ) - conc( 1, i - 1, n - 1 ) ) + D * dt &
* ( ( 1 / dx**2 ) * ( conc( 1, i + 1, n - 1 ) - 2 * conc( 1, i, n - 1 ) + conc( 1, i - 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( 2, i, n - 1 ) - 2 * conc( 1, 1, n - 1 ) ) )

! Top side
conc( ndy, i, n ) = conc( ndy, i, n - 1 ) - u * dt / dx * ( conc( ndy, i, n - 1 ) - conc( ndy, i - 1, n - 1 ) ) + D * dt &
* ( ( 1 / dx**2 ) * ( conc( ndy, i + 1, n - 1 ) - 2 * conc( ndy, i, n - 1 ) + conc( ndy, i - 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( 2 * conc( ndy - 1, i, n - 1 ) - 2 * conc( ndy, i, n - 1 ) ) )

end do

do j = 2, ndy - 1

! Left side ! is conc supposed to be i + 1 near boundary?
conc( j, 1, n ) = conc( j, 1, n - 1 ) - u * dt / dx * ( conc( j, 1, n - 1 ) - conc( j, 2, n - 1 ) ) + D * dt &
* ( ( 1 / dx**2 ) * ( 2 * conc( j, 2, n - 1 ) - 2 * conc( j, 1, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( conc( j + 1, 1, n - 1 ) - 2 * conc( j, 1, n - 1 ) + conc( j - 1, i, n - 1 ) ) )

! Right side
conc( j, ndx, n ) = conc( j, ndx, n - 1 ) - u * dt / dx * ( conc( j, ndx, n - 1 ) - conc( j, ndx - 1, n - 1 ) ) + D * dt &
* ( ( 1 / dx**2 ) * ( 2 * conc( j, ndx - 1, n - 1 ) - 2 * conc( j, ndx, n - 1 ) ) &
+ ( 1 / dy**2 ) * ( conc( j + 1, ndx, n - 1 ) - 2 * conc( j, ndx, n - 1 ) + conc( j - 1, ndx, n - 1 ) ) )

end do

end do


open(10, file = 'C:\g95\fortran\CEE 573\Project\diffusion25.csv', status = 'unknown')


do j = 1, ndy
y( j ) = dy * ( j - 1 )

write(10,110) ( conc( j, i, 8 ), i = 1, ndx )
110 format( 10000( f12.4, ',' ) )

end do


end program main

[PaulLaidler]

See my response to the same article under "General".