Ошибка при использовании переменных в модуле Fortran - PullRequest
Новая
       104

Ошибка при использовании переменных в модуле Fortran

1 голос
/ 16 января 2020

Я сталкиваюсь с очень и очень странной ошибкой при использовании переменных в модуле Fortran. Это довольно простой код. Каждый столбец F рассчитывается по каждому столбцу X. X и F - все переменные в модуле. Ниже приведен код.

В global.f90: объявление модуля, включающего X и F. 

module global
implicit none
!  Parameters
real(8),parameter :: GravConst = 6.674184D-11
real(8),parameter :: PI = 3.141592653589793D0
integer :: N
real(8) :: X(3,2)
real(8) :: F(3,2)
real(8) :: omega
real(8) :: density
end module

В main.f90: вычисление F. 

program main
use global
implicit none

integer I,J,K
N = 2
F = 0.0D0
X = 0.0D0
omega = 1.0
density = 0.0

!  read positions
open (1,FILE="pos.txt")
do I = 1,N
    read(1,*) (X(K,I),K=1,3)
end do

write(*,*) 'X:',(X(J,2),J=1,3)
write(*,*) 'F:',(F(J,2),J=1,3)
write(*,*) 'F(1,2):',omega*omega*X(1,2)- 4.0D0/3.0D0*PI*GravConst*density*X(1,2)
write(*,*) 'F(2,2):',omega*omega*X(2,2)- 4.0D0/3.0D0*PI*GravConst*density*X(2,2)
do I = 1,N
    !yorp force
    do K = 1,2
        F(K,I) = F(K,I) + omega*omega*X(K,I)
    end do  
    !gravity      
    do K = 1,3
        F(K,I) = F(K,I) - 4.0D0/3.0D0*PI*GravConst*density*X(K,I)
    end do
end do
write(*,*) (F(J,2),J=1,3)
end

Значения F (1,2) и F (2,2), рассчитанные напрямую, печатаются до l oop. После l oop печатаются окончательные значения F (:, 2). Вот результат использования intel / 18.0.2. Очевидно, что результат, вычисленный как l oop, неверен. 

X: -4.97716000000000       0.897938000000000       -11.4228000000000
F: 0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
F(1,2): -4.97716000000000
F(2,2): 0.897938000000000
-4.69400000000000       -4.97716000000000    0.000000000000000E+000

Кажется, что-то странное происходит в l oop ... Но если я изменю порядок вычисления в l oop, результат будет правильным. 

do I = 1,N 
    !gravity      
    do K = 1,3
        F(K,I) = F(K,I) - 4.0D0/3.0D0*PI*GravConst*density*X(K,I)
    end do
    !yorp force
    do K = 1,2
        F(K,I) = F(K,I) + omega*omega*X(K,I)
    end do 
end do

X: -4.97716000000000       0.897938000000000       -11.4228000000000
F: 0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
F(1,2): -4.97716000000000
F(2,2): 0.897938000000000
-4.97716000000000       0.897938000000000    0.000000000000000E+000

Если я использую локальные переменные, а не переменные в модуле, результат также будет правильным.

Ниже приведены pos.txt и makefile, которые я использовал. 

2.64530000  -7.20517000     -4.69400000 
-4.97716000     0.89793800  -11.42280000

# Start of the makefile
COMPILER=ifort -O2
DEMBody: global.o main.o
    $(COMPILER)  -o DEMBody global.o main.o
global.mod: global.o global.f90
    $(COMPILER)  -c global.f90
global.o: global.f90
    $(COMPILER)  -c global.f90
main.o: global.mod main.f90
    $(COMPILER)  -c main.f90
clean:
    rm global.mod global.o main.o
# End of the makefile

Или вы мог бы попробовать этот вход. Это было бы более понятно. 

X(1,1) = 1.0
X(2,1) = 2.0
X(3,1) = 3.0
X(1,2) = 4.0
X(2,2) = 5.0
X(3,2) = 6.0

Вы получите запутанный результат для 'F': 

F(:,1): 1.000        2.000       0.000E+000
F(:,2): 3.000        4.000       0.000E+000

Правильный результат должен быть 

F(:,1): 1.000        2.000       0.000E+000
F(:,2): 4.000        5.000       0.000E+000

Это кажется, что элемент X / F был перемещен на одну клетку.

...