Allow setting the viscosity in the creeping flow model (#102)

examples/creepingFlowScript/lidDrivenCavity2DCreepingFlow/lidDrivenCavity2DCreepingFlow.js

@@ -14,7 +14,11 @@ console.log("FEAScript Version:", printVersion);
 
 // 1. Create the model and select the physics
 const model = new FEAScriptModel();
-model.setModelConfig("creepingFlowScript");
+model.setModelConfig("creepingFlowScript", {
+  coefficientFunctions: {
+    viscosity: 1, // μ
+  },
+});
 
 // 2. Define the mesh
 model.setMeshConfig({

src/FEAScript.js

@@ -364,7 +364,11 @@ export class FEAScriptModel {
       }
     } else if (this.solverConfig === "creepingFlowScript") {
       // Use regular linear solver methods for steady creeping flow
-      const creepingFlowResult = assembleCreepingFlowMatrix(meshData, this.boundaryConditions);
+      const creepingFlowResult = assembleCreepingFlowMatrix(
+        meshData,
+        this.boundaryConditions,
+        this.coefficientFunctions,
+      );
       jacobianMatrix = creepingFlowResult.jacobianMatrix;
       residualVector = creepingFlowResult.residualVector;
 

src/physics/creepingFlow.js

@@ -25,6 +25,8 @@ import { basicLog, debugLog, errorLog } from "../utilities/logging.js";
  *
  * @param {object} meshData - Object containing prepared mesh data (must use quadratic elements)
  * @param {object} boundaryConditions - Object containing boundary conditions for the finite element analysis
+ * @param {object} [coefficientFunctions] - Coefficients of the creeping flow model:
+ *  - viscosity: Dynamic viscosity μ, a number or a function of (x, y) (optional, defaults to 1)
  * @returns {object} An object containing:
  *  - jacobianMatrix: The assembled (sparse) Jacobian matrix
  *  - residualVector: The assembled residual vector
@@ -40,7 +42,7 @@ import { basicLog, debugLog, errorLog } from "../utilities/logging.js";
  * classic stiffness/conductivity matrix and `residualVector`
  * corresponds to the traditional load (RHS) vector.
  */
-export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
+export function assembleCreepingFlowMatrix(meshData, boundaryConditions, coefficientFunctions) {
   basicLog("Starting creeping flow matrix assembly...");
 
   // Extract mesh data
@@ -144,8 +146,12 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
     pressureBasisFunctions.getBasisFunctions(ksi, eta),
   );
 
-  // Viscosity coefficient
-  const mu = 1.0;
+  // Dynamic viscosity, constant or a function of position
+  const viscosity = coefficientFunctions?.viscosity ?? 1;
+  const viscosityAt = typeof viscosity === "function" ? viscosity : () => viscosity;
+  if (typeof viscosity !== "function" && !(viscosity > 0)) {
+    errorLog(`Creeping flow viscosity must be a positive number or a function of (x, y); got ${viscosity}.`);
+  }
 
   // Local DOF ordering in the element matrix: [u_0 … u_{n_2−1}, v_0 … v_{n_2−1}, p_0 … p_{n_1−1}]
   const localDOFs = 2 * nodesPerVelocityElement + nodesPerPressureElement;
@@ -188,7 +194,7 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
       const pressureBasisFunction = pressureReferenceBasisFunctions[pointIndex].basisFunction;
 
       // Perform isoparametric mapping using Q2/P2 velocity basis functions
-      const { detJacobian } = performIsoparametricMapping2D({
+      const { xCoordinates, yCoordinates, detJacobian } = performIsoparametricMapping2D({
         basisFunction,
         basisFunctionDerivKsi,
         basisFunctionDerivEta,
@@ -202,6 +208,7 @@ export function assembleCreepingFlowMatrix(meshData, boundaryConditions) {
 
       // Gauss integration weight factors
       const weightFactor = integrationWeight * detJacobian;
+      const mu = viscosityAt(xCoordinates, yCoordinates); // Viscosity at this integration point
       const viscousFactor = -weightFactor * mu;
 
       // Assemble viscous stiffness terms (K block)

tests/regression/lidDrivenCavity2DCreepingFlow/README.md

@@ -47,6 +47,10 @@ The minimum x-velocity is the strongest return flow of the recirculation below t
 pressure nodes are the element corners (13 × 7 for the structured mesh) or the triangle vertices.
 Tolerance used in the assertions: `1e-4`.
 
+The test also solves each case with a viscosity of 2.5, given as a number and as a function of
+(x, y). The Stokes equations with velocity boundary conditions are linear in μ, so the velocity
+must be unchanged and the pressure scaled by 2.5 (relative tolerance `1e-10`).
+
 ## How to run
 
 From the repository root:

tests/regression/lidDrivenCavity2DCreepingFlow/regression.test.js

@@ -13,7 +13,9 @@
  * Replicates the physical setup from lidDrivenCavity2DCreepingFlow.js with Q2-Q1 quadrilaterals
  * (structured and Gmsh meshes) and P2-P1 triangles (Gmsh mesh), and asserts the strongest return
  * flow (minimum x-velocity), the nodal mean x-velocity and y-velocity and the number of pressure
- * nodes against known good values
+ * nodes against known good values. It also checks that a viscosity μ (constant or a function of
+ * position) leaves the velocity unchanged and scales the pressure by μ, as the Stokes equations
+ * with velocity boundary conditions require
  *
  * Run: node tests/regression/lidDrivenCavity2DCreepingFlow/regression.test.js (or npm test)
  */
@@ -74,10 +76,10 @@ async function getMeshConfig({ mesh, meshFile }) {
   return { parsedMesh, meshDimension: "2D", elementOrder: "quadratic" };
 }
 
-function runSimulation(meshConfig) {
+function runSimulation(meshConfig, modelOptions = {}) {
   const model = new FEAScriptModel();
 
-  model.setModelConfig("creepingFlowScript");
+  model.setModelConfig("creepingFlowScript", modelOptions);
   model.setMeshConfig(meshConfig);
   ["0", "1", "2", "3"].forEach((boundaryKey) => {
     model.addBoundaryCondition(
@@ -152,6 +154,24 @@ for (const testCase of cases) {
       pressureNodesInField >= testCase.expectedPressureNodes,
     `${label}: field p has a finite value at every node and holds the solved corner pressures`,
   );
+
+  // A viscosity μ leaves the velocity unchanged and scales the pressure by μ
+  const meshConfig = await getMeshConfig(testCase);
+  const maxPressure = Math.max(...fields.p.map(Math.abs));
+  for (const [viscosityLabel, viscosity] of [
+    ["viscosity 2.5", 2.5],
+    ["viscosity (x, y) => 2.5", (x, y) => 2.5],
+  ]) {
+    const { fields: scaledFields } = runSimulation(meshConfig, { coefficientFunctions: { viscosity } });
+    const velocityUnchanged = scaledFields.u.every((u, i) => Math.abs(u - fields.u[i]) < 1e-10);
+    const pressureScaled = scaledFields.p.every(
+      (p, i) => Math.abs(p - 2.5 * fields.p[i]) < 1e-10 * 2.5 * maxPressure,
+    );
+    assert(
+      velocityUnchanged && pressureScaled,
+      `${label}: ${viscosityLabel} leaves the velocity unchanged and scales the pressure by 2.5`,
+    );
+  }
 }
 
 basicLog("");

tests/unit/model.test.js

@@ -482,6 +482,29 @@ function assert(condition, message) {
         3: ["constantTemperature", 200],
       },
     },
+    {
+      label: "Creeping flow (varying viscosity)",
+      modelConfig: "creepingFlowScript",
+      options: {
+        coefficientFunctions: {
+          viscosity: (x, y) => 1 + x * y,
+        },
+      },
+      meshConfig: {
+        meshDimension: "2D",
+        elementOrder: "quadratic",
+        numElementsX: 4,
+        numElementsY: 2,
+        maxX: 4,
+        maxY: 2,
+      },
+      boundaryConditions: {
+        0: ["constantVelocity", 0, 0],
+        1: ["constantVelocity", 0, 0],
+        2: ["constantVelocity", 1, 0],
+        3: ["constantVelocity", 0, 0],
+      },
+    },
   ];
 
   for (const testCase of testCases) {