Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -17,3 +17,4 @@

# coverage artifact
.coverage
.idea
60 changes: 49 additions & 11 deletions heatrapy/dimension_1/solvers/_latent_heat.py
Original file line number Diff line number Diff line change
@@ -1,49 +1,87 @@
"""Shared latent heat computation for 1D solvers."""

from __future__ import annotations

import copy
from typing import TYPE_CHECKING

if TYPE_CHECKING:
from ..objects.object import Object


def apply_latent_heat(nx, obj):
def apply_latent_heat(
nx: list[float], obj: Object
) -> list[list[list[float]]]:
"""Apply latent heat corrections to computed temperatures.

Handles phase transitions by absorbing/releasing energy when the
temperature crosses a transition threshold.
Handles phase transitions by absorbing/releasing energy when
the temperature crosses a transition threshold.

Parameters
----------
nx : list
nx : list[float]
New temperatures for each grid point (modified in place).
obj : Object
Thermal object with current state.

Returns
-------
lheat : list
lheat : list[list[list[float]]]
Updated latent heat accumulation state.

Raises
------
TypeError
If ``nx`` is not a list or ``obj`` lacks required
thermal attributes.

"""
if not isinstance(nx, list):
raise TypeError(
f"nx must be a list, got {type(nx).__name__}"
)
if not hasattr(obj, 'num_points'):
raise TypeError(
f"obj must be a thermal Object, "
f"got {type(obj).__name__}"
)

lheat = copy.copy(obj.lheat)
for i in range(1, obj.num_points - 1):
j = 0
for lh in obj.latent_heat[i]:
temper = obj.temperature[i][0]
# heating: crossing transition from below
if nx[i] > lh[0] and temper <= lh[0] and lheat[i][j][1] != lh[1]:
if (
nx[i] > lh[0]
and temper <= lh[0]
and lheat[i][j][1] != lh[1]
):
en = obj.Cp[i] * obj.rho[i] * (nx[i] - temper)
if en + lheat[i][j][1] >= lh[1]:
lheat[i][j][1] = lh[1]
energy_temp = lheat[i][j][1] + en - lh[1]
nx[i] = temper + energy_temp / (obj.Cp[i] * obj.rho[i])
energy_temp = (
lheat[i][j][1] + en - lh[1]
)
nx[i] = temper + energy_temp / (
obj.Cp[i] * obj.rho[i]
)
else:
lheat[i][j][1] += en
nx[i] = temper
# cooling: crossing transition from above
if nx[i] < lh[0] and temper >= lh[0] and lheat[i][j][1] != 0:
if (
nx[i] < lh[0]
and temper >= lh[0]
and lheat[i][j][1] != 0
):
en = obj.Cp[i] * obj.rho[i] * (nx[i] - temper)
if en + lheat[i][j][1] <= 0.:
lheat[i][j][1] = 0.
energy_temp = (en + lheat[i][j][1])
nx[i] = temper + energy_temp / (obj.Cp[i] * obj.rho[i])
energy_temp = en + lheat[i][j][1]
nx[i] = temper + energy_temp / (
obj.Cp[i] * obj.rho[i]
)
else:
lheat[i][j][1] += en
nx[i] = temper
Expand Down
103 changes: 78 additions & 25 deletions heatrapy/dimension_1/solvers/explicit_general.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,41 +4,90 @@

"""

from __future__ import annotations

from typing import TYPE_CHECKING

import numpy as np

from ._latent_heat import apply_latent_heat

if TYPE_CHECKING:
from ..objects.object import Object

def explicit_general(obj):

def explicit_general(
obj: Object,
) -> tuple[list[list[float]], list[list[list[float]]]]:
"""explicit_general solver.

Used to compute one time step of 1D systems with fixed thermal
conductivity.

Parameters
----------
obj : Object
Thermal object with current state.

Returns
-------
y : list[list[float]]
Updated temperatures as ``[[T, T], ...]`` pairs.
lheat : list[list[list[float]]]
Updated latent heat state.

