Coverage for pyrc \ core \ solver \ stationary.py: 93%

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 numpy as np 

9import sympy as sp 

10from scipy.sparse import csr_matrix, spmatrix, sparray, issparse 

11from scipy.sparse.linalg import spsolve 

12from sympy import MutableDenseMatrix, MutableSparseMatrix, ImmutableSparseMatrix 

13 

14 

15def _to_sympy_matrix(matrix: csr_matrix | sp.MatrixBase) -> sp.MatrixBase: 

16 """Convert a sparse matrix to a sympy Matrix if not already one.""" 

17 if isinstance(matrix, sp.MatrixBase): 

18 return matrix 

19 return sp.Matrix(matrix.toarray()) 

20 

21 

22def _to_sparse_matrix(matrix: csr_matrix | sp.MatrixBase) -> csr_matrix: 

23 """Convert a sympy Matrix to a scipy csr_matrix if not already one.""" 

24 if issparse(matrix): 

25 return matrix 

26 return csr_matrix(np.array(matrix.tolist(), dtype=float)) 

27 

28 

29def _has_free_symbols(*matrices: csr_matrix | sp.MatrixBase | np.ndarray | None) -> bool: 

30 """Return True if any of the given matrices contains free sympy symbols.""" 

31 return any(isinstance(matrix, sp.MatrixBase) and matrix.free_symbols for matrix in matrices if matrix is not None) 

32 

33 

34def _solve_symbolic( 

35 system_matrix: csr_matrix | sp.MatrixBase, 

36 input_matrix: csr_matrix | sp.MatrixBase | None, 

37 input_vector: np.ndarray | sp.MatrixBase | None, 

38) -> sp.MatrixBase: 

39 """Solve the stationary system symbolically using sympy.""" 

40 system_matrix = _to_sympy_matrix(system_matrix) 

41 number_of_nodes = system_matrix.shape[0] 

42 if input_matrix is None: 

43 right_hand_side = sp.zeros(number_of_nodes, 1) 

44 else: 

45 input_matrix = _to_sympy_matrix(input_matrix) 

46 input_vector = sp.Matrix(input_vector) if not isinstance(input_vector, sp.MatrixBase) else input_vector 

47 if input_vector.shape[1] != 1: 

48 input_vector = input_vector.reshape(len(input_vector), 1) 

49 right_hand_side = -(input_matrix * input_vector) 

50 return system_matrix.solve(right_hand_side) 

51 

52 

53def _solve_numeric( 

54 system_matrix: csr_matrix | sp.MatrixBase, 

55 input_matrix: csr_matrix | sp.MatrixBase | None, 

56 input_vector: np.ndarray | sp.MatrixBase | None, 

57) -> np.ndarray: 

58 """Solve the stationary system numerically using scipy sparse.""" 

59 system_matrix = _to_sparse_matrix(system_matrix) 

60 number_of_nodes = system_matrix.shape[0] 

61 if input_matrix is None: 

62 right_hand_side = np.zeros(number_of_nodes) 

63 else: 

64 input_matrix = _to_sparse_matrix(input_matrix) 

65 input_vector_numeric = np.asarray(input_vector, dtype=float).ravel() 

66 right_hand_side = -(input_matrix @ input_vector_numeric) 

67 if hasattr(right_hand_side, "toarray"): 

68 right_hand_side = right_hand_side.toarray().ravel() 

69 stationary_temperatures = spsolve(system_matrix.tocsc(), right_hand_side) 

70 if np.any(np.isnan(stationary_temperatures)): 

71 raise ValueError("System matrix is singular; no unique stationary solution exists.") 

72 return stationary_temperatures 

73 

74 

75def solve_stationary( 

76 system_matrix: csr_matrix | sp.Matrix | MutableSparseMatrix | ImmutableSparseMatrix | spmatrix | sparray, 

77 input_matrix: csr_matrix | sp.Matrix | MutableSparseMatrix | ImmutableSparseMatrix | spmatrix | sparray = None, 

78 input_vector: np.ndarray | sp.Matrix | sparray = None, 

79) -> np.ndarray | sp.Matrix | MutableDenseMatrix: 

80 """ 

81 Solve the stationary solution of a linear RC network: x = -A⁻¹ · B · u. 

82 

83 Dispatches to a numeric (scipy sparse) or symbolic (sympy) solver based on whether free symbols are detected in 

84 the inputs. A purely numeric sympy matrix is converted and solved numerically. For homogeneous systems, 

85 ``input_matrix`` and ``input_vector`` can be omitted, returning the zero vector. 

86 

87 Parameters 

88 ---------- 

89 system_matrix : csr_matrix | sp.Matrix | MutableSparseMatrix | ImmutableSparseMatrix | spmatrix | sparray 

90 System matrix (Jacobian) of shape (n, n). 

91 input_matrix : csr_matrix | sp.Matrix | MutableSparseMatrix | ImmutableSparseMatrix | spmatrix | sparray, optional 

92 Input matrix of shape (n, m). Must be provided together with ``input_vector``. 

93 input_vector : np.ndarray | sp.Matrix, optional 

94 Input vector of shape (m,) or (m, 1). Must be provided together with ``input_matrix``. 

95 

96 Returns 

97 ------- 

98 np.ndarray | sp.Matrix | MutableDenseMatrix : 

99 Stationary temperature vector of shape (n,) for numeric, (n, 1) for symbolic. 

100 

101 Raises 

102 ------ 

103 ValueError : 

104 If only one of ``input_matrix`` and ``input_vector`` is provided. 

105 ValueError : 

106 If the system matrix is singular. 

107 """ 

108 if (input_matrix is None) != (input_vector is None): 

109 raise ValueError("Both input_matrix and input_vector must be provided together, or both omitted.") 

110 

111 if _has_free_symbols(system_matrix, input_matrix, input_vector): 

112 return _solve_symbolic(system_matrix, input_matrix, input_vector) 

113 else: 

114 return _solve_numeric(system_matrix, input_matrix, input_vector)