# Curious result in a simple program

This weekend I ran into a rather unexpected issue and I hope you can point out my error. I simply fail to see it and that annoys me.

The program is fairly straightforward: an attempt to solve an advection-diffusion problem with a compact implementation. In the absence of advection, the result should be symmetric.It is if I use method 2 in the code. But with method 1, the result is not and I fail to see why.

Here is the code:

``````! advdiff.f90 --
!     Module implementing a straightforeward solver for an advection-diffusion
!     equation in two dimensions
!
use iso_c_binding

implicit none

!
! Type holding the various parameters
! (alternatives: put them in separate arguments
! store them in module variables)
!
type, bind(C) :: param_def
real(kind=c_double) :: diffusion      ! Diffusion coefficient (m2/s)
real(kind=c_double) :: flow_velocity  ! Flow velocity (left to right; m/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 ) bind(c, name = 'update')
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, veldt

real(kind=c_double) :: tmp(n-2,n-2)

!
! We want a compact calculation
!
conc_out = conc_in

diffdt   = params%diffusion / params%deltax**2  * params%deltat
veldt    = params%flow_velocity / params%deltax * params%deltat

!
! Using associate we can formulate the difference equation in a way that avoids
! explicit loops and boundary conditions.
! Note:
! It is not necessarily very efficient ;)
!
! By first copying the initial values and then using the output array, we
! avoid intermediate copying.
!
associate( 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 + veldt  * ( cleft - ccentre ) + &
!!                    diffdt * ( cleft + cright - 2.0_c_double * ccentre )
!
! Method 2: produces the expected symmetric result
!
tmp     =           veldt  * ( cleft - ccentre ) + &
diffdt * ( cleft + cright - 2.0_c_double * ccentre )
ccentre  = ccentre + tmp
!write(20,'(10f8.4)') tmp
write(20,'(10f8.4)') ccentre
enddo
end associate

end subroutine update

!     Test program for the advection-diffusion solver
!
implicit none

type(param_def) :: params

integer, parameter  :: n = 12
integer             :: i, j
real(kind=c_double) :: conc_in(n,n), conc_out(n,n)

open( 20, file = 'test_advdiff.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%flow_velocity = 0.0_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)' ) 'Time: ', i
do j = 1,n
write( 20, '(100f8.4)' ) conc_out(:,j)
enddo

conc_in = conc_out
enddo
``````

The asymmetry occurs right at the first time step.

Also: the results are the same with gfortran and Intel Fortran oneAPI, so I suspect some problem in my code, but I do not see it 1 Like
``````program test_advdiff
implicit none
integer :: cout(3)
cout=[0,1,1]
associate( ccentre  => cout(2:3) , &
cleft    => cout(1:2)  &
)
![1,1] + [0,1] - [1,1]
ccentre=ccentre+(cleft-ccentre)
! [1,0]
write(*,*)ccentre
end associate
``````

I simplify the code, And I think it is the bug for gfortran and ifort. nvfortran shows the right result.

``````> nvfortran
0 1
> gfortran
0 0
> ifort
0 0
``````

and if change it as ` cleft => (cout(1:2))`,the result is right.

1 Like

Here is another minimal program, which gives a different result only for Pattern 1 (among the five patterns tried). Because I never use `associate` for array sections or expressions, I am not sure whether the result for Pattern 1 is correct or not…

``````program main
implicit none
integer a( 3 ), b( 3 ), c( 3 ), d( 3 ), tmp( 2 )
integer, target :: e( 3 )
integer, pointer :: pinp(:), pout(:)

a = [1,2,3]
b = a
c = a
d = a
e = a

!! Pat-1: associate by reference (?)
associate( inp => a( 1 : 2 ), &
out => a( 2 : 3 ) )
out = inp
end associate

!! Pat-2: associate by value (?)
associate( inp => ( b( 1 : 2 ) ), &
out => b( 2 : 3 ) )
out = inp
end associate

!! Pat-3: use a temporary array
tmp(:) = c( 1 : 2 )
c( 2 : 3 ) = tmp(:)

!! Pat-4: direct assignment of array sections
d( 2 : 3 ) = d( 1 : 2 )

!! Pat-5: use array pointers
pinp => e( 1 : 2 )
pout => e( 2 : 3 )
pout = pinp

print *, "a = ", a
print *, "b = ", b
print *, "c = ", c
print *, "d = ", d
print *, "e = ", e
end
``````

gfortran-10 test.f90

`````` a =            1           1           1
b =            1           1           2
c =            1           1           2
d =            1           1           2
e =            1           1           2
``````

Not a fan of associate, but are you sure the assignment is legal? Perhaps someone could comment on whether this applies ( I personally avoid ASSOCIATE) but perhaps this phrase from the standard implies the assignment result is not defined?

• If selector is not a variable or is a variable that has a vector
subscript, associate-name shall not appear in a variable definition
context(16.6.7).

What you are doing appears to fit the description of assigning to something with a vector subscript.

A vector subscript is an array of individual indices, but I use array sections. IIRC the “shall” bit in the above statement means that the compiler is supposed not to accept it. But I immediately admit that my knowledge of standardese is extremely limited.

The alternative programs are intriguing though.

Hm, remarkably enough, both gfortran and Intel Fortran oneAPI accept the following program - and produce an executable that does the right thing:

``````program chk_associate_vector
implicit none

real    :: array(10)
integer :: i

array = [ (1.0*i, i = 1,size(array)) ]

associate( a => array([2,4,6,8,10]) )
write(*,*) a
end associate
end program chk_associate_vector

``````

(A more arbitrary vector of indices is also accepted and the resulting program does the right thing…)

While the selector is a variable that has a vector subscript, the associate name (`a`) does not appear in a variable definition context, and thus this program is valid. However, adding an assignment to `a`, as in the following, makes it invalid.

``````program chk_associate_vector
implicit none

real    :: array(10)
integer :: i

array = [ (1.0*i, i = 1,size(array)) ]

associate( a => array([2,4,6,8,10]) )
a = 0.0
write(*,*) a
end associate
end program chk_associate_vector
``````

Looking at the subsequent examples and then the original code, I notice your associate names are to overlapping sections of the array, thus creating aliasing between elements. I’m guessing this is confusing the compiler, and that it is doing evaluation and assignment one element at a time, when it really, in this case, must create the entire array temporary before doing the assignment.

1 Like

Yes, as one gets into serious calculations and simulations involving the advection-diffusion problem, employing objects declared explicitly with the `POINTER` attribute will prove effective across the model, `ASSOCIATE` is more of a local solution:

``````    integer, target :: cout(3)
integer, pointer :: ccentre(:), cleft(:)
cout=[0,1,1]
ccentre => cout(2:3)
cleft => cout(1:2)
ccentre = ccentre + (cleft-ccentre)
print *, ccentre
print *, cout
end
``````

@Arjen, this reminds of a support request re: a possible bug I had submitted a long time ago with Intel Support Center, it is likely not fixed yet.

Meaning, can you try a workaround with Intel oneAPI `IFORT`?

And that is to retry with `IFORT` after applying the `TARGET` attribute on your `conc_out` dummy argument in your `update` subprogram.

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?

implicit none
integer,parameter:: c_double = kind(1d0)
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

! Test program for the advection-diffusion solver
!
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
``````