fd2d.f90

是一个接收永磁体磁场计算的程序

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

 ! 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
  call sparse ( 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 ) 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.'  
!  call timestamp ( )

  call i4vec2_sort_a_plus2 ( nz_num, b_i, b_j, b_u, b_v )

!  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
!
!  About 10 times during the iteration, print out the current time step.