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

## 44 statements

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

1import numpy as np

2from scipy.integrate import solve_ivp

4from pySDC.core.problem import Problem

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

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:

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

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

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

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

25 .. math::

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

28 using an initial condition.

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.

47 Attributes:

48 A : np.2darray

49 Coefficient matrix of the linear ODE system.

50 """

52 dtype_u = mesh

53 dtype_f = imex_mesh

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

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 )

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

75 def eval_f(self, u, t):

76 """

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

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.

86 Returns

87 -------

88 f : dtype_f

89 The right-hand side of the problem.

90 """

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

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

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

112 Returns

113 -------

114 me : dtype_u

115 The solution as mesh.

116 """

118 me = self.dtype_u(self.init)

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

120 return me

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.

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.

135 Returns

136 -------

137 me : dtype_u

138 The reference solution.

139 """

141 me = self.dtype_u(self.init)

143 # fill initial conditions

144 me[0] = 0.0 # v1

145 me[1] = 0.0 # v2

146 me[2] = 0.0 # p3

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

159 def rhs(t, u):

160 f = self.eval_f(u, t)

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

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

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

167 return me