Heat on an L-shaped domain
Heat on the L-shaped domain with fixed boundary temperatures, discretised by linear finite elements.
- Name:
heat_l_shape - Field:
physics/heat - Formalism: ordinary differential equations
- References: Strang and Fix 1973; Gustafsson and McBain 2020
Factory
@register_model(
field='physics/heat',
tags=('field', 'diffusion', 'finite element', 'mesh', 'pde'),
references=('Strang and Fix 1973', 'Gustafsson and McBain 2020'),
defaults={'until': 2.0, 'dt': 0.1, 'solver': IMEXEulerSolver(dt=0.05)},
)
def heat_l_shape(
refine: int = 3,
*,
D: float = 1.0,
temperature: Callable[[Array, Array], Array] | float = harmonic,
) -> System:
"""Return heat on the L-shaped domain with fixed boundary temperatures, discretised by linear finite elements.
The plate starts cold inside, and its temperature relaxes to the solution of Laplace's equation with the
boundary temperatures. The default boundary temperature is harmonic, so the steady state is known.
Args:
refine:
The number of uniform refinements of the mesh.
D:
The diffusion coefficient.
temperature:
The boundary temperature, a number or a function of the coordinates.
Returns:
system:
The model, with the temperature `plate.u` at the nodes of the mesh.
Examples:
>>> from skcomplex.models.physics.heat import heat_l_shape
>>> heat_l_shape(1).elements[0].check().prepare().size_
21
"""
grid = skfem.MeshTri.init_lshaped().refined(refine)
mesh = MeshSpace(grid.p, grid.t)
return System(
elements=[
FieldElement(
'plate',
element=HeatElement,
space=mesh,
boundary={'u': ('fixed', temperature)},
values={'D': D},
discretization='finite-element',
),
],
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))