Coverage for pyrc\core\nodes.py: 28%

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 

8from __future__ import annotations 

9 

10import warnings 

11from abc import ABC 

12from collections import deque 

13from typing import TYPE_CHECKING 

14 

15import numpy as np 

16 

17from pyrc.core.components.capacitor import Capacitor 

18from pyrc.core.components.templates import ( 

19 Cell, 

20 ConnectedFlowObject, 

21 Fluid, 

22 Geometric, 

23 Material, 

24 ObjectWithPorts, 

25 RCObjects, 

26 RCSolution, 

27 Solid, 

28 calculate_balance_for_resistors, 

29 initial_rc_objects, 

30 solution_object, 

31) 

32from pyrc.core.inputs import BoundaryCondition, FlowBoundaryCondition, InternalHeatSource 

33from pyrc.core.materials import Air 

34 

35if TYPE_CHECKING: 

36 from pyrc.core.components.resistor import Resistor 

37 from pyrc.core.resistors import MassTransport 

38 from pyrc.core.settings import Settings 

39 

40 

41class ConnectedFlowCapacitorObject(Capacitor, ConnectedFlowObject, ABC): 

42 pass 

43 

44 

45class Node(Capacitor, Cell): 

46 def __init__( 

47 self, 

48 material: Material, 

49 temperature: float | int, 

50 position: np.ndarray | tuple, 

51 temperature_derivative: float | int = 0, 

52 internal_heat_source: bool | int | float = False, 

53 internal_heat_source_type=InternalHeatSource, 

54 delta: np.ndarray | tuple | list = None, 

55 rc_objects=initial_rc_objects, 

56 rc_solution: RCSolution = solution_object, 

57 **kwargs: float | int | np.ndarray | Settings, 

58 ): 

59 """ 

60 

61 Parameters 

62 ---------- 

63 material : Material 

64 Container with all material properties stored in the node is made up from. 

65 temperature : float | int 

66 The temperature of the node. 

67 It is recommended to use the SI unit Kelvin instead of degrees Celsius. 

68 position : np.ndarray | tuple 

69 The position of the node in 2D/3D space. 

70 If 2D, a zero is added for the z coordinate. 

71 temperature_derivative : float | int, default=0 

72 The temperature derivative of the node. 

73 internal_heat_source : bool | int | float, default=False 

74 If True, an internal heat source term is added. 

75 It can be accessed using Node.internal_heat_source 

76 If a number, the number is used to set as volumetric heat source power in W/m³. 

77 internal_heat_source_type : type, default=InternalHeatSource 

78 The internal heat source type. 

79 Can be each :class:`InternalHeatSource` type. 

80 delta : np.ndarray | tuple, optional 

81 Delta vector [delta_x, delta_y, delta_z]. 

82 kwargs : dict 

83 Parameters passed to `InternalHeatSource`/its type. 

84 

85 Notes 

86 ----- 

87 Be aware of the fact that the **Node's properties can change during network creation!** 

88 This means, that if you add connections to a Node it's volume/mass can change. So every object, that needs 

89 the Node's volume/mass should be initialized **after** the whole network was build up. This affects, 

90 among other things, :class:`Input`\\s like :class:`InternalHeatSource`\\s, especially :class:`Radiation`\\. 

91 To avoid this, you must pass the whole object (node/BC) to the Input and program it in a way that the values 

92 are created earliest when values are requested (and therefore the network must be set up (hopefully)). 

93 """ 

94 Capacitor.__init__( 

95 self, 

96 capacity=None, 

97 temperature=temperature, 

98 temperature_derivative=temperature_derivative, 

99 rc_objects=rc_objects, 

100 rc_solution=rc_solution, 

101 ) 

102 Cell.__init__(self, position=position, delta=delta) 

103 self.material: Material | Fluid | Solid = material 

104 self._volume = None 

105 self._kwargs = kwargs # saved for __copy__ 

106 

107 # This must be at the end of init. It uses some attributes of node already (like node.id) 

108 if internal_heat_source: 

109 if not isinstance(internal_heat_source, bool): 

110 internal_heat_source = internal_heat_source_type( 

111 node=self, specific_power_in_w_per_cubic_meter=internal_heat_source, **kwargs 

112 ) 

113 elif internal_heat_source_type is not None: 

114 internal_heat_source = internal_heat_source_type(node=self, **kwargs) 

