heat1d-conduction-solver

1	#include "heat1d/solver.hpp"
2	
3	#include <cmath>
4	#include <cstddef>
5	#include <stdexcept>
6	#include <vector>
7	
8	#include "heat1d/types.hpp"
9	
10	namespace heat1d {
11	
12	// ============================================================================
13	// TODO(candidate): Implement the implicit ADI (Peaceman--Rachford) 2D
14	// heat-conduction time-stepper.
15	//
16	// The two member functions below are STUBS that currently return the field
17	// UNCHANGED (no heat actually flows), so every transient is wrong and the test
18	// suite FAILS. Replace the stub bodies of
19	//
20	//     std::vector<double> HeatSolver::step(const std::vector<double>& field,
21	//                                          double dt) const;
22	//     HeatResult          HeatSolver::solve(double dt, std::size_t n_steps) const;
23	//
24	// with correct implementations. Do NOT change the public signatures declared
25	// in include/heat1d/solver.hpp; reuse the Plate accessors and the diagnostics
26	// (fourier_x, fourier_y, total_energy) which are already implemented below and
27	// must not be altered.
28	//
29	// ---------------------------------------------------------------------------
30	// Governing equation (2D transient heat conduction):
31	//
32	//     dT/dt = alpha ( d2T/dx2 + d2T/dy2 ),  alpha > 0,
33	//
34	// on the rectangle [0, Lx] x [0, Ly], uniform grid (Nx+1)x(Ny+1), hx = Lx/Nx,
35	// hy = Ly/Ny. Fields are flattened ROW-MAJOR in x: node (i, j) at index
36	// i + (Nx+1)*j.
37	//
38	// ---------------------------------------------------------------------------
39	// ADI (Peaceman--Rachford), one full step of size dt = two half-steps of dt/2.
40	// With rx = alpha*(dt/2)/hx^2, ry = alpha*(dt/2)/hy^2:
41	//
42	//   Half-step 1 (implicit in x): for each row j solve the tridiagonal system
43	//     -rx T*_{i-1,j} + (1+2rx) T*_{i,j} - rx T*_{i+1,j}
44	//          = T^n_{i,j} + ry ( T^n_{i,j-1} - 2 T^n_{i,j} + T^n_{i,j+1} ).
45	//
46	//   Half-step 2 (implicit in y): for each column i solve the tridiagonal system
47	//     -ry T^{n+1}_{i,j-1} + (1+2ry) T^{n+1}_{i,j} - ry T^{n+1}_{i,j+1}
48	//          = T*_{i,j} + rx ( T*_{i-1,j} - 2 T*_{i,j} + T*_{i+1,j} ).
49	//
50	// Solve each tridiagonal system with the Thomas algorithm (implement it
51	// yourself). The scheme is UNCONDITIONALLY STABLE: there is NO Fourier-number
52	// stability cap. Do NOT throw on a large dt.
53	//
54	// ---------------------------------------------------------------------------
55	// Boundary conditions (per edge: left x=0, right x=Lx, bottom y=0, top y=Ly):
56	//   * Dirichlet (T = value): pin the edge nodes to `value` in both sweeps.
57	//   * Neumann (dT/dn = value outward; 0 == insulated): central ghost node,
58	//     e.g. at x=0  T_{-1,j} = T_{1,j} + 2 hx*value, folded into the end row of
59	//     the implicit sweep when the edge is normal to it, else into the RHS.
60	//   * Robin (dT/dn = -B (T - T_inf), B = bc.h, T_inf = bc.T_inf): central ghost
61	//     e.g. at x=0  T_{-1,j} = T_{1,j} - 2 hx*B*(T_{0,j} - T_inf).
62	//   * Corners: if either incident edge is Dirichlet the corner takes that
63	//     value (left/right take precedence over bottom/top); otherwise the corner
64	//     is an unknown carrying ghosts in both directions.
65	//
66	// ---------------------------------------------------------------------------
67	// Contract:
68	//   step(): throw std::invalid_argument on a field-size mismatch or dt <= 0.
69	//   solve(): call Plate::validate() first (propagate its std::runtime_error),
70	//            throw std::invalid_argument on dt <= 0, start from the plate's
71	//            initial() field, apply step() exactly n_steps times, and return
72	//            HeatResult{ temperature, steps = n_steps, time = n_steps*dt }.
73	//
74	// The required behaviour is locked down by tests/: the 2D analytic Fourier mode
75	// decay, unconditional stability at a LARGE dt (an explicit FTCS port diverges
76	// and fails here), mixed-BC steady states (Dirichlet ramp, Robin convective
77	// balance), the energy-conservation invariant for a fully insulated plate, the
78	// discrete maximum principle, a symmetry invariant, and the exception cases.
79	// ============================================================================
80	
81	std::vector<double> HeatSolver::step(const std::vector<double>& field,
82	                                     double dt) const {
83	    // TODO(candidate): perform one full ADI step (two implicit half-sweeps with
84	    // the proper boundary conditions) and return the next-level field.
85	    (void)dt;
86	    if (field.size() != plate_.num_nodes()) {
87	        throw std::invalid_argument("HeatSolver::step: field size mismatch");
88	    }
89	    return field; // unchanged -> wrong (no diffusion)
90	}
91	
92	HeatResult HeatSolver::solve(double dt, std::size_t n_steps) const {
93	    // TODO(candidate): validate, step n_steps times, return the final field.
94	    HeatResult result;
95	    result.temperature = plate_.initial(); // never evolved -> wrong
96	    result.steps = 0;
97	    result.time = 0.0;
98	    (void)dt;
99	    (void)n_steps;
100	    return result;
101	}
102	
103	// ----------------------------------------------------------------------------
104	// Diagnostics (already implemented; do NOT modify).
105	
106	double HeatSolver::fourier_x(double dt) const {
107	    const double hx = plate_.hx();
108	    return plate_.alpha() * (0.5 * dt) / (hx * hx);
109	}
110	
111	double HeatSolver::fourier_y(double dt) const {
112	    const double hy = plate_.hy();
113	    return plate_.alpha() * (0.5 * dt) / (hy * hy);
114	}
115	
116	double HeatSolver::total_energy(const std::vector<double>& field) const {
117	    return plate_.total_energy(field);
118	}
119	
120	} // namespace heat1d
121

