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

1import numpy as np 

2 

3from pySDC.core.errors import ParameterError 

4from pySDC.core.problem import Problem 

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

6 

7 

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 

15 

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

18 

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 

21 

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, 

24 

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

26 

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. 

37 

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. 

46 

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 """ 

53 

54 dtype_u = particles 

55 dtype_f = acceleration 

56 

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

58 """Initialization routine""" 

59 

60 if energy_modes is None: 

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

62 

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) 

66 

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) 

70 

71 def eval_f(self, u, t): 

72 """ 

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

74 

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. 

81 

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) 

88 

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

99 

100 return me 

101 

102 def u_exact(self, t): 

103 r""" 

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

105 

106 Parameters 

107 ---------- 

108 t : float 

109 Time of the exact solution. 

110 

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' 

117 

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

119 

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

121 me.vel[:] = 0.0 

122 

123 return me 

124 

125 def eval_hamiltonian(self, u): 

126 """ 

127 Routine to compute the Hamiltonian. 

128 

129 Parameters 

130 ---------- 

131 u : dtype_u 

132 The particles. 

133 

134 Returns 

135 ------- 

136 ham : float 

137 The Hamiltonian. 

138 """ 

139 

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 

148 

149 def eval_mode_energy(self, u): 

150 r""" 

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

152 

153 Parameters 

154 ---------- 

155 u : dtype_u 

156 Particles. 

157 

158 Returns 

159 ------- 

160 energy : dict 

161 Energies. 

162 """ 

163 

164 energy = {} 

165 

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

171 

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) 

175 

176 return energy