Forward coefficient functions to solveAsync, add varying-coefficient regression tests

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(

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.             |

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);

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.                        |