115 else: 

116 internal_heat_source = InternalHeatSource(node=self, **kwargs) 

117 self.add_internal_heat_source(internal_heat_source) 

118 

119 def __copy__(self): 

120 # TODO: This was not checked well. Have a closer look on it. 

121 cls = self.__class__ 

122 internal_heat_source = False if self.internal_heat_source is None else True 

123 internal_heat_source_type = InternalHeatSource 

124 if internal_heat_source: 

125 internal_heat_source_type = type(self.internal_heat_source) 

126 if isinstance(self.internal_heat_source, InternalHeatSource): 

127 Warning( 

128 "A Node with an InternalHeatSource is copied. The attribute\n" 

129 "\tspecific_power_in_w_per_meter_squared\n" 

130 "might be lost." 

131 ) 

132 if hasattr(self.internal_heat_source, "volume_specific_power"): 

133 internal_heat_source = self.internal_heat_source.volume_specific_power 

134 return cls( 

135 material=self.material, 

136 temperature=self.temperature, 

137 position=self.position, 

138 temperature_derivative=self.temperature_derivative, 

139 internal_heat_source=internal_heat_source, 

140 internal_heat_source_type=internal_heat_source_type, 

141 delta=self.delta, 

142 rc_objects=self.rc_objects, 

143 rc_solution=self.solutions, 

144 **self._kwargs, 

145 ) 

146 

147 def __str__(self): 

148 return self.__repr__() 

149 

150 def __repr__(self): 

151 return f"{self.__class__.__name__} {self.id}: {self.material.__class__.__name__}; ϑ={self.temperature}" 

152 

153 @property 

154 def capacity(self) -> float | int: 

155 """ 

156 Returns the capacity (not specific) of the node in J/K. 

157 

158 Returns 

159 ------- 

160 float | int : 

161 The capacity of the node in J/K. 

162 """ 

163 if Capacitor.capacity.fget(self) is None: 

164 self.calculate_capacity() 

165 return Capacitor.capacity.fget(self) 

166 

167 @property 

168 def volume(self) -> float | int: 

169 if self._volume is None: 

170 self.calculate_volume() 

171 return self._volume 

172 

173 @property 

174 def mass(self) -> float | int | np.number: 

175 """ 

176 Returns the mass flow of the node in kg. 

177 

178 Returns 

179 ------- 

180 float | int : 

181 The mass flow of the node in kg. 

182 If the density isn't defined np.nan is returned. 

183 """ 

184 if self.material.density: 

185 return self.volume * self.material.density 

186 return np.nan 

187 

188 def update_color(self, t_min: float = 263.15, t_max: float = 413.15, colormap="managua") -> None: 

189 """ 

190 Update the color of the vbox for visualization. 

191 

192 Parameters 

193 ---------- 

194 t_min : float | int, default=263.15 

195 The minimal temperature for the color code. 

196 t_max : float | int, default=413.15 

197 The maximal temperature for the color code. 

198 colormap : str, default="managua" 

199 The colormap to use. See pyrc.core.visualization.color.color.py or the txt files in there respectively. 

200 """ 

201 Cell.update_color(self.temperature, t_min=t_min, t_max=t_max, colormap=colormap) 

202 

203 def area(self, normal: np.ndarray) -> float: 

204 """ 

205 Returns the area of the desired plane in m**2. 

206 

207 Parameters 

208 ---------- 

209 normal : np.ndarray 

210 The direction to the wanted plane. 

211 E.g.: if (1,0,0) the y-z-plane is wanted so the area of the yz-plane is returned. 

212 

213 Returns 

214 ------- 

215 float : 

216 The area in m^2. 

217 """ 

218 direction_abs = np.abs(normal) 

219 

220 multiplication_vector = np.array([1, 1, 1]) - direction_abs 

221 area = np.prod(self.delta, where=multiplication_vector.astype(bool)) 

222 return float(area.item()) 

223 

224 def calculate_capacity(self): 

225 Capacitor.capacity.fset(self, self.material.heat_capacity * self.mass) # J/K 

226 

227 def calculate_volume(self): 

228 self._volume = self.delta_x * self.delta_y * self.delta_z # m^3 

229 

230 # add volume for each channel node 

231 channel_nodes = self.get_connected_nodes(ChannelNode) 

