Another piece of evidence that associating overlapping sections of the array may have confused the compiler: I removed the associate construct and C stuff. The revised version of Arjen’s program that used to give asymmetric results now gives symmetric ones agreeing with his other method.
(Sorry about putting the whole program here - I don’t know how to make fortran_lang pull in a file to appear in a separate “window”.)
! advdiff3.f90 – by JFH modifying Arjen’s advdiff.f90 to remove C stuff, flow_velocity and associate,
! and make n = 4 instead of 12 for simplicity.
! Module implementing a straightforward solver for an advection-diffusion
! equation in two dimensions. Uses Arjen’s “symmetric” method 1.
! RESULT: symmetric result. Is there a compiler bug in associate? Or are overlapping arrays illegal?
module advdiff
implicit none
integer,parameter:: c_double = kind(1d0)
character(*),parameter:: progname = ‘advdiff3’
character(12) :: fmta = ‘(NNNNf8.4,A)’ ! NNNN is changed in update
!
! Type holding the various parameters
! (alternatives: put them in separate arguments
! store them in module variables)
!
type :: param_def
real(kind=c_double) :: diffusion ! Diffusion coefficient (m2/s)
real(kind=c_double) :: deltax ! Grid cell size (m)
real(kind=c_double) :: deltat ! Time step (s)
integer :: nreport ! Number of steps to set before returning
end type param_def
contains
subroutine update( n, conc_in, conc_out, params )
integer, intent(in) :: n
real(kind=c_double), intent(in) :: conc_in(n,n)
real(kind=c_double), intent(out) :: conc_out(n,n)
type(param_def), intent(in) :: params
integer :: i
real(kind=c_double) :: diffdt
real(kind=c_double),dimension(n-2,n-2) :: tmp,ccentre,cleft,cright
if(index(fmta,'N')>0) write(fmta(2:5),'(I4)') size(tmp)
! We want a compact calculation
!
conc_out = conc_in
diffdt = params%diffusion / params%deltax**2 * params%deltat
ccentre = conc_out(2:n-1,2:n-1)
cleft = conc_out(1:n-2,2:n-1)
cright = conc_out(3:n ,2:n-1)
do i = 1,params%nreport
!
! Method 1: produces asymmetric result
!
ccentre = ccentre + &
diffdt * ( cleft + cright - 2.0_c_double * ccentre )
!
! Method 2: produces the expected symmetric result
!
!!tmp = diffdt * ( cleft + cright - 2.0_c_double * ccentre )
!!ccentre = ccentre + tmp
!write(20,'(10f8.4)') tmp
conc_out(2:n-1,2:n-1) = ccentre
write(20,fmta) ccentre ,' in update'
enddo
end subroutine update
end module advdiff
! test_advdiff.f90 –
! Test program for the advection-diffusion solver
!
program test_advdiff
use advdiff
implicit none
type(param_def) :: params
integer, parameter :: n = 4 ! was 12
integer :: i, j
real(kind=c_double) :: conc_in(n,n), conc_out(n,n)
open( 20, file = 'test_'//progname//'.out' )
conc_in = 0.0_c_double
conc_in(2:n-1,2:n-1) = 1.0_c_double
params%diffusion = 0.1_c_double
params%deltax = 0.1_c_double
params%deltat = 0.02_c_double
params%nreport = 1
do i = 1,10
call update( n, conc_in, conc_out, params )
write( 20, '(a,i0)' ) 'Program '//progname//' Time: ', i
do j = 1,n
write( 20, fmta) conc_out(:,j),' after update'
enddo
conc_in = conc_out
enddo
end program test_advdiff