As with met_rocket, calls to subroutines have been inserted automatically into the code to capture the results of all numerical computations. These are the routines trace_r8_data, trace_i4_data etc. shown below. On the first run, the results are captured to a binary file, the trace file. On subsequent runs, the trace file is read by the trace routines and the newly computed results are compared with the stored results. In the case of real and complex numbers, if the results differ by more than a criterion percentage the statement is reported as a possible breakdown in precision.
Note that the routines which capture the results have a second integer argument which is a unique identifier for the call-site. This is used to identify the statements where computation has broken down, and also to verify that the program flow has not changed.
!H!**************************************************************************** !H! File: ../fpt_output\ss_9b.f90 !H! Output by FPT 3.8-j Win32 On 30-SEP-15 At 16:45:10 Input files: !H! Main: j:\fpt\emulat~1\ss\fpt\ss_9b_precision_rtt.fsp !H! Current: j:\fpt\emulat~1\ss\original_source\ss_9b.f90 !H! Licensee: SimCon: Development version. !H!**************************************************************************** ! ss_9b.f90 16-Sep-13 John Collins ! 16-Sep-13 John Collins ! PROGRAM ss_9_body ! ! USE module_emulate_real_arithmetic ! IMPLICIT NONE ! INTEGER,PARAMETER :: kr = KIND(1.0D0) ! Real kind INTEGER,PARAMETER :: frames_per_day = 32 ! INTEGER,PARAMETER :: sample_interval = 30*frames_per_day ! 30 days INTEGER,PARAMETER :: end_frame = (10*365+2)*frames_per_day ! ! Note: We measure time in days, distance in AU, mass in kg ! TYPE(em_real_k8),PARAMETER :: au = em_real_k8(149597870700.0D0) TYPE(em_real_k8),PARAMETER :: gu = em_real_k8(-6.67300D-11) ! m**3/kg/s**2 TYPE(em_real_k8),PARAMETER :: g = em_real_k8(((60.0D0*60.0D0*24.0D0)**2/au%value**3)*gu%value) ! au**3/kg/day**2 TYPE(em_real_k8) & time & ,p_time & ! Previous time ,dt & ,end_time & ,mass(9) & ,s(1:3,1:9) & ! Position (axis,body) ,ds(1:3,1:9,1:9) & ! Difference in position, (axis,body1,body2) ,s2(1:3) & ! Square of difference in position, body pair ,r2(1:9,1:9) & ! Square of radial distance, body to body ,r(1:9,1:9) & ! Radial distances, body to body ,gf(1:9,1:9) & ! Total gravitational field, body to body ,proj(1:3,1:9,1:9) & ! Cosine projections on axes (axis,body1,body2) ,vf(1:3,1:9,1:9) & ! Vector field components (axis,body1,body2) ,a(1:3,1:9) & ! Accelerations (axis,body) ,pv(1:3,1:9) & ! Previous velocities (axis,body) ,v(1:3,1:9) & ! Velocities (axis,body) ,pa(1:3,1:9) ! Previous accelerations (axis,body) ! INTEGER & axis & ,body1 & ,body2 & ,frame ! LOGICAL end_run ! DATA & mass /em_real_k8(1.988544D+30) & ! Sun ,em_real_k8(3.302D+23) & ! Mercury ,em_real_k8(48.685D+23) & ! Venus ,em_real_k8(59.7219D+23) & ! Earth ,em_real_k8(6.4185D+23) & ! Mars ,em_real_k8(18981.3D+23) & ! Jupiter ,em_real_k8(5683.19D+23) & ! Saturn ,em_real_k8(868.103D+23) & ! Uranus ,em_real_k8(1024.1D+23)/ & ! Neptune ! ! Initialisation data from JPL, California Institute of Technology, http://ssd.jpl.nasa.gov/horizons.cgi#top ! Positions in AU ,s /em_real_k8(-3.747147775505799D-03),em_real_k8(2.926587035303590D-03),em_real_k8(4.446807296978120D-06) & ! Sun ,em_real_k8(4.715831111175716D-02),em_real_k8(3.056015767171350D-01),em_real_k8(2.006272440161816D-02) & ! Mercury ,em_real_k8(4.959853421725741D-02),em_real_k8(-7.222408775972484D-01),em_real_k8(-1.300479358383397D-02) & ! Venus ,em_real_k8(-1.797649390463135D-01),em_real_k8(9.703470886653921D-01),em_real_k8(-1.743826995752289D-05) & ! Earth ,em_real_k8(-7.334184274376777D-01),em_real_k8(1.457212408029754D+00),em_real_k8(4.839432410544958D-02) & ! Mars ,em_real_k8(4.505316309281244D+00),em_real_k8(-2.163555523180046D+00),em_real_k8(-9.190499900147728D-02) & ! Jupiter ,em_real_k8(-9.468003745121990D+00),em_real_k8(2.616109582176896D-01),em_real_k8(3.722149420390660D-01) & ! Saturn ,em_real_k8(2.003287446208430D+01),em_real_k8(-1.529803159589830D+00),em_real_k8(-2.652256724942619D-01) & ! Uranus ,em_real_k8(2.481340360454546D+01),em_real_k8(-1.689465149553225D+01),em_real_k8(-2.239358924631250D-01)/ & ! Neptune ! ! Velocities in AU per day ,v/em_real_k8(-2.990258994506073D-06),em_real_k8(-5.530878190433255D-06),em_real_k8(6.981801439481427D-08) & ! Sun ,em_real_k8(-3.338934398882634D-02),em_real_k8(5.702520017360675D-03),em_real_k8(3.529925953990557D-03) & ! Mercury ,em_real_k8(2.003227503747324D-02),em_real_k8(1.405250683890825D-03),em_real_k8(-1.136887572032341D-03) & ! Venus ,em_real_k8(-1.720173260060367D-02),em_real_k8(-3.148295197112937D-03),em_real_k8(8.961124784916037D-07) & ! Earth ,em_real_k8(-1.198027375283752D-02),em_real_k8(-5.092533276958321D-03),em_real_k8(1.875893458494313D-04) & ! Mars ,em_real_k8(3.174273016845083D-03),em_real_k8(7.161028006264430D-03),em_real_k8(-1.007937605937214D-04) & ! Jupiter ,em_real_k8(-4.520953258457782D-04),em_real_k8(-5.588245332485017D-03),em_real_k8(1.148915369554184D-04) & ! Saturn ,em_real_k8(2.706976434750568D-04),em_real_k8(3.738322243478300D-03),em_real_k8(1.038310207754998D-05) & ! Uranus ,em_real_k8(1.746288049471184D-03),em_real_k8(2.613667805270097D-03),em_real_k8(-9.397429335776478D-05)/ ! Naptune CALL trace_start_sub_program(1) ! frame = 0 CALL trace_i4_data(frame,181) time = 0.0D0 CALL trace_r8_data(time,182) p_time = 0.0D0 CALL trace_r8_data(p_time,183) ! 1000 CONTINUE ! IF(MOD(frame,sample_interval)==0)THEN WRITE(*,'(F8.2,18E21.12)')time,((s(axis,body1),axis = 1,2),body1 = 1,9) ENDIF frame = frame+1 CALL trace_i4_data(frame,184) ! Note that time will jitter with finite word length. Therefore: time = REAL(frame)/REAL(frames_per_day) CALL trace_r8_data(time,185) dt = time-p_time CALL trace_r8_data(dt,186) p_time = time CALL trace_r8_data(p_time,187) ! DO body1 = 1,9 CALL trace_i4_data(body1,188) DO body2 = body1+1,9 CALL trace_i4_data(body2,189) DO axis = 1,3 CALL trace_i4_data(axis,190) ! Difference in position ds(axis,body1,body2) = s(axis,body1)-s(axis,body2) CALL trace_r8_data(ds(axis,body1,body2),191) ! Square of difference in position s2(axis) = ds(axis,body1,body2)*ds(axis,body1,body2) CALL trace_r8_data(s2(axis),192) ENDDO ! Square of radial distance r2(body1,body2) = s2(1)+s2(2)+s2(3) CALL trace_r8_data(r2(body1,body2),193) ! Radial distance r(body1,body2) = SQRT(r2(body1,body2)) CALL trace_r8_data(r(body1,body2),194) ! Total gravitation field from body2 gf(body1,body2) = g*mass(body2)/r2(body1,body2) CALL trace_r8_data(gf(body1,body2),195) ! Total gravitational field from body1 gf(body2,body1) = g*mass(body1)/r2(body1,body2) CALL trace_r8_data(gf(body2,body1),196) ! Cosines and vector fields to and from each body DO axis = 1,3 CALL trace_i4_data(axis,197) proj(axis,body1,body2) = ds(axis,body1,body2)/r(body1,body2) CALL trace_r8_data(proj(axis,body1,body2),198) proj(axis,body2,body1) = -proj(axis,body1,body2) CALL trace_r8_data(proj(axis,body2,body1),199) vf(axis,body1,body2) = proj(axis,body1,body2)*gf(body1,body2) CALL trace_r8_data(vf(axis,body1,body2),200) vf(axis,body2,body1) = proj(axis,body2,body1)*gf(body2,body1) CALL trace_r8_data(vf(axis,body2,body1),201) ENDDO ENDDO ENDDO ! ! Sum accelerations and integrate DO axis = 1,3 CALL trace_i4_data(axis,202) ! Total accelerations of each body a(axis,1) = vf(axis,1,2)+vf(axis,1,3)+vf(axis,1,4)+vf(axis,1,5)+vf(axis,1,6)+vf(axis,1,7)+vf(axis,1,8)+vf(axis,1,9) CALL trace_r8_data(a(axis,1),203) a(axis,2) = vf(axis,2,1)+vf(axis,2,3)+vf(axis,2,4)+vf(axis,2,5)+vf(axis,2,6)+vf(axis,2,7)+vf(axis,2,8)+vf(axis,2,9) CALL trace_r8_data(a(axis,2),204) a(axis,3) = vf(axis,3,1)+vf(axis,3,2)+vf(axis,3,4)+vf(axis,3,5)+vf(axis,3,6)+vf(axis,3,7)+vf(axis,3,8)+vf(axis,3,9) CALL trace_r8_data(a(axis,3),205) a(axis,4) = vf(axis,4,1)+vf(axis,4,2)+vf(axis,4,3)+vf(axis,4,5)+vf(axis,4,6)+vf(axis,4,7)+vf(axis,4,8)+vf(axis,4,9) CALL trace_r8_data(a(axis,4),206) a(axis,5) = vf(axis,5,1)+vf(axis,5,2)+vf(axis,5,3)+vf(axis,5,4)+vf(axis,5,6)+vf(axis,5,7)+vf(axis,5,8)+vf(axis,5,9) CALL trace_r8_data(a(axis,5),207) a(axis,6) = vf(axis,6,1)+vf(axis,6,2)+vf(axis,6,3)+vf(axis,6,4)+vf(axis,6,5)+vf(axis,6,7)+vf(axis,6,8)+vf(axis,6,9) CALL trace_r8_data(a(axis,6),208) a(axis,7) = vf(axis,7,1)+vf(axis,7,2)+vf(axis,7,3)+vf(axis,7,4)+vf(axis,7,5)+vf(axis,7,6)+vf(axis,7,8)+vf(axis,7,9) CALL trace_r8_data(a(axis,7),209) a(axis,8) = vf(axis,8,1)+vf(axis,8,2)+vf(axis,8,3)+vf(axis,8,4)+vf(axis,8,5)+vf(axis,8,6)+vf(axis,8,7)+vf(axis,8,9) CALL trace_r8_data(a(axis,8),210) a(axis,9) = vf(axis,9,1)+vf(axis,9,2)+vf(axis,9,3)+vf(axis,9,4)+vf(axis,9,5)+vf(axis,9,6)+vf(axis,9,7)+vf(axis,9,8) CALL trace_r8_data(a(axis,9),211) ! ! Integrate for each body DO body1 = 1,9 CALL trace_i4_data(body1,212) IF(frame==1)THEN ! For the first frame we make the previous acceleration equal to the current. This causes an Euler step for ! the velocity integration in this frame. pa(axis,body1) = a(axis,body1) CALL trace_r8_data(pa(axis,body1),213) ENDIF ! Keep previous velocities for Tustin's method pv(axis,body1) = v(axis,body1) CALL trace_r8_data(pv(axis,body1),214) ! Trapezoidal integration for velocities v(axis,body1) = v(axis,body1)+(3*a(axis,body1)-pa(axis,body1))*dt*0.5D0 CALL trace_r8_data(v(axis,body1),215) ! Tustin's method for positions s(axis,body1) = s(axis,body1)+0.5D0*(v(axis,body1)+pv(axis,body1))*dt CALL trace_r8_data(s(axis,body1),216) ! Keep previous accelerations for next frame pa(axis,body1) = a(axis,body1) CALL trace_r8_data(pa(axis,body1),217) ENDDO ENDDO ! end_run = frame>=end_frame CALL trace_l4_data(end_run,218) IF(.NOT.end_run)THEN GOTO 1000 ENDIF ! WRITE(*,'(/,F8.2,18E21.12)')time,((s(axis,body1),axis = 1,2),body1 = 1,9) CALL trace_terminate(2) ! ! END PROGRAM