232 channel_node: ChannelNode 

233 for channel_node in channel_nodes: 

234 raw_volume = channel_node.delta_x * channel_node.delta_y * channel_node.delta_z 

235 solid_volume = raw_volume - channel_node.volume 

236 self._volume += solid_volume / 4 

237 

238 def get_area(self, asking_node: Node | Cell | BoundaryCondition | Capacitor) -> float | int | np.number: 

239 """ 

240 Returns the area that the asking node sees. Through this area travels the heat between both nodes. 

241 

242 Example: If the asking node has the same position in x and y direction but is shifted in z direction, 

243 the area of the cell face X-Y is returned. 

244 

245 A non-Cell class is only valid if its direction to/from self was set manually in 

246 self.manual_directions[asking_node]. Otherwise, it will throw an error. 

247 

248 Parameters 

249 ---------- 

250 asking_node : Node | Cell | BoundaryCondition | Capacitor 

251 The Node that asks for the area of this Node. It is used to calculate the shift in position and to return 

252 the associated surface area accordingly. 

253 If it doesn't inherit from Cell, it is assumed that the direction is set manually in manual_directions dict 

254 using ``self`` as key. 

255 

256 Returns 

257 ------- 

258 float | int | np.number : 

259 The area of the surface of the cell that has the normal vector of the line you get connecting both `self` 

260 and `asking_node`. 

261 """ 

262 # determine the shift in position between self and asking_node 

263 direction_to_asking_node = self.get_direction(asking_node=asking_node) 

264 

265 return self.area(direction_to_asking_node) 

266 

267 def get_conduction_length(self, asking_node: Node | ConnectedFlowObject | Geometric | Capacitor): 

268 """ 

269 Returns the length to the asking node. It's half the delta value in direction to the asking node. 

270 

271 Parameters 

272 ---------- 

273 asking_node : Geometric 

274 The Node asking for the (half) length between itself and self. 

275 

276 Returns 

277 ------- 

278 float : 

279 The half-length of self. Or: Half the delta value of self pointing towards asking_node. 

280 """ 

281 direction = self.get_direction(asking_node) # returns the direction to (not from) the asking node. 

282 

283 return np.linalg.norm(direction * self.delta) * 0.5 

284 

285 def get_both_perpendicular_lengths(self, asking_node: Geometric | Node | BoundaryCondition) -> tuple[float, float]: 

286 """ 

287 Returns both perpendicular lengths between the self and the asking node. 

288 

289 Parameters 

290 ---------- 

291 asking_node : Geometric | Node | BoundaryCondition 

292 The node asking. 

293 

294 Returns 

295 ------- 

296 tuple[float, float] : 

297 Both length in no particular order. 

298 """ 

299 # get the direction to the asking node 

300 direction = self.get_direction(asking_node) 

301 

302 # Define the mapping of axes 

303 axis_map = { 

304 (1, 0, 0): ["delta_y", "delta_z"], # If direction is X, take Y and Z 

305 (0, 1, 0): ["delta_x", "delta_z"], # If direction is Y, take X and Z 

306 (0, 0, 1): ["delta_x", "delta_y"], # If direction is Z, take X and Y 

307 } 

308 

309 # Convert direction to tuple 

310 direction_tuple = tuple(np.abs(direction).astype(int)) 

311 

312 if direction_tuple not in axis_map: 

313 return np.nan, np.nan # Shouldn't happen if input is correct 

314 

315 # Extract the two perpendicular delta values 

316 delta_1 = getattr(self, axis_map[direction_tuple][0]) 

317 delta_2 = getattr(self, axis_map[direction_tuple][1]) 

318 

319 return delta_1, delta_2 

320 

321 def get_smaller_edge(self, asking_node: Geometric | Node | BoundaryCondition | Capacitor): 

322 """ 

323 Returns the smaller edge of the cell using the asking node to set the direction. 

324 

325 Used to determine the diameter in `MassFlowNode` s. 

326 

327 Parameters 

328 ---------- 

329 asking_node : Geometric | Node | BoundaryCondition 

330 The node that asks for the diameter of this `Node`. It determines the direction to get the right dimension. 

331 

332 Returns 

333 ------- 

334 float | int | np.number : 

335 The diameter (smallest length of the cell side pointing towards `asking_node`). 

336 """ 