Raises
------
TypeError
If ``obj`` is not a thermal Object.

"""
n = obj.num_points
s = slice(1, n - 1)

# extract current temperatures and material properties as arrays
T = np.array([obj.temperature[i][0] for i in range(n)], dtype=float)
k = np.array([obj.k[i] if obj.k[i] is not None else 0.0 for i in range(n)])
rho = np.array([obj.rho[i] if obj.rho[i] is not None else 1.0
for i in range(n)])
Cp = np.array([obj.Cp[i] if obj.Cp[i] is not None else 1.0
for i in range(n)])
Q = np.array([obj.Q[i] if obj.Q[i] is not None else 0.0
for i in range(n)])
Q0 = np.array([obj.Q0[i] if obj.Q0[i] is not None else 0.0
for i in range(n)])
if not hasattr(obj, 'num_points'):
raise TypeError(
f"obj must be a thermal Object, "
f"got {type(obj).__name__}"
)

n: int = obj.num_points
s: slice = slice(1, n - 1)

# extract current temperatures and material properties
T: np.ndarray = np.array(
[obj.temperature[i][0] for i in range(n)], dtype=float
)
k: np.ndarray = np.array(
[obj.k[i] if obj.k[i] is not None else 0.0
for i in range(n)]
)
rho: np.ndarray = np.array(
[obj.rho[i] if obj.rho[i] is not None else 1.0
for i in range(n)]
)
Cp: np.ndarray = np.array(
[obj.Cp[i] if obj.Cp[i] is not None else 1.0
for i in range(n)]
)
Q: np.ndarray = np.array(
[obj.Q[i] if obj.Q[i] is not None else 0.0
for i in range(n)]
)
Q0: np.ndarray = np.array(
[obj.Q0[i] if obj.Q0[i] is not None else 0.0
for i in range(n)]
)

# vectorized FDM stencil for interior points
alpha = obj.dt * k[s] / (rho[s] * Cp[s] * obj.dx * obj.dx)
beta = obj.dt / (rho[s] * Cp[s])

nx = T.copy()
nx[s] = ((1 + beta * Q[s]) * T[s] +
alpha * (T[0:n-2] - 2 * T[s] + T[2:n]) +
beta * (Q0[s] - Q[s] * obj.amb_temperature))
alpha: np.ndarray = (
obj.dt * k[s] / (rho[s] * Cp[s] * obj.dx * obj.dx)
)
beta: np.ndarray = obj.dt / (rho[s] * Cp[s])

nx: np.ndarray = T.copy()
nx[s] = (
(1 + beta * Q[s]) * T[s]
+ alpha * (T[0:n-2] - 2 * T[s] + T[2:n])
+ beta * (Q0[s] - Q[s] * obj.amb_temperature)
)

# boundaries
if obj.boundaries[0] == 0:
Expand All @@ -52,10 +101,14 @@ def explicit_general(obj):
nx[n - 1] = obj.boundaries[1]

# latent heat (per-element, branching logic)
nx_list = nx.tolist()
lheat = apply_latent_heat(nx_list, obj)
nx_list: list[float] = nx.tolist()
lheat: list[list[list[float]]] = apply_latent_heat(
nx_list, obj
)

# pack into [current, next] pairs expected by the caller
y = [[nx_list[i], nx_list[i]] for i in range(n)]
y: list[list[float]] = [
[nx_list[i], nx_list[i]] for i in range(n)
]

return y, lheat
117 changes: 86 additions & 31 deletions heatrapy/dimension_1/solvers/explicit_k.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,43 +4,94 @@

"""

from __future__ import annotations

from typing import TYPE_CHECKING

import numpy as np

from ._latent_heat import apply_latent_heat

if TYPE_CHECKING:
from ..objects.object import Object


