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

## 58 statements

, created at 2024-09-20 17:10 +0000

1import numpy as np

3from pySDC.core.problem import Problem

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

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

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),

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

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).

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.

25 Parameters

26 ----------

27 sun_only : bool, optional

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

30 Attributes

31 ----------

32 G : float

33 Gravitational constant.

34 """

36 dtype_u = particles

37 dtype_f = acceleration

39 G = 2.95912208286e-4

41 def __init__(self, sun_only=False):

42 """Initialization routine"""

44 # invoke super init, passing nparts

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

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

48 def eval_f(self, u, t):

49 """

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

52 Parameters

53 ----------

54 u : dtype_u

55 The particles.

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

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)

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

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

84 return me

86 def u_exact(self, t):

87 r"""

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

90 Parameters

91 ----------

92 t : float

93 Time of the exact/initial trajectory.

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)

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]

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]

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

124 return me

126 def eval_hamiltonian(self, u):

127 """

128 Routine to compute the Hamiltonian.

130 Parameters

131 ----------

132 u : dtype_u

133 The particles.

135 Returns

136 -------

137 ham : float

138 The Hamiltonian.

139 """

141 ham = 0

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

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

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

152 return ham