Fortran NBody sim code won't run for high year values

Question

The program is as follows.

The issue occurs when I try to run the code for >~80 years, at which point the code apparently 'runs' instantly, generating an empty text file. The code runs fine for smaller timescales.

PROGRAM NBody
IMPLICIT NONE

DOUBLE PRECISION:: m(1:10), deltaR(1:3)
DOUBLE PRECISION:: G, r
DOUBLE PRECISION, DIMENSION(10,3):: pos, v, a0, a1 !x, y, z
INTEGER:: n,i,j,k,stepsize, year, zero, length
CHARACTER(len=13):: fname !xxxyrxxpl.txt
zero = 0

m(1) = 1988500e24    !sun
m(2) = 0.33e24    !mercury
m(3) = 4.87e24    !venus
m(4) = 5.97e24    !earth
m(5) = 0.642e24    !mars
m(6) = 1898e24    !jupiter
m(7) = 568e24    !saturn
m(8) = 86.8e24    !uranus
m(9) = 102e24    !!neptune
m(10) = 0.0146e24    !pluto

!Initial POS
pos = zero
pos(2,1) = 57.9e9
pos(3,1) = 108e9
pos(4,1) = 149e9
pos(5,1) = 227e9
pos(6,1) = 778e9
pos(7,1) = 1352.6e9
pos(8,1) = 2741.3e9
pos(9,1) = 4444.5e9
pos(10,1) = 4436.8e9
!FORTRAN works column,row: (particle,x/y/z)
!Momentum is initially non-zero due to planet and velocity starting points. Figure out a solution.

!Initial velocity
v = zero
v(2,2) = 47.4e3
v(3,2) = 35e3
v(4,2) = 29.8e3
v(5,2) = 24.1e3
v(6,2) = 13.1e3
v(7,2) = 9.7e3
v(8,2) = 6.8e3
v(9,2) = 5.4e3
v(10,2) = 4.7e3

g = 6.67e-11
stepsize = 1800 !3600 = 1 hour
year = 3.154e+7

!Calculate initial values
a0 = 0
a1 = 0
do i = 1,10
    do j = 1,10
    if(i==j) cycle
    deltaR(:) = (pos(i,:)-pos(j,:))
    r = -sqrt((pos(i,1)-pos(j,1))**2+(pos(i,2)-pos(j,2))**2+(pos(i,3)-pos(j,3))**2)
    a0(i,:) = a0(i,:) + g*M(j)*deltaR*r**(-3)
    END DO
END DO

write(6,*) "Specify length in years"
read (*,*) length
write(6,*) "Specify file name (xxxYRzzPL.txt)"
read(*,*) fname

!Just above is where I call for a length in the terminal, values of 40 will work, much higher do not. I don't know the exact cut-off.


open (unit = 2, file = fname)

!Do loop over time, planet and partners to step positions
do k=0, length*year,stepsize
    write(2,*) pos
    pos = pos + v*stepsize + 0.5*a0*stepsize**2
    do i = 1,10
        do j = 1,10
        if(i==j) cycle
        deltaR(:) = (pos(i,:)-pos(j,:))
        r = -sqrt((pos(i,1)-pos(j,1))**2+(pos(i,2)-pos(j,2))**2+(pos(i,3)-pos(j,3))**2)
        a1(i,:) = a1(i,:) + G*M(j)*deltaR/r**3
        END DO
    END DO
    v = v + 0.5*(a0+a1)*stepsize
    a0=a1
    a1=0

END DO

close (2)

END PROGRAM

I suspected it could be an issue with variable storage but I can't see any problems there.

So it is an issue with the number being too big? How do I fix this problem (I'm still new to fortran) — Craig Gardner, Feb 26 '18 at 16:50
If that is the only problem, then you'll find discussion of integer overflow in questions like [this one](https://stackoverflow.com/q/24460466/3157076). — francescalus, Feb 26 '18 at 17:20

score 3 · Accepted Answer · answered Feb 26 '18 at 17:29

Using an iterator like this can be dubious. Even an 8 byte integer will overflow if you go long enough. Considering how this code is set up, I would do something like this:

do iYear = 1, length
   do k = 0, year, stepsize
   ....
   enddo
enddo

Inner do loop loops over one year. Outer do loop loops over the years. Could go Gigayears like this with just 4 byte integers if you want to wait that long.

I would likely rename your variables too to make more sense. This would look better:

do iYear = 1, nYears
   do k = 0, YearLength, stepsize
   ....
   enddo
enddo

DavidH · Answer 2 · 2018-02-26T17:12:22.163

2

Expanding on @francescalus, you may need to specify your integers as 8-bytes rather than the default 4:

integer, parameter :: c_int8 = selected_int_kind (10)
integer(kind = c_int8) :: n,i,j,k,stepsize, year, zero, length

EDIT I added a parameter to determine the correct value for a 64-bit integer intrinsically.

edited Feb 26 '18 at 17:12

answered Feb 26 '18 at 16:57

DavidH

415
4
21

1

*"Note that .. the 8 in kind = 8 may need to be something different.."* That's why using a straight literal number is highly not-recommended and considered ugly. See https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4 for better alternatives. Especially considering that here we DO know the required supported value. – Vladimir F Героям слава Feb 26 '18 at 17:02
I usually specify the different `kind`s in a module (e.g. `int4 = 4`, then `kind = int4`). I didn't know you could do this with intrinsic functions, so thanks for the link! – DavidH Feb 26 '18 at 17:06
This didn't appear to work, however @Dan_Sp solution worked perfectly. – Craig Gardner Feb 26 '18 at 18:10
2

edit: I understand why it didn't work, it was my fault: I didn't change my count variable k to a int8, I changed the variable storing the length, which was never the issue. – Craig Gardner Feb 26 '18 at 18:17

Fortran NBody sim code won't run for high year values

2 Answers2