1	#ifndef HEAT1D_SOLVER_HPP
2	#define HEAT1D_SOLVER_HPP
3	
4	#include <cstddef>
5	#include <vector>
6	
7	#include "heat1d/rod.hpp"
8	
9	namespace heat1d {
10	
11	/// Result of a transient 2D heat-conduction simulation.
12	struct HeatResult {
13	    /// Temperature field at the final time, length == plate.num_nodes(),
14	    /// flattened ROW-MAJOR in x: node (i, j) at index i + (Nx+1)*j.
15	    std::vector<double> temperature;
16	
17	    /// Number of time steps actually taken.
18	    std::size_t steps{0};
19	
20	    /// Final simulated time [s] = steps * dt.
21	    double time{0.0};
22	};
23	
24	/// Implicit Alternating-Direction-Implicit (ADI, Peaceman--Rachford) finite-
25	/// difference solver for the 2D transient heat-conduction equation
26	///
27	///     dT/dt = alpha * ( d2T/dx2 + d2T/dy2 ),  alpha > 0,
28	///
29	/// on the rectangle [0, Lx] x [0, Ly] discretized on a uniform grid
30	/// (Nx+1) x (Ny+1), spacings hx = Lx/Nx, hy = Ly/Ny. The scheme is 2nd-order
31	/// accurate in both space and time and UNCONDITIONALLY STABLE: there is no
32	/// Fourier-number (stability) cap on the time step.
33	///
34	/// ---------------------------------------------------------------------------
35	/// ADI (Peaceman--Rachford) scheme
36	/// ---------------------------------------------------------------------------
37	/// One full step of size dt advances T^n -> T^{n+1} in two half-steps of dt/2,
38	/// each implicit in ONE direction and explicit in the other. With the central
39	/// second differences
40	///     dxx T_{ij} = (T_{i-1,j} - 2 T_{ij} + T_{i+1,j}) / hx^2,
41	///     dyy T_{ij} = (T_{i,j-1} - 2 T_{ij} + T_{i,j+1}) / hy^2,
42	/// and rx = alpha*(dt/2)/hx^2, ry = alpha*(dt/2)/hy^2:
43	///
44	///   Half-step 1 (implicit in x, explicit in y) -> intermediate field T*:
45	///       (I - rx*Dxx) T*  =  (I + ry*Dyy) T^n
46	///     i.e. for each row j, the unknowns T*_{.,j} satisfy a TRIDIAGONAL system
47	///       -rx T*_{i-1,j} + (1+2rx) T*_{i,j} - rx T*_{i+1,j}
48	///            = T^n_{i,j} + ry ( T^n_{i,j-1} - 2 T^n_{i,j} + T^n_{i,j+1} ).
49	///
50	///   Half-step 2 (implicit in y, explicit in x) -> new field T^{n+1}:
51	///       (I - ry*Dyy) T^{n+1}  =  (I + rx*Dxx) T*
52	///     i.e. for each column i, the unknowns T^{n+1}_{i,.} satisfy a TRIDIAGONAL
53	///     system
54	///       -ry T^{n+1}_{i,j-1} + (1+2ry) T^{n+1}_{i,j} - ry T^{n+1}_{i,j+1}
55	///            = T*_{i,j} + rx ( T*_{i-1,j} - 2 T*_{i,j} + T*_{i+1,j} ).
56	///
57	/// Each tridiagonal system (one per row / per column) is solved DIRECTLY by the
58	/// Thomas algorithm (below). Both half-steps' operators are SPD / diagonally
59	/// dominant, which is the source of the unconditional stability.
60	///
61	/// ---------------------------------------------------------------------------
62	/// Thomas algorithm (tridiagonal solve)
63	/// ---------------------------------------------------------------------------
64	/// For a system with sub-diagonal a[k], diagonal b[k], super-diagonal c[k] and
65	/// right-hand side d[k] (k = 0..m-1):
66	///   forward sweep:  c'_0 = c_0/b_0, d'_0 = d_0/b_0;  for k>=1
67	///       w = b_k - a_k c'_{k-1};  c'_k = c_k / w;  d'_k = (d_k - a_k d'_{k-1})/w
68	///   back substitution: x_{m-1} = d'_{m-1};  x_k = d'_k - c'_k x_{k+1}.
69	/// O(m) per line; stable for diagonally dominant systems.
70	///
71	/// ---------------------------------------------------------------------------
72	/// Boundary conditions (per edge, applied within the sweeps)
73	/// ---------------------------------------------------------------------------
74	///   * Dirichlet (T = value): the edge nodes are pinned to `value` in BOTH
75	///     sweeps (the trivial 1-row equation T = value), which keeps the boundary
76	///     fixed across the step for time-constant data.
77	///   * Neumann (dT/dn = g, outward; g = 0 == insulated): central ghost node,
78	///     e.g. at x = 0:  T_{-1,j} = T_{1,j} + 2 hx g  (since dT/dn = -dT/dx),
79	///     folded into the tridiagonal end row when the edge is normal to the sweep
80	///     direction, and into the explicit RHS otherwise. Reflective (g=0) ghosts
81	///     make the all-insulated problem exactly energy-conserving.
82	///   * Robin (-k dT/dn = h (T - T_inf)  =>  dT/dn = -B (T - T_inf), B = h/k):
83	///     central ghost node, e.g. at x = 0:
84	///       T_{-1,j} = T_{1,j} - 2 hx B (T_{0,j} - T_inf),
85	///     which adds 2 rx hx B to the end-row diagonal and 2 rx hx B T_inf to the
86	///     RHS (1st/2nd-order accurate and unconditionally stable).
87	///   * Corners (a node shared by two edges): if either incident edge is
88	///     Dirichlet the corner takes that Dirichlet value (left/right edges take
89	///     precedence over bottom/top); otherwise the corner is a genuine unknown
90	///     carrying ghost nodes in both directions.
91	class HeatSolver {
92	public:
93	    explicit HeatSolver(const Plate& plate) : plate_(plate) {}
94	
95	    /// Advance an existing temperature field by ONE full ADI step of size `dt`
96	    /// (the two implicit half-sweeps above, including all boundary handling).
97	    /// `field` is the current state (length num_nodes(), row-major); the
98	    /// returned vector is the state at the next time level. Throws
99	    /// std::invalid_argument on a field-size mismatch or if dt <= 0. The scheme
100	    /// is unconditionally stable, so there is NO Fourier-number guard.
101	    std::vector<double> step(const std::vector<double>& field, double dt) const;
102	
103	    /// Run the full simulation: start from the plate's initial field and take
104	    /// `n_steps` steps of size `dt`. Calls Plate::validate() first and
105	    /// propagates any std::runtime_error it throws. Throws
106	    /// std::invalid_argument if dt <= 0. Returns
107	    /// HeatResult{ temperature, steps = n_steps, time = n_steps*dt }; zero
108	    /// steps returns the initial field at time 0.
109	    HeatResult solve(double dt, std::size_t n_steps) const;
110	
111	    // ---- Diagnostics (fully implemented; do not depend on step()/solve()).
112	
113	    /// Fourier numbers per direction for a HALF step dt/2 used by ADI:
114	    ///   rx = alpha*(dt/2)/hx^2, ry = alpha*(dt/2)/hy^2.
115	    /// Reported purely for diagnostics; they impose NO stability constraint.
116	    double fourier_x(double dt) const;
117	    double fourier_y(double dt) const;
118	
119	    /// Total thermal-energy proxy (2D trapezoidal integral of T), forwarded to
120	    /// Plate::total_energy. Used to check the energy-conservation invariant for
121	    /// a fully insulated plate. Throws std::invalid_argument on size mismatch.
122	    double total_energy(const std::vector<double>& field) const;
123	
124	private:
125	    const Plate& plate_;
126	};
127	
128	} // namespace heat1d
129	
130	#endif // HEAT1D_SOLVER_HPP
131

1	#ifndef HEAT1D_TYPES_HPP
2	#define HEAT1D_TYPES_HPP
3	
4	#include <cstddef>
5	
6	namespace heat1d {
7	
8	/// Boundary-condition kind on one edge of the plate.
9	///
10	///  - Dirichlet: the edge temperature is prescribed (held fixed at `value`).
11	///  - Neumann:   a prescribed outward normal gradient, dT/dn = `value`
12	///               (`value == 0` is the insulated / zero-flux case). Implemented
13	///               with a reflective (central) ghost node. Insulating all four
14	///               edges conserves the plate's total thermal energy.
15	///  - Robin:     a convective (mixed) condition, -k dT/dn = h (T - T_inf),
16	///               i.e.  dT/dn = -(h/k) (T - T_inf). Here the model is written in
17	///               terms of the thermal diffusivity, so `h` is carried as a
18	///               NON-DIMENSIONAL Biot-like coefficient B = h/k with units of
19	///               1/length, and `T_inf` is the ambient temperature. The
20	///               condition is dT/dn = -B (T - T_inf), discretized with a
21	///               central ghost node (2nd-order, unconditionally stable inside
22	///               the implicit ADI sweeps).
23	enum class BCKind {
24	    Dirichlet,
25	    Neumann,
26	    Robin,
27	};
28	
29	/// A boundary condition on one edge of the plate.
30	///
31	/// Field meaning by kind:
32	///   - Dirichlet: `value` is the fixed edge temperature [K]; `h`, `T_inf`
33	///                ignored.
34	///   - Neumann:   `value` is the prescribed OUTWARD normal gradient dT/dn
35	///                [K/m] (0 == insulated); `h`, `T_inf` ignored.
36	///   - Robin:     `h` is the Biot-like coefficient B = h_conv/k [1/m] (>= 0)
37	///                and `T_inf` is the ambient temperature [K]; `value` ignored.
38	struct BoundaryCondition {
39	    BCKind kind{BCKind::Dirichlet};
40	    double value{0.0};  ///< Dirichlet temperature or Neumann outward gradient.
41	    double h{0.0};      ///< Robin Biot-like coefficient B = h/k [1/m].
42	    double T_inf{0.0};  ///< Robin ambient temperature [K].
43	};
44	
45	/// Convenience builders for the three boundary-condition kinds.
46	inline BoundaryCondition dirichlet(double value) {
47	    return BoundaryCondition{BCKind::Dirichlet, value, 0.0, 0.0};
48	}
49	inline BoundaryCondition neumann(double gradient = 0.0) {
50	    return BoundaryCondition{BCKind::Neumann, gradient, 0.0, 0.0};
51	}
52	inline BoundaryCondition insulated() {
53	    return BoundaryCondition{BCKind::Neumann, 0.0, 0.0, 0.0};
54	}
55	inline BoundaryCondition robin(double biot, double T_inf) {
56	    return BoundaryCondition{BCKind::Robin, 0.0, biot, T_inf};
57	}
58	
59	} // namespace heat1d
60	
61	#endif // HEAT1D_TYPES_HPP
62

1	#ifndef HEAT1D_ROD_HPP
2	#define HEAT1D_ROD_HPP
3	
4	#include <cstddef>
5	#include <vector>
6	
7	#include "heat1d/types.hpp"
8	
9	namespace heat1d {
10	
11	/// A rectangular plate for 2D transient heat conduction, discretized on a
12	/// uniform structured grid.
13	///
14	/// This is a plain data container plus light validation / derived-quantity
15	/// helpers. It is fully implemented; the time-stepping numerical scheme (the
16	/// implicit ADI solver) lives in HeatSolver (see solver.hpp).
17	///
18	/// Geometry / grid:
19	///   - The plate occupies the rectangle [0, Lx] x [0, Ly].
20	///   - It is divided into Nx equal cells in x and Ny equal cells in y, giving
21	///     (Nx+1)*(Ny+1) grid points at (x_i, y_j) = (i*hx, j*hy) with
22	///     hx = Lx/Nx, hy = Ly/Ny, i = 0..Nx, j = 0..Ny.
23	///   - Fields are stored FLATTENED, ROW-MAJOR in x: the value at node (i, j)
24	///     lives at linear index  idx(i, j) = i + (Nx+1)*j. Use index(i, j).
25	///
26	/// Physics:
27	///   - alpha is the (constant, positive) thermal diffusivity [m^2/s].
28	///   - The governing PDE is the 2D heat equation
29	///         dT/dt = alpha * ( d2T/dx2 + d2T/dy2 ).
30	///   - The four edges carry independent boundary conditions:
31	///         left   : x = 0      (outward normal -x)
32	///         right  : x = Lx     (outward normal +x)
33	///         bottom : y = 0      (outward normal -y)
34	///         top    : y = Ly     (outward normal +y)
35	///   - initial() is the initial temperature field T(x, y, 0), one value per
36	///     node, flattened row-major.
37	class Plate {
38	public:
39	    /// Construct a plate of size `Lx` x `Ly` [m] with thermal diffusivity
40	    /// `alpha` [m^2/s], divided into `Nx` x `Ny` equal cells, with the given
41	    /// edge boundary conditions. The initial field is zero until set with
42	    /// set_initial(). Throws std::invalid_argument if Lx, Ly, or alpha is
43	    /// non-positive or if Nx < 2 or Ny < 2.
44	    Plate(double Lx, double Ly, double alpha, std::size_t Nx, std::size_t Ny,
45	          BoundaryCondition left, BoundaryCondition right,
46	          BoundaryCondition bottom, BoundaryCondition top);
47	
48	    /// Set the initial temperature field (one value per node, flattened
49	    /// row-major). Throws std::invalid_argument if field.size() != num_nodes().
50	    void set_initial(const std::vector<double>& field);
51	
52	    /// Set a SPATIALLY VARYING thermal diffusivity field (one positive value per
53	    /// node, flattened row-major). When set, the medium is heterogeneous and the
54	    /// solver must use alpha_at(i, j) per node (the scalar alpha() is then only a
55	    /// nominal reference). Throws std::invalid_argument if the size does not
56	    /// match the grid or any value is non-positive.
57	    void set_diffusivity(const std::vector<double>& field);
58	
59	    double Lx() const { return Lx_; }
60	    double Ly() const { return Ly_; }
61	    /// Nominal (scalar) thermal diffusivity. With a variable-diffusivity field
62	    /// set this is just a reference value; use alpha_at(i, j) for the local one.
63	    double alpha() const { return alpha_; }
64	
65	    /// Local thermal diffusivity at node (i, j): the variable field value if one
66	    /// was set with set_diffusivity(), otherwise the scalar alpha(). Always
67	    /// positive. Throws std::out_of_range on a bad index.
68	    double alpha_at(std::size_t i, std::size_t j) const;
69	
70	    /// True iff a spatially varying diffusivity field has been set.
71	    bool has_variable_diffusivity() const { return !diffusivity_.empty(); }
72	    std::size_t Nx() const { return Nx_; }
73	    std::size_t Ny() const { return Ny_; }
74	
75	    /// Nodes per direction (Nx+1 in x, Ny+1 in y).
76	    std::size_t nx_nodes() const { return Nx_ + 1; }
77	    std::size_t ny_nodes() const { return Ny_ + 1; }
78	    /// Total number of grid points = (Nx+1)*(Ny+1).
79	    std::size_t num_nodes() const { return (Nx_ + 1) * (Ny_ + 1); }
80	
81	    /// Grid spacings.
82	    double hx() const { return Lx_ / static_cast<double>(Nx_); }
83	    double hy() const { return Ly_ / static_cast<double>(Ny_); }
84	
85	    /// Flattened (row-major in x) linear index of node (i, j).
86	    /// Throws std::out_of_range if i > Nx or j > Ny.
87	    std::size_t index(std::size_t i, std::size_t j) const;
88	
89	    /// Coordinates of node (i, j). Throw std::out_of_range on a bad index.
90	    double node_x(std::size_t i) const;
91	    double node_y(std::size_t j) const;
92	
93	    const BoundaryCondition& left() const { return left_; }
94	    const BoundaryCondition& right() const { return right_; }
95	    const BoundaryCondition& bottom() const { return bottom_; }
96	    const BoundaryCondition& top() const { return top_; }
97	    const std::vector<double>& initial() const { return initial_; }
98	
99	    /// True iff all four edges are insulated (Neumann zero-flux). In this case
100	    /// the plate's total thermal energy is conserved by the ADI scheme.
101	    bool is_insulated() const;
102	
103	    /// Throws std::runtime_error if the model is ill-formed for a solve:
104	    ///   - the initial field length does not match the grid, or
105	    ///   - a Robin edge carries a NEGATIVE Biot coefficient h (which would be a
106	    ///     surface heat source rather than convective loss -- unphysical and
107	    ///     not guaranteed stable).
108	    void validate() const;
109	
110	    /// Total thermal-energy proxy E = double integral of T over the plate,
111	    /// evaluated from `field` by the 2D trapezoidal (tensor-product) rule.
112	    /// `field` must have length num_nodes(). Used to check the
113	    /// energy-conservation invariant for a fully insulated plate. Throws
114	    /// std::invalid_argument on a size mismatch.
115	    double total_energy(const std::vector<double>& field) const;
116	
117	private:
118	    double Lx_{0.0};
119	    double Ly_{0.0};
120	    double alpha_{0.0};
121	    std::size_t Nx_{0};
122	    std::size_t Ny_{0};
123	    BoundaryCondition left_{};
124	    BoundaryCondition right_{};
125	    BoundaryCondition bottom_{};
126	    BoundaryCondition top_{};
127	    std::vector<double> initial_;
128	    std::vector<double> diffusivity_; // empty => uniform scalar alpha_
129	};
130	
131	} // namespace heat1d
132	
133	#endif // HEAT1D_ROD_HPP
134

