Struct Solid3D 

pub struct Solid3D { /* private fields */ }

the arithmetic one, and the difference is not a refinement away.

Implementations

impl Solid3D

pub fn new( name: impl Into<String>, substance: Substance, counts: (usize, usize, usize), dx: Length, initial: Temperature, ) -> Solid3D

A block of counts cubic cells of side dx, all starting at initial.

Each count is forced to at least one. A block one cell thick in two directions is a legitimate thing to ask for and reduces exactly to a bar, which is how the closed-form tests check the three axes against each other.

pub fn melted_fraction_at(&self, i: usize, j: usize, k: usize) -> f64

How much of one cell has melted: 0 entirely solid, 1 entirely liquid.

Always zero for a substance with no FusionProps, which is every entry in the catalogue but Substance::ice.

pub fn set_melted_fraction( &mut self, i: usize, j: usize, k: usize, fraction: f64, )

Declare how much of a cell has melted, for an initial condition a temperature cannot express.

The counterpart to set_temperature, and it exists because a temperature is not a state for a substance that melts: 0 °C is ice, water, or any mixture, and Stefan’s problem starts from liquid at exactly the melting point. Without this the only way to say that is to start a hair above it, which puts sensible heat in the initial condition that the closed form does not have.

Ignored for a substance that does not melt, and out of range is ignored, matching set_temperature.

Supercooling is not representable, and this keeps it that way

The state is one monotone number, so a fraction and a temperature cannot disagree: asking for liquid raises the temperature to at least the melting point and asking for solid lowers it to at most. Liquid below freezing and solid above it are real states of real matter and this model does not have them — the sharp-interface problem assumes the interface is at the melting point, which is what makes Neumann’s solution its solution.

pub fn melted_volume(&self) -> Volume

The volume that has melted, summed over every cell.

The integral quantity, and the one worth reading rather than a front position: a front is only a position if the problem is one-dimensional, and this is the same number in three. For a column of cells it gives the position anyway — melted_volume / area is how far in the front has reached, including the partial cell it is currently inside, which is what makes a front measurable to better than one cell.

pub fn fill( &mut self, substance: Substance, which: impl Fn(usize, usize, usize) -> bool, ) -> &mut Solid3D

Put a different substance in every cell the predicate accepts.

Composable: call it once per layer, per coating, per inclusion. Cells nobody claims keep the constructor’s substance, so a block is never partly undefined.

The harmonic mean, and why it is not a detail

Heat crosses a face, and a face has two materials touching it. The conductance of the half cell either side is in series, so the conductivity that governs the face is

  k_face = 2 k_L k_R / (k_L + k_R)

the harmonic mean, and the arithmetic mean (k_L + k_R)/2 is not an alternative convention — it is wrong. With aluminium against borosilicate the two differ by 38× (2.21 against 84.1 W/m/K), and the arithmetic one short-circuits the interface: a wall it models has 4.2% less resistance than its own layers add up to at 24 cells, and reaching the harmonic answer to 0.1% would take about a thousand.

The harmonic mean is not merely better. With the material interface on a cell face it makes the discrete series resistance exactly Σ Lᵢ/(kᵢA) at every resolution, which is why a_layered_wall.rs is an equality and not a tolerance.

“On a cell face” is the whole condition, and it is one this method cannot break: cells are whole, so an interface is always on a face. A scheme that placed a layer boundary partway through a cell would be first order there whichever mean it used, because the cell would be a mixture and no single conductivity describes one.

What it costs

Five f64 and a u32 per cell beyond the temperatures: three face conductivities, a capacity, its reciprocal scaled by dx, and which substance the cell is. All rebuilt here and none of it in the sweep, which is the trade — the alternative is six harmonic means and six divisions per cell per step.

A uniform block pays the same. That is deliberate: a fast path for one material would be a second implementation of the same physics, exercised by the tests that came before this method and by nothing after it.

pub fn substance_at(&self, i: usize, j: usize, k: usize) -> &Substance

The substance in one cell. Out of range reads the nearest one in range.

pub fn substances(&self) -> usize

How many distinct substances are in the block. One until something has filled it.

pub fn face_conductance( &self, a: (usize, usize, usize), b: (usize, usize, usize), ) -> Option<Conductance>

The conductance of the face between two cells, or None if they are not face neighbours.

k_face · A / dx, which for cubic cells is k_face · dx. This is the number the sweep actually uses, so a caller checking a joint against Σ Lᵢ/(kᵢA) by hand is checking the same arithmetic the march does rather than a restatement of it.

