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 Pi-line model of a transmission line as three linear ODEs, IMEX with the constant source term explicit. 

12 

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

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

15 

16 .. math:: 

17 \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}, 

18 

19 .. math:: 

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

21 

22 .. math:: 

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

24 

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

26 

27 .. math:: 

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

29 

30 using an initial condition. 

31 

32 Parameters 

33 ---------- 

34 Vs : float, optional 

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

36 Rs : float, optional 

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

38 C1 : float, optional 

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

40 Rpi : float, optional 

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

42 Lpi : float, optional 

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

44 C2 : float, optional 

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

46 Rl : float, optional 

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

48 

49 Attributes 

50 ---------- 

51 A : np.2darray 

52 Coefficient matrix of the linear ODE system. 

53 """ 

54 

55 dtype_u = mesh 

56 dtype_f = imex_mesh 

57 

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

59 """Initialization routine""" 

60 

61 nvars = 3 

62 # invoke super init, passing number of dofs 

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

64 self._makeAttributeAndRegister( 

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

66 ) 

67 

68 # compute dx and get discretization matrix A 

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

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

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

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

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

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

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

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

77 

78 def eval_f(self, u, t): 

79 """ 

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

81 

82 Parameters 

83 ---------- 

84 u : dtype_u 

85 Current values of the numerical solution. 

86 t : float 

87 Current time of the numerical solution is computed. 

88 

89 Returns 

90 ------- 

91 f : dtype_f 

92 The right-hand side of the problem. 

93 """ 

94 

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

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

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

98 return f 

99 

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

101 r""" 

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

103 

104 Parameters 

105 ---------- 

106 rhs : dtype_f 

107 Right-hand side for the linear system. 

108 factor : float 

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

110 u0 : dtype_u 

111 Initial guess for the iterative solver. 

112 t : float 

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

114 

115 Returns 

116 ------- 

117 me : dtype_u 

118 The solution as mesh. 

119 """ 

120 

121 me = self.dtype_u(self.init) 

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

123 return me 

124 

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

126 r""" 

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

128 

129 Parameters 

130 ---------- 

131 t : float 

132 Time of the exact solution. 

133 u_init : pySDC.problem.Piline.dtype_u 

134 Initial conditions for getting the exact solution. 

135 t_init : float 

136 The starting time. 

137 

138 Returns 

139 ------- 

140 me : dtype_u 

141 The reference solution. 

142 """ 

143 

144 me = self.dtype_u(self.init) 

145 

146 # fill initial conditions 

147 me[0] = 0.0 # v1 

148 me[1] = 0.0 # v2 

149 me[2] = 0.0 # p3 

150 

151 if t > 0.0: 

152 if u_init is not None: 

153 if t_init is None: 

154 raise ValueError('Please supply `t_init` when you want to get the exact solution from a point that \ 

155is not 0!') 

156 me = u_init 

157 else: 

158 t_init = 0.0 

159 

160 def rhs(t, u): 

161 f = self.eval_f(u, t) 

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

163 

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

165 

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

167 

168 return me