/opt/heat1d/include:
heat1d

/opt/heat1d/include/heat1d:
io.hpp
rod.hpp
solver.hpp
types.hpp

/opt/heat1d/src:
io.cpp
main.cpp
rod.cpp
solver.cpp

/opt/heat1d/tests:
CMakeLists.txt
test_dev.cpp
test_util.hpp

[stdout]
/opt/heat1d/include:
heat1d

/opt/heat1d/include/heat1d:
io.hpp
rod.hpp
solver.hpp
types.hpp

/opt/heat1d/src:
io.cpp
main.cpp
rod.cpp
solver.cpp

/opt/heat1d/tests:
CMakeLists.txt
test_dev.cpp
test_util.hpp

1	#include "heat1d/rod.hpp"
2	
3	#include <stdexcept>
4	#include <string>
5	
6	namespace heat1d {
7	
8	Plate::Plate(double Lx, double Ly, double alpha, std::size_t Nx, std::size_t Ny,
9	             BoundaryCondition left, BoundaryCondition right,
10	             BoundaryCondition bottom, BoundaryCondition top)
11	    : Lx_(Lx),
12	      Ly_(Ly),
13	      alpha_(alpha),
14	      Nx_(Nx),
15	      Ny_(Ny),
16	      left_(left),
17	      right_(right),
18	      bottom_(bottom),
19	      top_(top),
20	      initial_((Nx + 1) * (Ny + 1), 0.0) {
21	    if (Lx_ <= 0.0 || Ly_ <= 0.0) {
22	        throw std::invalid_argument("Plate: Lx and Ly must be positive");
23	    }
24	    if (alpha_ <= 0.0) {
25	        throw std::invalid_argument("Plate: alpha must be positive");
26	    }
27	    if (Nx_ < 2 || Ny_ < 2) {
28	        throw std::invalid_argument("Plate: need at least 2 cells per direction");
29	    }
30	}
31	
32	void Plate::set_initial(const std::vector<double>& field) {
33	    if (field.size() != num_nodes()) {
34	        throw std::invalid_argument("Plate::set_initial: field size mismatch");
35	    }
36	    initial_ = field;
37	}
38	
39	void Plate::set_diffusivity(const std::vector<double>& field) {
40	    if (field.size() != num_nodes()) {
41	        throw std::invalid_argument("Plate::set_diffusivity: field size mismatch");
42	    }
43	    for (double a : field) {
44	        if (a <= 0.0) {
45	            throw std::invalid_argument(
46	                "Plate::set_diffusivity: diffusivity must be positive");
47	        }
48	    }
49	    diffusivity_ = field;
50	}
51	
52	double Plate::alpha_at(std::size_t i, std::size_t j) const {
53	    if (i > Nx_ || j > Ny_) {
54	        throw std::out_of_range("Plate::alpha_at: node index out of range");
55	    }
56	    if (diffusivity_.empty()) {
57	        return alpha_;
58	    }
59	    return diffusivity_[i + (Nx_ + 1) * j];
60	}
61	
62	std::size_t Plate::index(std::size_t i, std::size_t j) const {
63	    if (i > Nx_ || j > Ny_) {
64	        throw std::out_of_range("Plate::index: node index out of range");
65	    }
66	    return i + (Nx_ + 1) * j;
67	}
68	
69	double Plate::node_x(std::size_t i) const {
70	    if (i > Nx_) {
71	        throw std::out_of_range("Plate::node_x: node index out of range");
72	    }
73	    return static_cast<double>(i) * hx();
74	}
75	
76	double Plate::node_y(std::size_t j) const {
77	    if (j > Ny_) {
78	        throw std::out_of_range("Plate::node_y: node index out of range");
79	    }
80	    return static_cast<double>(j) * hy();
81	}
82	
83	bool Plate::is_insulated() const {
84	    auto insul = [](const BoundaryCondition& bc) {
85	        return bc.kind == BCKind::Neumann && bc.value == 0.0;
86	    };
87	    return insul(left_) && insul(right_) && insul(bottom_) && insul(top_);
88	}
89	
90	void Plate::validate() const {
91	    if (initial_.size() != num_nodes()) {
92	        throw std::runtime_error("Plate::validate: initial field size mismatch");
93	    }
94	    auto check_robin = [](const BoundaryCondition& bc, const char* edge) {
95	        if (bc.kind == BCKind::Robin && bc.h < 0.0) {
96	            throw std::runtime_error(std::string("Plate::validate: ") + edge +
97	                                     " Robin edge has negative Biot coefficient");
98	        }
99	    };
100	    check_robin(left_, "left");
101	    check_robin(right_, "right");
102	    check_robin(bottom_, "bottom");
103	    check_robin(top_, "top");
104	}
105	
106	double Plate::total_energy(const std::vector<double>& field) const {
107	    if (field.size() != num_nodes()) {
108	        throw std::invalid_argument("Plate::total_energy: field size mismatch");
109	    }
110	    // Tensor-product trapezoidal rule: weight 1/2 on edge nodes, 1/4 at the
111	    // four corners, 1 in the interior, times hx*hy.
112	    const std::size_t nx = nx_nodes();
113	    const std::size_t ny = ny_nodes();
114	    double e = 0.0;
115	    for (std::size_t j = 0; j < ny; ++j) {
116	        const double wy = (j == 0 || j + 1 == ny) ? 0.5 : 1.0;
117	        for (std::size_t i = 0; i < nx; ++i) {
118	            const double wx = (i == 0 || i + 1 == nx) ? 0.5 : 1.0;
119	            e += wx * wy * field[i + nx * j];
120	        }
121	    }
122	    return e * hx() * hy();
123	}
124	
125	} // namespace heat1d
126

1	// Development test subset (NON-AUTHORITATIVE).
2	//
3	// These few checks give you a basic signal while implementing the scheme:
4	// Dirichlet edges get pinned, an insulated uniform field is steady, the
5	// fundamental Fourier mode decays roughly like the analytic solution, an
6	// insulated plate conserves energy, and the exception contract fires. They are
7	// deliberately coarse and incomplete.
8	//
9	// Passing them is NECESSARY but NOT SUFFICIENT: grading uses a stricter hidden
10	// suite (full Dirichlet / nonzero-Neumann / Robin / mixed boundary matrix,
11	// anisotropic grids, much tighter tolerances). Implement the scheme correctly
12	// from the discretization rather than tuning to these tests.
13	
14	#include <algorithm>
15	#include <cmath>
16	#include <vector>
17	
18	#include "heat1d/rod.hpp"
19	#include "heat1d/solver.hpp"
20	#include "test_util.hpp"
21	
22	using heat1d::dirichlet;
23	using heat1d::HeatResult;
24	using heat1d::HeatSolver;
25	using heat1d::insulated;
26	using heat1d::Plate;
27	
28	namespace {
29	double max_of(const std::vector<double>& v) {
30	    return *std::max_element(v.begin(), v.end());
31	}
32	} // namespace
33	
34	HEAT1D_TEST("dev_dirichlet_edges_are_pinned") {
35	    const double Lx = 2.0, Ly = 1.0, alpha = 0.5;
36	    const std::size_t Nx = 12, Ny = 8;
37	    Plate p(Lx, Ly, alpha, Nx, Ny, dirichlet(100.0), dirichlet(300.0),
38	            dirichlet(50.0), dirichlet(75.0));
39	    const std::size_t nx = p.nx_nodes(), ny = p.ny_nodes();
40	    HeatSolver s(p);
41	    std::vector<double> U = s.step(std::vector<double>(p.num_nodes(), 0.0), 0.01);
42	    for (std::size_t j = 1; j + 1 < ny; ++j) {
43	        heat1d_test::expect_near(U[0 + nx * j], 100.0, 1e-9, "left Dirichlet");
44	        heat1d_test::expect_near(U[(nx - 1) + nx * j], 300.0, 1e-9, "right Dirichlet");
45	    }
46	    for (std::size_t i = 1; i + 1 < nx; ++i) {
47	        heat1d_test::expect_near(U[i + nx * 0], 50.0, 1e-9, "bottom Dirichlet");
48	        heat1d_test::expect_near(U[i + nx * (ny - 1)], 75.0, 1e-9, "top Dirichlet");
49	    }
50	}
51	
52	HEAT1D_TEST("dev_insulated_uniform_field_is_steady") {
53	    Plate p(1.3, 0.8, 0.9, 10, 14, insulated(), insulated(), insulated(),
54	            insulated());
55	    HeatSolver s(p);
56	    std::vector<double> U = s.step(std::vector<double>(p.num_nodes(), 42.0), 0.05);
57	    double dev = 0.0;
58	    for (double v : U) dev = std::max(dev, std::fabs(v - 42.0));
59	    heat1d_test::expect(dev < 1e-9, "uniform insulated field stays uniform");
60	}
61	
62	HEAT1D_TEST("dev_fourier_mode_decays_roughly") {
63	    const double Lx = 1.0, Ly = 1.0, alpha = 0.1;
64	    const std::size_t Nx = 40, Ny = 40;
65	    Plate p(Lx, Ly, alpha, Nx, Ny, dirichlet(0.0), dirichlet(0.0),
66	            dirichlet(0.0), dirichlet(0.0));
67	    const std::size_t nx = p.nx_nodes(), ny = p.ny_nodes();
68	    std::vector<double> T0(p.num_nodes());
69	    for (std::size_t j = 0; j < ny; ++j)
70	        for (std::size_t i = 0; i < nx; ++i)
71	            T0[i + nx * j] = std::sin(M_PI * p.node_x(i) / Lx) *
72	                             std::sin(M_PI * p.node_y(j) / Ly);
73	    p.set_initial(T0);
74	    HeatSolver s(p);
75	    const std::size_t steps = 80;
76	    const double dt = 0.4 / steps;
77	    HeatResult res = s.solve(dt, steps);
78	    // Compare to the continuous analytic decay with a LOOSE tolerance.
79	    const double k2 = (M_PI / Lx) * (M_PI / Lx) + (M_PI / Ly) * (M_PI / Ly);
80	    const double decay = std::exp(-alpha * k2 * res.time);
81	    double err = 0.0;
82	    for (std::size_t j = 0; j < ny; ++j)
83	        for (std::size_t i = 0; i < nx; ++i) {
84	            const double exact = decay * std::sin(M_PI * p.node_x(i) / Lx) *
85	                                 std::sin(M_PI * p.node_y(j) / Ly);
86	            err = std::max(err, std::fabs(res.temperature[i + nx * j] - exact));
87	        }
88	    heat1d_test::expect(err < 5e-3, "mode roughly follows analytic decay");
89	    heat1d_test::expect(max_of(res.temperature) < 0.7, "mode measurably decayed");
90	}
91	
92	HEAT1D_TEST("dev_insulated_energy_roughly_conserved") {
93	    const double Lx = 1.0, Ly = 1.0, alpha = 0.01;
94	    const std::size_t Nx = 40, Ny = 40;
95	    Plate p(Lx, Ly, alpha, Nx, Ny, insulated(), insulated(), insulated(),
96	            insulated());
97	    const std::size_t nx = p.nx_nodes(), ny = p.ny_nodes();
98	    std::vector<double> T0(p.num_nodes());
99	    for (std::size_t j = 0; j < ny; ++j)
100	        for (std::size_t i = 0; i < nx; ++i) {
101	            const double sx = p.node_x(i) - 0.5, sy = p.node_y(j) - 0.5;
102	            T0[i + nx * j] = 20.0 + 50.0 * std::exp(-50.0 * (sx * sx + sy * sy));
103	        }
104	    p.set_initial(T0);
105	    HeatSolver s(p);
106	    const double E0 = s.total_energy(T0);
107	    HeatResult res = s.solve(2.0 * p.hx() * p.hx() / alpha, 400);
108	    heat1d_test::expect_rel(s.total_energy(res.temperature), E0, 1e-6,
109	                            "insulated plate roughly conserves energy");
110	}
111	
112	HEAT1D_TEST("dev_variable_diffusivity_layered_steady_basic") {
113	    // A two-material (layered in x) plate, hot/cold Dirichlet left/right,
114	    // insulated top/bottom. Coarse sanity only: the steady field must span the
115	    // wall temperatures, be transverse-uniform, and rise monotonically in x with
116	    // a visibly steeper gradient in the low-diffusivity (left) half. (The exact
117	    // profile is graded by the hidden suite.)
118	    const double L = 1.0, TL = 0.0, TR = 100.0;
119	    const std::size_t N = 16;
120	    Plate p(L, L, 1.0, N, N, dirichlet(TL), dirichlet(TR), insulated(),
121	            insulated());
122	    const std::size_t nx = p.nx_nodes(), ny = p.ny_nodes();
123	    std::vector<double> af(p.num_nodes());
124	    for (std::size_t j = 0; j < ny; ++j)
125	        for (std::size_t i = 0; i < nx; ++i)
126	            af[i + nx * j] = (p.node_x(i) < 0.5 - 1e-12) ? 1.0 : 5.0;
127	    p.set_diffusivity(af);
128	    p.set_initial(std::vector<double>(p.num_nodes(), 50.0));
129	    HeatSolver s(p);
130	    HeatResult res = s.solve(2.0 * p.hx() * p.hx() / 5.0, 8000);
131	    const std::size_t jmid = N / 2;
132	    double prev = -1e9;
133	    bool monotonic = true;
134	    for (std::size_t i = 0; i < nx; ++i) {
135	        const double v = res.temperature[i + nx * jmid];
136	        if (v < prev - 1e-9) monotonic = false;
137	        prev = v;
138	    }
139	    heat1d_test::expect(monotonic, "steady field rises monotonically across layers");
140	    heat1d_test::expect_near(res.temperature[0 + nx * jmid], TL, 1e-3, "cold wall held");
141	    heat1d_test::expect_near(res.temperature[(nx - 1) + nx * jmid], TR, 1e-3, "hot wall held");
142	    // Low-diffusivity (left) half carries the larger temperature drop.
143	    const double Tmid = res.temperature[(N / 2) + nx * jmid];
144	    heat1d_test::expect(Tmid - TL > TR - Tmid,
145	                        "steeper gradient in the low-diffusivity half");
146	}
147	
148	HEAT1D_TEST("dev_exception_contract") {
149	    Plate p(1.0, 1.0, 0.01, 8, 8, dirichlet(0.0), dirichlet(0.0), dirichlet(0.0),
150	            dirichlet(0.0));
151	    p.set_initial(std::vector<double>(p.num_nodes(), 0.0));
152	    HeatSolver s(p);
153	    bool a = false, b = false;
154	    try { s.step(std::vector<double>(p.num_nodes(), 0.0), 0.0); }
155	    catch (const std::invalid_argument&) { a = true; }
156	    try { s.solve(-1.0, 3); } catch (const std::invalid_argument&) { b = true; }
157	    heat1d_test::expect(a, "step(dt<=0) throws invalid_argument");
158	    heat1d_test::expect(b, "solve(dt<=0) throws invalid_argument");
159	}
160	
161	HEAT1D_TEST_MAIN()
162

1	#ifndef HEAT1D_TEST_UTIL_HPP
2	#define HEAT1D_TEST_UTIL_HPP
3	
4	// Minimal dependency-free test harness for the heat1d suite.
5	//
6	// Each test executable defines test functions and registers them with
7	// HEAT1D_TEST(name) { ... }. main() (provided by HEAT1D_TEST_MAIN) runs every
8	// registered test and returns non-zero if any assertion fails, which CTest
9	// interprets as a failure.
10	
11	#include <cmath>
12	#include <cstddef>
13	#include <functional>
14	#include <iostream>
15	#include <string>
16	#include <vector>
17	
18	namespace heat1d_test {
19	
20	struct Case {
21	    std::string name;
22	    std::function<void()> fn;
23	};
24	
25	inline std::vector<Case>& registry() {
26	    static std::vector<Case> cases;
27	    return cases;
28	}
29	
30	struct Registrar {
31	    Registrar(const std::string& name, std::function<void()> fn) {
32	        registry().push_back(Case{name, std::move(fn)});
33	    }
34	};
35	
36	// Thrown by an assertion failure to abort the current test case.
37	struct AssertionError {
38	    std::string message;
39	};
40	
41	inline void expect(bool cond, const std::string& what) {
42	    if (!cond) {
43	        throw AssertionError{what};
44	    }
45	}
46	
47	inline void expect_near(double got, double want, double tol, const std::string& what) {
48	    if (std::fabs(got - want) > tol) {
49	        throw AssertionError{what + " (got " + std::to_string(got) +
50	                             ", want " + std::to_string(want) +
51	                             ", tol " + std::to_string(tol) + ")"};
52	    }
53	}
54	
55	// Relative-error check (with an absolute floor) for quantities whose scale
56	// varies widely across cases.
57	inline void expect_rel(double got, double want, double rel_tol,
58	                       const std::string& what) {
59	    const double denom = std::fabs(want) > 1e-300 ? std::fabs(want) : 1.0;
60	    if (std::fabs(got - want) / denom > rel_tol) {
61	        throw AssertionError{what + " (got " + std::to_string(got) +
62	                             ", want " + std::to_string(want) +
63	                             ", rel_tol " + std::to_string(rel_tol) + ")"};
64	    }
65	}
66	
67	inline int run_all() {
68	    int failures = 0;
69	    for (const auto& c : registry()) {
70	        try {
71	            c.fn();
72	            std::cout << "[ PASS ] " << c.name << "\n";
73	        } catch (const AssertionError& e) {
74	            std::cout << "[ FAIL ] " << c.name << ": " << e.message << "\n";
75	            ++failures;
76	        } catch (const std::exception& e) {
77	            std::cout << "[ FAIL ] " << c.name
78	                      << ": unexpected exception: " << e.what() << "\n";
79	            ++failures;
80	        } catch (...) {
81	            std::cout << "[ FAIL ] " << c.name << ": unknown exception\n";
82	            ++failures;
83	        }
84	    }
85	    std::cout << "----\n"
86	              << (registry().size() - failures) << "/" << registry().size()
87	              << " tests passed\n";
88	    return failures == 0 ? 0 : 1;
89	}
90	
91	} // namespace heat1d_test
92	
93	#define HEAT1D_CONCAT_INNER(a, b) a##b
94	#define HEAT1D_CONCAT(a, b) HEAT1D_CONCAT_INNER(a, b)
95	
96	#define HEAT1D_TEST(NAME)                                                      \
97	    static void HEAT1D_CONCAT(heat1d_test_fn_, __LINE__)();                    \
98	    static ::heat1d_test::Registrar HEAT1D_CONCAT(heat1d_test_reg_,           \
99	                                                  __LINE__){                   \
100	        NAME, &HEAT1D_CONCAT(heat1d_test_fn_, __LINE__)};                      \
101	    static void HEAT1D_CONCAT(heat1d_test_fn_, __LINE__)()
102	
103	#define HEAT1D_TEST_MAIN()                                                     \
104	    int main() { return ::heat1d_test::run_all(); }
105	
106	#endif // HEAT1D_TEST_UTIL_HPP
107

BUILD OK
Test project /opt/heat1d/build
    Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.15 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.15 sec

[stdout]
BUILD OK
Test project /opt/heat1d/build
    Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.15 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.15 sec

1) discrete eigenmode 1-step max err = 4.406e-16 (factor g=0.092208)
2) nonzero-Neumann steady max err = 2.576e-14
3) Robin steady max err = 3.268e-13 (slope=50.00000)
4) high-contrast interface T: got 99.89754 want 99.90010 (err 2.556e-03)
5) diagonal symmetry max err = 6.939e-18
6) insulated energy: E0=0.3000732422 Ef=0.3000732422 reldiff=5.918e-13

[stdout]
1) discrete eigenmode 1-step max err = 4.406e-16 (factor g=0.092208)
2) nonzero-Neumann steady max err = 2.576e-14
3) Robin steady max err = 3.268e-13 (slope=50.00000)
4) high-contrast interface T: got 99.89754 want 99.90010 (err 2.556e-03)
5) diagonal symmetry max err = 6.939e-18
6) insulated energy: E0=0.3000732422 Ef=0.3000732422 reldiff=5.918e-13

A) high-contrast steady flux range [204.912784,204.912784] spread=2.103e-08
B) extreme dt=1e6: finite=1 maxabs=1.000e+00 (should be tiny, decayed)
C) Dirichlet decay: max-monotone=1 stayed-bounded=1
D) all-Robin steady -> [25.000000, 25.000000] (want ~25.0)

[stdout]
A) high-contrast steady flux range [204.912784,204.912784] spread=2.103e-08
B) extreme dt=1e6: finite=1 maxabs=1.000e+00 (should be tiny, decayed)
C) Dirichlet decay: max-monotone=1 stayed-bounded=1
D) all-Robin steady -> [25.000000, 25.000000] (want ~25.0)

Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.15 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.15 sec

[stdout]
    Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.15 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.15 sec