Skip to content

Commit 0ae6e24

Browse files
author
Berry Wang
committed
refactored rkf45() since there are some problems with the original one
1 parent 2bf6422 commit 0ae6e24

7 files changed

+382
-25
lines changed

‎solver.f90

+1-9
Original file line numberDiff line numberDiff line change
@@ -44,36 +44,28 @@ module solver
4444
!------------------------------------------------------------
4545

4646
subroutine rkf45_step(t,h,w,err)
47-
4847
real*16, intent(in) :: t, h
4948
type(triBody) :: w, w1, z1
5049
real*16 :: err
5150
type(triBody), dimension(6) :: s
5251
integer :: i, j
53-
54-
! w1 currently used as a temporary variable
5552
do i=1,6
5653
w1 = w
5754
do j=1,i - 1
5855
w1 = w1 + ( rkf_csu(i - 1,j)/rkf_csl(i - 1) )*h*s(j)
5956
end do
60-
!print *, w1%p(1)%m
6157
s(i) = fdt( rkf_ct(i) * h,w1 )
6258
end do
63-
6459
w1 = w
6560
do j=1,6
6661
w1 = w1 + ( rkf_cw(j) )*h*s(j)
6762
end do
68-
6963
z1 = w
7064
do j=1,6
7165
z1 = z1 + ( rkf_cz(j) )*h*s(j)
7266
end do
73-
7467
err = modulus(z1-w1)/modulus(z1)
7568
w = z1
76-
7769
end subroutine rkf45_step
7870

7971
!--------------------------------------------------------
@@ -102,7 +94,7 @@ subroutine rkf45(h,t,delta_t,state,err,tol)
10294
call rkf45_step(tt,dt,tstate,terr)
10395

10496
do ! main loop
105-
dt = 0.8*( (tol/terr)**(1/dble(p+1)) )*dt
97+
dt = 0.84q0*( (tol/terr)**(1/dble(p)) )*dt
10698
its = 0
10799
do ! subloop
108100
tt = ct

‎solver_c.f90