337 deltas = self.get_both_perpendicular_lengths(asking_node=asking_node) 

338 

339 # Return the smaller of the two 

340 return min(deltas) 

341 

342 

343class MassFlowNode(Node, ConnectedFlowObject): 

344 def __init__( 

345 self, 

346 material: Material, 

347 temperature, 

348 position, 

349 flow_area: float | int | str = None, 

350 flow_length: float | int | str = None, 

351 temperature_derivative=0, 

352 flow_velocity=None, 

353 mass_flow=None, 

354 volume_flow=None, 

355 rc_objects: RCObjects = initial_rc_objects, 

356 rc_solution: RCSolution = solution_object, 

357 **kwargs, 

358 ): 

359 """ 

360 

361 Parameters 

362 ---------- 

363 temperature : float | int | np.number 

364 The (initial) temperature of the node. 

365 position : np.ndarray 

366 The position of the node in 2D/3D space. 

367 If 2D, a zero is added for the z coordinate. 

368 flow_area : float | int | str 

369 Defining the area that is used for flow calculations (determining the velocity and the Courant number). 

370 Can be a string: 

371 xy: xy-plane, xz: xz-plane, yz: yz-plane 

372 If not given the smallest flow area is used that points to a `ConnectedFlowObject`. 

373 flow_length : float | int | str 

374 Defining the length that the flow takes. 

375 Can be a string: 

376 x: delta_x, y: delta_y, z: delta_z 

377 If not given the matching one for the flow area is used. 

378 temperature_derivative : float | int | np.number 

379 The temperature derivative of the node. 

380 flow_velocity : float | int | np.number 

381 The flow velocity of the node. 

382 It doesn't matter which sign it has. Only the absolute value is used. 

383 """ 

384 ConnectedFlowObject.__init__(self) 

385 Node.__init__( 

386 self, 

387 material=material, 

388 temperature=temperature, 

389 position=position, 

390 temperature_derivative=temperature_derivative, 

391 rc_objects=rc_objects, 

392 rc_solution=rc_solution, 

393 **kwargs, 

394 ) 

395 # Compute area if it's given as a plane 

396 self._flow_length = None 

397 self._flow_area = None 

398 if isinstance(flow_length, str): 

399 self._flow_length = getattr(self, f"delta_{flow_length}") 

400 else: 

401 self._flow_length = flow_length 

402 if isinstance(flow_area, str): 

403 match flow_area.lower(): 

404 case "xy" | "yx": 

405 self._flow_area = self.delta_x * self.delta_y 

406 self._flow_length = self.delta_z 

407 case "xz" | "zx": 

408 self._flow_area = self.delta_x * self.delta_z 

409 self._flow_length = self.delta_y 

410 case "yz" | "zy": 

411 self._flow_area = self.delta_y * self.delta_z 

412 self._flow_length = self.delta_x 

413 case _: 

414 raise ValueError("Invalid plane. Use 'xy', 'xz', or 'yz'.") 

415 else: 

416 self._flow_area = flow_area # Assume it's given as a direct value in m². 

417 # TODO: set flow_area and flow_direction depending on each other. 

418 

419 self._flow_propagation = False 

420 

421 if sum(val is not None for val in [flow_velocity, mass_flow, volume_flow]) != 1: 

422 self._flow_propagation = True 

423 

424 self._passed_flow_velocity = flow_velocity 

425 

426 if volume_flow is not None: 

427 self._volume_flow = volume_flow 

428 elif flow_velocity is not None: 

429 # calculate later because the area isn't known until all connections were made 

430 self._volume_flow = None 

431 elif mass_flow is not None: 

432 self._volume_flow = mass_flow / self.material.density 

433 

434 # Initialize cache 

435 self._velocity = None 

436 self._mass_flow = None 

437 self._flow_direction = None 

438 

439 def courant_number(self, time_step: float): 

440 """ 

441 Calculates the courant number of the node. 

442 

443 Parameters 

444 ---------- 

445 time_step : float | int | np.number 

446 The (maximum) time step used for calculation. 

447 

448 Returns 

449 ------- 

450 float | np.number : 

451 The courant number. 

452 """ 

453 return self.velocity * time_step / self.flow_length 

454 

455 @property 

456 def flow_length(self) -> float: 