pub fn gap_patches(&self) -> Vec<GapPatch>

The clearances in this block, as sheets of facing area rather than as cell pairs.

A pair is what the exchange is computed on; a patch is what a view factor is a statement about, and the two are not the same object. Grouped by axis and separation, then by connectivity across the gap, so two unrelated clearances at the same width do not average into one patch that describes neither.

Empty for a block with no void, which is every block that existed before there was one.

pub fn empty(self, which: impl Fn(usize, usize, usize) -> bool) -> Solid3D

Mark cells as nothing — not a substance, not part of the block.

A grid had no void until now, so the cells a part did not occupy were some other material, and an assembly of two parts in air was two parts buried in whatever the block was made of. ARCHITECTURE.md names it: “insulating it is a substance with a low conductivity, which is not the same thing”. A low conductivity still conducts, still stores heat, and still sets a stability limit; nothing does none of those.

A void cell holds no heat, conducts to nothing — every face it touches carries zero, which the harmonic mean already gives for a zero conductivity — takes no share of what arrives on the bus, and is left out of every average. Its temperature is not a number, because there is nothing there to have one, and temperature_at says so rather than returning a zero somebody would plot.

What crosses a gap, and what does not

Radiation does. Two solid cells facing each other along a grid line across a run of void exchange σA(T₁⁴ − T₂⁴)/(1/ε₁ + 1/ε₂ − 1), the parallel-plate series, with an exchange factor of one.

That factor was written down here as a known approximation — a wide gap has a view factor well under one, so charging it as one was said to couple a wide gap too hard. Measuring it says otherwise, and the correction is worth more than the caveat was. The sides of a gap in this model are the block’s own outer faces, and an insulated boundary is implemented as a mirror; a mirror puts an image of each surface beyond it and the images tile the plane, so the pair is two infinite parallel plates and one is exact at every width. What the old note described was a different geometry from the one the model has.

The geometry it described is a real one, though — two parts floating in vacuum, open to space, where most of what leaves one surface does miss the other. Which of the two a scene means is a statement about its boundary that a grid cannot infer, so it is not chosen here. It is reported: gap_patches groups a clearance into the sheets a view factor is a statement about and GapPatch::view_factor gives the open-gap number, so the difference between the two readings is a factor somebody can see. For a 32 mm square part 16 mm under a lid it is 2.4x.

Still not here: a radiative-exchange solver. Side walls that are real material at a real temperature do not radiate into the gap at all, and that is a domain rather than a boundary condition.

Convection does not. A gap full of air carries heat by moving that air, which needs a Rayleigh number and a correlation, and a correlation is not a closed form. So a gap in air is coupled less here than it really is, by however much the convection would have carried — and for a millimetre-scale gap at modest temperatures that is the same order as the radiation, so the answer is a lower bound rather than an estimate.

Conduction does not, and that is right: there is nothing there to conduct through.

It is still not a fluid. A part in air also loses heat to the room, and that is Solid3D::losing_from, which is about the block’s outer faces.

pub fn is_void(&self, i: usize, j: usize, k: usize) -> bool

Whether this cell is void — nothing rather than a substance.

pub fn void_cells(&self) -> usize

How many cells hold nothing.

pub fn losing_from(self, face: Face, environment: Environment) -> Solid3D

Expose a face to an environment, so the block can lose heat through it.

Until this existed a Solid3D was adiabatic on all six faces, which means no three-dimensional thermal scene could reach a steady state: every one of them warmed for as long as it ran. That is honest for a pulse and useless for the question a designer actually asks — what temperature does this run at — which is the question a chip package, a magnet busbar, a motor and a factory cell are all asking.

The loss is Environment::loss_from’s, convective and radiative, applied per boundary cell at that cell’s own temperature and its own emissivity. Per cell rather than to a mean, because that is what makes a gradient: the middle of a face runs hotter than its edge, and a lumped loss says it does not.

The environment’s area is the whole face’s area; each boundary cell is charged its share, area / cells_on_the_face. Stating the face’s area rather than a cell’s keeps the number a caller writes independent of the grid they chose, which is what lets the same scene be refined without becoming a different problem — the property pantometry-world’s verify sweep depends on.

Exposing a face tightens the stability limit, because a cell that can also lose to air has more conductance leaving it. Solid3D::max_stable_dt accounts for it; a caller stepping by hand past the returned limit is refused as ever.

