Help Solve FORTRAN Post Newtonian Gravity Simulation

Click For Summary
SUMMARY

The discussion centers on implementing a Post Newtonian Expansion for General Relativity simulations using FORTRAN. The user has successfully created a Newtonian gravity simulation with four different integrators but is encountering a significant error in their calculations, specifically being off by two orders of magnitude. They seek assistance in identifying the error in their code, which includes various calculations for gravitational forces and accelerations among multiple bodies. The user acknowledges that optimization is not yet a priority, focusing first on achieving correct functionality.

PREREQUISITES
  • FORTRAN programming language proficiency
  • Understanding of Newtonian gravity principles
  • Familiarity with General Relativity concepts
  • Knowledge of numerical integration techniques
NEXT STEPS
  • Research FORTRAN best practices for numerical simulations
  • Study the Post Newtonian Expansion in detail
  • Learn about debugging techniques specific to FORTRAN
  • Explore optimization strategies for FORTRAN code
USEFUL FOR

Researchers, physicists, and developers interested in gravitational simulations, particularly those using FORTRAN for numerical modeling and simulations of physical phenomena.

goalie39
Messages
2
Reaction score
0
Hi all:

In my free time, I've been playing with creating a code to help do toy simulations of gravity in my personal code of choice, FORTRAN. The first step was to get Newtonian gravity up and running and this has been pretty successful so far and I've been able to implement about 4 different integrators to play and test with. My next step is to try to get General Relativity going via the Post Newtonian Expansion.

I'm all self taught in this area, so I've been leaning on what I could find by scouring the internet. I found, what seemed to be, a fairly straight forward way to implement this in my code (see the attached image), however, somewhere I've gone wrong and I'm about 2 orders of magnitude off in my calculation. The source of this is here, a website that seems to be dead these days.

I've been staring at this quite a bit and I'm thinking that there has to be something simple I'm just not catching. My hunch is that someone who knows what they're doing might be able to look at it and understand my error fairly quickly. If anyone has the time to read through the attached code, I'd truly appreciate the assistance.

PS - I realize that things are definitely not optimized yet. I'd like to get it working and then focus on making it fast.

Thanks everyone.

Code: 
DYDT%DDY = 0D0
DYDT%DDZ = 0D0
C1(:)    = 0D0
C2(:)    = 0D0
C3(:)    = 0D0
A10(:)   = 0D0
A1       = 0D0
X_I      = BODIES(BODY)%MYPOSITION(1)
Y_I      = BODIES(BODY)%MYPOSITION(2)
Z_I      = BODIES(BODY)%MYPOSITION(3)
VX_I     = BODIES(BODY)%VELOCITY(1)
VY_I     = BODIES(BODY)%VELOCITY(2)
VZ_I     = BODIES(BODY)%VELOCITY(3)
AX_I     = BODIES(BODY)%ACCELERATION(1)
AY_I     = BODIES(BODY)%ACCELERATION(2)
AZ_I     = BODIES(BODY)%ACCELERATION(3)
M        = BODIES(BODY)%MASS
BETA  = 1D0
GAMMA = 1D0

A2 = 0D0
DO K = 1,NBODIES
    IF(K.EQ.BODY)CYCLE
    DX2 = BODIES(K)%MYPOSITION(1) - X_I
    DY2 = BODIES(K)%MYPOSITION(2) - Y_I
    DZ2 = BODIES(K)%MYPOSITION(3) - Z_I
    A2 = A2 + (G_C*BODIES(K)%MASS) / ((DX2*DX2+DY2*DY2+DZ2*DZ2)**(0.5D0))
ENDDO
A2 = 2D0*(BETA+GAMMA)*A2

