Curious result in a simple program

Help
1

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
!
module advdiff
    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

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 = 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
end program test_advdiff

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 :frowning:

1 Like
2 
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
end program test_advdiff

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
3

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
4

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.

5

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.

6

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…)

7

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
8

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
9

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
10

@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.

11

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