pub fn dissipating( self, watts: f64, where_: impl Fn(usize, usize, usize) -> bool, ) -> Solid3D

Generate watts spread evenly over the cells where_ selects, replacing whatever those cells generated before.

The total is the watts given, not watts per cell: a die dissipating 50 W dissipates 50 W whether the grid gives it eight cells or eight thousand, so a scene’s answer does not move when its grid refines. That is what makes a source stated this way survive verify’s resolution sweep, and it is the opposite of the choice a per-cell figure would force.

Void cells are skipped and do not count toward the spread — nothing generates nothing — so a box drawn around a part and its clearance heats the part at the full rate rather than losing a share of it to the gap.

Selecting no solid cell is not an error here, because a caller building an assembly cell by cell passes through that state; the watts simply have nowhere to go and none are generated. A scene refuses it, because a scene saying 50 W and meaning none is a different mistake.

pub fn cells_on_where( &self, where_: &dyn Fn(usize, usize, usize) -> bool, ) -> usize

How many solid cells a predicate selects.

What a caller needs in order to refuse a source that would be generated nowhere: the block itself allows that state, because assembling cell by cell passes through it, and a scene does not.

pub fn generated_power(&self) -> Power

Total watts generated, summed over the cells.

What a caller gets back is what it asked for, which is worth being able to check: a source spread over cells and then read back through a different route is exactly where a factor of the cell count hides.

pub fn generated_energy(&self) -> Energy

Joules generated over the run.

pub fn lost_energy(&self) -> Energy

Heat given up to the exposed faces’ environments over the run.

pub fn exposed_faces(&self) -> impl Iterator<Item = (Face, &Environment)> + '_

Which faces are exposed, and to what.

pub fn counts(&self) -> (usize, usize, usize)

How many cells along each axis.

pub fn spacing(&self) -> Length

The cell side.

pub fn size(&self) -> LengthVec

The block’s extent, which is counts × dx — the outer faces, not the cell centres.

pub fn index(&self, i: usize, j: usize, k: usize) -> Option<usize>

The flat index of a cell, or None if any component is out of range.

Returned rather than panicking because the natural way to write a stencil is to ask for a neighbour that may not exist, and a boundary is exactly where that happens.

pub fn centre_of(&self, i: usize, j: usize, k: usize) -> LengthVec

Where the centre of cell (i, j, k) is, in the block’s own coordinates.

pub fn temperature_at(&self, i: usize, j: usize, k: usize) -> Temperature

The temperature of one cell. Out of range reads the nearest one in range.

pub fn set_temperature(&mut self, i: usize, j: usize, k: usize, t: Temperature)

Set one cell, for an initial condition a constructor cannot express.

This does not change what the block has absorbed. It is a statement about the initial state, not a delivery of heat, so stored_heat moves and absorbed_energy does not — and a simulation started this way and then audited will show the difference as its opening balance rather than as a leak. Use deposit for heat that arrived.

Out of range is ignored rather than a panic: a caller writing a hot spot in a loop over a radius is the expected use, and clipping at the edge is what they mean.

pub fn deposit(&mut self, i: usize, j: usize, k: usize, joules: Energy)

Put joules into one cell, as heat that arrived there.

Counts toward absorbed_energy, so the books balance. Out of range is ignored, which would silently lose energy — so it does not: the joules are refused, and nothing is added to either total.

pub fn mean_temperature(&self) -> Temperature

Mean over every cell. Every cell has the same volume, so this is the volume average.

pub fn peak_temperature(&self) -> Temperature

The hottest cell — the number a hot spot exists to produce, and the one a lumped model reports as the mean.

pub fn coldest_temperature(&self) -> Temperature

The coldest cell.

pub fn absorbed_energy(&self) -> Energy

Heat taken from the bus over the run.

pub fn fourier_number(&self, dt: Time) -> f64

α·dt/dx² for the constructor’s substance — the classic Fourier number.

This is the quantity mode_amplification is written in terms of, so it stays the textbook one and does not become shape- or fill-aware. For a block of one material with at least two cells on every axis it is the stability number, and fourier_number(max_stable_dt) == 1/6 exactly.

It is not the stability number for a thin block or a filled one, and those are the two cases where a caller sizing a step by hand needs the other one — see stability_ratio, which is what step actually enforces.

pub fn stability_ratio(&self, dt: Time) -> f64

dt as a fraction of the largest step this block is stable at — exactly one at the limit.

