From bb6badab5a880467d044a23399172a45c5f313bb Mon Sep 17 00:00:00 2001 From: Ed Date: Tue, 18 Aug 2026 13:02:44 +0100 Subject: [PATCH 1/2] Forward coefficient functions to the asynchronous assembly path `solve()` passes `this.coefficientFunctions` to `assembleHeatConductionMat`, but `solveAsync()` called the same assembler with only the mesh and the boundary conditions. Any model configured with spatially varying coefficients silently fell back to a uniform conductivity of 1 and a heat source of 0 when solved through the asynchronous path. Pass the coefficients at that call site as well, so both paths assemble the same system. Refs #82 --- src/FEAScript.js | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/FEAScript.js b/src/FEAScript.js index 4977dba..e8f07b6 100644 --- a/src/FEAScript.js +++ b/src/FEAScript.js @@ -288,7 +288,11 @@ export class FEAScriptModel { basicLog(`Using solver: ${this.solverConfig}`); if (this.solverConfig === "heatConductionScript") { - ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(meshData, this.boundaryConditions)); + ({ jacobianMatrix, residualVector } = assembleHeatConductionMat( + meshData, + this.boundaryConditions, + this.coefficientFunctions, + )); if (this.solverMethod === "jacobi-gpu") { const { solutionVector: x } = await solveLinearSystemAsync( From 65e51e000007b5ad87e109cb965f4468a78c8aa0 Mon Sep 17 00:00:00 2001 From: Ed Date: Tue, 18 Aug 2026 15:48:10 +0100 Subject: [PATCH 2/2] Add regression tests for spatially varying thermal coefficients The existing regression tests pass no `coefficientFunctions`, so the coefficients resolve to a uniform conductivity of 1 and a heat source of 0. Both are then invisible to the result: multiplying by 1 and adding 0 leaves the assembled system identical to one that never read them. Nothing in the suite could distinguish correct coefficient handling from none at all. These two tests assert closed-form solutions rather than stored reference values, which for a new feature would only record whatever the code produced when the test was written. Each case is chosen so the finite element solution is exact at the nodes, giving a tolerance of 1e-10 instead of 1e-4 and expected values that never need re-deriving when the mesh or element order changes. 1D covers a uniform source against T = x(1 - x)/2, whose exact solution lies outside the finite element space and so pins the quadrature of the source term; the manufactured solution T = x under k = 1 + x and Q = -1 for both element orders, which pins the Gauss point as the evaluation point; the frontal assembler against the matrix assembler; and the coefficient forwarding in `solveAsync`. 2D covers T = x and its rotation T = y, since the 2D assembler is a separate implementation and the 1D path calls the coefficients with x alone. Confirmed to have teeth by mutation: swapping x and y in the 2D assembler moves both 2D cases from 1e-15 to 1e-1, and dropping the coefficients from the `solveAsync` call site fails the 1D suite. Refs #82 --- .../REGRESSION.md | 59 +++++ .../regression.test.js | 211 ++++++++++++++++++ .../REGRESSION.md | 52 +++++ .../regression.test.js | 155 +++++++++++++ 4 files changed, 477 insertions(+) create mode 100644 tests/regression/HeatConduction1DVaryingCoefficients/REGRESSION.md create mode 100644 tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js create mode 100644 tests/regression/HeatConduction2DVaryingCoefficients/REGRESSION.md create mode 100644 tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js diff --git a/tests/regression/HeatConduction1DVaryingCoefficients/REGRESSION.md b/tests/regression/HeatConduction1DVaryingCoefficients/REGRESSION.md new file mode 100644 index 0000000..edbfe0c --- /dev/null +++ b/tests/regression/HeatConduction1DVaryingCoefficients/REGRESSION.md @@ -0,0 +1,59 @@ +# Regression Test — HeatConduction1DVaryingCoefficients + +## Purpose + +This test guards the spatially varying `thermalConductivity` and `heatSource` coefficients of +`heatConductionScript` in 1D, covering the matrix assembler, the frontal assembler, and the +coefficient forwarding performed by `solveAsync`. + +It differs from the other regression tests in that no stored reference value is used. A stored +value can only record whatever the code produced when the test was written; here each case has a +closed-form solution of the underlying PDE, and the setups are chosen so the finite element +solution is exact at the nodes. That allows a tolerance of `1e-10` instead of `1e-4`, and the +expected values never need re-deriving when the mesh or element order changes. + +## Problem setup + +Common to every case: domain `x ∈ [0, 1]`, 8 elements, Dirichlet boundaries. The 1D `convection` +condition is avoided so that the expected solution is unambiguous. + +| Case | k(x) | Q(x) | Boundaries | Exact solution | Elements | +| ---- | ----- | ---- | ---------------- | ---------------- | ----------------- | +| 1 | 1 | 1 | T(0) = T(1) = 0 | x (1 − x) / 2 | linear | +| 2 | 1 + x | −1 | T(0) = 0, T(1) = 1 | x | linear, quadratic | +| 3 | 1 + x | 5x | T(0) = 0, T(1) = 1 | frontal vs `lusolve` | linear | +| 4 | counter | counter | T(0) = 0, T(1) = 1 | coefficients reach `solveAsync` | linear | + +Case 1 is the only one whose exact solution lies outside the finite element space, so it is what +pins the quadrature of the source term; case 2's `T = x` would still be reproduced by an +under-integrated source. Case 2 is also what pins the evaluation point, as it fails if the +conductivity is sampled anywhere other than the Gauss points. + +Case 4 cannot be driven end to end, since `jacobi-gpu` requires a WebGPU compute engine. +Assembly happens before the solver method is branched on, so coefficients that count their own +invocations are enough to prove they reach the assembler. + +## Expected values + +Every nodal temperature must match the closed-form solution to within `1e-10`. Observed largest +deviations are of order `1e-15`. + +## How to run + +From the repository root: + +```bash +node tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js +``` + +A passing run prints five `PASS:` lines and `5 passed, 0 failed.`; a failing run prints `FAIL:` +with the largest deviation and its node, and exits with code 1. + +## After modifying the code + +| Situation | Action | +| ---------------------------------------------------- | ----------------------------------------------------------------------------------------------- | +| Bug fix that should not change results | Run the test — it must still pass. | +| Change to the coefficient API | Update the cases; the analytical solutions themselves stay valid. | +| Intentional change to quadrature or element mapping | The expected values do not move. If a case now fails, the change altered the physics, not the reference. | +| New assembler or solver path reading the coefficients | Add a case for it here, as cases 3 and 4 do for the frontal and asynchronous paths. | diff --git a/tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js b/tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js new file mode 100644 index 0000000..f4a3d7f --- /dev/null +++ b/tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js @@ -0,0 +1,211 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Lightweight Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +/** + * Regression test for HeatConduction1DVaryingCoefficients + * + * Guards the spatially varying `thermalConductivity` and `heatSource` coefficients of + * heatConductionScript in 1D + * + * Unlike the other regression tests the expected values are not stored reference numbers + * but closed-form solutions of the underlying PDE. Each case is chosen so the finite + * element solution is exact at the nodes, which allows a tolerance of 1e-10 rather than + * the 1e-4 used where a stored value is compared. + * + * Run: node tests/regression/HeatConduction1DVaryingCoefficients/regression.test.js (or npm test) + */ + +import * as mathjs from "mathjs"; +import { FEAScriptModel } from "../../../src/FEAScript.js"; +import { basicLog, errorLog } from "../../../src/utilities/logging.js"; + +// FEAScript.js references `math` as a global (loaded via CDN in browser). +// Set it here before any solve() call. +globalThis.math = mathjs; + +const TOLERANCE = 1e-10; + +function runSimulation(coefficientFunctions, elementOrder, boundaryConditions, solverMethod = "lusolve") { + const model = new FEAScriptModel(); + + model.setModelConfig("heatConductionScript", { coefficientFunctions }); + model.setMeshConfig({ + meshDimension: "1D", + elementOrder, + numElementsX: 8, + maxX: 1, + }); + + Object.entries(boundaryConditions).forEach(([boundaryKey, condition]) => { + model.addBoundaryCondition(boundaryKey, condition); + }); + model.setSolverMethod(solverMethod); + + const { solutionVector, nodesCoordinates } = model.solve(); + + // solutionVector from math.lusolve is a nested array: [[T0], [T1], ...] + return { + temperatures: solutionVector.map((value) => (Array.isArray(value) ? value[0] : value)), + nodesXCoordinates: nodesCoordinates.nodesXCoordinates, + }; +} + +let passed = 0; +let failed = 0; + +function assert(condition, message) { + if (!condition) { + errorLog(`FAIL: ${message}`); + failed++; + } else { + basicLog(`PASS: ${message}`); + passed++; + } +} + +/** + * Function to assert that every nodal temperature matches an analytical solution + * @param {string} label - Description of the case under test + * @param {object} result - Object containing the computed temperatures and node coordinates + * @param {function} analyticalSolution - Function returning the exact temperature at a coordinate + */ +function assertMatchesAnalyticalSolution(label, result, analyticalSolution) { + const { temperatures, nodesXCoordinates } = result; + + let maxError = 0; + let maxErrorNodeIndex = 0; + for (let nodeIndex = 0; nodeIndex < temperatures.length; nodeIndex++) { + const error = Math.abs(temperatures[nodeIndex] - analyticalSolution(nodesXCoordinates[nodeIndex])); + if (error > maxError) { + maxError = error; + maxErrorNodeIndex = nodeIndex; + } + } + + assert( + maxError < TOLERANCE, + `${label}: largest nodal deviation ${maxError.toExponential(3)} at ` + + `x = ${nodesXCoordinates[maxErrorNodeIndex]} (tolerance ${TOLERANCE})`, + ); +} + +basicLog(""); +basicLog("================================"); +basicLog("Starting regression test for solid heat transfer in 1D with varying coefficients..."); + +/** + * Case 1 - uniform heat source + * + * With k = 1 and Q = 1 on [0, 1] and T = 0 at both ends, div(k * grad(T)) + Q = 0 reduces + * to T'' = -1, so T(x) = x * (1 - x) / 2. The solution is quadratic while the elements are + * linear, so this is the case that pins the quadrature of the source term. + */ +assertMatchesAnalyticalSolution( + "Uniform heat source, linear elements", + runSimulation({ heatSource: 1 }, "linear", { 0: ["constantTemp", 0], 1: ["constantTemp", 0] }), + (x) => (x * (1 - x)) / 2, +); + +/** + * Case 2 - conductivity and heat source together (method of manufactured solutions) + * + * Picking k(x) = 1 + x and Q = -1 makes T(x) = x an exact solution, since + * div(k * grad(T)) + Q = d(1 + x)/dx - 1 = 0. Imposing T = 0 and T = 1 at the two ends + * therefore has to reproduce the identity function. + * + * This case fails if the conductivity is evaluated anywhere other than the Gauss points, + * so it pins down the evaluation point as well as the coefficient itself. + */ +for (const elementOrder of ["linear", "quadratic"]) { + assertMatchesAnalyticalSolution( + `Manufactured solution T = x, ${elementOrder} elements`, + runSimulation({ thermalConductivity: (x) => 1 + x, heatSource: -1 }, elementOrder, { + 0: ["constantTemp", 0], + 1: ["constantTemp", 1], + }), + (x) => x, + ); +} + +/** + * Case 3 - the frontal assembler must agree with the matrix assembler + * + * `assembleHeatConductionFront` carries its own copy of the coefficient handling, so it is + * compared against `lusolve` on a problem where both coefficients vary. + */ +{ + const coefficientFunctions = { thermalConductivity: (x) => 1 + x, heatSource: (x) => 5 * x }; + const boundaryConditions = { 0: ["constantTemp", 0], 1: ["constantTemp", 1] }; + + const luResult = runSimulation(coefficientFunctions, "linear", boundaryConditions); + const frontalResult = runSimulation(coefficientFunctions, "linear", boundaryConditions, "frontal"); + + let maxDifference = 0; + for (let nodeIndex = 0; nodeIndex < luResult.temperatures.length; nodeIndex++) { + maxDifference = Math.max( + maxDifference, + Math.abs(luResult.temperatures[nodeIndex] - frontalResult.temperatures[nodeIndex]), + ); + } + + assert( + maxDifference < TOLERANCE, + `Frontal assembler matches lusolve: largest difference ${maxDifference.toExponential(3)} ` + + `(tolerance ${TOLERANCE})`, + ); +} + +/** + * Case 4 - the asynchronous path forwards the coefficients + * + * `solveAsync` holds a second call into `assembleHeatConductionMat`, which is easy to miss + * when the signature changes. It cannot be driven end to end here because `jacobi-gpu` + * needs a WebGPU compute engine, but assembly happens before the solver method is branched + * on, so a coefficient that counts its own invocations is enough to prove the coefficients + * reach the assembler. + */ +{ + let thermalConductivityCalls = 0; + let heatSourceCalls = 0; + + const model = new FEAScriptModel(); + model.setModelConfig("heatConductionScript", { + coefficientFunctions: { + thermalConductivity: () => { + thermalConductivityCalls++; + return 1; + }, + heatSource: () => { + heatSourceCalls++; + return 0; + }, + }, + }); + model.setMeshConfig({ meshDimension: "1D", elementOrder: "linear", numElementsX: 8, maxX: 1 }); + model.addBoundaryCondition("0", ["constantTemp", 0]); + model.addBoundaryCondition("1", ["constantTemp", 1]); + model.setSolverMethod("lusolve"); + + await model.solveAsync(null); + + assert( + thermalConductivityCalls > 0 && heatSourceCalls > 0, + `solveAsync forwards the coefficients to the assembler: thermalConductivity evaluated ` + + `${thermalConductivityCalls} times, heatSource ${heatSourceCalls} times`, + ); +} + +basicLog(""); +if (failed > 0) { + errorLog(`${passed} passed, ${failed} failed.`); +} else { + basicLog(`${passed} passed, ${failed} failed.`); +} +basicLog("================================"); +if (failed > 0) process.exit(1); diff --git a/tests/regression/HeatConduction2DVaryingCoefficients/REGRESSION.md b/tests/regression/HeatConduction2DVaryingCoefficients/REGRESSION.md new file mode 100644 index 0000000..f39c65b --- /dev/null +++ b/tests/regression/HeatConduction2DVaryingCoefficients/REGRESSION.md @@ -0,0 +1,52 @@ +# Regression Test — HeatConduction2DVaryingCoefficients + +## Purpose + +This test guards the spatially varying `thermalConductivity` and `heatSource` coefficients of +`heatConductionScript` in 2D, where the coefficients are evaluated at the physical coordinates +produced by the 2D isoparametric mapping. + +The 2D assembly path is a separate implementation from the 1D one, with its own Gauss loop and +its own mapping, so the 1D test does not cover it. As there, the expected values are closed-form +solutions rather than stored reference numbers, with the setups chosen so the finite element +solution is exact at the nodes and the tolerance can be `1e-10`. + +## Problem setup + +Common to both cases: domain `x ∈ [0, 1]`, `y ∈ [0, 1]`, 4 × 3 quadratic elements, `lusolve`. +Boundaries left unspecified are natural (zero flux), which both exact solutions satisfy. + +| Case | k(x, y) | Q | Boundaries | Exact solution | +| ---- | ------- | --- | -------------------------------------- | -------------- | +| 1 | 1 + x | −1 | left (1) T = 0, right (3) T = 1 | T = x | +| 2 | 1 + y | −1 | bottom (0) T = 0, top (2) T = 1 | T = y | + +Case 2 is case 1 rotated onto the other axis. It is what confirms the y-coordinate reaches the +coefficients rather than being dropped or swapped with x — a mutation swapping the two arguments +is invisible to case 1 alone and to the whole of the 1D test, where the coefficient is called +with x only. + +## Expected values + +Every nodal temperature must match the closed-form solution to within `1e-10`. Observed largest +deviations are of order `1e-15`; swapping x and y in the 2D assembler moves them to `1e-1`. + +## How to run + +From the repository root: + +```bash +node tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js +``` + +A passing run prints two `PASS:` lines and `2 passed, 0 failed.`; a failing run prints `FAIL:` +with the largest deviation and the node it occurred at, and exits with code 1. + +## After modifying the code + +| Situation | Action | +| --------------------------------------------------- | ----------------------------------------------------------------------------------------------- | +| Bug fix that should not change results | Run the test — it must still pass. | +| Change to the coefficient API | Update the cases; the analytical solutions themselves stay valid. | +| Intentional change to quadrature or element mapping | The expected values do not move. If a case now fails, the change altered the physics, not the reference. | +| Adding the frontal solver to the 2D coverage | Add a case comparing it against `lusolve`, as case 3 of the 1D test does. | diff --git a/tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js b/tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js new file mode 100644 index 0000000..a9e2a23 --- /dev/null +++ b/tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js @@ -0,0 +1,155 @@ +/** + * ════════════════════════════════════════════════════════════════ + * FEAScript Core Library + * Lightweight Finite Element Simulation in JavaScript + * Version: 0.3.0 (RC) | https://feascript.com + * MIT License © 2023–2026 FEAScript + * ════════════════════════════════════════════════════════════════ + */ + +/** + * Regression test for HeatConduction2DVaryingCoefficients + * + * Guards the spatially varying `thermalConductivity` and `heatSource` coefficients of + * heatConductionScript in 2D, where the coefficients are evaluated at the physical + * coordinates produced by the 2D isoparametric mapping + * + * As in the 1D counterpart the expected values are closed-form solutions rather than stored + * reference numbers, with the cases chosen so the finite element solution is exact at the + * nodes and the tolerance can be 1e-10. + * + * Run: node tests/regression/HeatConduction2DVaryingCoefficients/regression.test.js (or npm test) + */ + +import * as mathjs from "mathjs"; +import { FEAScriptModel } from "../../../src/FEAScript.js"; +import { basicLog, errorLog } from "../../../src/utilities/logging.js"; + +// FEAScript.js references `math` as a global (loaded via CDN in browser). +// Set it here before any solve() call. +globalThis.math = mathjs; + +const TOLERANCE = 1e-10; + +function runSimulation(coefficientFunctions, boundaryConditions) { + const model = new FEAScriptModel(); + + model.setModelConfig("heatConductionScript", { coefficientFunctions }); + model.setMeshConfig({ + meshDimension: "2D", + elementOrder: "quadratic", + numElementsX: 4, + numElementsY: 3, + maxX: 1, + maxY: 1, + }); + + Object.entries(boundaryConditions).forEach(([boundaryKey, condition]) => { + model.addBoundaryCondition(boundaryKey, condition); + }); + model.setSolverMethod("lusolve"); + + const { solutionVector, nodesCoordinates } = model.solve(); + + // solutionVector from math.lusolve is a nested array: [[T0], [T1], ...] + return { + temperatures: solutionVector.map((value) => (Array.isArray(value) ? value[0] : value)), + nodesXCoordinates: nodesCoordinates.nodesXCoordinates, + nodesYCoordinates: nodesCoordinates.nodesYCoordinates, + }; +} + +let passed = 0; +let failed = 0; + +function assert(condition, message) { + if (!condition) { + errorLog(`FAIL: ${message}`); + failed++; + } else { + basicLog(`PASS: ${message}`); + passed++; + } +} + +/** + * Function to assert that every nodal temperature matches an analytical solution + * @param {string} label - Description of the case under test + * @param {object} result - Object containing the computed temperatures and node coordinates + * @param {function} analyticalSolution - Function returning the exact temperature at (x, y) + */ +function assertMatchesAnalyticalSolution(label, result, analyticalSolution) { + const { temperatures, nodesXCoordinates, nodesYCoordinates } = result; + + let maxError = 0; + let maxErrorNodeIndex = 0; + for (let nodeIndex = 0; nodeIndex < temperatures.length; nodeIndex++) { + const error = Math.abs( + temperatures[nodeIndex] - + analyticalSolution(nodesXCoordinates[nodeIndex], nodesYCoordinates[nodeIndex]), + ); + if (error > maxError) { + maxError = error; + maxErrorNodeIndex = nodeIndex; + } + } + + assert( + maxError < TOLERANCE, + `${label}: largest nodal deviation ${maxError.toExponential(3)} at ` + + `(x = ${nodesXCoordinates[maxErrorNodeIndex]}, y = ${nodesYCoordinates[maxErrorNodeIndex]}) ` + + `(tolerance ${TOLERANCE})`, + ); +} + +basicLog(""); +basicLog("================================"); +basicLog("Starting regression test for solid heat transfer in 2D with varying coefficients..."); + +/** + * Case 1 - conductivity varying along x, with a matching heat source + * + * With k(x, y) = 1 + x and Q = -1 the field T = x satisfies div(k * grad(T)) + Q = 0. The + * left and right boundaries are held at 0 and 1, while the bottom and top are left + * unspecified and so are natural (zero flux), which T = x also satisfies as it has no + * y-dependence. + */ +assertMatchesAnalyticalSolution( + "Manufactured solution T = x", + runSimulation( + { thermalConductivity: (x, y) => 1 + x, heatSource: -1 }, + { + 1: ["constantTemp", 0], // Left boundary (x = 0) + 3: ["constantTemp", 1], // Right boundary (x = 1) + }, + ), + (x, y) => x, +); + +/** + * Case 2 - the same problem rotated onto the y-axis + * + * With k(x, y) = 1 + y and Q = -1 the exact solution is T = y, held by the bottom and top + * boundaries while the sides stay natural. Rotating the previous case confirms that the + * y-coordinate reaches the coefficients rather than being dropped or swapped with x. + */ +assertMatchesAnalyticalSolution( + "Manufactured solution T = y", + runSimulation( + { thermalConductivity: (x, y) => 1 + y, heatSource: -1 }, + { + 0: ["constantTemp", 0], // Bottom boundary (y = 0) + 2: ["constantTemp", 1], // Top boundary (y = 1) + }, + ), + (x, y) => y, +); + +basicLog(""); +if (failed > 0) { + errorLog(`${passed} passed, ${failed} failed.`); +} else { + basicLog(`${passed} passed, ${failed} failed.`); +} +basicLog("================================"); +if (failed > 0) process.exit(1);