Coverage for pyrc\core\solver\handler.py: 56%
248 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 os.path
11import time
12import warnings
13from collections import deque
14from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
15from datetime import datetime
16from multiprocessing import cpu_count
17from typing import Any, Callable
19import numpy as np
20from scipy.integrate import solve_ivp
21from scipy.sparse import sparray, spmatrix
23from pyrc.core.components.templates import RCSolution
24from pyrc.core.settings import Settings, SolveSettings, initial_settings
25from pyrc.core.solver.symbolic import ArraySymbolicEvaluator, SparseSymbolicEvaluator
26from pyrc.tools.science import get_free_ram_gb
29class HomogeneousSystemHandler:
30 def __init__(
31 self,
32 system_matrix: SparseSymbolicEvaluator | spmatrix | sparray,
33 rc_solution: RCSolution,
34 settings: Settings,
35 print_progress=True,
36 print_points: list | np.ndarray = None,
37 batch_end=None,
38 time_dependent_function=None,
39 _initialize_temperature_function=True,
40 **kwargs,
41 ):
42 """
43 System Handler that is called ("as function") by ``solve_ivp``.
45 Parameters
46 ----------
47 system_matrix : spmatrix | sparray | SparseSymbolicEvaluator
48 The system matrix.
49 If constant, the type must be a sparse scipy matrix.
50 If time dependent symbols are contained, the system matrix must be given in a `SparseSymbolicEvaluator`
51 object.
52 rc_solution
53 settings
54 print_progress
55 print_points
56 batch_end
57 time_dependent_function : Callable, optional
58 A function that calculates all time dependent variables within the time step and returns them in the same order as
59 RCNetwork.get_time_dependent_symbols().
60 It gets parameters like this:\n
61 ``time_dependent_function(time, temperature_vector, input_vector)``\\.\n
62 This function is required if time dependent symbols exist.
63 It must return an iterable (e.g. list).
64 _initialize_temperature_function : bool, optional
65 If ``False``, the ``dtemperature_dt_function`` attribute is not set using the method `_init_temperature_equation()`
66 Used in subclasses to prevent an early call of attributes that are created in the subclasses later on.
67 kwargs : dict, optional
68 Not used: Just to not raise an error if too much information is passed.
69 """
70 self.system_matrix: spmatrix | sparray | SparseSymbolicEvaluator = system_matrix
71 self.rc_solution = rc_solution
72 self.settings = settings
74 self.time_dependent_function: Callable | None = time_dependent_function # must return the value(s) as iterable
75 self.time_dependent_active: bool = True
76 if time_dependent_function is None:
77 self.time_dependent_active = False
78 self.time_dependent_function = lambda *args, **keyword_args: []
80 if print_progress and batch_end is None:
81 import warnings
83 warnings.warn("Print progress might not work as expected because batch_end value is not given.")
84 batch_end = np.inf
85 self.batch_end = batch_end
87 # print progress
88 if print_points is None:
89 print_points = [0]
90 if isinstance(print_points, np.ndarray):
91 print_points = print_points.tolist()
92 self.print_points: deque = deque(print_points)
93 self.next_printed_time_step = self.print_points.popleft()
94 self.print_progress = print_progress
96 # initialize the temperature equation using a lambda
97 self.dtemperature_dt_function = lambda *args, **keyword_args: None
98 if _initialize_temperature_function:
99 # this can be switched off because it always have to run at the end of all subclasses inits
100 self.dtemperature_dt_function = self._init_temperature_equation()
102 def _init_temperature_equation(self) -> Callable:
103 if isinstance(self.system_matrix, SparseSymbolicEvaluator):
104 assert self.time_dependent_active
105 def system_matrix(v):
106 return self.system_matrix.evaluate(v)
107 else:
108 if self.time_dependent_active:
109 warnings.warn(
110 "Function to calculate time dependent values passed, but no time dependent system matrix detected.\n"
111 "The passed function has no effect."
112 )
113 return lambda temperature: self.system_matrix @ temperature
114 return lambda temperature, v: system_matrix(v) @ temperature
116 def __call__(self, t, temperature):
117 """
118 The function the solve_ivp is going to solve.
120 Parameters
121 ----------
122 t
123 temperature
125 Notes
126 -----
127 The input vector that is saved during the iteration of the solver at the t_eval
129 Returns
130 -------
131 np.ndarray :
132 The resulting vector of the temperature derivative.
133 """
134 if self.print_progress:
135 self._update_progress_print(t)
137 if self.time_dependent_active:
138 time_dependent_values = self.time_dependent_function(t, temperature)
139 dtemperature_dt = self.dtemperature_dt_function(temperature.reshape(-1, 1), time_dependent_values)
140 else:
141 dtemperature_dt = self.dtemperature_dt_function(temperature.reshape(-1, 1))
143 return dtemperature_dt.flatten()
145 def set_new_t_eval(self, new_t_eval):
146 """
147 Change the current t_eval.
149 Parameters
150 ----------
151 new_t_eval : array_like
152 The new t_eval.
153 """
154 # Currently only used in the child class
155 pass
157 def _update_progress_print(self, t):
158 if t >= self.next_printed_time_step:
159 # prevent initialization print out (is not fail save, but okay)
160 if t != self.batch_end or (self.print_points and t == self.print_points[0]):
161 if len(self.print_points) > 0:
162 self.next_printed_time_step = self.print_points.popleft()
163 self.print_out_progress(t)
164 else:
165 self.next_printed_time_step = np.inf
166 # deactivate printing for better performance
167 self.print_progress = False
169 def print_out_progress(self, t):
170 """
171 Prints a formatted time stamp.
173 Parameters
174 ----------
175 t : float | int
176 The time to format/print.
177 """
178 print(f"Progress: t = {self.format_t(t)}")
180 @staticmethod
181 def format_t(t):
182 days = int(t) // 86400
183 hours = (int(t) % 86400) // 3600
184 minutes = (int(t) % 3600) // 60
185 seconds = int(t) % 60
186 return (
187 f"{days:>4} days, {hours:02}:{minutes:02}:{seconds:02} - current time: "
188 f"{datetime.now().strftime('%d %H:%M:%S')}"
189 )
192class InhomogeneousSystemHandler(HomogeneousSystemHandler):
193 def __init__(
194 self,
195 system_matrix,
196 input_matrix,
197 input_vector : np.ndarray | ArraySymbolicEvaluator,
198 rc_solution,
199 t_eval=None,
200 time_dependent_function=None,
201 time_dependent_function_input: Callable | None=None,
202 print_progress=True,
203 print_points=3600,
204 settings: Settings = initial_settings,
205 first_time: float | int = 0,
206 **kwargs,
207 ):
208 """
210 Parameters
211 ----------
212 system_matrix
213 input_matrix
214 input_vector : np.ndarray | ArraySymbolicEvaluator
215 rc_solution
216 t_eval
217 time_dependent_function : Callable, optional
218 A function that calculates all time dependent variables within the time step and returns them in the same order as
219 RCNetwork.get_time_dependent_symbols().
220 It gets parameters like this:\n
221 ``time_dependent_function(time, temperature_vector, input_vector)``\\.\n
222 This function is required if time dependent symbols exist.
223 It must return an iterable (e.g. list).
224 To not run into Errors just use ``*args``\\, ``**kwargs`` at the end in case more values are passed then
225 needed.
226 time_dependent_function_input : Callable, optional
227 A function that calculates all time dependent variables within the time step and returns them in the same order as
228 RCNetwork.variable_input_vector_symbols.
229 It gets parameters like this:\n
230 ``time_dependent_function(time, temperature_vector)`` \n
231 This function is required if time dependent symbols exist in the input vector.
232 It must return an iterable (e.g. list).
233 print_progress
234 print_points
235 settings
236 first_time
237 kwargs
238 """
239 super().__init__(
240 system_matrix,
241 rc_solution,
242 settings,
243 time_dependent_function=time_dependent_function,
244 print_progress=print_progress,
245 print_points=print_points,
246 batch_end=t_eval[-1],
247 _initialize_temperature_function=False,
248 )
250 self.time_dependent_function_input: Callable | None = time_dependent_function_input # must return the value(s) as iterable
251 self.time_dependent_input_active: bool = True
252 if time_dependent_function_input is None:
253 self.time_dependent_input_active = False
254 self.time_dependent_function_input = lambda *args, **keyword_args: []
256 self.input_matrix = input_matrix
257 self.input_vector_function = None
258 if isinstance(input_vector, np.ndarray):
259 self.input_vector: np.ndarray = input_vector.reshape(-1, 1)
260 elif isinstance(input_vector, ArraySymbolicEvaluator):
261 self.input_vector_function: Callable = lambda v: input_vector.evaluate(v)
262 self.input_vector: np.ndarray = np.array([])
263 else:
264 raise ValueError("input vector must be a numpy array or an ArraySymbolicEvaluator")
266 self.after_iteration = np.inf
267 self.next_eval_time = 0
268 self.t_eval_iter = iter(())
269 self.set_new_t_eval(t_eval)
271 self.last_t = first_time
272 self.last_eval_time = 0
273 self.current_eval_time = 0
274 self.use_current_eval_time = False
276 # initialize the temperature equation using a lambda
277 self.dtemperature_dt_function = self._init_temperature_equation()
279 def _init_temperature_equation(self) -> Callable:
280 if isinstance(self.system_matrix, SparseSymbolicEvaluator):
281 def system_matrix_fun(v):
282 return self.system_matrix.evaluate(v)
283 else:
284 def system_matrix_fun(*args, **kwargs):
285 return self.system_matrix
286 if isinstance(self.input_matrix, SparseSymbolicEvaluator):
287 def input_matrix_fun(v):
288 return self.input_matrix.evaluate(v)
289 else:
290 def input_matrix_fun(*args, **kwargs):
291 return self.input_matrix
292 return lambda temperature, input_v, v: (system_matrix_fun(v) @ temperature + input_matrix_fun(v) @ input_v)
294 def set_new_t_eval(self, new_t_eval):
295 """
296 Change the current t_eval.
298 Parameters
299 ----------
300 new_t_eval : array_like
301 The new t_eval.
302 """
303 if new_t_eval is None:
304 self.t_eval_iter = iter([-np.inf])
305 self.after_iteration = -np.inf
306 else:
307 self.t_eval_iter = iter(new_t_eval)
308 self.after_iteration = np.inf
309 self.next_eval_time = next(self.t_eval_iter, self.after_iteration)
310 self.use_current_eval_time = False
312 @staticmethod
313 def _worker_call(args):
314 f, t, temperature, input_vector, kwargs = args
315 return f(t, temperature, input_vector, **kwargs)
317 def __call__(self, t, temperature: np.ndarray):
318 """
319 The function the solve_ivp is going to solve.
321 Parameters
322 ----------
323 t
324 temperature
326 Notes
327 -----
328 The input vector that is saved during the iteration of the solver at the t_eval
330 Returns
331 -------
332 np.ndarray :
333 The resulting vector of the temperature derivative.
334 """
335 if t < self.last_t:
336 if t == self.batch_end or (self.last_t >= self.last_eval_time > t):
337 # Delete the last result because the solver is in initialization or one time step was saved to early.
338 # Revert the next_eval_time to the last value.
339 self.next_eval_time = self.last_eval_time
340 self.use_current_eval_time = True
341 self.rc_solution.delete_last_input()
342 if self.print_progress:
343 self._update_progress_print(t)
345 if self.time_dependent_input_active:
346 self.input_vector = self.input_vector_function(self.time_dependent_function_input(t, temperature))
348 time_dependent_values = self.time_dependent_function(
349 t,
350 temperature,
351 self.input_vector.reshape(-1,),
352 )
354 # Check if the input vector has to be saved.
355 if t >= self.next_eval_time:
356 # NOTE: This is not perfect, because the input vector of the next time step is saved for the result of
357 # the last time step. However, the solver will perform such small iteration steps that this will not
358 # become a problem.
359 self.rc_solution.append_to_input(self.input_vector)
360 self.last_eval_time = self.next_eval_time
361 if self.use_current_eval_time:
362 self.next_eval_time = self.current_eval_time
363 self.use_current_eval_time = False
364 else:
365 self.next_eval_time = next(self.t_eval_iter, self.after_iteration)
367 dtemperature_dt = self.dtemperature_dt_function(
368 temperature.reshape(-1, 1), self.input_vector, time_dependent_values
369 )
371 self.last_t = t
372 return dtemperature_dt.flatten()
375class SolveIVPHandler:
376 def __init__(
377 self,
378 system_handler: HomogeneousSystemHandler | InhomogeneousSystemHandler,
379 max_saved_steps=None,
380 method=None,
381 max_step=None,
382 rtol=None,
383 atol=None,
384 save_interval=None,
385 save_path=None,
386 save_prefix="",
387 minimize_ram_usage=None,
388 **kwargs,
389 ):
390 """
391 Handler of the solving process with solve_ivp.
393 Parameters
394 ----------
395 system_handler : HomogeneousSystemHandler | InhomogeneousSystemHandler
396 The system handler that holds the function to solve in __call__().
397 max_saved_steps : int, optional
398 The maximum number of used seconds during one solve_ivp call.
399 It defines the batch size in seconds.
400 Using this prevents long solving time because the matrices become very big.
401 method : str, optional
402 The method to use to solve the system (see scipy.solve_ivp).
403 max_step
404 rtol
405 atol
406 save_interval
407 save_path
408 save_prefix : str, optional
409 This is the beginning of the name of the pickle file that is saved during the solving.
410 minimize_ram_usage : bool, optional
411 If True, the solution is deleted if saved to minimize the RAM usage during runtime.
412 kwargs
413 """
414 self.solve_settings = system_handler.settings.solve_settings
415 max_saved_steps, method, max_step, rtol, atol, save_interval, minimize_ram_usage = self.set_initial_values(
416 max_saved_steps=max_saved_steps,
417 method=method,
418 max_step=max_step,
419 rtol=rtol,
420 atol=atol,
421 save_interval=save_interval,
422 minimize_ram_usage=minimize_ram_usage,
423 )
424 self.system_handler: HomogeneousSystemHandler | InhomogeneousSystemHandler = system_handler
425 self.max_saved_steps = int(max_saved_steps)
426 self.method = method
427 self.max_step = max_step
428 self.rtol = rtol
429 self.atol = atol
431 self.save_interval = save_interval
432 self.save_counter = 0
433 if save_path is None:
434 save_path = self.system_handler.settings.save_folder_path
435 if save_path is None:
436 self.save_path = None
437 else:
438 self.save_path = os.path.normpath(save_path)
439 self.save_prefix = save_prefix
441 if self.save_path is None and minimize_ram_usage:
442 print("Minimize RAM usage deactivated because no save_path is declared in Settings.")
443 print("This also deactivates every save during solving.")
444 minimize_ram_usage = False
445 self.minimize_ram_usage = minimize_ram_usage
447 self.kwargs = kwargs
449 def set_initial_values(self, **kwargs):
450 result = []
451 settings_dict: dict = self.solve_settings.dict
452 for key, arg in zip(kwargs.keys(), kwargs.values()):
453 if arg is None:
454 result.append(settings_dict[key])
455 else:
456 result.append(arg)
457 print(f"Solve setting {key} is overwritten by manual/coded value: {arg}")
458 return result
460 def get_batches(self, t_span):
461 start, end = t_span
462 splits = np.arange(start, end, self.max_saved_steps)
463 if splits[-1] != end:
464 splits = np.append(splits, end)
465 return splits
467 def solve(
468 self,
469 t_span,
470 y0,
471 t_eval=None,
472 continued_simulation: bool = False,
473 expected_solution_size_mb=5000,
474 ):
475 """
476 Runs the solve_ivp in batches.
478 Parameters
479 ----------
480 t_span : tuple[int | float, int| float] :
481 The start and end of the simulation in seconds.
482 Usually it starts at 0: (0, end_seconds)
483 See also: scipy.integrate.solve_ivp()
484 y0 : np.ndarray | Iterable | float | int
485 The initial state of the system.
486 See also: scipy.integrate.solve_ivp()
487 t_eval : np.ndarray | Iterable | float | int, optional
488 Times at which to store the computed solution, must be sorted and lie within `t_span`\\.
489 If None (default), use points selected by the solver.
490 See also: scipy.integrate.solve_ivp()
491 continued_simulation : bool, optional
492 If True, the first value of the first batch is not kept, because the simulation is continued.
493 expected_solution_size_mb : int | float, optional
494 The expected solution size in Megabytes. Is used for RAM-Management: if not enough free memory is available,
495 the method waits for 10 seconds and tries again for a total of 360 times before raising an error.
496 The solution size depends on the size of the RC network and number of saved time steps.
497 """
498 batches = self.get_batches(t_span)
500 list_t, list_y = [], []
501 current_y0 = y0
503 if self.save_path is not None:
504 intermediate_save_prefix = os.path.join(self.save_path, self.save_prefix + f"_{int(t_span[-1])}_")
505 else:
506 intermediate_save_prefix = "dummy_path"
508 for i in range(len(batches) - 1):
509 batch_start, batch_end = float(batches[i]), float(batches[i + 1])
510 if t_eval is not None:
511 if i == 0 and not continued_simulation:
512 mask = (t_eval >= batch_start) & (t_eval <= batch_end)
513 else:
514 # exclude the batch_start for every batch except the first one to prevent duplicates
515 mask = (t_eval > batch_start) & (t_eval <= batch_end)
516 t_eval_batch = t_eval[mask]
517 else:
518 t_eval_batch = None
520 self.system_handler.set_new_t_eval(t_eval_batch)
522 last_step_not_in_solution = False
523 if t_eval_batch.size == 0 or t_eval_batch[-1] != batch_end:
524 # The last step must be returned by solve_ivp anyway to pass it as new y0 in the next batch
525 t_eval_batch = np.append(t_eval_batch, batch_end)
526 last_step_not_in_solution = True
528 # update end of batch (used for correct print out and saves)
529 self.system_handler.batch_end = batch_end
531 sol: Any = solve_ivp(
532 fun=self.system_handler,
533 t_span=(batch_start, batch_end),
534 y0=current_y0,
535 method=self.method,
536 t_eval=t_eval_batch,
537 max_step=self.max_step,
538 rtol=self.rtol,
539 atol=self.atol,
540 **self.kwargs,
541 )
542 current_y0 = sol.y[:, -1]
544 if last_step_not_in_solution:
545 # delete last solution if not requested in settings (but it was needed for y0)
546 self.system_handler.rc_solution.delete_last_input()
547 new_time_steps = sol.t[:-1]
548 new_y_values = sol.y[:, :-1]
549 else:
550 new_time_steps = sol.t
551 new_y_values = sol.y
553 if new_time_steps.size != 0:
554 list_t.append(new_time_steps)
555 list_y.append(new_y_values)
556 # increase counter only if new solution was added. Otherwise it saves the same solution / an empty one
557 self.save_counter += 1
559 if self.save_counter // self.save_interval > 0 or (
560 batch_end == t_span[-1]
561 and (self.system_handler.rc_solution.last_saved_timestep_index > 0 or self.minimize_ram_usage)
562 ):
563 save_path = f"{intermediate_save_prefix}{float(batch_end):09.0f}_s.pickle"
564 if self.minimize_ram_usage:
565 # save the solution and delete it afterward
566 assert self.save_path is not None
567 self.system_handler.rc_solution.save_to_file_only(
568 t=np.concatenate(list_t),
569 y=np.concatenate(list_y, axis=1).T,
570 path_with_name_and_ending=save_path,
571 )
572 else:
573 self.system_handler.rc_solution.add_to_solution(list_t, list_y)
574 if self.save_path is not None:
575 self.system_handler.rc_solution.save_new_solution(save_path)
577 list_y = []
578 list_t = []
579 self.save_counter = 0
581 # Make print out for the last time step that else will be missed
582 if self.system_handler.print_progress:
583 if self.system_handler.next_printed_time_step <= t_span[-1]:
584 self.system_handler.print_out_progress(t_span[-1])
586 save_path = f"{intermediate_save_prefix}result.pickle"
587 if self.minimize_ram_usage:
588 # load solution in and save it as one big solution
589 # it will load all intermediate saves because it will not find the requested path with "_result.pickle"
590 # first check, if enough RAM is available
591 assert self.save_path is not None
592 loop_counter = 0
593 max_loops = 360
594 while get_free_ram_gb() < expected_solution_size_mb / 1000 and loop_counter < max_loops:
595 loop_counter += 1
596 time.sleep(10) # waiting for more free RAM
597 if not loop_counter >= max_loops:
598 self.system_handler.rc_solution.load_solution(save_path)
599 # Now save the accumulated result in one pickle with the "_result" ending
600 self.system_handler.rc_solution.save_solution(save_path)
601 else:
602 print("Not enough RAM for over one hour. Incremental solutions are not combined.")
603 else:
604 self.system_handler.rc_solution.add_to_solution(list_t, list_y)
605 if self.save_path is not None:
606 self.system_handler.rc_solution.save_solution(save_path)