Coverage for pySDC/implementations/problem_classes/OuterSolarSystem.py: 100%

1import numpy as np 

2 

3from pySDC.core.problem import Problem 

4from pySDC.implementations.datatype_classes.particles import particles, acceleration 

5 

6 

7# noinspection PyUnusedLocal 

8class outer_solar_system(Problem): 

9 r""" 

10 Gravitational N-body problem of the outer solar system: the sun, Jupiter, Saturn, Uranus, Neptune and Pluto. 

11 

12 The :math:`N`-body problem describes the mutual influence of the motion of :math:`N` bodies. Formulation of the problem is 

13 based on Newton's second law. Therefore, the :math:`N`-body problem is formulated as 

14 

15 .. math:: 

16 m_i \frac{d^2 {\bf r}_i}{d t^2} = \sum_{j=1, i\neq j}^N G \frac{m_i m_j}{|{\bf r}_i - {\bf r}_j|^3}({\bf r}_j - {\bf r}_i), 

17 

18 where :math:`m_i` is the :math:`i`-th mass point with position described by the vector :math:`{\bf r}_i`, and :math:`G` 

19 is the gravitational constant. If only the sun (body 1) influences the motion of the bodies gravitationally, the 

20 equations become 

21 

22 .. math:: 

23 \frac{d^2 {\bf r}_i}{d t^2} = G \frac{m_1}{|{\bf r}_i - {\bf r}_1|^3}({\bf r}_1 - {\bf r}_i), \quad i \geq 2, 

24 

25 and the sun is not accelerated. 

26 

27 This class implements the outer solar system consisting of the sun and five outer bodies: Jupiter, Saturn, Uranus, 

28 Neptune, and Pluto, i.e., :math:`N=6`. The masses of the inner planets are added to the sun. 

29 

30 Parameters 

31 ---------- 

32 sun_only : bool, optional 

33 If True, only the gravity of the sun acts on the other bodies, see above. If False (default), all bodies 

34 attract each other. 

35 

36 Attributes 

37 ---------- 

38 G : float 

39 Gravitational constant. 

40 """ 

41 

42 dtype_u = particles 

43 dtype_f = acceleration 

44 

45 G = 2.95912208286e-4 

46 

47 def __init__(self, sun_only=False): 

48 """Initialization routine""" 

49 

50 # invoke super init, passing nparts 

51 super().__init__(((3, 6), None, np.dtype('float64'))) 

52 self._makeAttributeAndRegister('sun_only', localVars=locals()) 

53 

54 def eval_f(self, u, t): 

55 """ 

56 Routine to compute the right-hand side of the problem. 

57 

58 Parameters 

59 ---------- 

60 u : dtype_u 

61 The particles. 

62 t : float 

63 Current time at which the particles are computed (not used here). 

64 

65 Returns 

66 ------- 

67 me : dtype_f 

68 The right-hand side of the problem. 

69 """ 

70 me = self.dtype_f(self.init, val=0.0) 

71 

72 # compute the acceleration due to gravitational forces 

73 # ... only with respect to the sun 

74 if self.sun_only: 

75 for i in range(1, self.init[0][-1]): 

76 dx = u.pos[:, i] - u.pos[:, 0] 

77 r = np.sqrt(np.dot(dx, dx)) 

78 df = self.G * dx / (r**3) 

79 me[:, i] -= u.m[0] * df 

80 

81 # ... or with all planets involved 

82 else: 

83 for i in range(self.init[0][-1]): 

84 for j in range(i): 

85 dx = u.pos[:, i] - u.pos[:, j] 

86 r = np.sqrt(np.dot(dx, dx)) 

87 df = self.G * dx / (r**3) 

88 me[:, i] -= u.m[j] * df 

89 me[:, j] += u.m[i] * df 

90 

91 return me 

92 

93 def u_exact(self, t): 

94 r""" 

95 Routine to compute the exact/initial trajectory at time :math:`t`. 

96 

97 Parameters 

98 ---------- 

99 t : float 

100 Time of the exact/initial trajectory. 

101 

102 Returns 

103 ------- 

104 me : dtype_u 

105 The exact/initial position and velocity. 

106 """ 

107 assert t == 0.0, 'error, u_exact only works for the initial time t0=0' 

108 me = self.dtype_u(self.init) 

109 

110 me.pos[:, 0] = [0.0, 0.0, 0.0] 

111 me.pos[:, 1] = [-3.5025653, -3.8169847, -1.5507963] 

112 me.pos[:, 2] = [9.0755314, -3.0458353, -1.6483708] 

113 me.pos[:, 3] = [8.3101420, -16.2901086, -7.2521278] 

114 me.pos[:, 4] = [11.4707666, -25.7294829, -10.8169456] 

115 me.pos[:, 5] = [-15.5387357, -25.2225594, -3.1902382] 

116 

117 me.vel[:, 0] = [0.0, 0.0, 0.0] 

118 me.vel[:, 1] = [0.00565429, -0.00412490, -0.00190589] 

119 me.vel[:, 2] = [0.00168318, 0.00483525, 0.00192462] 

120 me.vel[:, 3] = [0.00354178, 0.00137102, 0.00055029] 

121 me.vel[:, 4] = [0.00288930, 0.00114527, 0.00039677] 

122 me.vel[:, 5] = [0.00276725, -0.0017072, -0.00136504] 

123 

124 me.m[0] = 1.00000597682 # Sun 

125 me.m[1] = 0.000954786104043 # Jupiter 

126 me.m[2] = 0.000285583733151 # Saturn 

127 me.m[3] = 0.0000437273164546 # Uranus 

128 me.m[4] = 0.0000517759138449 # Neptune 

129 me.m[5] = 1.0 / 130000000.0 # Pluto 

130 

131 return me 

132 

133 def eval_hamiltonian(self, u): 

134 """ 

135 Routine to compute the Hamiltonian. 

136 

137 Parameters 

138 ---------- 

139 u : dtype_u 

140 The particles. 

141 

142 Returns 

143 ------- 

144 ham : float 

145 The Hamiltonian. 

146 """ 

147 

148 ham = 0 

149 

150 for i in range(self.init[0][-1]): 

151 ham += 0.5 * u.m[i] * np.dot(u.vel[:, i], u.vel[:, i]) 

152 

153 for i in range(self.init[0][-1]): 

154 for j in range(i): 

155 dx = u.pos[:, i] - u.pos[:, j] 

156 r = np.sqrt(np.dot(dx, dx)) 

157 ham -= self.G * u.m[i] * u.m[j] / r 

158 

159 return ham