📄 fd2d.f
字号:
program main !*****************************************************************************80!!! MAIN is the main program for FD2D.!! Discussion:!! FD2D simulates the dynamics of a 2D predator-prey system.!! Modified:!! 29 November 2004!! Author:!! Marcus Garvie (MATLAB version)! John Burkardt (FORTRAN90 translation)!! Reference:!! M R Garvie,! Computational Algorithms for Spatially Extended Predator-Prey! Systems with a Holling Type II Functional Response,! To appear.!! Parameters:!! User input, real ( kind = 8 ) A, the left endpoints of the interval.!! User input, real ( kind = 8 ) ALPHA, the value of the parameter ALPHA.!! User input, real ( kind = 8 ) B, the right endpoints of the interval.!! User input, real ( kind = 8 ) BETA, the value of the parameter BETA.!! Local parameter, integer BIGJ, the number of intervals in either! spatial coordinate.!! User input, real ( kind = 8 ) DELT, the time step to use.!! User input, real ( kind = 8 ) DELTA, the value of the parameter DELTA.!! Local parameter, integer DIMJ, the number of nodes in either! ! User input, real ( kind = 8 ) GAMMA, the value of the parameter GAMMA.!! User input, real ( kind = 8 ) H, the "space step", the desired spacing! between nodes in [A,B].!! Local parameter, integer N, the total number of spatial nodes.!! Local parameter, integer SOLVE, 0 for GMRES routine, 1 for Jacobi.!! User input, real ( kind = 8 ) T, the final time. The problem is to be! integrated from 0 to T.! implicit none integer, parameter :: maxl = 10 integer, parameter :: ligw = 20 real ( kind = 8 ) a real ( kind = 8 ) alpha real ( kind = 8 ) b integer, allocatable, dimension ( : ) :: b_i integer, allocatable, dimension ( : ) :: b_j integer, allocatable, dimension ( : ) :: b_ppm real ( kind = 8 ), allocatable, dimension ( : ) :: b_u real ( kind = 8 ), allocatable, dimension ( : ) :: b_v real ( kind = 8 ) beta integer bigj real ( kind = 8 ) delt real ( kind = 8 ) delta real ( kind = 8 ) diff integer dimj real ( kind = 8 ) dummy(1) real ( kind = 8 ) err
real ( kind = 8 ), allocatable, dimension ( : ) :: f
real ( kind = 8 ), allocatable, dimension ( : ) :: g
integer, allocatable, dimension ( : ) :: g_ppm
real ( kind = 8 ) gamma
real ( kind = 8 ) h
real ( kind = 8 ), allocatable, dimension ( : ) :: hhat
integer i
integer i2
integer idummy(1)
integer ierr
integer ierror
integer igwk(ligw)
integer isym
integer it
integer it_max
integer itol
integer iunit
integer j
integer job
integer k
integer lrgw
external matvec_triad
integer mode
external msolve_identity
real ( kind = 8 ) mu
integer n
integer nz
integer nz_num
integer, allocatable, dimension ( : ) :: r_ppm
real ( kind = 8 ), allocatable, dimension ( : ) :: rgwk
integer solve
real ( kind = 8 ) t
integer time_step
integer time_steps
real ( kind = 8 ) tol
real ( kind = 8 ), allocatable, dimension ( : ) :: u
real ( kind = 8 ), allocatable, dimension ( :, : ) :: u_grid
real ( kind = 8 ) u_max
real ( kind = 8 ), allocatable, dimension ( :, : ) :: u0
real ( kind = 8 ), allocatable, dimension ( : ) :: v
real ( kind = 8 ), allocatable, dimension ( :, : ) :: v_grid
real ( kind = 8 ) v_max
real ( kind = 8 ), allocatable, dimension ( :, : ) :: v0
real ( kind = 8 ), allocatable, dimension ( : ) :: x
real ( kind = 8 ), allocatable, dimension ( :, : ) :: x_grid
real ( kind = 8 ), allocatable, dimension ( : ) :: y
real ( kind = 8 ), allocatable, dimension ( :, : ) :: y_grid
real ( kind = 8 ), allocatable, dimension ( : ) :: y1
real ( kind = 8 ), allocatable, dimension ( : ) :: y2
c call timestamp ( )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'FD2D'
write ( *, '(a)' ) ' FORTRAN90 version'
write ( *, '(a)' ) ' A two dimensional finite difference algorithm'
write ( *, '(a)' ) ' for a predator-prey system.'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Enter parameter alpha:'
read ( *, * ) alpha
write ( *, '(a)' ) ' Enter parameter beta:'
read ( *, * ) beta
write ( *, '(a)' ) ' Enter parameter gamma:'
read ( *, * ) gamma
write ( *, '(a)' ) ' Enter parameter delta:'
read ( *, * ) delta
write ( *, '(a)' ) ' Enter a in [a,b]:'
read ( *, * ) a
write ( *, '(a)' ) ' Enter b in [a,b]:'
read ( *, * ) b
write ( *, '(a)' ) ' Enter space-step h:'
read ( *, * ) h
write ( *, '(a)' ) ' Enter maximum time t:'
read ( *, * ) t
write ( *, '(a)' ) ' Enter time-step Delta t:'
read ( *, * ) delt
write ( *, '(a)' ) ' Enter SOLVE ( 0 = GMRES, 1 = Jacobi ):'
read ( *, * ) solve
if ( solve /= 0 ) then
solve = 1
end if
time_steps = nint ( t / delt )
mu = delt / ( h * h )
bigj = ( b - a ) / h
dimj = 1 + int ( ( b - a ) / h )
n = dimj * dimj
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Parameters:'
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' ALPHA = ', alpha
write ( *, '(a,g14.6)' ) ' BETA = ', beta
write ( *, '(a,g14.6)' ) ' GAMMA = ', gamma
write ( *, '(a,g14.6)' ) ' DELTA = ', delta
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' A = ', a
write ( *, '(a,g14.6)' ) ' B = ', b
write ( *, '(a,g14.6)' ) ' H = ', h
write ( *, '(a,i6)' ) ' BIGJ = ', bigj
write ( *, '(a)' ) ' (number of intervals on one side)'
write ( *, '(a,i6)' ) ' DIMJ = ', dimj
write ( *, '(a)' ) ' (number of nodes on one side)'
write ( *, '(a,i12)' ) ' N = ', n
write ( *, '(a)' ) ' (total number of nodes)'
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' T = ', t
write ( *, '(a)' ) ' (final time)'
write ( *, '(a,g14.6)' ) ' DELT = ', delt
write ( *, '(a)' ) ' (time step)'
write ( *, '(a,i12)' ) ' TIME_STEPS = ', time_steps
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' MU = ', mu
write ( *, '(a,g14.6)' ) ' SOLVE = ', mu
if ( solve == 0 ) then
write ( *, '(a)' ) ' = GMRES iteration.'
else
write ( *, '(a)' ) ' = Jacobi iteration.'
end if
allocate ( f(1:n) )
allocate ( g(1:n) )
allocate ( hhat(1:n) )
allocate ( u(1:n) )
allocate ( u_grid(1:dimj,1:dimj) )
allocate ( u0(1:dimj,1:dimj) )
allocate ( v(1:n) )
allocate ( v_grid(1:dimj,1:dimj) )
allocate ( v0(1:dimj,1:dimj) )
allocate ( x(1:dimj) )
allocate ( x_grid(1:dimj,1:dimj) )
allocate ( y(1:dimj) )
allocate ( y_grid(1:dimj,1:dimj) )
allocate ( y1(1:n) )
allocate ( y2(1:n) )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Arrays allocated.'
do i = 1, dimj
x(i) = ( i - 1 ) * h
end do
do i = 1, dimj
y(i) = ( i - 1 ) * h
end do
do j = 1, dimj
x_grid(1:dimj,j) = x(j)
end do
do i = 1, dimj
y_grid(i,1:dimj) = y(i)
end do
call u_init ( alpha, beta, gamma, delta, dimj, dimj, x_grid,
& y_grid, u0 )
call v_init ( alpha, beta, gamma, delta, dimj, dimj, x_grid,
& y_grid, v0 )
k = 0
do j = 1, dimj
do i = 1, dimj
k = k + 1
u(k) = u0(i,j)
v(k) = v0(i,j)
end do
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Arrays initialized.'
!
! MODE = 0, initialize SPARSE.
!
allocate ( b_i(1) )
allocate ( b_j(1) )
allocate ( b_v(1) )
mode = 0
callsparse( 0, 0, 0.0D+00, n, n, mode, nz_num, nz, b_i, b_j, b_v )
!
! MODE = 1, count entries.
! MODE = 2, store entries into B_V.
!
do mode = 1, 2
call sparse(1,1,3.0D+00, n, n, mode, nz_num, nz, b_i, b_j, b_v )
call sparse(1,2,-1.5D+00, n, n, mode, nz_num, nz, b_i, b_j, b_v )
call sparse ( bigj+1, bigj+1, 6.0D+00, n, n, mode, nz_num, nz,
& b_i, b_j, b_v )
call sparse ( bigj+1, bigj, -3.0D+00, n, n, mode, nz_num, nz,
& b_i, b_j, b_v )
do i = 2, bigj
j = i+1
call sparse(i,j,-1.0D+00,n, n,mode,nz_num,nz, b_i, b_j, b_v )
end do
do i = 2, bigj
j = i
call sparse (i,j,4.0D+00,n,n,mode,nz_num, nz, b_i, b_j, b_v )
end do
do i = 2, bigj
j = i-1
call sparse (i,j,-1.0D+00,n,n,mode, nz_num, nz,b_i, b_j, b_v )
end do
!
! Construct block T.
!
i = 1
j = bigj + 2
call sparse (i,j,-1.5D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
i = bigj + 1
j = 2 * bigj + 2
call sparse (i,j,-3.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
do i = 2, bigj
j = bigj + 1 + i
call sparse (i,j,-2.0D+00,n,n, mode, nz_num, nz, b_i, b_j, b_v )
end do
!
! Construct block Z.
!
i = n - bigj
j = n - bigj
call sparse (i,j, 6.0D+00, n, n, mode, nz_num, nz, b_i, b_j, b_v )
i = n - bigj
j = n - bigj + 1
call sparse (i,j,-3.0D+00, n, n, mode, nz_num, nz, b_i, b_j, b_v )
i = n
j = n
call sparse (i,j,3.0D+00,n, n, mode, nz_num, nz, b_i, b_j, b_v )
i = n
j = n - 1
call sparse (i,j, -1.5D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
do i = n-bigj+1, n-1
j = i + 1
call sparse (i,j,-1.0D+00,n, n, mode,nz_num,nz,b_i, b_j, b_v )
end do
do i = n-bigj+1, n-1
j = i
call sparse(i,j,4.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
do i = n-bigj+1, n-1
j = i - 1
call sparse(i,j,-1.0D+00, n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
!
! Construct block Y.
!
i = n - bigj
j = n - ( 2 * bigj + 1 )
call sparse (i,j,-3.0D+00,n,n, mode, nz_num, nz, b_i, b_j, b_v )
i = n
j = n - ( bigj + 1 )
call sparse (i,j,-1.5D+00,n, n, mode, nz_num, nz, b_i, b_j, b_v )
do i = n-bigj+1, n-1
j = i - bigj - 1
call sparse (i,j,-2.0D+00, n,n,mode,nz_num, nz, b_i, b_j, b_v )
end do
!
! Upper W blocks.
!
do i = bigj+2, n-bigj-1
j = i + bigj + 1
call sparse(i,j,-1.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
!
! Lower W blocks.
!
do i = bigj+2, n-bigj-1
j = i - bigj - 1
call sparse (i,j,-1.0D+00,n,n,mode,nz_num,nz, b_i, b_j, b_v )
end do
do i = bigj+2, n-bigj-1
j = i
call sparse (i,j,4.0D+00,n, n,mode, nz_num, nz, b_i, b_j, b_v )
end do
!
! Upper diagonals of X blocks.
!
do i = bigj+2, n-bigj-2
j = i + 1
call sparse (i,j,-1.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
do i = bigj+2, n-2*bigj-1, bigj+1
j = i + 1
call sparse (i,j,-1.0D+00,n, n,mode,nz_num,nz, b_i, b_j, b_v )
end do
do i = 2*bigj+2, n-2*bigj-2, bigj+1
j = i + 1
call sparse (i,j,+1.0D+00, n, n, mode,nz_num, nz,b_i,b_j, b_v )
end do
!
! Lower diagonals of X blocks.
!
do i = bigj+3, n-bigj-1
j = i - 1
call sparse (i,j,-1.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
do i = 2*bigj+2, n-bigj-1, bigj+1
j = i - 1
call sparse (i,j,-1.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
do i = 2*bigj+3, n-2*bigj-1, bigj+1
j = i - 1
call sparse (i,j,+1.0D+00,n,n,mode, nz_num, nz, b_i, b_j, b_v )
end do
!
! Once the entries are counted, you can allocate space.
!
if ( mode == 1 )
c then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'PRED_PREY1:'
write ( *, '(a,i12)' ) ' Maximum number of nonzeros = ', nz_num
deallocate ( b_i )
deallocate ( b_j )
deallocate ( b_v )
allocate ( b_i(1:nz_num) )
allocate ( b_j(1:nz_num) )
allocate ( b_u(1:nz_num) )
allocate ( b_v(1:nz_num) )
end if
end do
!
! Construct B_U = I + mu * L.
!
do k = 1, nz_num
b_u(k) = mu * b_v(k)
i = b_i(k)
j = b_j(k)
if ( i == j ) then
b_u(k) = b_u(k) + 1.0D+00
end if
end do
!
! Construct B_V = I + delta * mu * L.
!
do k = 1, nz_num
b_v(k) = delta * mu * b_v(k)
i = b_i(k)
j = b_j(k)
if ( i == j ) then
b_v(k) = b_v(k) + 1.0D+00
end if
end do
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Sparse matrix values assigned.'
!
! Now sort the (I,J) data, and adjust the B_U and B_V data accordingly.
!
write ( *, * ) ' '
write ( *, '(a)' ) ' Begin sorting.'
c call timestamp ( )
call i4vec2_sort_a_plus2 ( nz_num, b_i, b_j, b_u, b_v )
c call timestamp ( )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Sorting done.'
if ( .false. ) then
do i = 1, nz_num
write ( *, '(2x,i4,2x,i4,2x,g14.6)' ) b_i(i), b_j(i), b_v(i)
end do
end if
!
! Merge any multiple occurrences of (I,J) information.
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Compress the matrix data.'
write ( *, '(a,i12)' ) 'Initial number of matrix entries,
& NZ_NUM = ', nz_num
i2 = 1
do i = 2, nz_num
if ( b_i(i2) == b_i(i) .and. b_j(i2) == b_j(i) ) then
b_u(i2) = b_u(i2) + b_u(i)
b_v(i2) = b_v(i2) + b_v(i)
else
i2 = i2 + 1
b_i(i2) = b_i(i)
b_j(i2) = b_j(i)
b_u(i2) = b_u(i)
b_v(i2) = b_v(i)
end if
end do
nz_num = i2
write ( *, '(a,i12)' ) ' After compression, NZ_NUM = ', nz_num
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Begin time loop.'
!
! Time loop.
!
if ( solve == 0 ) then
isym = 0
itol = 0
tol = 0.00001D+00
it_max = 10
iunit = 0
lrgw = 1 + n * ( maxl + 6 ) + maxl * ( maxl + 3 )
allocate ( rgwk(1:lrgw) )
igwk(1) = maxl
igwk(2) = maxl
igwk(3) = 0
igwk(4) = 0
igwk(5) = 10
else
tol = 0.00001D+00
isym = 0
it_max = 10
job = 0
call ds3_diagonal2 ( n, nz_num, isym, b_i, b_j, b_u, b_v )
end if
do time_step = 1, time_steps
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -