Refactor solver functions to accept options object for parameters and improve documentation

Author: nikoscham

CONTRIBUTING.md

@@ -32,7 +32,7 @@ Thank you for your interest in contributing! FEAScript is in early development,
 
 All files in the FEAScript-core codebase should follow this structure:
 
-1.  Banner: All files start with the FEAScript ASCII art banner
+1.  Banner: All files start with the FEAScript banner
 2.  Imports:
     - External imports first, alphabetically ordered
     - Internal imports next, grouped by module/folder

src/FEAScript.js

@@ -46,10 +46,17 @@ export class FEAScriptModel {
     this.solverConfig = solverConfig;
 
     // Store coefficient functions if provided
-    if (options && options.coefficientFunctions) {
+    if (options?.coefficientFunctions) {
       this.coefficientFunctions = options.coefficientFunctions;
       debugLog("Coefficient functions set");
     }
+    // Only update if a value is provided, otherwise keep the default
+    if (options?.maxIterations !== undefined) {
+      this.maxIterations = options.maxIterations;
+    }
+    if (options?.tolerance !== undefined) {
+      this.tolerance = options.tolerance;
+    }
 
     debugLog(`Solver config set to: ${solverConfig}`);
   }
@@ -69,67 +76,15 @@ export class FEAScriptModel {
     debugLog(`Solver method set to: ${solverMethod}`);
   }
 
-  async solveWithWebgpu(computeEngine) {
-    if (!this.solverConfig || !this.meshConfig || !this.boundaryConditions) {
-      errorLog("Solver config, mesh config, and boundary conditions must be set before solving.");
-    }
-
-    let jacobianMatrix = [];
-    let residualVector = [];
-    let solutionVector = [];
-    let nodesCoordinates = {};
-
-    // Prepare the mesh
-    basicLog("Preparing mesh...");
-    const meshData = prepareMesh(this.meshConfig);
-    basicLog("Mesh preparation completed");
-
-    // Extract node coordinates from meshData
-    nodesCoordinates = {
-      nodesXCoordinates: meshData.nodesXCoordinates,
-      nodesYCoordinates: meshData.nodesYCoordinates,
-    };
-
-    // Assembly matrices
-    basicLog("Beginning matrix assembly...");
-    console.time("assemblyMatrices");
-    basicLog(`Using solver: ${this.solverConfig}`);
-    if (this.solverConfig === "solidHeatTransferScript") {
-      ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(
-        meshData,
-        this.boundaryConditions
-      ));
-    }
-    console.timeEnd("assemblyMatrices");
-    basicLog("Matrix assembly completed");
-
-    // System solving with WebGPU Jacobi
-    basicLog("Solving system using WebGPU Jacobi...");
-    console.time("systemSolving");
-
-    // Convert matrices to arrays for WebGPU
-    const A = Array.isArray(jacobianMatrix) ? jacobianMatrix : jacobianMatrix.toArray();
-    const b = Array.isArray(residualVector) ? residualVector : residualVector.toArray();
-
-    // For heat conduction FEM, the matrix might be negative definite
-    console.log("Matrix diagonal sample:", A.slice(0, 5).map((row, i) => row[i]));
-    console.log("RHS sample:", b.slice(0, 5));
-
-    // Use WebGPU Jacobi method
-    const initialGuess = new Array(b.length).fill(0);
-    solutionVector = await computeEngine.webgpuJacobiSolver(A, b, initialGuess, 10000, 1e-3);
-
-    console.timeEnd("systemSolving");
-    basicLog("System solved successfully with WebGPU Jacobi");
-
-    return { solutionVector, nodesCoordinates };
-  }
-
-  solve() {
+  /**
+   * Function to solve the finite element problem synchronously
+   * @param {object} [options] - Additional parameters for the solver, such as `maxIterations` and `tolerance`
+   * @returns {object} An object containing the solution vector and the coordinates of the mesh nodes
+   */
+  solve(options = {}) {
     if (!this.solverConfig || !this.meshConfig || !this.boundaryConditions) {
       errorLog("Solver config, mesh config, and boundary conditions must be set before solving.");
     }
-
     /**
      * For consistency across both linear and nonlinear formulations,
      * this project always refers to the assembled right-hand side vector
@@ -172,7 +127,10 @@ export class FEAScriptModel {
       } else {
         // Use regular linear solver methods
         ({ jacobianMatrix, residualVector } = assembleHeatConductionMat(meshData, this.boundaryConditions));
-        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector);
+        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
+          maxIterations: options.maxIterations ?? this.maxIterations,
+          tolerance: options.tolerance ?? this.tolerance,
+        });
         solutionVector = linearSystemResult.solutionVector;
       }
     } else if (this.solverConfig === "frontPropagationScript") {
@@ -187,6 +145,9 @@ export class FEAScriptModel {
         eikonalActivationFlag: eikonalActivationFlag,
         solverMethod: this.solverMethod,
         initialSolution,
+        // TODO: Consider using different maxIterations/tolerance for Newton-Raphson and linear solver
+        maxIterations: options.maxIterations ?? this.maxIterations,
+        tolerance: options.tolerance ?? this.tolerance,
       };
 
       while (eikonalActivationFlag <= 1) {
@@ -199,7 +160,7 @@ export class FEAScriptModel {
         }
 
         // Solve the assembled non-linear system
-        const newtonRaphsonResult = newtonRaphson(assembleFrontPropagationMat, context, 100, 1e-4);
+        const newtonRaphsonResult = newtonRaphson(assembleFrontPropagationMat, context);
 
         // Extract results
         jacobianMatrix = newtonRaphsonResult.jacobianMatrix;
@@ -223,7 +184,10 @@ export class FEAScriptModel {
           this.coefficientFunctions
         ));
 
-        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector);
+        const linearSystemResult = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector, {
+          maxIterations: options.maxIterations ?? this.maxIterations,
+          tolerance: options.tolerance ?? this.tolerance,
+        });
         solutionVector = linearSystemResult.solutionVector;
       }
     }
@@ -233,6 +197,12 @@ export class FEAScriptModel {
     return { solutionVector, nodesCoordinates };
   }
 
+  /**
+   * Function to solve the finite element problem asynchronously
+   * @param {object} computeEngine - The compute engine to use for the asynchronous solver (e.g., a worker or a WebGPU context)
+   * @param {object} [options] - Additional parameters for the solver, such as `maxIterations` and `tolerance`
+   * @returns {Promise<object>} A promise that resolves to an object containing the solution vector and the coordinates of the mesh nodes
+   */
   async solveAsync(computeEngine, options = {}) {
     if (!this.solverConfig || !this.meshConfig || !this.boundaryConditions) {
       errorLog("Solver config, mesh config, and boundary conditions must be set before solving.");
@@ -260,18 +230,18 @@ export class FEAScriptModel {
       if (this.solverMethod === "jacobi-gpu") {
         const { solutionVector: x } = await solveLinearSystemAsync("jacobi-gpu", jacobianMatrix, residualVector, {
           computeEngine,
-          maxIterations: options.maxIterations,
-          tolerance: options.tolerance,
+          maxIterations: options.maxIterations ?? this.maxIterations,
+          tolerance: options.tolerance ?? this.tolerance,
         });
         solutionVector = x;
       } else {
-        const { solutionVector: x } = solveLinearSystem(this.solverMethod, jacobianMatrix, residualVector);
-        solutionVector = x;
+        // Other async solver
       }
     }
     console.timeEnd("totalSolvingTime");
     basicLog("Solving process completed");
 
     return { solutionVector, nodesCoordinates };
   }
+
 }

src/methods/jacobiSolverScript.js

@@ -11,16 +11,16 @@
  * @param {array} A - The system matrix
  * @param {array} b - The right-hand side vector
  * @param {array} x0 - Initial guess for solution vector
- * @param {object} [options] - Additional options for the solver
- * @param {number} [options.maxIterations=1000] - Maximum number of iterations
- * @param {number} [options.tolerance=1e-6] - Convergence tolerance
+ * @param {object} [options] - Optional parameters for the solver, such as `maxIterations` and `tolerance`
  * @returns {object} An object containing:
  *  - solutionVector: The solution vector
  *  - iterations: The number of iterations performed
  *  - converged: Boolean indicating whether the method converged
  */
 export function jacobiSolver(A, b, x0, options = {}) {
-  const { maxIterations = 1000, tolerance = 1e-6 } = options;
+  // Extract options
+  const { maxIterations, tolerance } = options;
+
   const n = A.length;
   let x = [...x0];
   let xNew = new Array(n);

src/methods/linearSystemSolverScript.js

@@ -16,16 +16,16 @@ import * as Comlink from "../vendor/comlink.mjs";
  * @param {string} solverMethod - The solver method to use ("lusolve" or "jacobi")
  * @param {Array} jacobianMatrix - The coefficient matrix
  * @param {Array} residualVector - The right-hand side vector
- * @param {object} [options] - Additional options for the solver
- * @param {number} [options.maxIterations=10000] - Maximum iterations for iterative methods
- * @param {number} [options.tolerance=1e-3] - Convergence tolerance for iterative methods
+ * @param {object} [options] - Optional parameters for the solver, such as `maxIterations` and `tolerance`
  * @returns {object} An object containing:
  *  - solutionVector: The solution vector
  *  - converged: Boolean indicating whether the method converged (for iterative methods)
  *  - iterations: Number of iterations performed (for iterative methods)
  */
 export function solveLinearSystem(solverMethod, jacobianMatrix, residualVector, options = {}) {
-  const { maxIterations = 10000, tolerance = 1e-3 } = options;
+
+  // Extract options
+  const { maxIterations = 10000, tolerance = 1e-4 } = options;
 
   let solutionVector = [];
   let converged = true;
@@ -85,16 +85,16 @@ async function createDefaultComputeEngine() {
  * @param {string} solverMethod - The solver method to use (e.g., "jacobi-gpu")
  * @param {array} jacobianMatrix - The coefficient matrix
  * @param {array} residualVector - The right-hand side vector
- * @param {object} [options] - Additional options for the solver
- * @param {number} [options.maxIterations=10000] - Maximum iterations for iterative methods
- * @param {number} [options.tolerance=1e-3] - Convergence tolerance for iterative methods
+ * @param {object} [options] - Optional parameters for the solver, such as `maxIterations` and `tolerance`
  * @returns {Promise<object>} A promise that resolves to an object containing:
  *  - solutionVector: The solution vector
  *  - converged: Boolean indicating whether the method converged (for iterative methods)
  *  - iterations: Number of iterations performed (for iterative methods)
  */
 export async function solveLinearSystemAsync(solverMethod, jacobianMatrix, residualVector, options = {}) {
-  const { maxIterations = 10000, tolerance = 1e-3 } = options;
+  
+  // Extract options
+  const { maxIterations = 10000, tolerance = 1e-4 } = options;
 
   basicLog(`Solving system using ${solverMethod}...`);
   console.time("systemSolving");
@@ -137,7 +137,7 @@ export async function solveLinearSystemAsync(solverMethod, jacobianMatrix, resid
   basicLog(`System solved successfully (${solverMethod})`);
 
   if (created) {
-    await computeEngine?.destroy?.().catch(() => {});
+    await computeEngine?.destroy?.().catch(() => { });
     created.worker.terminate();
   }
 

src/methods/newtonRaphsonScript.js

@@ -17,15 +17,13 @@ import { assembleFrontPropagationFront } from "../models/frontPropagationScript.
  * Function to solve a system of non-linear equations using the Newton-Raphson method
  * @param {function} assembleMat - Matrix assembler based on the physical model
  * @param {object} context - Context object containing simulation data and options
- * @param {number} [maxIterations=100] - Maximum number of iterations
- * @param {number} [tolerance=1e-4] - Convergence tolerance
  * @returns {object} An object containing:
  *  - solutionVector: The solution vector
  *  - iterations: The number of iterations performed
  *  - converged: Boolean indicating whether the method converged
  */
 
-export function newtonRaphson(assembleMat, context, maxIterations = 100, tolerance = 1e-4) {
+export function newtonRaphson(assembleMat, context = {}) {
   let errorNorm = 0;
   let converged = false;
   let iterations = 0;
@@ -34,6 +32,9 @@ export function newtonRaphson(assembleMat, context, maxIterations = 100, toleran
   let jacobianMatrix = [];
   let residualVector = [];
 
+  // Extract context
+  const { maxIterations = 100, tolerance = 1e-4 } = context;
+
   // Calculate system size
   let totalNodes = context.meshData.nodesXCoordinates.length;
 

src/workers/webgpuWorkerScript.js

@@ -49,9 +49,7 @@ class webgpuFEAScriptWorker {
    * @param {array} A - The system matrix
    * @param {array} b - The right-hand side vector
    * @param {array} x0 - The initial guess
-   * @param {object} [options] - Additional options for the solver
-   * @param {number} [options.maxIterations=1000] - The maximum number of iterations
-   * @param {number} [options.tolerance=1e-6] - The convergence tolerance
+   * @param {object} [options] - Optional parameters for the solver, such as `maxIterations` and `tolerance`
    * @returns {Promise<object>} An object containing the solution vector, iterations, and convergence status
    */
   async webgpuJacobiSolver(A, b, x0, options = {}) {