The number step refuses on, and the honest one for a block that is thin or filled. It is dt · maxᵢ (Σ_f k_f / Cᵢ) · dx: a maximum over cells, because stability is a statement about a row of the update matrix and each cell has its own.

pub fn heat_capacity(&self) -> HeatCapacity

The heat capacity of the whole block, which for a filled one is a sum and not a product.

pub fn mode_amplification(&self, mode: (usize, usize, usize), dt: Time) -> f64

The exact per-step amplification of one separable cosine mode.

(a, b, c) are half-wave counts along the three axes: (1, 0, 0) is the longest mode along x with the other two flat. Returns the factor the mode’s amplitude is multiplied by in one step of dt, which is 1 + F·Σ(−4 sin²(mπ/2n)) — exact, not a linearisation, because that mode is an eigenvector of the discrete operator this domain steps with.

Public because it is what makes this domain checkable without a reference implementation: a caller can predict an amplitude arbitrarily far ahead and compare. It is also what a grid designer wants, since a factor outside (−1, 1] is the instability itself.

pub fn release_mode( &mut self, mode: (usize, usize, usize), mean: Temperature, amplitude: f64, )

Fill the block with one separable cosine mode about a mean.

The initial condition the closed-form tests use, and a genuinely useful one for anybody checking a grid: it is the only shape whose future this domain can state exactly.

pub fn mode_amplitude(&self, mode: (usize, usize, usize)) -> f64

The amplitude of one mode currently present, by projection.

The counterpart to release_mode, and what makes a decay measurable rather than merely visible. Cosine modes on this grid are orthogonal, so this is exact for a block holding one of them and is the correct coefficient for a block holding several.

pub fn volume(&self) -> Volume

The volume of the whole block.

Trait Implementations

impl Clone for Solid3D

fn clone(&self) -> Solid3D

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for Solid3D

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Domain for Solid3D

fn max_stable_dt(&self, _now: Time) -> Time

minᵢ Cᵢ / (dx · Σ_f k_f) — the tightest row of the update matrix, and dx²/(6α) when the block is one material with at least two cells on every axis.

Why a maximum over cells, and what summing the faces buys

Stability is Gershgorin’s condition on one row: the update is I + dt·D⁻¹L, row i has diagonal −θᵢ and off-diagonals summing to +θᵢ for θᵢ = dt Σ_f G_f / Cᵢ, so every eigenvalue is in [1 − 2θ_max, 1] and θ_max ≤ 1 is the limit.

On a uniform grid that disc is tight rather than cautious — the sharpest representable mode saturates it — and this is where 1/2, 1/4 and 1/6 come from. They are not three dimensions; they are however many axes have more than one cell, and the block is stepped at whichever applies. This used to report dx²/(6α) for every shape, which cost a bar-shaped block three times the steps for nothing.

On a filled grid the value is the opposite one: the limit is usually far looser than dx²/(6·α_max), because k_f ≤ 2·min(k_L, k_R) means a cell cannot be heated through a face faster than its worse side allows. Heat does not reach a fast material at that material’s own rate; it reaches it at the rate the neighbour delivers. Measured on one aluminium cell embedded in borosilicate, the honest limit is 75× the one aluminium’s diffusivity would name, and that factor is wall-clock.

It can in principle go the other way — up to a factor of two, when a cell’s neighbours conduct better and store more — but that needs volumetric heat capacity to vary as widely as conductivity, and across solids it varies by one order of magnitude where conductivity varies by four. So the tightening is a bound this catalogue cannot reach, and the loosening is a saving any coating or inclusion gets.

A third of the bar’s, for a cube

A third of what Bar1D reports for the same spacing and material. That factor is the reason Schedule::Multirate exists: a block and a lumped mass in one scene differ by five orders of magnitude in the step they can take, and a single global step would make the cheap domain pay the expensive one’s bill.

Stable is not accurate, and this is only the first

At exactly this step the sharpest mode the grid can hold has an amplification factor of −1. It flips sign every step and never decays. That is what marginal stability means and it is not a defect — the scheme does not diverge there, which is the whole of what a stability limit claims.

It does mean sharp initial data is carried badly. A point source excites that mode as hard as anything can, and the peak comes out 1.96× the exact answer while the conservation audit stays exact to the last bit. At half this step it is 1.005×.

So take a fraction of it when the initial condition is sharp. Schedule::Multirate divides by ceil(dt / limit) and so usually lands comfortably inside, but a caller stepping by hand can sit exactly on it. cargo run --example heat_in_three_dimensions is the measurement.

