Temperature Computation

Temperature Physics

The temperature system simulates heat diffusion across the world using an implicit finite-difference solver. Each game tick, it assembles and solves a sparse linear system to advance temperature by one timestep — this allows large timesteps without numerical instability, unlike explicit methods.

---

Overview

The system is built on three layers of abstraction:

• Data layers — store per-voxel physical properties compactly in 16×16×16 sections.
• Physics solver (AbstractMatrixPhysicsSolver) — handles all matrix infrastructure: building, updating, and solving the linear system using Conjugate Gradient.
• Temperature solver (TemperatureSolver) — implements the heat equation physics on top of the base solver.

---

The Heat Equation

Each tick solves the following discretised heat equation implicitly:

C * (T_next - T_current) / dt = k∇²T_next + β * resilience * (T_default - T_next)

Rearranged into the linear system A * T_next = T_current + b_source, where:

• C (CAPACITY = 3e4) is the thermal capacity.
• dt (DT = 1/20) is the timestep (one game tick).
• k is the local thermal conductivity (from ConductionDataLayer).
• resilience is how strongly a voxel is pulled toward its default temperature (from ResilienceDataLayer).
• T_default is the voxel's natural equilibrium temperature.

The matrix A is symmetric and diagonally dominant, which makes it well-suited for the Conjugate Gradient solver.

Per-voxel matrix row

For each voxel at global index i, the diagonal entry is:

A[i][i] = 1 + resilience * β + Σ k_eff(i, j)  (sum over all 6 face neighbors)

For each neighbor j in the matrix:

A[i][j] = -k_eff(i, j)

where k_eff is the harmonic mean of the two voxels' conductivities:

k_eff = 2 * k_self * k_neighbor / (k_self + k_neighbor)

Boundary neighbors (outside the active matrix) contribute their temperature directly to the RHS instead of as off-diagonal entries.

The RHS for voxel i is:

rhs[i] = T_current[i] + resilience * β * T_default[i] + Σ k_eff * T_boundary

---

Data Layers

Each section (16×16×16 voxels) has four associated data layers. All layers extend AbstractDataLayer, which provides 3D indexing via:

index = x + z * 16 + y * 256

TemperatureDataLayer

Stores current and default temperatures as 16-bit fixed-point values.

• Encoding: stored = (short)(temperature * 10 + Short.MIN_VALUE)
• Range: 0 K to ~6553 K
• Precision: 0.1 K

ConductionDataLayer

Stores thermal conductivity as a compact 8-bit mini-float.

• Format: 2-bit mantissa + 6-bit signed exponent (bias −16), giving values from 2^-16 to 1.75 × 2^47.
• Decoding: value = (1 + mantissa/4) * 2^exponent
• A cache of decoded float[] values is maintained for fast read access.

ResilienceDataLayer

Stores resilience (how strongly a voxel resists temperature change from diffusion) as an 8-bit unsigned value.

• Range: 0.0 (freely diffuses) to 1.0 (fully locked to default temperature)
• Decoding: value = (stored + 128) / 255

---

Matrix Infrastructure (AbstractMatrixPhysicsSolver)

The base solver manages the full lifecycle of the sparse linear system.

Matrix Storage

The system uses a PaddedCSRMatrix — a flat-array sparse matrix with a fixed number of slots per row (7 for the temperature solver: 1 diagonal + 6 face neighbors). Rows are written in-place; no HashMap or CSR compilation step is needed.

Section Lifecycle

Active sections are tracked as a LongSet. Each tick, the set of ticking sections is compared to the cached matrix:

Situation	Action
No change	Return cached matrix as-is
New sections added	extendMatrix: grow flat arrays, fill new rows, rebuild face-boundary rows of adjacent sections
Sections removed	shrinkMatrix: swap removed sections with the last block, rewrite affected column indices, rebuild boundary rows
Both	Shrink first, then extend

Solver

The linear system is solved with Conjugate Gradient (ConjugateGradient2). Working arrays (cgR, cgP, cgAp) are allocated once and stored in the matrix object to avoid per-tick heap allocation (which would otherwise produce ~240 MB/s of garbage at 20 ticks/s).

The previous tick's solution is used as the warm start for the next solve, which significantly reduces the number of iterations needed in steady state.

---

Dynamic Sources (Block Entities)

Block entities that emit or absorb heat (e.g. furnaces, heaters) register as dynamic sources. Each tick, before the matrix is solved:

1. Their voxel's temperature, defaultTemperature, conductivity, and resilience are overwritten in the relevant data layers.
2. The surrounding voxels are "stamped" to mark those rows as dirty so the matrix is rebuilt for that area.

After the solve, each dynamic source receives the temperature delta T_next - T_current at its voxel via addTemperature(dT), so it can track its own thermal state.

---

Ticking Propagation

After each solve, if any voxel's temperature changed by more than ε = 0.1, the section is marked as needing to tick next step. Additionally, if the changed voxel lies on a section boundary, the adjacent section is also marked for ticking. This ensures thermal fronts propagate naturally across section borders without simulating quiescent regions.