Heat diffusion
Heat diffusing on a plate from a hot square in its middle, du/dt = D laplacian(u).
- Name:
heat_diffusion - Field:
physics/heat - Formalism: ordinary differential equations
- References: Fourier 1822; Datseris and Parlitz 2022
Factory
@register_model(
field='physics/heat',
tags=('field', 'diffusion', 'pde'),
references=('Fourier 1822', 'Datseris and Parlitz 2022'),
defaults={'until': 50.0, 'dt': 1.0},
)
def heat_diffusion(
shape: Sequence[int] = (64, 64),
*,
D: float = 1.0,
spacing: float = 1.0,
boundary: str = 'zero-flux',
) -> System:
"""Return heat diffusing on a plate from a hot square in its middle, `du/dt = D laplacian(u)`.
Args:
shape:
The number of sites along every dimension.
D:
The diffusion coefficient.
spacing:
The distance between sites.
boundary:
`'periodic'`, `'zero-flux'` for an insulated plate, or `'fixed'` for edges held at zero.
Returns:
system:
The model, with the temperature `plate.u`.
Examples:
>>> from skcomplex.models.physics.heat import heat_diffusion
>>> heat_diffusion((10, 10)).elements[0].check().prepare().size_
100
"""
shape = tuple(shape)
u = np.zeros(shape)
middle = tuple(slice(size // 2 - max(size // 8, 1), size // 2 + max(size // 8, 1)) for size in shape)
u[middle] = 1.0
condition = ('fixed', 0.0) if boundary == 'fixed' else boundary
space = LatticeSpace(shape, spacing=spacing)
return System(
elements=[
FieldElement('plate', element=HeatElement, space=space, boundary={'u': condition}, values={'u': u, 'D': D}),
],
interactions=[DiffusionInteraction('conduction', target='plate', state='u', coefficient=lambda s: s.target.D)],
)
Components
class HeatElement(Element):
"""Site of a plate, with the temperature `u` and the diffusion coefficient `D` of heat.
Args:
label:
The label of the site.
visual:
How the site is drawn, or `None` for none.
u:
The initial temperature.
D:
The diffusion coefficient.
Examples:
>>> from skcomplex.models.physics.heat import HeatElement
>>> HeatElement('site', D=0.5).D
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.D = D
class DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import DiffusionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
Equations
plate.u[i]: diff(u_plate(t), t) == D_plate*laplacian(u_plate(t))