DO J = 1,NBODIES
    IF(J.EQ.BODY)CYCLE

    X_J  = BODIES(J)%MYPOSITION(1)
    Y_J  = BODIES(J)%MYPOSITION(2)
    Z_J  = BODIES(J)%MYPOSITION(3)
    VX_J = BODIES(J)%VELOCITY(1)
    VY_J = BODIES(J)%VELOCITY(2)
    VZ_J = BODIES(J)%VELOCITY(3)
    AX_J = BODIES(J)%ACCELERATION(1)
    AY_J = BODIES(J)%ACCELERATION(2)
    AZ_J = BODIES(J)%ACCELERATION(3)
    M_J  = BODIES(J)%MASS
   
    !...Distance between bodies
    DX = X_J - X_I
    DY = Y_J - Y_I
    DZ = Z_J - Z_I

    A1 = A1 + (G_C*M_J) / ((DX*DX+DY*DY+DZ*DZ)**(1.5D0))

    A3 = 0D0
    DO K = 1,NBODIES
        IF(K.EQ.J)CYCLE
        DX2 = BODIES(K)%MYPOSITION(1) - X_J
        DY2 = BODIES(K)%MYPOSITION(2) - Y_J
        DZ2 = BODIES(K)%MYPOSITION(3) - Z_J
        A3 = A3 + (G_C*BODIES(K)%MASS) / ((DX2*DX2+DY2*DY2+DZ2*DZ2)**(0.5D0))
    ENDDO
    A3 = (2D0*BETA-1D0)*A3

    !...Sum of square velocities
    A5 = GAMMA*(VX_I*VX_I+VY_I*VY_I+VZ_I*VZ_I)
   
    !...Sum of square velocities - Sum of product of velocities
    A6 = (1D0+GAMMA)*((VX_J*VX_J+VY_J*VY_J+VZ_J*VZ_J) - & 
                 2D0*(VX_I*VX_J-VY_I*VY_J-VZ_I*VZ_J))
   
    !...DX*Velocity Term / Square distance
    A7 = 1.5D0*((-DX*VX_J-DY*VY_J-DZ*VX_J)**2D0  &
             /  (DX*DX+DY*DY+DZ*DZ))
         
    !...DX*Acceleration Term * DX
    A8 = 0.5D0*(DX*AX_J+DY*AY_J+DZ*AZ_J)

    !...DX * GAMMA * Velocity
    A9 = (DX*(2D0+2D0*GAMMA)*VX_I - (1D0+2D0*GAMMA)*VX_J) + &
         (DY*(2D0+2D0*GAMMA)*VY_I - (1D0+2D0*GAMMA)*VY_J) + &
         (DZ*(2D0+2D0*GAMMA)*VZ_I - (1D0+2D0*GAMMA)*VZ_J) 

    A10(1) = A10(1) + (G_C*M_J*AX_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))
    A10(2) = A10(2) + (G_C*M_J*AY_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))
    A10(3) = A10(3) + (G_C*M_J*AZ_J)/((DX*DX+DY*DY+DZ*DZ)**(0.5D0))

    !...Combination of terms
    C1(1)  = C1(1) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DX ) + A9 * (VX_I - VX_J)
    C1(2)  = C1(2) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DY ) + A9 * (VY_I - VY_J)
    C1(3)  = C1(3) + ( ( C*C - A2 - A3 + A5 + A6 - A7 + A8 )* DZ ) + A9 * (VZ_I - VZ_J)

ENDDO

DYDT%DDX = ((A1*C1(1)) + 0.5D0*(3D0+4D0*GAMMA)*A10(1))/(C*C)
DYDT%DDY = ((A1*C1(2)) + 0.5D0*(3D0+4D0*GAMMA)*A10(2))/(C*C)
DYDT%DDZ = ((A1*C1(3)) + 0.5D0*(3D0+4D0*GAMMA)*A10(3))/(C*C)
 

Attachments

  • gr.png
    gr.png
    19 KB · Views: 543
Physics news on Phys.org
My apologies, I guess this should have gone in the homework area. If someone could move it, I'd appreciate it. Sorry for the confusion.
 
goalie39 said:
I guess this should have gone in the homework area.

Not really, since it's not a homework problem, it's just something you're doing in your free time.

However, I'm not sure how many FORTRAN coders there are here, so you might not get much useful input.