457 """ 

458 The effective length of the channel. Use only if no asking_node is given! 

459 

460 This assumes that each ChannelNode is connected to at least one solid node at its sides. This node is used to 

461 get the perpendicular direction and define which sides defining the length. 

462 

463 Returns 

464 ------- 

465 float : 

466 The length, but it can be wrong in few cases. 

467 """ 

468 if self._flow_length is None: 

469 # request flow_area and calculate the value in it 

470 _ = self.flow_area 

471 return self._flow_length 

472 

473 @property 

474 def flow_area(self): 

475 if self._flow_area is None or self._flow_length is None: 

476 # get the areas between self and ConnectedFlowObjects and take the smallest one. 

477 result = np.inf 

478 flow_length = 0 

479 flow_obj: MassFlowNode | BoundaryCondition 

480 for flow_obj in self.mass_flow_neighbours: 

481 area = self.get_area(flow_obj) 

482 if result > area: 

483 result = area 

484 flow_length = self.get_direction(flow_obj) * self.delta 

485 self._flow_area = result 

486 self._flow_length = np.linalg.norm(flow_length) 

487 return self._flow_area 

488 

489 @property 

490 def volume_flow(self): 

491 """Returns the cached volume flow rate (m³/s).""" 

492 if self._flow_propagation and self._volume_flow is None: 

493 self.propagate_flow() 

494 # volume flow should be set now 

495 self._flow_propagation = False 

496 if self._volume_flow is None: 

497 if self._passed_flow_velocity: 

498 self._volume_flow = self._passed_flow_velocity * self.flow_area 

499 else: 

500 raise ValueError("This shouldn't have happened!") 

501 return self._volume_flow 

502 

503 @property 

504 def velocity(self): 

505 """Returns the cached flow velocity (m/s) or calculates it if needed.""" 

506 if self._velocity is None: 

507 self._velocity = self.volume_flow / self.flow_area 

508 return self._velocity 

509 

510 @property 

511 def mass_flow(self): 

512 """Returns the cached mass flow rate (kg/s) or calculates it if needed.""" 

513 if self._mass_flow is None: 

514 self._mass_flow = self.volume_flow * self.material.density 

515 return self._mass_flow 

516 

517 def reset_properties(self): 

518 """Clears cached flow velocity and mass flow to force recalculation.""" 

519 super().reset_properties() 

520 self._velocity = None 

521 self._mass_flow = None 

522 

523 @property 

524 def sources(self) -> list[Capacitor | Node]: 

525 """ 

526 Returns a list with all nodes where volume flow is flowing into the node. 

527 

528 No check is performed which node type is put into the list, but it should only be `MassFlowNode` s and its 

529 descendants or `FlowBoundaryCondition` s. 

530 

531 Returns 

532 ------- 

533 list : 

534 All nodes serving as volume flow source. 

535 """ 

536 result = [] 

537 seen = set() 

538 resistor: MassTransport 

539 for resistor in self.connected_mass_transport_resistors(): 

540 if resistor.sink == self: 

541 connected = resistor.get_connected_node(self) 

542 if connected not in seen: 

543 seen.add(connected) 

544 result.append(connected) 

545 return result 

546 

547 @property 

548 def sinks(self) -> list[MassFlowNode | FlowBoundaryCondition]: 

549 """ 

550 Returns a list of all `MassFlowNode` s into which the flow is flowing. 

551 

552 No check is performed which node type is put into the list, but it should only be `MassFlowNode` s and its 

553 descendants or `FlowBoundaryCondition` s. 

554 

555 Returns 

556 ------- 

557 list : 

558 All nodes serving as volume flow sink. 

559 """ 

560 result = [] 

561 seen = set() 

562 resistor: MassTransport 

563 for resistor in self.connected_mass_transport_resistors(): 

564 if resistor.source == self: 

565 connected = resistor.get_connected_node(self) 

566 if connected not in seen: 

567 seen.add(connected) 

568 result.append(connected) 

569 return result 

570 

571 @property 

572 def diameter(self) -> float: 

573 """ 

574 Returns the hydraulic diameter. 

575 

576 Returns 

577 ------- 

578 float : 

579 The hydraulic diameter of the node. 

580 """ 

581 len1, len2 = self.get_minimal_perpendicular_lengths_to_flow() 

582 len1 = float(len1) 

583 len2 = float(len2) 

584 # return hydraulic diameter 

