Coverage for pyrc\core\network.py: 67%
500 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-29 15:57 +0200
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# ------------------------------------------------------------------------------
8from __future__ import annotations
10import operator
11import os
12import pickle
13import threading
14from abc import abstractmethod
15from copy import copy
16from typing import Any, TYPE_CHECKING, Callable, Iterable
18import numpy as np
19from sympy import (
20 lambdify,
21 latex,
22 SparseMatrix,
23 Matrix,
24 diag,
25 Symbol,
26 Basic,
27 ImmutableSparseMatrix,
28 MutableSparseMatrix,
29 Expr,
30)
31import sympy as sym
32import hashlib
33import networkx as nx
34from scipy.sparse import coo_matrix
36from pyrc.core.solver.stationary import solve_stationary
37from pyrc.paths import run_folder
38from pyrc.core.components.resistor import Resistor
39from pyrc.core.inputs import InternalHeatSource
40from pyrc.core.settings import initial_settings, Settings
41from pyrc.core.solver.symbolic import SparseSymbolicEvaluator, ArraySymbolicEvaluator
42from pyrc.core.solver.handler import InhomogeneousSystemHandler, SolveIVPHandler, HomogeneousSystemHandler
43from pyrc.tools.errors import HighCourantNumberError
44from pyrc.tools.functions import add_leading_underscore
46from pyrc.tools.science import build_jacobian
47from pyrc.core.components.templates import RCObjects, RCSolution, EquationItem, initial_rc_objects, Cell
48from pyrc.core.components.templates import solution_object
49from pyrc.core.components.capacitor import Capacitor
51if TYPE_CHECKING:
52 from pyrc.core.nodes import MassFlowNode
53 from pyrc.core.components.templates import EquationItemInput
54 from scipy.sparse import spmatrix, sparray
57class RCNetwork:
58 def __init__(
59 self,
60 save_to_pickle: bool = True,
61 load_from_pickle: bool = True,
62 load_solution: bool = True,
63 num_cores_jacobian: int = 1,
64 # TODO: Parallelization only works on linux currently (because it forks instead
65 # of spawning it). Problem is: Sympy symbols are not picklable.
66 settings: Settings = initial_settings,
67 rc_objects: RCObjects = initial_rc_objects,
68 rc_solution: RCSolution = solution_object,
69 # strategy: str = "lambdify",
70 continue_canceled: bool = True,
71 skip_infinity_check: bool = False,
72 ):
73 """
75 Parameters
76 ----------
77 save_to_pickle
78 load_from_pickle
79 load_solution
80 num_cores_jacobian
81 settings
82 rc_objects
83 rc_solution
84 continue_canceled : bool, default=True
85 If True, continue canceled simulations.
86 Overwrites load_solution
87 skip_infinity_check : bool, default=False
88 If True the infinity check of the created matrices is skipped.
89 This is recommended for complicated dynamic (time-dependent) systems.
90 """
91 self.rc_objects: RCObjects = rc_objects
92 self.settings: Settings = settings
93 self._system_matrix_symbol: MutableSparseMatrix = None
94 self._input_matrix_symbol: MutableSparseMatrix = None
95 self._temperature_vector_symbol = None
96 self._input_vector_symbol = None
97 self._input_matrix = None
98 self._system_matrix = None
99 self.__hash = None
100 self.rc_solution: RCSolution = rc_solution
101 # link rc objects to rc solution
102 self.rc_solution.rc_objects = rc_objects
103 self.skip_infinity_check = skip_infinity_check
105 # For each time dependent symbol the following values are not allowed because they would make an entry in the
106 # system matrix be infinity
107 self.forbidden_values: dict[Symbol, set] = {}
109 self.save_to_pickle: bool = save_to_pickle
110 self.load_from_pickle: bool = load_from_pickle
111 self.load_solution: bool = load_solution
112 self.continue_canceled: bool = continue_canceled
114 if num_cores_jacobian is None:
115 num_cores_jacobian = os.cpu_count()
116 self.num_cores_jacobian = min(max(num_cores_jacobian, 1), os.cpu_count())
118 self.groups: dict = {}
120 def __getattr__(self, item):
121 if getattr(self.rc_objects, item, None):
122 return self.rc_objects.__getattribute__(item)
123 return AttributeError
125 def create_network(self, *args, **kwargs):
126 EquationItem.reset_counter()
127 Resistor.reset_counter()
128 self.rc_objects.wipe_all()
129 self._create_network(*args, **kwargs)
130 assert self.nodes is not None
132 @abstractmethod
133 def _create_network(self, *args, **kwargs):
134 pass
136 def copy_dict(self):
137 return {
138 "save_to_pickle": self.save_to_pickle,
139 "load_from_pickle": self.load_from_pickle,
140 "load_solution": self.load_solution,
141 "num_cores_jacobian": self.num_cores_jacobian,
142 "settings": copy(self.settings),
143 }
145 def __copy__(self):
146 return RCNetwork(**self.copy_dict())
148 @property
149 def nodes(self):
150 return self.rc_objects.nodes
152 @property
153 def hash(self):
154 """
155 Compute a deterministic hash for the :class:`RCNetwork` where node identity (class name + obj.id) and topology matter.
157 For :class:`Capacitor` instances, an optional internal_heat_source (identified by its id) is encoded as a node attribute.
159 For :class:`Cell` subclass instances, position and delta arrays are encoded.
161 Returns
162 -------
163 str :
164 SHA-256 hex digest of the canonical edge and node representation.
165 """
166 if self.__hash is None:
168 def node_key(_obj) -> str:
169 return f"{type(_obj).__name__}:{_obj.id}"
171 def array_str(arr: np.ndarray) -> str:
172 return ",".join(f"{v:.14g}" for v in arr.ravel())
174 graph = nx.Graph()
175 for obj in self.rc_objects.all:
176 attrs = {}
177 if isinstance(obj, Capacitor):
178 attrs["heat_source_id"] = (
179 str(obj.internal_heat_source.id) if obj.internal_heat_source is not None else None
180 )
181 attrs["heat_source_symbol_name"] = (
182 obj.internal_heat_source.power.name
183 if obj.internal_heat_source is not None and isinstance(obj.internal_heat_source.power, Basic)
184 else None
185 )
186 attrs["capacitor_symbol_name"] = obj.capacity_symbol.name
187 attrs["capacity_variable"] = "yes" if isinstance(obj.capacity, Basic) else "no"
188 if isinstance(obj, Cell):
189 attrs["position"] = array_str(obj.position)
190 attrs["delta"] = array_str(obj.delta)
191 if isinstance(obj, Resistor):
192 attrs["resistance_symbol"] = obj.resistance_symbol.name
193 attrs["resistance_variable"] = "yes" if isinstance(obj.resistance, Basic) else "no"
194 graph.add_node(node_key(obj), **attrs)
196 # loop over all linked/connected network objects
197 for obj in self.rc_objects.all:
198 for neighbour in obj.neighbours:
199 graph.add_edge(node_key(obj), node_key(neighbour))
201 node_str = "|".join(
202 f"{nid}#{data}"
203 for nid, data in sorted((nid, str(sorted(data.items()))) for nid, data in graph.nodes(data=True))
204 )
205 edge_str = "|".join(sorted(f"{min(u, v)}-{max(u, v)}" for u, v in graph.edges()))
206 canonical = f"nodes:{node_str};edges:{edge_str}"
207 self.__hash = hashlib.sha256(canonical.encode()).hexdigest()
208 return self.__hash
210 @property
211 def time_steps(self) -> list:
212 """
213 Returns a list with all time steps.
215 Returns
216 -------
217 list :
218 The time steps.
219 """
220 return self.rc_solution.time_steps.tolist()
222 def reset_properties(self):
223 self._temperature_vector_symbol = None
225 @property
226 def inputs(self) -> list[EquationItemInput]:
227 return self.rc_objects.inputs
229 def change_static_dynamic(self, calculate_static: bool | str):
230 if isinstance(calculate_static, str):
231 match calculate_static.lower():
232 case "static" | "statisch" | "s":
233 calculate_static = True
234 case _:
235 calculate_static = False
236 self.settings.calculate_static = calculate_static
238 def get_node_by_id(self, _id: str | int):
239 """
240 Returns the node with the given id.
242 Parameters
243 ----------
244 _id : str | int
245 The id of the node.
247 Returns
248 -------
249 TemperatureNode | InternalHeatSource :
250 """
251 if isinstance(_id, str):
252 _id = int(_id)
253 return next((x for x in self.rc_objects.all_equation_objects if x.id == _id), None)
255 def no_infinity(self, sp_matrix: SparseMatrix | spmatrix | sparray) -> bool:
256 """
257 Checks the given sparse matrix for infinity entries and symbolic singularities.
259 Notes
260 -----
261 Saves a dict, which maps each free sympy symbol to a set of values at which at least one symbolic entry
262 diverges. Empty if no symbolic entries exist.
264 Parameters
265 ----------
266 sp_matrix : SparseMatrix | spmatrix | sparray
267 A sparse matrix/vector to check, may contain numeric values or sympy
268 expressions.
270 Returns
271 -------
272 tuple[bool, dict] :
273 bool : True if no numeric infinity was found.
274 """
275 sp_matrix: SparseMatrix | spmatrix | sparray
276 if isinstance(sp_matrix, SparseMatrix) or isinstance(sp_matrix, ImmutableSparseMatrix):
277 entries = sp_matrix.todok().values()
278 else:
279 sp_matrix: spmatrix | sparray
280 entries = sp_matrix.data
282 def _substitute_thetas(expression: Expr) -> Expr:
283 theta_symbols = sorted(
284 {s for s in expression.free_symbols if s.name.startswith("theta_")}, key=lambda s: s.name
285 )
286 num_thetas = len(theta_symbols)
287 substitution = {symbol: 270 + (index / num_thetas) * 60 for index, symbol in enumerate(theta_symbols)}
288 return expression.subs(substitution)
290 def _solve_with_timeout(expression: Expr, symbol: sym.Symbol, timeout: float = 0.5) -> list:
291 result = []
292 timed_out = threading.Event()
294 def target():
295 result.extend(sym.solve(expression, symbol))
297 thread = threading.Thread(target=target)
298 thread.start()
299 thread.join(timeout)
300 if thread.is_alive():
301 timed_out.set()
303 if timed_out.is_set():
304 raise TimeoutError(f"sym.solve timed out after {timeout}s for symbol {symbol}")
306 return result
308 for entry in entries:
309 if isinstance(entry, Expr):
310 if not entry.free_symbols:
311 if entry in (sym.oo, -sym.oo, sym.zoo, sym.nan):
312 return False
313 else:
314 _, denominator = sym.fraction(sym.cancel(entry))
315 if denominator == sym.Integer(0):
316 return False
317 if denominator != sym.Integer(1):
318 substituted_denominator = _substitute_thetas(denominator)
319 float_atoms = {a for a in substituted_denominator.atoms(sym.Float)}
320 dummy_map = {f: sym.Symbol(f"staticVALUE_{i}") for i, f in enumerate(float_atoms)}
321 inverse_dummy_map = {v: k for k, v in dummy_map.items()}
322 cleaned_denominator = substituted_denominator.subs(dummy_map)
323 try:
324 for s in cleaned_denominator.free_symbols - set(dummy_map.values()):
325 solution = _solve_with_timeout(cleaned_denominator, s, timeout=1)
326 solution = [sol.subs(inverse_dummy_map) for sol in solution]
327 self.forbidden_values.setdefault(s, set()).update(solution)
328 except (NotImplementedError, TimeoutError):
329 print(f"No infinity check for eq: {cleaned_denominator}")
330 else:
331 if np.isinf(entry):
332 return False
334 return True
336 def _get_matrix_symbol(self, key: str) -> SparseMatrix | ImmutableSparseMatrix:
337 """
338 Returns the desired matrix in symbolic notation.
340 Parameters
341 ----------
342 key : str
343 Either "symbol" or "input" to determine the matrix type.
345 Returns
346 -------
347 SparseMatrix | ImmutableSparseMatrix :
348 The (sympy) symbolic notation of the matrix.
349 """
350 attr = f"_{key}_matrix_symbol"
351 if getattr(self, attr) is None:
352 self.make_system_matrices()
353 return getattr(self, attr)
355 def _build_matrix_function(self, key: str) -> Callable:
356 """
357 Builds the desired matrix lambda function using the symbols matrix, respectively.
359 Parameters
360 ----------
361 key : str
362 Either "symbol" or "input" to determine the matrix type.
364 Returns
365 -------
366 Callable :
367 The function of the matrix.
368 """
369 matrix_symbol = self._get_matrix_symbol(key)
370 all_symbols = self.system_symbols
371 if self.time_dependent_system_symbols:
372 free_symbols_in_matrix = matrix_symbol.free_symbols
373 ordered_symbols = [
374 next((s for s in free_symbols_in_matrix if s.name == symbol.name), symbol) for symbol in all_symbols
375 ]
376 return lambdify(ordered_symbols, matrix_symbol, "sympy")
377 return lambdify(all_symbols, matrix_symbol, "scipy")
379 def _build_matrix_direct(self, key: str) -> spmatrix:
380 """
381 Bypasses lambdify by substituting all symbol values via name-based xreplace.
383 Parameters
384 ----------
385 key : str
386 Either "symbol" or "input" to determine the matrix type.
388 Notes
389 -----
390 The replace must be done using the name of the symbols and not the objects itself because they might be loaded
391 from disk and might not match, but the names do.
393 Returns
394 -------
395 scipy.spmatrix :
396 A scipy sparse matrix with all values substituted.
397 """
398 matrix_symbol = self._get_matrix_symbol(key)
399 name_to_value = {s.name: v for s, v in zip(self.system_symbols, self.system_values)}
400 free = matrix_symbol.free_symbols
401 substitution_map = {s: name_to_value[s.name] for s in free if s.name in name_to_value}
403 dok = matrix_symbol.todok()
404 rows, cols, data = [], [], []
405 for (i, j), expr in dok.items():
406 rows.append(i)
407 cols.append(j)
408 data.append(float(expr.xreplace(substitution_map)))
410 return coo_matrix((data, (rows, cols)), shape=matrix_symbol.shape).tocsr()
412 # import warnings
413 # import logging
414 # import pickle
415 # from concurrent.futures import ProcessPoolExecutor
416 # from pathlib import Path
417 # from typing import Callable
418 #
419 # import numba
420 # import numpy as np
421 # from scipy.sparse import coo_matrix, csr_matrix
422 # from sympy import lambdify
423 # from sympy.printing import ccode
424 #
425 # def _xreplace_element(args: tuple) -> object:
426 # """Apply xreplace to a single expression."""
427 # expr, static_map = args
428 # return expr.xreplace(static_map)
429 #
430 # def _build_c_source(reduced_exprs: list, td_ordered: list, func_name: str) -> str:
431 # """
432 # Generates a C source file with a single function iterating over all nonzero
433 # expressions.
434 #
435 # Parameters
436 # ----------
437 # reduced_exprs : list
438 # Sympy expressions with static symbols already substituted.
439 # td_ordered : list
440 # Ordered time-dependent sympy symbols.
441 # func_name : str
442 # Name of the generated C function.
443 #
444 # Returns
445 # -------
446 # str :
447 # C source code as a string.
448 # """
449 # args = ", ".join(f"double {s.name}" for s in td_ordered)
450 # lines = [
451 # "#include <stddef.h>",
452 # f"void {func_name}({args}, double *out, size_t n) {{",
453 # ]
454 # for i, expr in enumerate(reduced_exprs):
455 # lines.append(f" out[{i}] = {ccode(expr)};")
456 # lines.append("}")
457 # return "\n".join(lines)
458 #
459 # def _compile_c_function(
460 # source: str,
461 # func_name: str,
462 # n_out: int,
463 # cache_path: Path,
464 # ) -> Callable:
465 # """
466 # Compiles C source to a shared library and returns a callable via ctypes.
467 #
468 # Parameters
469 # ----------
470 # source : str
471 # C source code string.
472 # func_name : str
473 # Name of the C function to load.
474 # n_out : int
475 # Number of output values (nnz).
476 # cache_path : Path
477 # Path to write the compiled .so / .dll.
478 #
479 # Returns
480 # -------
481 # Callable :
482 # Python callable wrapping the compiled C function.
483 # """
484 # import ctypes
485 # import subprocess
486 # import tempfile
487 #
488 # src_path = cache_path.with_suffix(".c")
489 # src_path.write_text(source)
490 #
491 # result = subprocess.run(
492 # ["gcc", "-O3", "-shared", "-fPIC", "-o", str(cache_path), str(src_path)],
493 # capture_output=True,
494 # text=True,
495 # )
496 # if result.returncode != 0:
497 # raise RuntimeError(result.stderr)
498 #
499 # lib = ctypes.CDLL(str(cache_path))
500 # fn = getattr(lib, func_name)
501 # fn.restype = None
502 #
503 # out_type = ctypes.c_double * n_out
504 #
505 # def c_callable(*td_vals: float) -> np.ndarray:
506 # out = out_type()
507 # fn(*[ctypes.c_double(v) for v in td_vals], out, ctypes.c_size_t(n_out))
508 # return np.frombuffer(out, dtype=np.float64).copy()
509 #
510 # return c_callable
512 # def _prepare_matrix_ingredients(
513 # self, key: str, parallel: bool
514 # ) -> tuple[list, list, np.ndarray, np.ndarray, tuple, list]:
515 # """
516 # Extracts and reduces nonzero expressions, substituting all static symbols.
517 #
518 # Parameters
519 # ----------
520 # key : str
521 # Either "symbol" or "input" to determine the matrix type.
522 # parallel : bool
523 # If True, parallelizes static substitution via ProcessPoolExecutor.
524 #
525 # Returns
526 # -------
527 # tuple :
528 # td_ordered, reduced_exprs, rows, cols, shape, raw exprs (unreduced)
529 # """
530 # matrix_symbol = self._get_matrix_symbol(key)
531 # all_symbols = self.get_symbols()
532 # all_values = self.get_values()
533 # td_symbols = self.get_time_dependent_symbols()
534 # td_names = {s.name for s in td_symbols}
535 #
536 # static_map = {
537 # s_mat: val
538 # for sym, val in zip(all_symbols, all_values)
539 # if sym.name not in td_names
540 # for s_mat in matrix_symbol.free_symbols
541 # if s_mat.name == sym.name
542 # }
543 #
544 # td_ordered = [
545 # next(s for s in matrix_symbol.free_symbols if s.name == td.name)
546 # for td in td_symbols
547 # if any(s.name == td.name for s in matrix_symbol.free_symbols)
548 # ]
549 #
550 # dok = matrix_symbol.todok()
551 # shape = matrix_symbol.shape
552 # keys, exprs = zip(*dok.items()) if dok else ([], [])
553 # rows = np.array([k[0] for k in keys], dtype=np.int32)
554 # cols = np.array([k[1] for k in keys], dtype=np.int32)
555 #
556 # TODO: do not do it like that, it is worse than the solution that already exists.
557 # if parallel:
558 # tasks = [(expr, static_map) for expr in exprs]
559 # with ProcessPoolExecutor() as executor:
560 # reduced_exprs = list(executor.map(_xreplace_element, tasks))
561 # else:
562 # reduced_exprs = [expr.xreplace(static_map) for expr in exprs]
563 #
564 # return td_ordered, reduced_exprs, rows, cols, shape
565 #
566 # TODO: This is worse than the current solution, but the other version (numby/c...) should be implemented to speed
567 # things up. But attention: Currently, it is just vibe coded...
568 # def _build_matrix_function_lambdify(
569 # self, key: str, parallel: bool = False
570 # ) -> Callable:
571 # """
572 # Builds the matrix function using lambdify with numpy backend.
573 #
574 # Parameters
575 # ----------
576 # key : str
577 # Either "symbol" or "input" to determine the matrix type.
578 # parallel : bool
579 # If True, parallelizes static symbol substitution.
580 #
581 # Returns
582 # -------
583 # Callable :
584 # Function accepting td symbol values, returning a csr_matrix.
585 # """
586 # td_ordered, reduced_exprs, rows, cols, shape = (
587 # self._prepare_matrix_ingredients(key, parallel)
588 # )
589 #
590 # if not len(rows):
591 # empty = csr_matrix(shape, dtype=float)
592 #
593 # def assembled(*_) -> csr_matrix:
594 # return empty
595 #
596 # return assembled
597 #
598 # data_func = lambdify(td_ordered, reduced_exprs, "numpy")
599 # template = coo_matrix(
600 # (np.ones(len(rows), dtype=float), (rows, cols)), shape=shape
601 # ).tocsr()
602 #
603 # def assembled(*td_vals: float) -> csr_matrix:
604 # template.data[:] = np.asarray(data_func(*td_vals), dtype=float)
605 # return template
606 #
607 # return assembled
608 #
609 # def _build_matrix_function_numba(self, key: str, parallel: bool = False) -> Callable:
610 # """
611 # Builds the matrix function using a numba-JIT compiled callable, with
612 # caching to self.save_folder_path keyed by self.hash.
613 #
614 # Parameters
615 # ----------
616 # key : str
617 # Either "symbol" or "input" to determine the matrix type.
618 # parallel : bool
619 # If True, parallelizes static symbol substitution.
620 #
621 # Returns
622 # -------
623 # Callable :
624 # Function accepting td symbol values, returning a csr_matrix.
625 # """
626 # td_ordered, reduced_exprs, rows, cols, shape = (
627 # self._prepare_matrix_ingredients(key, parallel)
628 # )
629 #
630 # if not len(rows):
631 # empty = csr_matrix(shape, dtype=float)
632 #
633 # def assembled(*_) -> csr_matrix:
634 # return empty
635 #
636 # return assembled
637 #
638 # cache_path = Path(self.save_folder_path) / f"{self.hash}_{key}_numba.pkl"
639 #
640 # if cache_path.exists():
641 # with open(cache_path, "rb") as f:
642 # data_func = pickle.load(f)
643 # else:
644 # data_func = lambdify(td_ordered, reduced_exprs, "numpy")
645 # data_func_nb = numba.njit(data_func, cache=True)
646 # with open(cache_path, "wb") as f:
647 # pickle.dump(data_func_nb, f)
648 # data_func = data_func_nb
649 #
650 # template = coo_matrix(
651 # (np.ones(len(rows), dtype=float), (rows, cols)), shape=shape
652 # ).tocsr()
653 #
654 # def assembled(*td_vals: float) -> csr_matrix:
655 # template.data[:] = np.asarray(data_func(*td_vals), dtype=float)
656 # return template
657 #
658 # return assembled
659 #
660 # def _build_matrix_function_c(self, key: str, parallel: bool = False) -> Callable:
661 # """
662 # Builds the matrix function using a compiled C shared library via ctypes,
663 # falling back to numba if compilation fails.
664 #
665 # Parameters
666 # ----------
667 # key : str
668 # Either "symbol" or "input" to determine the matrix type.
669 # parallel : bool
670 # If True, parallelizes static symbol substitution.
671 #
672 # Returns
673 # -------
674 # Callable :
675 # Function accepting td symbol values, returning a csr_matrix.
676 # """
677 # td_ordered, reduced_exprs, rows, cols, shape = (
678 # self._prepare_matrix_ingredients(key, parallel)
679 # )
680 #
681 # if not len(rows):
682 # empty = csr_matrix(shape, dtype=float)
683 #
684 # def assembled(*_) -> csr_matrix:
685 # return empty
686 #
687 # return assembled
688 #
689 # func_name = f"matrix_func_{self.hash}_{key}"
690 # so_path = Path(self.save_folder_path) / f"{func_name}.so"
691 #
692 # try:
693 # if so_path.exists():
694 # import ctypes
695 # lib = ctypes.CDLL(str(so_path))
696 # fn = getattr(lib, func_name)
697 # n_out = len(rows)
698 # out_type = ctypes.c_double * n_out
699 #
700 # def data_func(*td_vals: float) -> np.ndarray:
701 # out = out_type()
702 # fn(
703 # *[ctypes.c_double(v) for v in td_vals],
704 # out,
705 # ctypes.c_size_t(n_out),
706 # )
707 # return np.frombuffer(out, dtype=np.float64).copy()
708 # else:
709 # source = _build_c_source(reduced_exprs, td_ordered, func_name)
710 # data_func = _compile_c_function(
711 # source, func_name, len(rows), so_path
712 # )
713 # except Exception as e:
714 # print(
715 # f"C compilation failed for '{key}' matrix ({e}), falling back to numba."
716 # )
717 # return self._build_matrix_function_numba(key, parallel=parallel)
718 #
719 # template = coo_matrix(
720 # (np.ones(len(rows), dtype=float), (rows, cols)), shape=shape
721 # ).tocsr()
722 #
723 # def assembled(*td_vals: float) -> csr_matrix:
724 # template.data[:] = np.asarray(data_func(*td_vals), dtype=float)
725 # return template
726 #
727 # return assembled
728 #
729 # def _build_matrix_function(
730 # self, key: str, strategy: str | None = None
731 # ) -> Callable:
732 # """
733 # Dispatcher selecting the matrix function build strategy.
734 #
735 # Parameters
736 # ----------
737 # key : str
738 # Either "symbol" or "input" to determine the matrix type.
739 # strategy : str | None
740 # One of "lambdify", "parallel", "numba", "c". If None, falls back
741 # to self.strategy set at __init__.
742 #
743 # Returns
744 # -------
745 # Callable :
746 # Function accepting td symbol values, returning a csr_matrix.
747 # """
748 # s = strategy if strategy is not None else self.strategy
749 # if s == "lambdify":
750 # return self._build_matrix_function_lambdify(key, parallel=False)
751 # if s == "parallel":
752 # return self._build_matrix_function_lambdify(key, parallel=True)
753 # if s == "numba":
754 # return self._build_matrix_function_numba(key, parallel=False)
755 # if s == "c":
756 # return self._build_matrix_function_c(key, parallel=False)
757 # raise ValueError(f"Unknown strategy '{s}'. Choose from: lambdify, parallel, numba, c.")
758 #
759 # def _build_matrix(
760 # self, key: str, strategy: str | None = None
761 # ) -> csr_matrix:
762 # """
763 # Inputs all values into the desired matrix function.
764 #
765 # Parameters
766 # ----------
767 # key : str
768 # Either "symbol" or "input" to determine the matrix type.
769 # strategy : str | None
770 # Overrides self.strategy if provided.
771 #
772 # Returns
773 # -------
774 # csr_matrix :
775 # The desired matrix.
776 # """
777 # attr = f"_{key}_matrix"
778 # if getattr(self, attr) is None:
779 # if not self.get_time_dependent_symbols():
780 # result = self._build_matrix_direct(key)
781 # else:
782 # result = self._build_matrix_function(key, strategy=strategy)(
783 # *self.get_values()
784 # )
785 # assert self.no_infinity(result), (
786 # f"Infinity detected in {key} matrix. System cannot be solved."
787 # )
788 # setattr(self, attr, result)
789 # return getattr(self, attr)
791 def _build_matrix(self, key: str) -> SparseMatrix | ImmutableSparseMatrix | spmatrix | sparray:
792 """
793 Inputs all values (can be symbolic values) into the desired matrix function.
795 Parameters
796 ----------
797 key : str
798 Either "symbol" or "input" to determine the matrix type.
800 Returns
801 -------
802 SparseMatrix | ImmutableSparseMatrix | spmatrix | sparray :
803 The desired matrix. If no time dependent variables are used, it is a scipy sparse matrix, otherwise it's a
804 sympy SparseMatrix.
805 """
806 attr = f"_{key}_matrix"
807 if getattr(self, attr) is None:
808 if self.time_dependent_system_symbols:
809 result = self._build_matrix_function(key)(*self.system_values)
810 else:
811 result = self._build_matrix_direct(key)
812 if not self.skip_infinity_check:
813 assert self.no_infinity(result), f"Infinity detected in {key} matrix. System cannot be solved."
814 setattr(self, attr, result)
815 return getattr(self, attr)
817 @property
818 def system_matrix_symbol(self) -> SparseMatrix | ImmutableSparseMatrix:
819 return self._get_matrix_symbol("system")
821 @property
822 def input_matrix_symbol(self) -> SparseMatrix | ImmutableSparseMatrix:
823 return self._get_matrix_symbol("input")
825 @property
826 def system_matrix_function(self) -> Callable:
827 return self._build_matrix_function("system")
829 @property
830 def input_matrix_function(self) -> Callable:
831 return self._build_matrix_function("input")
833 @property
834 def system_matrix(self) -> SparseMatrix | ImmutableSparseMatrix | spmatrix | sparray:
835 return self._build_matrix("system")
837 @property
838 def input_matrix(self) -> SparseMatrix | ImmutableSparseMatrix | spmatrix | sparray:
839 return self._build_matrix("input")
841 @property
842 def temperature_vector_symbol(self) -> list:
843 if self._temperature_vector_symbol is None:
844 result = [node.temperature_symbol for node in self.nodes]
845 self._temperature_vector_symbol = result
846 return self._temperature_vector_symbol
848 @property
849 def temperature_vector(self) -> np.ndarray:
850 result = [node.temperature for node in self.nodes]
851 return np.array(result).flatten()
853 @property
854 def input_vector(self) -> np.ndarray:
855 result = [n for input_item in self.inputs for n in input_item.values]
856 return np.array(result).reshape(-1,1)
858 @property
859 def input_vector_symbol(self) -> np.ndarray:
860 if self._input_vector_symbol is None:
861 result = []
862 if self.inputs is not None and self.inputs != []:
863 result = [n for input_item in self.inputs for n in input_item.symbols]
864 self._input_vector_symbol = np.array(result).flatten()
865 return self._input_vector_symbol
867 # @property
868 # def seconds_of_day(self):
869 # """
870 # Returns the seconds that have already passed on the day of ``self.start_time``.
871 #
872 # Returns
873 # -------
874 # float | int :
875 # The passed seconds of the day.
876 # """
877 # return (self.settings.start_date - self.settings.start_date.astype("datetime64[D]")) / np.timedelta64(1, "s")
879 @property
880 def save_folder_path(self):
881 if self.settings.save_folder_path is None:
882 return run_folder
883 return self.settings.save_folder_path
885 @property
886 def pickle_path(self) -> str:
887 filename = f"{self.hash}.pickle"
888 return os.path.normpath(os.path.join(self.save_folder_path, filename))
890 def pickle_path_result(self, t_span, name_add_on=None):
891 name_add_on = add_leading_underscore(name_add_on)
892 filename = f"{self.hash}{name_add_on}_{int(t_span[1])}_result.pickle"
893 return os.path.normpath(os.path.join(self.save_folder_path, filename))
895 @property
896 def pickle_path_single_solution(self) -> str:
897 filename = f"{self.hash}_single_solution.pickle"
898 return os.path.normpath(os.path.join(self.save_folder_path, filename))
900 def pickle_path_solution(self, t_span, name_add_on=None) -> str:
901 name_add_on = add_leading_underscore(name_add_on)
902 filename = f"{self.hash}{name_add_on}_{int(t_span[0])}_{int(t_span[1])}_solution.pickle"
903 return os.path.normpath(os.path.join(self.save_folder_path, filename))
905 def save_matrices(self):
906 if self.save_to_pickle:
907 file_path = self.pickle_path
908 matrices = [
909 self._system_matrix_symbol,
910 self._input_matrix_symbol,
911 self._input_vector_symbol,
912 self._temperature_vector_symbol,
913 self.forbidden_values,
914 ]
915 with open(file_path, "wb") as f:
916 pickle.dump(matrices, f)
918 def load_matrices(self):
919 file_path = self.pickle_path
920 if os.path.exists(file_path):
921 with open(file_path, "rb") as f:
922 matrices = pickle.load(f)
923 self._system_matrix_symbol = matrices[0]
924 self._input_matrix_symbol = matrices[1]
925 self._input_vector_symbol = matrices[2]
926 self._temperature_vector_symbol = matrices[3]
927 self.forbidden_values = matrices[4]
928 print("matrices loaded")
929 return True
930 else:
931 print("No network found nor loaded.")
932 return False
934 def save_last_step_solution(self, pickle_path_single_solution: str = None):
935 """
936 Saves the last time step solution.
937 """
938 if pickle_path_single_solution is None:
939 pickle_path_single_solution = self.pickle_path_single_solution
940 file_path = pickle_path_single_solution
941 self.rc_solution.save_last_step(file_path)
943 def load_initial_values(self, return_bool: bool = False, pickle_path_single_solution: str = None) -> None | bool:
944 """
945 Loads the last time step solution (temperature vector) and sets it as initial values.
946 """
947 if pickle_path_single_solution is None:
948 pickle_path_single_solution = self.pickle_path_single_solution
949 file_path = pickle_path_single_solution
950 if os.path.exists(file_path):
951 with open(file_path, "rb") as f:
952 solution_tuple = pickle.load(f)
953 temp_vector = solution_tuple[0]
954 input_vector = solution_tuple[1]
955 for node, temp in zip(self.rc_objects.nodes, temp_vector):
956 node.initial_temperature = temp
957 for item, value in zip(self.inputs, input_vector):
958 item.initial_value = value
959 if return_bool:
960 return True
961 else:
962 print("No temperature vector found nor loaded.")
963 if return_bool:
964 return False
965 return None
967 @property
968 def inhomogeneous_system(self) -> bool:
969 """
970 Returns True if inputs are existing and the network is inhomogeneous. False otherwise.
972 Returns
973 -------
974 bool :
975 True if inputs are existing and the network is inhomogeneous. False otherwise.
976 """
977 if self.rc_objects.inputs is None or self.rc_objects.inputs == []:
978 return False
979 return True
981 def determine_max_step(self) -> float | None:
982 max_step = 1
983 loop_counter = 0
984 if self.rc_objects.mass_flow_nodes is not None or self.rc_objects.mass_flow_nodes != []:
985 while loop_counter < 1000:
986 try:
987 self.check_courant(max_step) # approximated maximum time step used during iterations.
988 # if no error
989 break
990 except HighCourantNumberError:
991 loop_counter += 1
992 max_step *= 0.99
993 if loop_counter >= 1000:
994 raise HighCourantNumberError
995 if loop_counter == 0:
996 return None
997 return max_step
999 def solve_stationary(self, time_step: float | int = 0):
1000 """
1001 Determines the analytic solution of the network (if it exists). It works for symbolic and numeric matrices.
1003 Parameters
1004 ----------
1005 time_step : float | int, default=0
1006 The time step as which the solution is saved in the rc_solution object.
1007 This parameter has no effect if the network is symbolic.
1009 Returns
1010 -------
1011 sp.Matrix | MutableDenseMatrix or None:
1012 Stationary symbolic temperature vector of shape (n, 1) (only if result is symbolic, otherwise nothing is
1013 returned and the result is written to the rc_solution object using time_step).
1014 """
1015 if self.inhomogeneous_system:
1016 result = solve_stationary(
1017 system_matrix=self.system_matrix,
1018 input_matrix=self.input_matrix,
1019 input_vector=self.input_vector,
1020 )
1021 else:
1022 result = solve_stationary(
1023 system_matrix=self.system_matrix,
1024 )
1026 if isinstance(result, Basic) and result.free_symbols is not None:
1027 # return symbolic result
1028 return result
1029 else:
1030 # save result in rc_solution
1031 self.rc_solution.add_to_solution(
1032 [
1033 np.array(time_step).reshape(
1034 -1,
1035 )
1036 ],
1037 [np.array(result).reshape(-1, 1)],
1038 )
1039 return None
1041 def solve_network(
1042 self,
1043 t_span: tuple,
1044 print_progress=False,
1045 name_add_on: str = "",
1046 check_courant: bool = True,
1047 time_dependent_tuple: tuple[Iterable, Callable] | None = None,
1048 time_dependent_tuple_input: tuple[Iterable, Callable] | None = None,
1049 hook_function: Callable = None,
1050 expected_solution_size_mb=5000,
1051 **kwargs: dict[str, int | Any],
1052 ):
1053 """
1054 Solves the network at the given times.
1056 If time_vector is None, one second is calculated.
1058 Parameters
1059 ----------
1060 t_span : tuple | Any
1061 Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf.
1062 Both t0 and tf must be floats or values interpretable by the float conversion function.
1063 print_progress : bool, default=False
1064 If True, some print-outs are made during solving to get the time step that is currently simulated.
1065 name_add_on : str, default=""
1066 Optional add-on to the name of in-between saves that is placed after the hash value (separated by "_").
1067 Example save name:
1068 42395d3d9f07f06ce9c9cb609b406aa54fc2dc5e8c4d9cc75987d9a02541ad2f_nameAddOn_zw_000001080_h.pickle
1069 check_courant : bool, default=True
1070 If True, self.check_courant(0.4) is run before simulating.
1071 time_dependent_tuple : tuple[Iterable, Callable] | None, optional
1072 A ordered list with the time dependent symbols and the function that calculates their values.
1073 The list represents the order of the output of the function.
1074 The function calculates the value of the time dependent symbols in the order of the list. It gets passed
1075 the time step, temperature vector and input vector (last one only if existing):
1076 value1, value2, ... = my_function(time_step, temperature_vector, input_vector)
1077 time_dependent_tuple_input : tuple[Iterable, Callable] | None, optional
1078 A ordered list with the time dependent input vector symbols and the function that calculates their values.
1079 The list represents the order of the output of the function.
1080 A function that calculates all time dependent variables within the time step and returns them in the same order as
1081 RCNetwork.variable_input_vector_symbols.
1082 It gets parameters like this:\n
1083 ``time_dependent_function(time, temperature_vector)`` \n
1084 This function is required if time dependent symbols exist in the input vector.
1085 It must return an iterable (e.g. list).
1086 hook_function : Callable, optional
1087 A function that is run before the network is solved using solve_ivp.
1088 Can be used to catch the time before the solving actually starts (e.g. used in `PerformanceTest`\\).
1089 expected_solution_size_mb : int | float, default=5000
1090 The expected solution size in Megabytes. Is used for RAM-Management in ``SystemHandler.solve()``\\.
1091 The solution size depends on the size of the RC network and number of saved time steps.
1092 kwargs : dict[str: int | Any], optional
1093 Passed to SystemHandler.
1094 """
1095 time_dependent_function = None
1096 if self.time_dependent_system_symbols:
1097 assert time_dependent_tuple is not None
1098 assert time_dependent_tuple[0] is not None
1099 assert time_dependent_tuple[1] is not None
1100 time_dependent_function: Callable = self.remap_symbol_function(*time_dependent_tuple)
1101 time_dependent_function_input = None
1102 if self.variable_input_vector_symbols:
1103 assert time_dependent_tuple_input is not None
1104 assert time_dependent_tuple_input[0] is not None
1105 assert time_dependent_tuple_input[1] is not None
1106 time_dependent_function_input = self.remap_symbol_function(
1107 *time_dependent_tuple_input, wanted_order=self.variable_input_vector_symbols
1108 )
1110 continued_simulation = False
1111 if self.load_solution or self.continue_canceled:
1112 success = self.rc_solution.load_solution(
1113 self.pickle_path_solution(t_span, name_add_on), last_time_step=t_span[-1]
1114 )
1115 if not success:
1116 success = self.rc_solution.load_solution(
1117 self.pickle_path_result(t_span, name_add_on), last_time_step=t_span[-1]
1118 )
1119 if success or (not isinstance(success, bool) and success == 0): # ==0 because 0 could be the last time
1120 # step (will never be the case.... -.-)
1121 if not isinstance(success, bool):
1122 if not self.continue_canceled:
1123 self.rc_solution.delete_solution_except_first()
1124 print(f"Partial solution detected but not used because continue_canceled is False.")
1125 else:
1126 print(
1127 f"Simulation is started at loaded time step {success} on {success / t_span[-1] * 100:.2f} %."
1128 )
1129 t_span = (success, t_span[-1])
1130 continued_simulation = True
1131 else:
1132 print("Solution was loaded from file.")
1133 return None
1135 if check_courant:
1136 new_max_step = self.determine_max_step()
1137 if new_max_step is not None:
1138 if self.settings.solve_settings.max_step > new_max_step:
1139 print(
1140 f"Because of Courant Number: New maximum step is set to {new_max_step} (previous:"
1141 f" {self.settings.solve_settings.max_step})"
1142 )
1143 self.settings.solve_settings.max_step = new_max_step
1145 temperature_vector = np.array(self.temperature_vector).flatten()
1146 system_matrix = self.system_matrix
1147 time_symbols = self.time_dependent_system_symbols
1148 if time_symbols:
1149 system_matrix = SparseSymbolicEvaluator(system_matrix, time_symbols)
1151 system_handler_kwargs = {
1152 "system_matrix": system_matrix,
1153 "rc_solution": self.rc_solution,
1154 "settings": self.settings,
1155 "print_progress": print_progress,
1156 "print_points": np.linspace(t_span[0], t_span[-1], 20 + 1),
1157 "batch_end": t_span[-1],
1158 "time_dependent_function": time_dependent_function,
1159 }
1161 # create the time_steps to solve the results (every X seconds a value)
1162 t_eval = np.linspace(
1163 t_span[0], t_span[1], max(1, int((t_span[1] - t_span[0]) / self.settings.save_all_x_seconds) + 1)
1164 )
1166 if self.inhomogeneous_system:
1167 input_matrix = self.input_matrix
1168 if time_symbols:
1169 input_matrix = SparseSymbolicEvaluator(input_matrix, time_symbols)
1170 input_vector_free_symbols = self.variable_input_vector_symbols
1171 if input_vector_free_symbols:
1172 input_vector = ArraySymbolicEvaluator(self.input_vector, input_vector_free_symbols)
1173 else:
1174 input_vector = self.input_vector
1176 system_handler_kwargs.update(
1177 {
1178 "input_matrix": input_matrix,
1179 "input_vector": input_vector,
1180 "t_eval": t_eval,
1181 "first_time": t_span[0],
1182 "time_dependent_function_input": time_dependent_function_input,
1183 }
1184 )
1186 # only update the keys that are in the predefined dictionary (it would also work if everything is passed anyway)
1187 system_handler_kwargs.update({k: kwargs[k] for k in system_handler_kwargs if k in kwargs})
1189 system_handler = self.system_handler_type(**system_handler_kwargs)
1190 save_prefix = f"{self.hash}"
1191 if name_add_on:
1192 save_prefix += f"_{name_add_on.removeprefix('_').removesuffix('_')}"
1193 solve_handler = SolveIVPHandler(
1194 system_handler=system_handler,
1195 save_prefix=save_prefix,
1196 )
1197 if hook_function is not None:
1198 hook_function()
1199 solve_handler.solve(
1200 t_span=t_span,
1201 y0=temperature_vector,
1202 t_eval=t_eval,
1203 continued_simulation=continued_simulation,
1204 expected_solution_size_mb=expected_solution_size_mb,
1205 )
1207 def remap_symbol_function(
1208 self, symbols: Iterable[Symbol] | Symbol, function: Callable, wanted_order: Iterable[Symbol] | Symbol = None
1209 ) -> Callable:
1210 """
1211 Remap the functions return so that it maps self.get_time_dependent_symbols().
1212 The symbols iterable shows the order of the current return.
1214 Parameters
1215 ----------
1216 symbols : Iterable[Symbol] | Symbol
1217 The order of the current return of the function.
1218 function : Callable
1219 The function that calculates the values for the symbols in the symbols iterable.
1220 wanted_order : Iterable[Symbol] | Symbol, optional
1221 A list with the expected order of the function output.
1222 If None, self.get_time_dependent_symbols() is used.
1224 Returns
1225 -------
1226 Callable :
1227 The function which return was reordered.
1228 """
1229 if wanted_order is None:
1230 wanted_order = self.time_dependent_system_symbols
1231 if isinstance(wanted_order, Symbol):
1232 wanted_order = [wanted_order]
1233 if isinstance(symbols, Symbol):
1234 symbols = [symbols]
1235 pos = {s.name: i for i, s in enumerate(symbols)}
1236 mapping = [pos[s.name] for s in wanted_order]
1238 if len(mapping) == 1:
1240 def reordered_function(*args, **kwargs):
1241 return (function(*args, **kwargs)[mapping[0]],)
1242 else:
1243 getter = operator.itemgetter(*mapping)
1245 def reordered_function(*args, **kwargs):
1246 return getter(function(*args, **kwargs))
1248 return reordered_function
1250 @property
1251 def system_handler_type(self):
1252 """
1253 Returns the SystemHandler type depending on which system (in/homogeneous) is needed.
1255 Returns
1256 -------
1257 type :
1258 The system handler type.
1259 """
1260 if self.inhomogeneous_system:
1261 return InhomogeneousSystemHandler
1262 return HomogeneousSystemHandler
1264 @property
1265 def all_objects(self) -> list:
1266 return [
1267 *(self.rc_objects.resistors or []),
1268 *(self.rc_objects.capacities or []),
1269 *(self.rc_objects.inputs or []),
1270 ]
1272 @property
1273 def resistors_filtered_equivalent(self) -> list[Resistor]:
1274 """
1275 Get all resistor objects which equivalent resistor symbols are unique in defined order.
1277 In other words: If two Capacitors are connected over multiple Resistors (parallel or serial) only one of
1278 these is returned.
1280 The order is kept like in self.rc_objects.resistors, except all objects that are filtered out.
1282 Returns
1283 -------
1284 list[Resistor] :
1285 All resistor objects without doubled equivalent symbols.
1286 """
1287 resistor_list: list = self.rc_objects.resistors or []
1288 all_resistors = set(resistor_list)
1289 to_exclude = set()
1291 # Collect all resistors to exclude in order of appearance
1292 for obj in resistor_list:
1293 if obj not in to_exclude:
1294 inbetween = set(obj.all_resistors_inbetween) - {obj}
1295 to_exclude.update(inbetween & all_resistors)
1297 # Filter maintaining original order
1298 return [obj for obj in resistor_list if obj not in to_exclude]
1300 @property
1301 def system_objects_unique(self) -> list[Resistor | Capacitor]:
1302 """
1303 Returns all Resistor and Capacitor objects that are needed for system and B-matrix (using
1304 resistors_filtered_equivalent instead of all resistors).
1306 Returns
1307 -------
1308 list[Resistor|Capacitor] :
1309 All (equivalent representative) Resistors and Capacitor objects in the network.
1310 """
1311 return [*self.resistors_filtered_equivalent, *(self.rc_objects.capacities or [])]
1313 @property
1314 def system_symbols(self) -> list:
1315 """
1316 Returns a list with all R and C symbols for A- and B-matrix.
1318 Returns
1319 -------
1320 list :
1321 All R and C symbols for A- and B-matrix.
1322 """
1323 return [symb for rc_object in self.system_objects_unique for symb in rc_object.symbols]
1325 @property
1326 def system_values(self) -> list:
1327 """
1328 Returns a list with all R and C values for A- and B-matrix.
1330 Returns
1331 -------
1332 list :
1333 All R and C values for A- and B-matrix (might be sympy Expressions).
1334 """
1335 return [
1336 val if isinstance(val, Expr) else np.float64(val)
1337 for rc_object in self.system_objects_unique
1338 for val in rc_object.values
1339 ]
1341 @property
1342 def time_dependent_system_symbols(self) -> list:
1343 """
1344 A list with all time dependent symbols of the system (A-matrix).
1346 Returns
1347 -------
1348 list :
1349 The time dependent symbols of the system (appearing in system (A) matrix).
1350 """
1351 symbols = []
1352 seen = set()
1353 for rc_object in self.system_objects_unique:
1354 for symb in rc_object.time_dependent_symbols:
1355 if symb not in seen:
1356 seen.add(symb)
1357 symbols.append(symb)
1358 return symbols
1360 @property
1361 def time_dependent_input_symbols(self):
1362 """
1363 A list with all time dependent input symbols of the system (u vector).
1365 Returns
1366 -------
1367 list :
1368 """
1369 return self.variable_input_vector_symbols
1371 @property
1372 def variable_input_vector_symbols(self) -> list:
1373 """
1374 A list of all symbols in the input vector.
1376 Returns
1377 -------
1378 list :
1379 The symbols in the input vector.
1380 """
1381 symbols = []
1382 seen = set()
1383 for value in self.input_vector.flatten():
1384 if isinstance(value, Expr):
1385 for symbol in value.free_symbols:
1386 if symbol not in seen:
1387 seen.add(symbol)
1388 symbols.append(symbol)
1389 return symbols
1391 def get_symbols_and_values(self) -> tuple[list, list]:
1392 """
1393 Return two lists of first: all symbols and second: all associated values of all rc objects.
1395 Returns
1396 -------
1397 tuple[list, list] :
1399 """
1400 return self.system_symbols, self.system_values
1402 def make_system_matrices(self):
1403 """
1404 Creates the Jacobian matrices of the RC Network.
1406 If self.load_from_pickle the whole network will be loaded from pickle file if it exists to prevent heavy
1407 calculations.
1409 Returns
1410 -------
1411 sympy expression :
1412 The system matrix as sympy expression.
1413 """
1414 success = False
1415 if self.load_from_pickle:
1416 try:
1417 success = self.load_matrices()
1418 except Exception:
1419 pass
1420 if not success:
1421 # fill system matrix
1422 terms = []
1423 involved_symbols: list[set] = []
1424 temperature_symbols: list[set] = []
1425 variables = []
1426 for node in self.nodes:
1427 # works only for one term
1428 term, t_symbols, all_symbols = node.temperature_derivative_term()
1429 terms.append(term)
1430 involved_symbols.append(all_symbols)
1431 temperature_symbols.append(t_symbols)
1432 variables.append(node.temperature_symbol)
1434 self._system_matrix_symbol = build_jacobian(
1435 terms, variables, temperature_symbols, num_cores=self.num_cores_jacobian
1436 )
1437 print("system matrix created")
1439 self.save_matrices()
1441 input_variables = list(self.input_vector_symbol)
1442 self._input_matrix_symbol = build_jacobian(
1443 terms, input_variables, involved_symbols, num_cores=self.num_cores_jacobian
1444 )
1445 print("input matrix created")
1447 self.save_matrices()
1449 def check_courant(self, time_step: float):
1450 """
1451 Calculates the Courant number for the whole network and raises an error if its larger than 1.
1453 Parameters
1454 ----------
1455 time_step : float
1456 The maximum time step in seconds.
1458 Raises
1459 ------
1460 HighCourantNumberError :
1461 If the courant number is greater than 1 the network will calculate shit.
1463 """
1464 # TODO: Do not use the courant number but return a critical maximum time_step to easy raise an error during
1465 # iterations. So: no iteration, but calculate the maximum time_step and return this so that it can be set.
1466 mass_flow_nodes = self.rc_objects.mass_flow_nodes
1467 mass_flow_node: MassFlowNode
1468 for mass_flow_node in mass_flow_nodes:
1469 if mass_flow_node.courant_number(time_step) > 1:
1470 raise HighCourantNumberError
1472 def matrix_to_latex_diag(self):
1473 matrix = self.system_matrix_symbol
1474 diag_elements = matrix.diagonal()
1475 matrix_no_diag = matrix - diag(*diag_elements)
1477 matrix_latex = latex(matrix_no_diag)
1478 diag_latex = latex(Matrix(diag_elements).reshape(len(diag_elements), 1))
1480 return f"{matrix_latex} + \\mathrm{{diag}}\\left({diag_latex}\\right)"