def explicit_k(obj):
def explicit_k(
obj: Object,
) -> tuple[list[list[float]], list[list[list[float]]]]:
"""explicit_k solver.

Used to compute one time step of 1D systems with x-dependent thermal
conductivity.
Used to compute one time step of 1D systems with x-dependent
thermal conductivity.

Parameters
----------
obj : Object
Thermal object with current state.

Returns
-------
y : list[list[float]]
Updated temperatures as ``[[T, T], ...]`` pairs.
lheat : list[list[list[float]]]
Updated latent heat state.

Raises
------
TypeError
If ``obj`` is not a thermal Object.

"""
n = obj.num_points
s = slice(1, n - 1)

# extract current temperatures and material properties as arrays
T = np.array([obj.temperature[i][0] for i in range(n)], dtype=float)
k = np.array([obj.k[i] if obj.k[i] is not None else 0.0 for i in range(n)])
rho = np.array([obj.rho[i] if obj.rho[i] is not None else 1.0
for i in range(n)])
Cp = np.array([obj.Cp[i] if obj.Cp[i] is not None else 1.0
for i in range(n)])
Q = np.array([obj.Q[i] if obj.Q[i] is not None else 0.0
for i in range(n)])
Q0 = np.array([obj.Q0[i] if obj.Q0[i] is not None else 0.0
for i in range(n)])

# vectorized FDM stencil for interior points (k varies with x)
eta = obj.dt / (2.0 * rho[s] * Cp[s] * obj.dx * obj.dx)
beta = obj.dt / (rho[s] * Cp[s])

nx = T.copy()
nx[s] = ((1 + beta * Q[s]) * T[s] +
eta * ((k[2:n] + k[s]) * T[2:n] -
(k[0:n-2] + k[2:n] + 2 * k[s]) * T[s] +
(k[0:n-2] + k[s]) * T[0:n-2]) +
beta * (Q0[s] - Q[s] * obj.amb_temperature))
if not hasattr(obj, 'num_points'):
raise TypeError(
f"obj must be a thermal Object, "
f"got {type(obj).__name__}"
)

n: int = obj.num_points
s: slice = slice(1, n - 1)

# extract current temperatures and material properties
T: np.ndarray = np.array(
[obj.temperature[i][0] for i in range(n)], dtype=float
)
k: np.ndarray = np.array(
[obj.k[i] if obj.k[i] is not None else 0.0
for i in range(n)]
)
rho: np.ndarray = np.array(
[obj.rho[i] if obj.rho[i] is not None else 1.0
for i in range(n)]
)
Cp: np.ndarray = np.array(
[obj.Cp[i] if obj.Cp[i] is not None else 1.0
for i in range(n)]
)
Q: np.ndarray = np.array(
[obj.Q[i] if obj.Q[i] is not None else 0.0
for i in range(n)]
)
Q0: np.ndarray = np.array(
[obj.Q0[i] if obj.Q0[i] is not None else 0.0
for i in range(n)]
)

# vectorized FDM stencil (k varies with x)
eta: np.ndarray = obj.dt / (
2.0 * rho[s] * Cp[s] * obj.dx * obj.dx
)
beta: np.ndarray = obj.dt / (rho[s] * Cp[s])

nx: np.ndarray = T.copy()
nx[s] = (
(1 + beta * Q[s]) * T[s]
+ eta * (
(k[2:n] + k[s]) * T[2:n]
- (k[0:n-2] + k[2:n] + 2 * k[s]) * T[s]
+ (k[0:n-2] + k[s]) * T[0:n-2]
)
+ beta * (Q0[s] - Q[s] * obj.amb_temperature)
)

# boundaries
if obj.boundaries[0] == 0:
Expand All @@ -54,10 +105,14 @@ def explicit_k(obj):
nx[n - 1] = obj.boundaries[1]

# latent heat (per-element, branching logic)
nx_list = nx.tolist()
lheat = apply_latent_heat(nx_list, obj)
nx_list: list[float] = nx.tolist()
lheat: list[list[list[float]]] = apply_latent_heat(
nx_list, obj
)

# pack into [current, next] pairs expected by the caller
y = [[nx_list[i], nx_list[i]] for i in range(n)]
y: list[list[float]] = [
[nx_list[i], nx_list[i]] for i in range(n)
]

return y, lheat
Loading
Loading