585 return 2 * len1 * len2 / (len1 + len2) 

586 

587 def get_resistor_volume_flow(self, resistor: MassTransport): 

588 """ 

589 Returns the volume flow for the given resistor using mass balance (assuming a constant density). 

590 

591 Parameters 

592 ---------- 

593 resistor : MassTransport 

594 The `Resistor` to determine the volume flow for. 

595 

596 Returns 

597 ------- 

598 float : 

599 The volume flow of the given `Resistor`. 

600 """ 

601 result = 0 

602 other_node_on_resistor = resistor.get_connected_node(self) 

603 for sink in self.sinks: 

604 if sink != other_node_on_resistor: 

605 result -= sink.volume_flow 

606 source: ConnectedFlowObject 

607 for source in self.sources: 

608 if source != other_node_on_resistor: 

609 result += source.volume_flow 

610 

611 # if other_node_on_resistor is in sinks, the volume flow is positive (if it flows from source to sink) and 

612 # vice versa: if in sources, the volume flow is negative 

613 if other_node_on_resistor in self.sources: 

614 result *= -1 

615 return result 

616 

617 @property 

618 def start_nodes(self) -> list: 

619 """ 

620 Returns all nodes that are defined as a start of the mass flow. 

621 

622 Used for the algorithm to set the direction of the mass flow (make_velocity_direction). 

623 

624 Returns 

625 ------- 

626 list : 

627 A list with all mass flow start BoundaryConditions. 

628 """ 

629 result = [] 

630 element: ObjectWithPorts 

631 for element in self.rc_objects.all: 

632 if isinstance(element, BoundaryCondition): 

633 if element.is_mass_flow_start: 

634 result.append(element) 

635 return result 

636 

637 @property 

638 def mass_flow_neighbours(self): 

639 """ 

640 Returns all neighbour `MassFlowNode` s. 

641 

642 Returns 

643 ------- 

644 list : 

645 The list with all neighbour `MassFlowNode` s 

646 """ 

647 return [*self.sinks, *self.sources] 

648 

649 @property 

650 def balance(self): 

651 return calculate_balance_for_resistors(self, self.connected_mass_transport_resistors()) 

652 

653 def check_balance(self) -> bool: 

654 """ 

655 Returns True if every parent `ConnectedFlowObject` has a balanced volume flow. 

656 

657 This implies that the algorithm to set the volume flow always walks from the flow start to its end and not 

658 the other way around. 

659 

660 Returns 

661 ------- 

662 bool : 

663 True if every parent `ConnectedFlowObject` has a balanced volume flow. False otherwise. 

664 """ 

665 if self.volume_flow_is_balanced: 

666 return True 

667 # walk the flow upwards and check every source if they are balanced. 

668 source: ConnectedFlowObject 

669 for source in self.sources: 

670 if source.check_balance(): 

671 continue 

672 else: 

673 return False 

674 # check if self is balanced 

675 balance = self.balance 

676 

677 if np.isclose(balance, 0): 

678 self.volume_flow_is_balanced = True 

679 return True 

680 return False 

681 

682 @property 

683 def flow_direction(self) -> list[np.ndarray]: 

684 result = [] 

685 if self._flow_direction is None: 

686 flow_direction = np.array((0, 0, 0)) 

687 for neighbour in self.mass_flow_neighbours: 

688 flow_direction += np.abs(neighbour.get_direction(self)) 

689 for i, value in enumerate(flow_direction): 

690 if value > 0: 

691 result.append(np.zeros(3)) 

692 result[-1][i] = 1 

693 self._flow_direction = flow_direction 

694 return self._flow_direction 

695 

696 def get_minimal_perpendicular_lengths_to_flow(self) -> tuple: 

697 flow_direction: list = self.flow_direction 

698 

699 if len(flow_direction) > 1: 

700 minima = np.partition(self.delta, 1)[0:2] 

701 return minima[0], minima[1] 

702 

703 mask_vector = np.array([1, 1, 1]) - np.abs(flow_direction[0]) 

704 

705 masked_array = np.ma.array(self.delta, mask=mask_vector) 

706 # invert mask of masked array 

707 masked_array = np.ma.array(masked_array.data, mask=~masked_array.mask) 

708 return tuple(masked_array.compressed()) 

709 

710 def get_perpendicular_length_parallel_to_flow(self, asking_node: Node): 

