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

1import numpy as np 

2from scipy.integrate import solve_ivp 

3 

4from pySDC.core.problem import Problem 

5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh 

6 

7 

8# noinspection PyUnusedLocal 

9class piline(Problem): 

10 r""" 

11 Example implementing the model of the piline. It serves as a transmission line in an energy grid. The problem of simulating the 

12 piline consists of three ordinary differential equations (ODEs) with nonhomogeneous part: 

13 

14 .. math:: 

15 \frac{d v_{C_1} (t)}{dt} = -\frac{1}{R_s C_1}v_{C_1} (t) - \frac{1}{C_1} i_{L_\pi} (t) + \frac{V_s}{R_s C_1}, 

16 

17 .. math:: 

18 \frac{d v_{C_2} (t)}{dt} = -\frac{1}{R_\ell C_2}v_{C_2} (t) + \frac{1}{C_2} i_{L_\pi} (t), 

19 

20 .. math:: 

21 \frac{d i_{L_\pi} (t)}{dt} = \frac{1}{L_\pi} v_{C_1} (t) - \frac{1}{L_\pi} v_{C_2} (t) - \frac{R_\pi}{L_\pi} i_{L_\pi} (t), 

22 

23 which can be expressed as a nonhomogeneous linear system of ODEs 

24 

25 .. math:: 

26 \frac{d u(t)}{dt} = A u(t) + f(t) 

27 

28 using an initial condition. 

29 

30 Parameters 

31 ---------- 

32 Vs : float, optional 

33 Voltage at the voltage source :math:`V_s`. 

34 Rs : float, optional 

35 Resistance of the resistor :math:`R_s` at the voltage source. 

36 C1 : float, optional 

37 Capacitance of the capacitor :math:`C_1`. 

38 Rpi : float, optional 

39 Resistance of the resistor :math:`R_\pi`. 

40 Lpi : float, optional 

41 Inductance of the inductor :math:`L_\pi`. 

42 C2 : float, optional 

43 Capacitance of the capacitor :math:`C_2`. 

44 Rl : float, optional 

45 Resistance of the resistive load :math:`R_\ell`. 

46 

47 Attributes: 

48 A : np.2darray 

49 Coefficient matrix of the linear ODE system. 

50 """ 

51 

52 dtype_u = mesh 

53 dtype_f = imex_mesh 

54 

55 def __init__(self, Vs=100.0, Rs=1.0, C1=1.0, Rpi=0.2, Lpi=1.0, C2=1.0, Rl=5.0): 

56 """Initialization routine""" 

57 

58 nvars = 3 

59 # invoke super init, passing number of dofs 

60 super().__init__(init=(nvars, None, np.dtype('float64'))) 

61 self._makeAttributeAndRegister( 

62 'nvars', 'Vs', 'Rs', 'C1', 'Rpi', 'Lpi', 'C2', 'Rl', localVars=locals(), readOnly=True 

63 ) 

64 

65 # compute dx and get discretization matrix A 

66 self.A = np.zeros((3, 3)) 

67 self.A[0, 0] = -1 / (self.Rs * self.C1) 

68 self.A[0, 2] = -1 / self.C1 

69 self.A[1, 1] = -1 / (self.Rl * self.C2) 

70 self.A[1, 2] = 1 / self.C2 

71 self.A[2, 0] = 1 / self.Lpi 

72 self.A[2, 1] = -1 / self.Lpi 

73 self.A[2, 2] = -self.Rpi / self.Lpi 

74 

75 def eval_f(self, u, t): 

76 """ 

77 Routine to evaluate the right-hand side of the problem. 

78 

79 Parameters 

80 ---------- 

81 u : dtype_u 

82 Current values of the numerical solution. 

83 t : float 

84 Current time of the numerical solution is computed. 

85 

86 Returns 

87 ------- 

88 f : dtype_f 

89 The right-hand side of the problem. 

90 """ 

91 

92 f = self.dtype_f(self.init, val=0.0) 

93 f.impl[:] = self.A.dot(u) 

94 f.expl[0] = self.Vs / (self.Rs * self.C1) 

95 return f 

96 

97 def solve_system(self, rhs, factor, u0, t): 

98 r""" 

99 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`. 

100 

101 Parameters 

102 ---------- 

103 rhs : dtype_f 

104 Right-hand side for the linear system. 

105 factor : float 

106 Abbrev. for the local stepsize (or any other factor required). 

107 u0 : dtype_u 

108 Initial guess for the iterative solver. 

109 t : float 

110 Current time (e.g. for time-dependent BCs). 

111 

112 Returns 

113 ------- 

114 me : dtype_u 

115 The solution as mesh. 

116 """ 

117 

118 me = self.dtype_u(self.init) 

119 me[:] = np.linalg.solve(np.eye(self.nvars) - factor * self.A, rhs) 

120 return me 

121 

122 def u_exact(self, t, u_init=None, t_init=None): 

123 r""" 

124 Routine to approximate the exact solution at time :math:`t` by ``SciPy`` as a reference. 

125 

126 Parameters 

127 ---------- 

128 t : float 

129 Time of the exact solution. 

130 u_init : pySDC.problem.Piline.dtype_u 

131 Initial conditions for getting the exact solution. 

132 t_init : float 

133 The starting time. 

134 

135 Returns 

136 ------- 

137 me : dtype_u 

138 The reference solution. 

139 """ 

140 

141 me = self.dtype_u(self.init) 

142 

143 # fill initial conditions 

144 me[0] = 0.0 # v1 

145 me[1] = 0.0 # v2 

146 me[2] = 0.0 # p3 

147 

148 if t > 0.0: 

149 if u_init is not None: 

150 if t_init is None: 

151 raise ValueError( 

152 'Please supply `t_init` when you want to get the exact solution from a point that \ 

153is not 0!' 

154 ) 

155 me = u_init 

156 else: 

157 t_init = 0.0 

158 

159 def rhs(t, u): 

160 f = self.eval_f(u, t) 

161 return f.impl + f.expl # evaluate only explicitly rather than IMEX 

162 

163 tol = 100 * np.finfo(float).eps 

164 

165 me[:] = solve_ivp(rhs, (t_init, t), me, rtol=tol, atol=tol).y[:, -1] 

166 

167 return me