Coverage for pyrc \ tests \ test_stationary.py: 100%

1# ------------------------------------------------------------------------------- 

2# Copyright (C) 2026 Joel Kimmich, Tim Jourdan 

3# ------------------------------------------------------------------------------ 

4# License 

5# This file is part of PyRC, distributed under GPL-3.0-or-later. 

6# ------------------------------------------------------------------------------ 

7 

8import unittest 

9import numpy as np 

10import sympy as sp 

11from scipy.sparse import csr_matrix 

12from pyrc.core.solver.stationary import solve_stationary 

13 

14 

15class TestSolveStationary(unittest.TestCase): 

16 def setUp(self): 

17 # 2-node network: A = [[-2, 1], [1, -2]], B = identity, u = [10, 20] 

18 # Stationary solution: x = -A⁻¹ u = [13.33, 16.67] 

19 self.system_matrix_dense = np.array([[-2.0, 1.0], [1.0, -2.0]]) 

20 self.input_matrix_dense = np.array([[1.0, 0.0], [0.0, 1.0]]) 

21 self.input_vector = np.array([10.0, 20.0]) 

22 self.expected = np.linalg.solve(-self.system_matrix_dense, self.input_matrix_dense @ self.input_vector) 

23 

24 self.resistance, self.capacity = sp.symbols("resistance capacity", positive=True) 

25 self.system_matrix_symbol = sp.Matrix( 

26 [ 

27 [-2 / (self.resistance * self.capacity), 1 / (self.resistance * self.capacity)], 

28 [1 / (self.resistance * self.capacity), -2 / (self.resistance * self.capacity)], 

29 ] 

30 ) 

31 

32 def test_numeric_sparse_returns_correct_temperatures(self): 

33 """Numeric sparse solve returns correct stationary temperatures for a 2-node network.""" 

34 result = solve_stationary( 

35 csr_matrix(self.system_matrix_dense), 

36 csr_matrix(self.input_matrix_dense), 

37 self.input_vector, 

38 ) 

39 np.testing.assert_allclose(result, self.expected, rtol=1e-10) 

40 

41 def test_numeric_result_is_one_dimensional(self): 

42 """Numeric solve returns a 1D array.""" 

43 result = solve_stationary( 

44 csr_matrix(self.system_matrix_dense), 

45 csr_matrix(self.input_matrix_dense), 

46 self.input_vector, 

47 ) 

48 self.assertEqual(result.ndim, 1) 

49 

50 def test_numeric_input_vector_as_column_vector(self): 

51 """Numeric solve handles input vector passed as a 2D column vector.""" 

52 result = solve_stationary( 

53 csr_matrix(self.system_matrix_dense), 

54 csr_matrix(self.input_matrix_dense), 

55 self.input_vector.reshape(-1, 1), 

56 ) 

57 np.testing.assert_allclose(result, self.expected, rtol=1e-10) 

58 

59 def test_numeric_homogeneous_system_returns_zero(self): 

60 """Homogeneous system with zero input vector returns zero temperatures.""" 

61 result = solve_stationary( 

62 csr_matrix(self.system_matrix_dense), 

63 csr_matrix(self.input_matrix_dense), 

64 np.zeros(2), 

65 ) 

66 np.testing.assert_allclose(result, np.zeros(2), atol=1e-15) 

67 

68 def test_numeric_non_identity_input_matrix(self): 

69 """Numeric solve handles a non-identity input matrix with a known analytical solution.""" 

70 system_matrix = csr_matrix(np.array([[-3.0, 1.0], [1.0, -3.0]])) 

71 input_matrix = csr_matrix(np.array([[2.0, 0.0], [0.0, 3.0]])) 

72 input_vector = np.array([5.0, 10.0]) 

73 expected = np.linalg.solve( 

74 np.array([[-3.0, 1.0], [1.0, -3.0]]), 

75 -(input_matrix.toarray() @ input_vector), 

76 ) 

77 result = solve_stationary(system_matrix, input_matrix, input_vector) 

78 np.testing.assert_allclose(result, expected, rtol=1e-10) 

79 

80 def test_numeric_single_node_network(self): 

81 """Single-node system returns the correct scalar stationary temperature.""" 

82 result = solve_stationary( 

83 csr_matrix(np.array([[-1.0]])), 

84 csr_matrix(np.array([[1.0]])), 

85 np.array([42.0]), 

86 ) 

87 np.testing.assert_allclose(result, [42.0], rtol=1e-10) 

88 

89 def test_numeric_large_tridiagonal_system(self): 

90 """Numeric solve handles a large tridiagonal sparse system correctly.""" 

91 number_of_nodes = 500 

92 diagonal = np.full(number_of_nodes, -2.0) 

93 off_diagonal = np.ones(number_of_nodes - 1) 