711 """ 

712 Returns the length that is perpendicular to the connection asking_node - self and is parallel to the flow. 

713 

714 Used for the effective area between a Solid and a MassFlowNode as well as to get the length of ChannelNodes. 

715 

716 Parameters 

717 ---------- 

718 asking_node : Node 

719 

720 Returns 

721 ------- 

722 float : 

723 The length. 

724 """ 

725 first_vector = self.mass_flow_neighbours[0].get_direction(self) 

726 

727 second_vector = self.get_direction(asking_node) 

728 

729 # make cross product to get the third vector 

730 third_vector = np.abs(np.cross(first_vector, second_vector)) 

731 

732 return float((third_vector.reshape(1, 3) @ self.delta.reshape(3, 1)).item()) 

733 

734 def connected_mass_transport_resistors(self) -> list[MassTransport]: 

735 """ 

736 Returns all connected mass transport resistors. 

737 

738 Returns 

739 ------- 

740 list[MassTransport] : 

741 The list. 

742 """ 

743 from pyrc.core.resistors import MassTransport 

744 

745 result = [] 

746 for neighbour in self.neighbours: 

747 if isinstance(neighbour, MassTransport): 

748 result.append(neighbour) 

749 return result 

750 

751 @property 

752 def mass_transport_resistors(self) -> list: 

753 """ 

754 Returns a list with all MassTransport resistors in the network. 

755 

756 Returns 

757 ------- 

758 list : 

759 The list with all MassTransport resistors in the network. 

760 """ 

761 from pyrc.core.resistors import MassTransport 

762 

763 return self.rc_objects.get_all_objects(MassTransport) 

764 

765 @property 

766 def mass_flow_nodes(self) -> list: 

767 """ 

768 Returns a list with all MassTransport resistors in the network. 

769 

770 Returns 

771 ------- 

772 list : 

773 The list with all MassTransport resistors in the network. 

774 """ 

775 from pyrc.core.resistors import MassFlowNode 

776 

777 return self.rc_objects.get_all_objects(MassFlowNode) 

778 

779 # def make_velocity_direction(self): 

780 # """ 

781 # Determines and sets the source and sink for each resistor in a network using iterative DFS. 

782 # """ 

783 # from pyrc.core.resistors import MassTransport 

784 # resistors: list[MassTransport] = self.mass_transport_resistors 

785 # start_nodes = self.start_nodes 

786 # 

787 # graph = {} 

788 # resistor_map = {} # Maps node pairs (unordered) to resistors 

789 # 

790 # # Step 1: Build adjacency list and resistor mapping 

791 # for resistor in resistors: 

792 # u, v = resistor.neighbours # Each resistor connects exactly two nodes 

793 # key = frozenset((u, v)) # Store as an unordered key 

794 # 

795 # graph[u].append(v) 

796 # graph[v].append(u) 

797 # resistor_map[key] = resistor # Single entry instead of two 

798 # 

799 # # Step 2: Start iterative DFS from each start node 

800 # for start_node in start_nodes: 

801 # stack = [(start_node, None)] 

802 # visited = set() 

803 # 

804 # while stack: 

805 # node, parent = stack.pop() 

806 # if node in visited: 

807 # continue 

808 # visited.add(node) 

809 # 

810 # for neighbor in graph[node]: 

811 # if neighbor == parent: 

812 # continue 

813 # key = frozenset((node, neighbor)) 

814 # resistor = resistor_map[key] 

815 # resistor.source = node 

816 # resistor.sink = neighbor 

817 # stack.append((neighbor, node)) 

818 # 

819 # # Ensure all resistors have source and sink set 

820 # for resistor in resistors: 

821 # if resistor.source is None or resistor.sink is None: 

822 # raise ValueError( 

823 # f"Resistor between {resistor.neighbours[0]} and {resistor.neighbours[1]} is not properly set.") 

824 

825 def propagate_flow(self): 

826 """ 

827 Propagates the volume flow through all connected `ConnectedFlowObject` s from start to end. 

828 

829 It considers already set volume flows on its way to the end, but until now there is now backflow allowed. So 

830 if the set volume flow of an sink exceeds the balance it will be overwritten. 

831 

832 Only works if the volume flow is set at the start `FlowBoundaryCondition`! It walks from the start to the end 

833 and not the other way round. 

834 """ 

