2. Recovering a heat source

Can a warm field reveal its source?

Consider a square conductor with a spatially varying internal heat source and all four edges held at one temperature. At steady state, the generated heat must leave through the boundary. Fourier conduction gives a heat flux $\boldsymbol q=-k\nabla T$ and steady balance gives $\nabla\cdot\boldsymbol q=q_v$ . Thus $-\nabla\cdot(k\nabla T)=q_v$ . The heat-transfer book develops these physical ingredients.

Choose a reference temperature $T_{\mathrm{ref}}$ , a positive temperature-difference scale $T_*$ , and a reference conductivity $k_*$ . Write $T=T_{\mathrm{ref}}+T_*u$ and $k=k_*a$ . For a source $q_v=k_*T_*s\sin(wx)\sin(wy)$ , our dimensionless temperature $u$ and positive dimensionless conductivity $a$ satisfy

-a\Delta u=s\sin(wx)\sin(wy),\qquad u=b\text{ on the boundary}.

The square has side $L=1\,\mathrm m$ , $w=\pi/L$ , source scale $s$ has units $\mathrm m^{-2}$ , and the boundary value $b$ is dimensionless. Constant $a$ and a sinusoidal source make a hand solution possible; these are explicit assumptions of this example.

Conceptual geometry. The source vanishes on the edges; the boundary fixes the temperature there regardless of the source strength.

Solve before fitting

Try $u=b+A\sin(wx)\sin(wy)$ . Each second derivative contributes $-w^2$ , so substitution gives $2aw^2A=s$ . Therefore

u(x,y)=b+\frac{s}{2aw^2}\sin(wx)\sin(wy).

With $a=1$ , $b=0$ , and $s=2\pi^2\,\mathrm m^{-2}$ , the centre value is 1. Doubling $s$ doubles every temperature excess above $b$ . Doubling $a$ halves it. These predictions follow directly from the equation.

Suppose $a$ and $b$ are known and point measurements are $d_i$ . Define $g_i=\sin(wx_i)\sin(wy_i)/(2aw^2)$ . The model predicts $b+sg_i$ . Minimizing the sum of squared differences gives

\Phi(s)=\frac12\sum_i(b+sg_i-d_i)^2,\quad \Phi'(s)=\sum_i g_i(b+sg_i-d_i),

and hence

\widehat s=\frac{\sum_i g_i(d_i-b)}{\sum_i g_i^2}.

The denominator must be positive. Sampling only boundary points makes every $g_i$ zero: those data cannot recover the source. This scalar calculation is a worked least-squares problem; the general formulation is in Boyd and Vandenberghe, chapter 12 of Introduction to Applied Linear Algebra.

Write the same balance in Eqiora

The component takes its geometry, material, source, and boundary values as inputs:

public component InversePoisson(
  support square: volume(ambient_dimension = 2),
  support x_lower: boundary(parent = square),
  support x_upper: boundary(parent = square),
  support y_lower: boundary(parent = square),
  support y_upper: boundary(parent = square),
  parameter diffusion: 1,
  parameter wave_number: 1 / m,
  parameter source_scale: 1 / m ^ 2,
  parameter boundary_offset: 1
) {

  variable potential: 1 on square;
  relation balance on square {
    -div(diffusion * grad(potential))
      - source_scale * math.sin(wave_number * coordinate(0))
        * math.sin(wave_number * coordinate(1)) = 0;
  }
  relation x_lower_value on x_lower { trace(potential) - boundary_offset = 0; }
  relation x_upper_value on x_upper { trace(potential) - boundary_offset = 0; }
  relation y_lower_value on y_lower { trace(potential) - boundary_offset = 0; }
  relation y_upper_value on y_upper { trace(potential) - boundary_offset = 0; }
}

Read -div(diffusion * grad(potential)) as the conduction balance, the two sine factors as the source shape, and each trace relation as an edge temperature. The caller supplies the square and its four edges. This separation lets one component serve several meshes and parameter choices.

Run the recovery experiment

Use the environment from Get started. Save inverse-poisson.eqi and learn_inverse_poisson.py beside .venv, then run:

uv run --no-project --python .venv/bin/python python learn_inverse_poisson.py inverse-poisson.eqi

The script chooses Q1 finite elements, computes a baseline response, and generates synthetic data with source $1.4(2\pi^2)$ . It applies the scalar formula above to the computed response vector. Expect a recovered source close to $27.6348923\,\mathrm m^{-2}$ . No iterative optimization is needed for this one-parameter linear problem.

"""Recover a Poisson source amplitude and inspect solved-field sensitivities."""

from __future__ import annotations

import argparse
from pathlib import Path