94 system_matrix_dense = np.diag(diagonal) + np.diag(off_diagonal, 1) + np.diag(off_diagonal, -1) 

95 input_vector = np.random.default_rng(0).uniform(0, 100, number_of_nodes) 

96 result = solve_stationary( 

97 csr_matrix(system_matrix_dense), 

98 csr_matrix(np.eye(number_of_nodes)), 

99 input_vector, 

100 ) 

101 np.testing.assert_allclose(result, np.linalg.solve(system_matrix_dense, -input_vector), rtol=1e-8) 

102 

103 def test_numeric_singular_system_matrix_raises(self): 

104 """Singular system matrix raises a ValueError in the numeric solver.""" 

105 singular_system_matrix = csr_matrix(np.array([[-1.0, 1.0], [1.0, -1.0]])) 

106 with self.assertRaises(ValueError): 

107 solve_stationary(singular_system_matrix, csr_matrix(self.input_matrix_dense), self.input_vector) 

108 

109 def test_sympy_matrix_without_free_symbols_uses_numeric_solver(self): 

110 """Sympy matrix with no free symbols dispatches to the numeric solver and returns an ndarray.""" 

111 result = solve_stationary( 

112 sp.Matrix(self.system_matrix_dense.tolist()), 

113 sp.Matrix(self.input_matrix_dense.tolist()), 

114 sp.Matrix(self.input_vector.tolist()), 

115 ) 

116 self.assertIsInstance(result, np.ndarray) 

117 np.testing.assert_allclose(result, self.expected, rtol=1e-10) 

118 

119 def test_symbolic_with_free_symbols_returns_sympy_matrix(self): 

120 """Symbolic solve returns a sympy Matrix when free symbols are present.""" 

121 result = solve_stationary(self.system_matrix_symbol, sp.Matrix([[1, 0], [0, 1]]), sp.Matrix([10, 20])) 

122 self.assertIsInstance(result, sp.MatrixBase) 

123 self.assertEqual(result.shape, (2, 1)) 

124 

125 def test_symbolic_and_numeric_solvers_agree_after_substitution(self): 

126 """Symbolic and numeric solvers return the same result after substituting numerical values.""" 

127 result_symbolic: sp.Matrix = solve_stationary( 

128 self.system_matrix_symbol, sp.Matrix([[1, 0], [0, 1]]), sp.Matrix([10, 20]) 

129 ) 

130 result_substituted = np.array( 

131 result_symbolic.subs({self.resistance: 1.0, self.capacity: 1.0}).evalf().tolist(), dtype=float 

132 ).ravel() 

133 np.testing.assert_allclose(result_substituted, self.expected, rtol=1e-10) 

134 

135 def test_symbolic_homogeneous_system_returns_zero(self): 

136 """Symbolic homogeneous system with zero input vector returns zero temperatures.""" 

137 result = solve_stationary(self.system_matrix_symbol, sp.Matrix([[1, 0], [0, 1]]), sp.Matrix([0, 0])) 

138 self.assertEqual(result, sp.zeros(2, 1)) 

139 

140 def test_symbolic_input_vector_as_column_matrix(self): 

141 """Symbolic solve gives the same result for a flat and a (m, 1) column input vector.""" 

142 input_matrix = sp.Matrix([[1, 0], [0, 1]]) 

143 result_flat = solve_stationary(self.system_matrix_symbol, input_matrix, sp.Matrix([10, 20])) 

144 result_column = solve_stationary(self.system_matrix_symbol, input_matrix, sp.Matrix([[10], [20]])) 

145 self.assertEqual(result_flat, result_column) 

146 

147 def test_homogeneous_numeric_returns_zero(self): 

148 """Numeric homogeneous system without input matrices returns zero temperatures.""" 

149 result = solve_stationary(csr_matrix(self.system_matrix_dense)) 

150 np.testing.assert_allclose(result, np.zeros(2), atol=1e-15) 

151 

152 def test_homogeneous_symbolic_returns_zero(self): 

153 """Symbolic homogeneous system without input matrices returns zero temperatures.""" 

154 result = solve_stationary(self.system_matrix_symbol) 

155 self.assertEqual(result, sp.zeros(2, 1)) 

156 

157 def test_only_input_matrix_raises(self): 

158 """Passing only input_matrix without input_vector raises a ValueError.""" 

159 with self.assertRaises(ValueError): 

160 solve_stationary(csr_matrix(self.system_matrix_dense), input_matrix=csr_matrix(self.input_matrix_dense)) 

161 

162 def test_only_input_vector_raises(self): 

163 """Passing only input_vector without input_matrix raises a ValueError.""" 

164 with self.assertRaises(ValueError): 

165 solve_stationary(csr_matrix(self.system_matrix_dense), input_vector=self.input_vector)