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 The :math:`N`-body problem describes the mutual influence of the motion of :math:`N` bodies. Formulation of the problem is 

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

12 

13 .. math:: 

14 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}_i - {\bf r}_j), 

15 

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

17 is the gravitational constant. If only the sun influences the motion of the bodies gravitationally, the equations become 

18 

19 .. math:: 

20 m_i \frac{d^2 {\bf r}_i}{d t^2} = G \frac{m_1}{|{\bf r}_i - {\bf r}_1|^3}({\bf r}_i - {\bf r}_1). 

21 

22 This class implements the outer solar system consisting of the six outer planets: the sun, Jupiter, Saturn, Uranus, 

23 Neptune, and Pluto, i.e., :math:`N=6`. 

24 

25 Parameters 

26 ---------- 

27 sun_only : bool, optional 

28 If False, only the sun is taken into account for the influence of the motion. 

29 

30 Attributes 

31 ---------- 

32 G : float 

33 Gravitational constant. 

34 """ 

35 

36 dtype_u = particles 

37 dtype_f = acceleration 

38 

39 G = 2.95912208286e-4 

40 

41 def __init__(self, sun_only=False): 

42 """Initialization routine""" 

43 

44 # invoke super init, passing nparts 

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

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

47 

48 def eval_f(self, u, t): 

49 """ 

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

51 

52 Parameters 

53 ---------- 

54 u : dtype_u 

55 The particles. 

56 t (float): Current time at which the particles are computed (not used here). 

57 

58 Returns 

59 ------- 

60 me : dtype_f 

61 The right-hand side of the problem. 

62 """ 

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

64 

65 # compute the acceleration due to gravitational forces 

66 # ... only with respect to the sun 

67 if self.sun_only: 

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

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

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

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

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

73 

74 # ... or with all planets involved 

75 else: 

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

77 for j in range(i): 

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

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

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

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

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

83 

84 return me 

85 

86 def u_exact(self, t): 

87 r""" 

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

89 

90 Parameters 

91 ---------- 

92 t : float 

93 Time of the exact/initial trajectory. 

94 

95 Returns 

96 ------- 

97 me : dtype_u 

98 The exact/initial position and velocity. 

99 """ 

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

101 me = self.dtype_u(self.init) 

102 

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

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

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

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

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

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

109 

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

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

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

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

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

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

116 

117 me.m[0] = 1.00000597682 # Sun 

118 me.m[1] = 0.000954786104043 # Jupiter 

119 me.m[2] = 0.000285583733151 # Saturn 

120 me.m[3] = 0.0000437273164546 # Uranus 

121 me.m[4] = 0.0000517759138449 # Neptune 

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

123 

124 return me 

125 

126 def eval_hamiltonian(self, u): 

127 """ 

128 Routine to compute the Hamiltonian. 

129 

130 Parameters 

131 ---------- 

132 u : dtype_u 

133 The particles. 

134 

135 Returns 

136 ------- 

137 ham : float 

138 The Hamiltonian. 

139 """ 

140 

141 ham = 0 

142 

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

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

145 

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

147 for j in range(i): 

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

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

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

151 

152 return ham