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

cmake_minimum_required(VERSION 3.16)
project(heat1d LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

if(NOT CMAKE_BUILD_TYPE)
  set(CMAKE_BUILD_TYPE Release)
endif()

add_compile_options(-Wall -Wextra)

# --- Core library -----------------------------------------------------------
add_library(heat1d
  src/rod.cpp
  src/solver.cpp
  src/io.cpp
)
target_include_directories(heat1d PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)

# --- Demo CLI ---------------------------------------------------------------
add_executable(heat1d_cli src/main.cpp)
target_link_libraries(heat1d_cli PRIVATE heat1d)

# --- Tests ------------------------------------------------------------------
enable_testing()
add_subdirectory(tests)
---TESTS---
# Development test subset (non-authoritative). Each test file is a standalone
# executable using the tiny header-only harness in test_util.hpp. CTest runs it
# and treats a non-zero exit as failure. Grading replaces this directory with a
# stricter hidden suite (see the task tests/ verifier bundle).

add_executable(test_dev test_dev.cpp)
target_link_libraries(test_dev PRIVATE heat1d)
target_include_directories(test_dev PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME dev COMMAND test_dev)

[stdout]
cmake_minimum_required(VERSION 3.16)
project(heat1d LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

if(NOT CMAKE_BUILD_TYPE)
  set(CMAKE_BUILD_TYPE Release)
endif()

add_compile_options(-Wall -Wextra)

# --- Core library -----------------------------------------------------------
add_library(heat1d
  src/rod.cpp
  src/solver.cpp
  src/io.cpp
)
target_include_directories(heat1d PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/include)

# --- Demo CLI ---------------------------------------------------------------
add_executable(heat1d_cli src/main.cpp)
target_link_libraries(heat1d_cli PRIVATE heat1d)

# --- Tests ------------------------------------------------------------------
enable_testing()
add_subdirectory(tests)
---TESTS---
# Development test subset (non-authoritative). Each test file is a standalone
# executable using the tiny header-only harness in test_util.hpp. CTest runs it
# and treats a non-zero exit as failure. Grading replaces this directory with a
# stricter hidden suite (see the task tests/ verifier bundle).

add_executable(test_dev test_dev.cpp)
target_link_libraries(test_dev PRIVATE heat1d)
target_include_directories(test_dev PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME dev COMMAND test_dev)

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

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.14 sec

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

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.14 sec

Exit code 1
  ok   discrete ADI mode decay matches g^n exactly
  ok   nonzero-Neumann linear steady state
  ok   Robin convective steady state
  FAIL heterogeneous interface temp (harmonic) (got 79.3814433 want 80)
  ok   diagonal symmetry preserved
  ok   energy conserved checkerboard huge dt
  ok   no blowup at huge dt
  ok   discrete maximum principle

FAILURES (1 failures)

[error] tool reported failure

max|num-discrete_steady| = 1.421e-13  (Tmid num=79.381443 exact=79.381443)
PASS: matches true discrete steady

[stdout]
max|num-discrete_steady| = 1.421e-13  (Tmid num=79.381443 exact=79.381443)
PASS: matches true discrete steady

Exit code 1
/tmp/v5.cpp: In function 'int main()':
/tmp/v5.cpp:41:41: error: conversion from 'std::vector<double>' to non-scalar type 'heat1d::HeatResult' requested
   41 |     HeatSolver s(p); HeatResult r=s.step(p.initial(),0.01);
      |                                   ~~~~~~^~~~~~~~~~~~~~~~~~
/tmp/v5.cpp:50:41: error: conversion from 'std::vector<double>' to non-scalar type 'heat1d::HeatResult' requested
   50 |     HeatSolver s(p); HeatResult r=s.step(p.initial(),0.01);
      |                                   ~~~~~~^~~~~~~~~~~~~~~~~~

[error] tool reported failure

Exit code 1
  ok   all-Robin equilibrium -> Tinf
  FAIL 2D Dirichlet Laplace steady ~ analytic (0.7181112619 vs 0)
  ok   corner left-precedence over bottom
  ok   corner takes bottom Dirichlet when left non-Dir
  FAIL hetero+huge dt decays to 0 steady (stable) (0.9999977216 vs 0)

FAIL (2)

[error] tool reported failure

interior discrete-Laplace residual = 5.912e-10
x-symmetry error = 1.137e-13
center temp = 25.000000 (continuous harmonic answer = 25.0)
PASS

[stdout]
interior discrete-Laplace residual = 5.912e-10
x-symmetry error = 1.137e-13
center temp = 25.000000 (continuous harmonic answer = 25.0)
PASS

no solver.cpp warnings
---
    Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.13 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.13 sec

[stdout]
no solver.cpp warnings
---
    Start 1: dev
1/1 Test #1: dev ..............................   Passed    0.13 sec

100% tests passed, 0 tests failed out of 1

Total Test time (real) =   0.13 sec