fn ledger(&self) -> Ledger

Heat gained since the start. The faces are insulated, so this is exactly what came in.

fn checkpoint(&mut self)

The temperatures and the phase, because a temperature alone is not a state.

A cell at 0 °C is ice, water or any mixture, so saving cells and not melted would restore a block that was half melted as one that was entirely solid at the same temperature — losing 306 mJ per cubic millimetre with nothing to say it had gone. Schedule::Iterative and the audit’s retry both restore, so this is a live path and not a precaution.

fn readings(&self) -> Vec<Reading>

Peak, mean and coldest, in celsius, and what it has absorbed.

All three ends of the distribution, because the whole reason to pay for a 3D grid is that they differ. A block reported by its mean alone is a LumpedMass that cost n³ times as much, and the gap between peak and mean is the number that says whether the reduction would have been honest.

fn as_field(&self) -> Option<&dyn ScalarField>

A temperature field, so nothing above has to know this is a block.

fn books_balance(&self) -> bool

Whether this domain’s books are exact: its ledger changes by precisely what it takes from the bus minus what it publishes, every step.

fn name(&self) -> &str

What this domain is called. Used to look it up and to name it in a violation.

fn step( &mut self, _t: Time, dt: Time, bus: &mut Exchange, ) -> Result<(), Violation>

Advance by dt from t, reading inputs from bus and publishing outputs to it. A quasi-static domain ignores dt.

fn restore(&mut self)

Restore the last Domain::checkpoint.

fn supports_restore(&self) -> bool

Whether Domain::checkpoint and Domain::restore actually do something.

fn as_any(&self) -> Option<&dyn Any>

This domain as Any, so a caller can get the concrete type back out of a Simulation — see Simulation::domain_as.

fn as_any_mut(&mut self) -> Option<&mut dyn Any>

The same, mutably, so a caller can write to a domain between steps.

fn kind(&self) -> Kind

Whether it has state to roll forward. Defaults to Kind::Evolving.

fn residual(&self) -> f64

How far this domain still is from agreeing with its neighbours, for Schedule::Iterative. Zero means converged.

fn as_bodies(&self) -> Option<&dyn Bodies>

This domain as a countable set of bodies, if that is what it is.

impl ScalarField for Solid3D

fn unit(&self) -> &'static str

Kelvin, because that is what the cells hold. See Bar1D.

fn at(&self, p: LengthVec, _t: Time) -> f64

Trilinear between cell centres, clamped at the faces, and masked over void.

Clamped rather than extrapolated: outside an insulated face the temperature is not defined, and continuing the gradient would draw a block hotter than any cell in it.

Masked because self.cells still holds whatever an emptied cell held when it was emptied, and reading it raw is how a gap came out of every exporter as a piece of the block sitting at ambient forever — a plausible, unchanging number, which is worse than no number. The weights are renormalised over the solid corners, so a sample at a cell centre is that cell exactly, a sample inside a clearance is NaN, and a sample straddling the two is the material’s own value rather than a blend with something that is not there.

fn gradient(&self, p: LengthVec, _t: Time, _h: Length) -> DVec3

Central differences on the cell grid, mirrored at the faces and at void.

A clearance is a boundary, and the block already knows what to do at one: mirror. Walking off the outer edge and walking into nothing are the same situation — there is no sample that way — so both return the asking cell and the difference becomes one-sided, which is the insulated condition rather than a slope towards a number that is not there.

NaN inside a clearance, for the same reason temperature_at gives one: a gradient there is a value somebody would plot.

fn laplacian(&self, p: LengthVec, _t: Time, _h: Length) -> f64

∇²T on the seven-point stencil.

The Laplacian of the temperature, which is what the trait asks for and is a statement about the field rather than about the material. It is the operator the sweep uses only when the block is one material; for a filled one the sweep uses ∇·(k∇T), and rate is where that appears.

fn rate(&self, p: LengthVec, _t: Time, _dt: Time) -> f64

∂T/∂t = (1/Cᵢ)·Σ_f G_f (T_f − Tᵢ), which is α∇²T where the block is one material.

The conductance form rather than α·laplacian, so that it is the sweep’s own operator at a point in a filled block too — read off the same face conductivities, not reconstructed from a diffusivity the block may not have a single value of.

Conduction only: heat arriving over the bus is a source this cannot see.