+8-8
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
! Deprecated
12
module glue
23
use tri_body_prob
34
use solver
@@ -9,14 +10,15 @@ subroutine init(h, t, delta_t, err, tol, p_num, p_m, p_x, p_v) bind(c, name='ini
910
integer(c_int), intent(in) :: p_num
1011
real(c_double), intent(in) :: h, t, delta_t, err, tol, p_m(p_num), p_x(p_num, 3), p_v(p_num, 3)
1112
type(planet) :: plist(p_num)
12-
integer :: i,j
13+
integer :: i
1314
gh = h
15+
!print *, h, t, delta_t, err, tol, p_num
1416
gt = t
1517
gdt = delta_t
1618
gerr = err
1719
gtol = tol
1820
do i=1, p_num
19-
!print *, i
21+
!print *, 'No.', i
2022
plist(i)%m = p_m(i)
2123
!print *, p_m(i)
2224
plist(i)%x(:) = p_x(i,:)
@@ -26,22 +28,20 @@ subroutine init(h, t, delta_t, err, tol, p_num, p_m, p_x, p_v) bind(c, name='ini
2628
end do
2729
state = triBody(plist)
2830
end
29-
subroutine iter_step(h, t, delta_t, err, tol, p_num, p_m, p_x, p_v) bind(c, name='iter_step')
31+
subroutine iter_step(t, err, energy, p_num, p_x, p_v) bind(c, name='iter_step')
3032
use, intrinsic :: ISO_C_BINDING
3133
integer(c_int), intent(in) :: p_num
32-
real(c_double), intent(inout) :: h, t, delta_t, err, tol, p_m(p_num), p_x(p_num, 3), p_v(p_num, 3)
34+
real(c_double), intent(out) :: t, err, energy, p_x(p_num, 3), p_v(p_num, 3)
3335
type(planet) :: plist(p_num)
34-
integer :: i,j
36+
integer :: i
3537
call rkf45(gh,gt,gdt,state,gerr,gtol)
3638
do i=1, p_num
3739
!p_m(i) = state%p(i)%m
3840
p_x(i,:) = state%p(i)%x(:)
3941
p_v(i,:) = state%p(i)%v(:)
4042
end do
41-
h = gh
4243
t = gt
43-
delta_t = gdt
4444
err = gerr
45-
tol = gtol
45+
energy = state%energy()
4646
end subroutine iter_step
4747
end module

‎solver_c_new.f90

+79
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
module new_glue
2+
use tri_body_prob
3+
use solver_new
4+
implicit none
5+
type(triBody) :: gstate
6+
real*16 :: gt, gdt, gerr, gtol
7+
contains
8+
subroutine init(t, delta_t, tol, p_num, p_m, p_x, p_v) bind(c, name='init')
9+
use, intrinsic :: ISO_C_BINDING
10+
integer(c_int), intent(in) :: p_num
11+
real(c_double), intent(in) :: t, delta_t, tol, p_m(p_num), p_x(p_num, 3), p_v(p_num, 3)
12+
type(planet) :: plist(p_num)
13+
integer :: i
14+
!print *, t, delta_t, tol, p_num
15+
gt = t
16+
gdt = delta_t
17+
gtol = tol
18+
do i=1, p_num
19+
!print *, 'No.', i
20+
plist(i)%m = p_m(i)
21+
!print *, p_m(i)
22+
plist(i)%x(:) = p_x(i,:)
23+
!print *, p_x(i,:)
24+
plist(i)%v(:) = p_v(i,:)
25+
!print *, p_v(i,:)
26+
end do
27+
gstate = triBody(plist)
28+
end
29+
subroutine iter_step(t, err, energy, p_num, p_x, p_v) bind(c, name='iter_step')
30+
use, intrinsic :: ISO_C_BINDING
31+
integer(c_int), intent(in) :: p_num
32+
real(c_double), intent(out) :: t, err, energy, p_x(p_num, 3), p_v(p_num, 3)
33+
type(planet) :: plist(p_num)
34+
integer :: i
35+
gerr = 0.q0
36+
!print *, gt, gdt, gerr, gtol
37+
call rkf45(gt,gdt,gstate,gerr,gtol)
38+
39+
do i=1, p_num
40+
!p_m(i) = state%p(i)%m
41+
p_x(i,:) = gstate%p(i)%x(:)
42+
p_v(i,:) = gstate%p(i)%v(:)
43+
end do
44+
t = gt
45+
err = gerr
46+
energy = gstate%energy()
47+
end subroutine iter_step
48+
subroutine iter_step_new(t, err, energy, p_num, p_x, p_v) bind(c, name='iter_step_new')
49+
use, intrinsic :: ISO_C_BINDING
50+
integer(c_int), intent(in) :: p_num
51+
real(c_double), intent(out) :: t, err, energy, p_x(p_num, 3), p_v(p_num, 3)
52+
type(planet) :: plist(p_num)
53+
integer :: i
54+
real*16 :: t_r, t_l, err_tmp
55+
type(triBody) :: state, state_r
56+
gerr = 0.q0
57+
!print *, gt, gdt, gerr, gtol
58+
t_l = gt
59+
t_r = t_l
60+
gt = gt + gdt
61+
state_r = gstate
62+
call rkf45_iter(t_r, gdt, state_r, gerr, gtol)
63+
do while(t_r < gt)
64+
t_l = t_r
65+
gstate = state_r
66+
call rkf45_iter(t_r, gt - t_l, state_r, err_tmp, gtol)
67+
gerr = gerr + err_tmp
68+
end do
69+
state = gstate + ((gt - t_l) / (t_r - t_l)) * (state_r - gstate)
70+
do i=1, p_num
71+
!p_m(i) = state%p(i)%m
72+
p_x(i,:) = state%p(i)%x(:)
73+
p_v(i,:) = state%p(i)%v(:)
74+
end do
75+
t = gt
76+
err = gerr
77+
energy = state%energy()
78+
end subroutine
79+
end module new_glue

‎solver_new.f90

+144
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
module solver_new
2+
use tri_body_prob
3+
implicit none
4+
real*16, parameter :: min_adj_dt = 1.Q-9
5+
real*16, parameter, dimension(6) :: rkf_tki = &
6+
[0.q0, 1.q0 / 4.q0, 3.q0 / 8.q0, 12.q0 / 13.q0, 1.q0, 1.q0 / 2.q0]
7+
real*16, parameter, dimension(6, 6) :: rkf_yki = transpose( reshape( [&
8+
0.q0, 0.q0, 0.q0, 0.q0, 0.q0, 0.q0, &
9+
1.q0 / 4.q0, 0.q0, 0.q0, 0.q0, 0.q0, 0.q0, &
10+
3.q0 / 32.q0, 9.q0 / 32.q0, 0.q0, 0.q0, 0.q0, 0.q0, &
11+
1932.q0 / 2197.q0, -7200.q0 / 2197.q0, 7296.q0 / 2197.q0, 0.q0, 0.q0, 0.q0, &
12+
439.q0 / 216.q0, -8.q0, 3680.q0 / 513.q0, -845.q0 / 4104.q0, 0.q0, 0.q0, &
13+
-8.q0 / 27.q0, 2.q0, -3544.q0 / 2565.q0, 1859.q0 / 4104.q0, -11.q0 / -40.q0, 0.q0 &
14+
], [6, 6]))
15+
real*16, parameter, dimension(6) :: rkf_yk = &
16+
[25.q0 / 216.q0, 0.q0, 1408.q0 / 2565.q0, 2197.q0 / 4101.q0, -1.q0 / 5.q0, 0.q0]
17+
real*16, parameter, dimension(6) :: rkf_zk = &
18+
[16.q0 / 135.q0, 0.q0, 6656.q0 / 12825.q0, 28561.q0 / 56430.q0, -9.q0 / 50.q0, 2.q0 / 55.q0]
19+
contains
20+
subroutine rk5(h, w, err)
21+
real*16, intent(in) :: h
22+
type(triBody) :: w, yl
23+
real*16 :: err
24+
type(triBody), dimension(6) :: k
25+
integer :: i, j
26+
27+
do i=1, 6
28+
k(i) = w
29+
do j=1, i-1
30+
k(i) = k(i) + rkf_yki(i, j) * h * k(j)
31+
end do
32+
k(i) = fdt(rkf_tki(i) * h, k(i))
33+
!call k(i)%show_triBody()
34+
end do
35+
yl = w
36+
do i=1, 5
37+
yl = yl + rkf_yk(i) * h * k(i)
38+
end do
39+
do i=1, 6
40+
w = w + rkf_zk(i) * h * k(i)
41+
end do
42+
err = modulus(w - yl) / modulus(w)
43+
!call w%show_triBody()
44+
!print *, 'err in rk5', err
45+
end
46+
recursive subroutine rkf45(t, delta_t, state, err, tol)
47+
real*16 :: adjusted_dt, err_tmp
48+
real*16, intent(inout) :: t, err
49+
real*16, intent(in) :: delta_t, tol
50+
type(triBody) :: state, test_state
51+
test_state = state
52+
call rk5(delta_t, test_state, err_tmp)
53+
if ((isnan(err_tmp) .eqv. .true.) .or. (err_tmp < tol)) then
54+
!print *, 'passed on first try'
55+
state = test_state
56+
t = t + delta_t
57+
err = err + err_tmp
58+
else
59+
!print *, 'F'
60+
adjusted_dt = delta_t * 0.84q0 * (tol / err_tmp)**(0.25q0)
61+
!print *, adjusted_dt / delta_t * 100.q0, '%'
62+
if (adjusted_dt < min_adj_dt) then
63+
!print *, 'adjusted length too short'
64+
call rk5(adjusted_dt, state, err_tmp)
65+
t = t + adjusted_dt
66+
err = err + err_tmp
67+
else
68+
!print *, 'L'
69+
call rkf45(t, adjusted_dt, state, err, tol)
70+
end if
71+
if (min_adj_dt > delta_t - adjusted_dt) then
72+
!print *, 'length left too short'
73+
call rk5(delta_t - adjusted_dt, state, err_tmp)
74+
t = t + delta_t - adjusted_dt
75+
err = err + err_tmp
76+
else
77+
!print *, 'R'
78+
call rkf45(t, delta_t - adjusted_dt, state, err, tol)
79+
end if
80+
end if
81+
!print *, 'EOF'
82+
end
83+
subroutine rkf45_loop(t, delta_t, state, err, tol)
84+
real*16 :: adjusted_dt, err_tmp, t_f
85+
real*16, intent(inout) :: t, err
86+
real*16, intent(in) :: delta_t, tol
87+
type(triBody) :: state, test_state
88+
test_state = state
89+
t_f = t + delta_t
90+
call rk5(delta_t, test_state, err_tmp)
91+
if ((isnan(err_tmp) .eqv. .true.) .or. (err_tmp < tol)) then
92+
!print *, 'passed on first try'
93+
state = test_state
94+
t = t + delta_t
95+
err = err + err_tmp
96+
return
97+
end if
98+
adjusted_dt = delta_t * 0.84q0 * (tol / err_tmp)**(0.25q0)
99+
do while (t + adjusted_dt < t_f .and. adjusted_dt > min_adj_dt)
100+
call rk5(adjusted_dt, test_state, err_tmp)
101+
t = t + adjusted_dt
102+
err = err + err_tmp
103+
adjusted_dt = adjusted_dt * 0.84q0 * (tol / err_tmp)**(0.25q0)
104+
end do
105+
call rk5(t_f - t, test_state, err_tmp)
106+
state = test_state
107+
err = err + err_tmp
108+
t = t_f
109+
end
110+
subroutine rkf45_iter(t, delta_t, state, err, tol)
111+
real*16 :: adjusted_dt
112+
real*16, intent(inout) :: t, err
113+
real*16, intent(in) :: delta_t, tol
114+
type(triBody) :: state, test_state
115+
test_state = state
116+
call rk5(delta_t, test_state, err)
117+
adjusted_dt = delta_t
118+
do while ((isnan(err) .eqv. .false.) .and. (err > tol) .and. (adjusted_dt > min_adj_dt))
119+
!print *, 'passed on first try'
120+
call rk5(adjusted_dt, test_state, err)
121+
adjusted_dt = adjusted_dt * 0.84q0 * (tol / err)**(0.25q0)
122+
end do
123+
state = test_state
124+
t = t + adjusted_dt
125+
end
126+
subroutine rkf45_interpolation(t, delta_t, state, err, tol)
127+
real*16, intent(inout) :: t, err
128+
real*16, intent(in) :: delta_t, tol
129+
real*16 :: t_r, t_l, err_tmp
130+
type(triBody) :: state, state_r
131+
t_l = t
132+
t_r = t_l
133+
t = t + delta_t
134+
state_r = state
135+
call rkf45_iter(t_r, delta_t, state_r, err, tol)
136+
do while(t_r < t)
137+
t_l = t_r
138+
state = state_r
139+
call rkf45_iter(t_r, t - t_l, state_r, err_tmp, tol)
140+
err = err + err_tmp
141+
end do
142+
state = state + ((t - t_l) / (t_r - t_l)) * (state_r - state)
143+
end
144+
end module solver_new

‎tri_body_prob.f90

+5-8
Original file line numberDiff line numberDiff line change
@@ -182,18 +182,15 @@ end function crush
182182
function modulus( triObj )
183183
class(triBody) :: triObj
184184
real*16 :: modulus
185-
real*16 :: m_m_square, m_x_square, m_v_square
186185
integer :: i, n
187186
n = triObj%p_num
188-
m_m_square = 0.q0
189-
m_x_square = 0.q0
190-
m_v_square = 0.q0
187+
modulus = 0.q0
191188
do i=1, n
192-
m_m_square = m_m_square + metric(1)*((triObj%p(i)%m)**2)
193-
m_x_square = m_x_square + metric(2)*(dot_product( triObj%p(i)%x, triObj%p(i)%x ))
194-
m_v_square = m_v_square + metric(3)*(dot_product( triObj%p(i)%v, triObj%p(i)%v ))
189+
modulus = modulus + metric(2)*(dot_product( triObj%p(i)%x, triObj%p(i)%x ))&
190+
+ metric(3)*(dot_product( triObj%p(i)%v, triObj%p(i)%v ))
191+
! + metric(1)*((triObj%p(i)%m)**2)
195192
end do
196-
modulus = sqrt( m_m_square + m_x_square + m_v_square )
193+
modulus = sqrt( modulus )
197194
end function modulus
198195

199196
subroutine show_triBody( triObj )

‎tribody.py

+5
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,8 @@
1+
"""
2+
Deprecated
3+
"""
4+
5+
16
from ctypes import *
27
import numpy as np
38
import time

0 commit comments

Comments
 (0)