Fortran 90 N-body problem code not working, help please!

[JoshuaJoseph23]

Hey guys,

as an assignment, we are coding up the N-body problem for the solar system. I'm running into some problems however.

We need to code up the forces from one body to another, F_ij, the force on i due to j. We have been told do not code it in matrix form, as for problems with a million bodies, this becomes a 3 * a million * a million matrix, which represents formiddable memory challenges. So we have to do it as a rank 1 array, i.e. F = [F_1x, F_1y, F_1z, F_2x, F_2y, F_2z, ...] where for example, F_1x is the total force on body 1 in the x-direction.

The force has the form F_ij = G * m_i * m_j * (r_i - r_j)/((abs(r_i - r_j)^3)) - just Newton's force law. Another condition to be satisfied is F_ij = - F_ji. Obviously F_ii = 0.

I'm running into a big problem here, as I have no idea how to do it. What I have done so far hasn't worked. My latest attempt is below. nDim = number of dimensions, so 3, and in this case we're only doing the sun and the planets up to Jupiter, so 6 bodies (nBody). r is a rank-2 array with the positions of each of the bodies, with dimension r(nDim,nBody), mass is a rank-1 array of the masses of the bodies, F is the force array that I am trying to solve for. I have also defined another variable, Force, which is an intermediate step which I don't believe I need but I'm really stuck. nVector is just nDim * nBody * 2, that was defined for a different purpose. Never mind it, shouldn't make a difference.

subroutine forceij(nDim,nBody,nVector,mass,r,F)
real(DP), dimension,, intent(in) :: r,
real(DP), dimension), intent(in) :: mass
integer, intent(in) :: nBody, nDim, nVector
real(DP), dimension), intent(out), allocatable :: F)

integer :: i, j, k ! counting variables
real(DP), dimension(nDim,nBody) :: Force

allocate(F(nVector/2))

do i = 1, nBody-1
do j = i+1, nBody

Force(i,j) = (G * mass(i) * mass(j) * (r(i,i) - r(i,j))/((abs(r(i,i) - r(i,j)))**3))

F(k) = sum(Force,j))
end do
end do
end subroutine forceij

Any help/suggestions/tips and whatnot is appreciated, thanks!

 

[mikrom]

I'm not sure if I understood you right, because you write once about [F_1x, F_1y, F_1z, F_2x, F_2y, F_2z, ...] and then about F_ij=... Does this mean that you have n*3 matrix F_ij where i=1,..,n and j=1,2,3 (or x,y,z) ?
If so, then you want to store this matrix as an 1-dimensional array?
In this case the relation between the elements of the Matrix M and the Vector v could be: M(i,j) = v(3*(i-1) + j)
so M(1,1)= v(1), M(1,2) = v(2), .., M(3,2)=v(8), M(3,3)= v(9), ..etc

In this case I would first compute all forces into the array, e.g:

subroutine set_force(force_array,...)
  real, dimesnion(3*n) :: force_array
  real :: f
  ...
  do i=1,n
    do j=1,3
      ! compute Force
      f = (G * mass(i) * mass(j) * (r(i,i) - r(i,j))/((abs(r(i,i) - r(i,j)))**3))
      ! store Force into array
      force_array(3*(i-1) + j) = f
    end do
  end do
end subroutine set_force

and for further computation with forces I would create a function which for given values of i and j returns the value of the force from the array

real function get_force(i,j,...)
  ...
  ! return force from array
  get_force = force_array(3*(i-1) + j)
end function get_force

 

[JoshuaJoseph23]

Okay great, thanks for your help!