835 start_nodes = self.start_nodes 

836 queue = deque(start_nodes) 

837 considered = set(start_nodes) 

838 

839 while queue: 

840 node: ConnectedFlowCapacitorObject = queue.popleft() 

841 

842 if node.guess_volume_flow is None: 

843 warnings.warn("That shouldn't have happened.") 

844 continue # Skip nodes without a defined volume flow 

845 node.volume_flow_is_balanced = True 

846 

847 sink_nodes: list[ConnectedFlowObject] = [n for n in node.sinks if n not in considered] 

848 if not sink_nodes: 

849 # dead end. 

850 continue 

851 

852 defined_sinks = [n for n in sink_nodes if n.volume_flow_is_balanced] 

853 undefined_sinks = [n for n in sink_nodes if not n.volume_flow_is_balanced] 

854 

855 if len(defined_sinks) > 0: 

856 total_defined_flow = sum(n.guess_volume_flow for n in defined_sinks) 

857 remaining_flow = node.guess_volume_flow - total_defined_flow 

858 if np.isclose(remaining_flow, 0): 

859 remaining_flow = 0 

860 if remaining_flow < 0: 

861 # overwrite all incorrect sink volume flows with equal distribution 

862 # TODO: Maybe allow negative volume flows (back-flows)? 

863 print("Warning: Set volume flow is overwritten to create no backflow.") 

864 undefined_sinks.extend(defined_sinks) 

865 defined_sinks = [] 

866 for sink in defined_sinks: 

867 if sink not in considered: 

868 # TODO: check if sources are all having a defined volume flow (this would save circle circuits) 

869 add_sink_to_queue = True 

870 for source in sink.sources: 

871 if source.volume_flow_is_balanced: 

872 continue 

873 else: 

874 add_sink_to_queue = False 

875 break 

876 if add_sink_to_queue: 

877 queue.append(sink) 

878 considered.add(sink) 

879 else: 

880 remaining_flow = node.guess_volume_flow 

881 

882 if len(undefined_sinks) > 0: 

883 equal_flow = remaining_flow / len(undefined_sinks) 

884 for sink in undefined_sinks: 

885 if sink not in considered: 

886 if sink.guess_volume_flow is None: 

887 sink._volume_flow = 0 # Initialize if undefined 

888 

889 sink._volume_flow += equal_flow 

890 # set volume flow in resistor in between 

891 resistor = node.get_mass_transport_to_node(sink) 

892 resistor._volume_flow = equal_flow 

893 

894 if all(s.volume_flow_is_balanced for s in sink.sources): 

895 queue.append(sink) 

896 considered.add(sink) 

897 

898 self.set_remaining_resistor_volume_flow() 

899 

900 def set_remaining_resistor_volume_flow(self): 

901 """ 

902 Sets all unset volume flows of all `MassTransport` `Resistors`. 

903 

904 The most of them should have been set already in `self.propagate_flow()`. 

905 

906 """ 

907 resistors = [r for r in self.mass_transport_resistors if r.guess_volume_flow is None] 

908 max_attempts = len(resistors) 

909 queue = deque(resistors) 

910 attempts = {resistor: 0 for resistor in resistors} # keep track how often the resistor is re-added to the queue 

911 

912 resistor: MassTransport 

913 while queue: 

914 resistor = queue.popleft() 

915 

916 node: ConnectedFlowCapacitorObject 

917 for node in resistor.neighbours: 

918 # First: check if node only has two resistors connected. If so, set volume_flow of node as own one. 

919 if len(node.sinks) == 1 and len(node.sources) == 1: 

920 resistor._volume_flow = node.volume_flow 

921 break # go to next resistor 

922 else: # Only executed if the first loop did not break 

923 node: ConnectedFlowCapacitorObject 

924 from pyrc.core.resistors import MassTransport 

925 

926 for node in resistor.neighbours: 

927 # Determine the volume_flow balance of the node using the resistor.guess_volume_flow 

928 # At least 'resistor' has no volume_flow set. If it is the only one, it can be calculated. 

929 # Otherwise, resistor is added back to queue. 

930 

931 defined_resistors = {"sink": [], "source": []} 

932 undefined_resistors = {"sink": [], "source": []} 

933 

934 connected_resistors = [n for n in node.neighbours if isinstance(n, MassTransport)]