import eqiora
import numpy as np


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("model", type=Path)
    parser.add_argument("--cells", type=int, default=12)
    args = parser.parse_args()
    if args.cells < 2:
        parser.error("--cells must be at least 2")

    graph = eqiora.geometry.GeometryGraph()
    rectangle = graph.rectangle(x_bounds=(0.0, 1.0), y_bounds=(0.0, 1.0))
    names = ("x_lower", "x_upper", "y_lower", "y_upper")
    geometry = graph.build(
        rectangle,
        named_topology={
            "square": rectangle.region,
            **dict(zip(names, rectangle.boundaries, strict=True)),
        },
    )
    mesh = eqiora.meshing.generate(
        eqiora.meshing.resolve(
            geometry, eqiora.meshing.CartesianMesher(cells=(args.cells, args.cells))
        )
    )
    model = eqiora.compile(
        path=args.model,
        geometry=geometry,
        entry="InversePoisson",
        bindings={
            "square": geometry.selection("square"),
            **{
                name: (geometry.selection(name), geometry.selection("square"))
                for name in names
            },
            "diffusion": 1.0,
            "wave_number": np.pi,
            "source_scale": 2.0 * np.pi**2,
            "boundary_offset": 0.0,
        },
    )
    plan = eqiora.resolve(
        model,
        mesh=mesh,
        spatial=eqiora.fem.Q1(),
        solve=eqiora.solve.Linear(
            algorithm=eqiora.solve.LinearSolver.BiConjugateGradientStabilized,
            preconditioner=eqiora.solve.Preconditioner.Identity,
            reduction=eqiora.solve.Reduction.Reproducible,
            provider=eqiora.solve.SolverProvider.reference(),
            relative_tolerance=1.0e-12,
            absolute_tolerance=1.0e-14,
            maximum_iterations=10_000,
        ),
    )
    program = eqiora.diff.compile(
        plan,
        inputs=[
            model.parameter("source_scale"),
            model.parameter("diffusion"),
            model.parameter("boundary_offset"),
        ],
        output=plan.capability.fields[0],
    )
    nominal = np.array([2.0 * np.pi**2, 1.0, 0.0], dtype=np.float64)
    response = program.primal().output.numpy()
    # Linear source response: g is the field produced per unit source scale.
    g = response / nominal[0]
    # A controlled synthetic recovery: use the same discretization for both solves.
    true_source = 1.4 * nominal[0]
    data = program.evaluate(
        np.array([true_source, 1.0, 0.0], dtype=np.float64)
    ).primal().output.numpy()
    recovered = float(np.dot(g, data) / np.dot(g, g))
    print(f"synthetic source: {true_source:.10f}")
    print(f"recovered source: {recovered:.10f}")
    print(f"relative recovery error: {abs(recovered / true_source - 1.0):.3e}")

    # Analytic continuum area mean: integral sin(pi*x) sin(pi*y) = 4/pi**2.
    # Trapezoidal weights integrate the uniform Q1 field, including its boundary.
    weights = np.ones((args.cells + 1, args.cells + 1), dtype=np.float64)
    weights[[0, -1], :] *= 0.5
    weights[:, [0, -1]] *= 0.5
    weights = weights.ravel() / args.cells**2
    print(f"Q1 area mean: {float(np.dot(weights, response)):.10f}")
    print(f"continuum area mean: {4.0 / np.pi**2:.10f}")

    # Differentiate the area mean with respect to [source, diffusion, boundary].
    gradient = program.vjp(weights).input_cotangent.numpy()
    print("area-mean gradient:", gradient)
    print("continuum gradient:", np.array([2.0 / np.pi**4, -4.0 / np.pi**2, 1.0]))
    direction = np.array([0.7, -0.2, 0.3], dtype=np.float64)
    tangent = program.jvp(direction).tangent.numpy()
    pairing_difference = abs(float(np.dot(weights, tangent) - np.dot(direction, gradient)))
    print(f"forward/reverse pairing difference: {pairing_difference:.3e}")
    step = 1.0e-4
    for _ in range(4):
        plus = program.evaluate(nominal + step * direction).primal().output.numpy()
        minus = program.evaluate(nominal - step * direction).primal().output.numpy()
        finite_difference = (plus - minus) / (2.0 * step)
        difference = np.max(np.abs(finite_difference - tangent))
        print(f"step={step:.0e}, maximum derivative difference={difference:.3e}")
        step *= 0.1


if __name__ == "__main__":
    main()

Read the script · Read the component.

The same-mesh synthetic experiment isolates parameter recovery: both the data and the fitted response use the same discrete model. A second printed comparison uses the independently integrated continuum area mean, $4/\pi^2$ . Repeat with --cells 24 to see how the computed mean changes with resolution. Recovering a synthetic parameter precisely and approximating the continuous field accurately are two different questions.

Reuse what you have derived

The Python calls to eqiora.diff reuse the compiled component to evaluate new parameter points and derivatives. There is no need to write a second Poisson solver or recopy the balance for each candidate source. The next chapters explain the additional derivative lines already in the script.

Exercises

A centre observation is 0.6 with $a=2$ and $b=0.1$ . Recover $s$ .
Derive the area mean $\bar u=b+2sL^2/(a\pi^4)$ for this one-metre square.
Replace $b=0$ by a nonzero known boundary value in the analytic fit. Why must you subtract it before fitting the source?
Repeat the script with 6, 12, and 24 cells per direction. Record recovery error and continuum-mean error separately. Explain why their sizes differ.

Check your reasoning: exercise 1 gives $s=2\pi^2\,\mathrm m^{-2}$ .

Reading

Stephen Boyd and Lieven Vandenberghe, Introduction to Applied Linear Algebra, Cambridge University Press, 2018, chapter 12. Authors’ book. For the physical balance, continue with heat transfer; for the spatial approximation, see numerical simulation.

Previous: Measurements to parameters · Book map · Next: Sensitivity and identifiability