Coverage for pySDC/implementations/problem_classes/FermiPastaUlamTsingou.py: 98%

40 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

3from pySDC.core.errors import ParameterError

4from pySDC.core.problem import Problem

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

8# noinspection PyUnusedLocal

9class fermi_pasta_ulam_tsingou(Problem):

10 r"""

11 The Fermi-Pasta-Ulam-Tsingou (FPUT) problem was one of the first computer experiments. E. Fermi, J. Pasta

12 and S. Ulam investigated the behavior of a vibrating spring with a weak correction term (which is quadratic

13 for the FPU-:math:\alpha model, and cubic for the FPU-:math:\beta model [1]_). This can be modelled by

14 the second-order problem

16 .. math::

17 \frac{d^2 u_j(t)}{d t^2} = (u_{j+1}(t) - 2 u_j(t) + u_{j-1}(t)) (1 + \alpha (u_{j+1}(t) - u_{j-1}(t))),

19 where :math:u_j(t) is the position of the :math:j-th particle. [2]_ is used as setup for this

20 implemented problem class. The Hamiltonian of this problem (needed for second-order SDC) is

22 .. math::

23 \sum_{i=1}^n \frac{1}{2}v^2_{i-1}(t) + \frac{1}{2}(u_{i+1}(t) - u_{i-1}(t))^2 + \frac{\alpha}{3}(u_{i+1}(t) - u_{i-1}(t))^3,

25 where :math:v_j(t) is the velocity of the :math:j-th particle.

27 Parameters

28 ----------

29 npart : int, optional

30 Number of particles.

31 alpha : float, optional

32 Factor of the nonlinear force :math:\alpha.

33 k : float, optional

34 Mode for initial conditions.

35 energy_modes : list, optional

36 Energy modes.

38 Attributes

39 ----------

40 dx : float

41 Mesh grid size.

42 xvalues : np.1darray

43 Spatial grid.

44 ones : np.1darray

45 Vector containing ones.

47 References

48 ----------

49 .. [1] E. Fermi, J. Pasta, S. Ulam. Studies of nonlinear problems (1955). I. Los Alamos report LA-1940.

50 Collected Papers of Enrico Fermi, E. Segré (Ed.), University of Chicago Press (1965)

51 .. [2] http://www.scholarpedia.org/article/Fermi-Pasta-Ulam_nonlinear_lattice_oscillations

52 """

54 dtype_u = particles

55 dtype_f = acceleration

57 def __init__(self, npart=2048, alpha=0.25, k=1.0, energy_modes=None):

58 """Initialization routine"""

60 if energy_modes is None:

61 energy_modes = [1, 2, 3, 4]

63 # invoke super init, passing nparts

64 super().__init__((npart, None, np.dtype('float64')))

65 self._makeAttributeAndRegister('npart', 'alpha', 'k', 'energy_modes', localVars=locals(), readOnly=True)

67 self.dx = (self.npart / 32) / (self.npart + 1)

68 self.xvalues = np.array([(i + 1) * self.dx for i in range(self.npart)])

69 self.ones = np.ones(self.npart - 2)

71 def eval_f(self, u, t):

72 """

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

75 Parameters

76 ----------

77 u : dtype_u

78 Current values of the numerical solution.

79 t : float

80 Current time of the numerical solution is computed.

82 Returns

83 -------

84 f : dtype_f

85 The right-hand side of the problem.

86 """

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

89 # me[1:-1] = u.pos[:-2] - 2.0 * u.pos[1:-1] + u.pos[2:] + \

90 # self.alpha * ((u.pos[2:] - u.pos[1:-1]) ** 2 -

91 # (u.pos[1:-1] - u.pos[:-2]) ** 2)

92 # me[0] = -2.0 * u.pos[0] + u.pos[1] + \

93 # self.alpha * ((u.pos[1] - u.pos[0]) ** 2 - (u.pos[0]) ** 2)

94 # me[-1] = u.pos[-2] - 2.0 * u.pos[-1] + \

95 # self.alpha * ((u.pos[-1]) ** 2 - (u.pos[-1] - u.pos[-2]) ** 2)

96 me[1:-1] = (u.pos[:-2] - 2.0 * u.pos[1:-1] + u.pos[2:]) * (self.ones + self.alpha * (u.pos[2:] - u.pos[:-2]))

97 me[0] = (-2.0 * u.pos[0] + u.pos[1]) * (1 + self.alpha * (u.pos[1]))

98 me[-1] = (u.pos[-2] - 2.0 * u.pos[-1]) * (1 + self.alpha * (-u.pos[-2]))

100 return me

102 def u_exact(self, t):

103 r"""

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

106 Parameters

107 ----------

108 t : float

109 Time of the exact solution.

111 Returns

112 -------

113 me : dtype_u

114 The exact/initial position and velocity.

115 """

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

118 me = self.dtype_u(self.init, val=0.0)

120 me.pos[:] = np.sin(self.k * np.pi * self.xvalues)

121 me.vel[:] = 0.0

123 return me

125 def eval_hamiltonian(self, u):

126 """

127 Routine to compute the Hamiltonian.

129 Parameters

130 ----------

131 u : dtype_u

132 The particles.

134 Returns

135 -------

136 ham : float

137 The Hamiltonian.

138 """

140 ham = sum(

141 0.5 * u.vel[:-1] ** 2

142 + 0.5 * (u.pos[1:] - u.pos[:-1]) ** 2

143 + self.alpha / 3.0 * (u.pos[1:] - u.pos[:-1]) ** 3

144 )

145 ham += 0.5 * u.vel[-1] ** 2 + 0.5 * (u.pos[-1]) ** 2 + self.alpha / 3.0 * (-u.pos[-1]) ** 3

146 ham += 0.5 * (u.pos[0]) ** 2 + self.alpha / 3.0 * (u.pos[0]) ** 3

147 return ham

149 def eval_mode_energy(self, u):

150 r"""

151 Routine to compute the energy following [1]_.

153 Parameters

154 ----------

155 u : dtype_u

156 Particles.

158 Returns

159 -------

160 energy : dict

161 Energies.

162 """

164 energy = {}

166 for k in self.energy_modes:

167 # Qk = np.sqrt(2.0 / (self.npart + 1)) * np.dot(u.pos, np.sin(np.pi * k * self.xvalues))

168 Qk = np.sqrt(2.0 * self.dx) * np.dot(u.pos, np.sin(np.pi * k * self.xvalues))

169 # Qkdot = np.sqrt(2.0 / (self.npart + 1)) * np.dot(u.vel, np.sin(np.pi * k * self.xvalues))

170 Qkdot = np.sqrt(2.0 * self.dx) * np.dot(u.vel, np.sin(np.pi * k * self.xvalues))

172 # omegak2 = 4.0 * np.sin(k * np.pi / (2.0 * (self.npart + 1))) ** 2

173 omegak2 = 4.0 * np.sin(k * np.pi * self.dx / 2.0) ** 2

174 energy[k] = 0.5 * (Qkdot**2 + omegak2 * Qk**2)

176 return energy