Document the creeping flow viscosity

Replace the "μ is fixed to 1" notes with the viscosity coefficient
(number or function of position), in the creeping flow page, the
FEAScriptModel options table and llms.txt.

Author: nikoscham

docs/api/feascript-model.md

@@ -44,12 +44,12 @@ Selects the physical model and sets its options.
 | `"creepingFlowScript"`       | [Creeping flow](../models/creeping-flow.md)               | 2D         |
 | `"eulerBernoulliBeamScript"` | [Euler–Bernoulli beam](../models/euler-bernoulli-beam.md) | 1D         |
 
-| Option                 | Type      | Used by                                                 | Description                                            |
-| ---------------------- | --------- | ------------------------------------------------------- | ------------------------------------------------------ |
-| `coefficientFunctions` | `object`  | Heat conduction, general form PDE, Euler–Bernoulli beam | Coefficients of the equation; see each model page      |
-| `nonlinear`            | `boolean` | General form PDE                                        | Solve with the Newton–Raphson method (default `false`) |
-| `maxIterations`        | `number`  | All                                                     | Maximum iterations of the iterative method in use      |
-| `tolerance`            | `number`  | All                                                     | Convergence tolerance of the iterative method in use   |
+| Option                 | Type      | Used by                                                                | Description                                            |
+| ---------------------- | --------- | ---------------------------------------------------------------------- | ------------------------------------------------------ |
+| `coefficientFunctions` | `object`  | Heat conduction, general form PDE, creeping flow, Euler–Bernoulli beam | Coefficients of the equation; see each model page      |
+| `nonlinear`            | `boolean` | General form PDE                                                       | Solve with the Newton–Raphson method (default `false`) |
+| `maxIterations`        | `number`  | All                                                                    | Maximum iterations of the iterative method in use      |
+| `tolerance`            | `number`  | All                                                                    | Convergence tolerance of the iterative method in use   |
 
 `maxIterations` and `tolerance` apply to the Jacobi solvers in linear models and to the Newton–Raphson method in
 nonlinear models; see [Solvers](../fundamentals/solvers.md) for their defaults.

docs/models/creeping-flow.md

@@ -21,9 +21,8 @@ where $\bar{u} = (u, v)$ is the velocity, $p$ the pressure and $\mu$ the dynamic
 forces are neglected, which is valid when the Reynolds number is much smaller than one.
 
 :::note
-The viscosity is fixed to $\mu = 1$ and cannot be changed yet. For a Newtonian fluid with a different viscosity
-and velocity boundary conditions only, the velocity field is the same and the pressure scales linearly with
-$\mu$.
+With velocity boundary conditions only, a constant viscosity does not change the velocity field: the pressure
+scales linearly with $\mu$.
 :::
 
 ## Taylor–Hood elements
@@ -74,7 +73,11 @@ import { FEAScriptModel } from "feascript";
 
 // 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({
@@ -97,7 +100,24 @@ model.setSolverMethod("lusolve");
 const { solutionVector, nodesCoordinates, fields } = model.solve();
 ```
 
-The model takes no coefficient functions.
+### Coefficient functions
+
+| Property    | Symbol | Default | Type                             |
+| ----------- | ------ | ------- | -------------------------------- |
+| `viscosity` | $\mu$  | `1`     | `number` or function of position |
+
+The dynamic viscosity can be a positive constant or a function of the physical coordinates `(x, y) => ...`,
+evaluated at every integration point, so regions of different viscosity are supported. `coefficientFunctions`
+can be omitted entirely to use $\mu = 1$.
+
+```javascript
+model.setModelConfig("creepingFlowScript", {
+  coefficientFunctions: {
+    // A ten times more viscous fluid in the right half of the cavity
+    viscosity: (x, y) => (x < 2 ? 1 : 10),
+  },
+});
+```
 
 ### Boundary conditions
 

static/llms.txt

@@ -61,7 +61,7 @@ const { solutionVector, nodesCoordinates } = model.solve();
 | `"heatConductionScript"`     | −∇·(k ∇T) = Q                                                | 1D or 2D, linear or quadratic |
 | `"generalFormPDEScript"`     | d/dx(A du/dx) + B du/dx + C u = D                            | 1D only                      |
 | `"frontPropagationScript"`   | Eikonal equation \|∇u\| = 1 (vanishing viscosity)            | 1D or 2D                     |
-| `"creepingFlowScript"`       | Stokes flow: −μ∇²u + ∇p = 0, ∇·u = 0, with μ = 1 (fixed)     | 2D, quadratic only           |
+| `"creepingFlowScript"`       | Stokes flow: −μ∇²u + ∇p = 0, ∇·u = 0                         | 2D, quadratic only           |
 | `"eulerBernoulliBeamScript"` | d²/dx²(EI d²w/dx²) + c0 w = q                                | 1D, `elementOrder: "linear"` only |
 
 ### heatConductionScript
@@ -109,7 +109,7 @@ const result = model.solve({ maxIterations: 100, tolerance: 1e-5 });
 
 ### creepingFlowScript (2D, quadratic only)
 
-- No options. The viscosity μ is fixed to 1.
+- Options: `coefficientFunctions: { viscosity }`, the dynamic viscosity μ, a positive number or a function `(x, y) => ...`. Default `viscosity = 1`; `coefficientFunctions` can be omitted. With velocity boundary conditions only, a constant μ leaves the velocity unchanged and scales the pressure by μ.
 - Mesh: `meshDimension: "2D"` and `elementOrder: "quadratic"` (Taylor–Hood Q2–Q1 on quadrilaterals, P2–P1 on triangles from Gmsh).
 - Boundary conditions: `["constantVelocity", u, v]` (Dirichlet), `["stressFree"]` (natural outflow). If no boundary is `"stressFree"`, the pressure is set to 0 at the first pressure node.
 - Solvers: `"lusolve"` only.
@@ -237,8 +237,8 @@ worker.terminate();
 - 3D meshes, and 1D or 3D Gmsh meshes
 - 2D general form PDE
 - Transient (time-dependent) problems
-- Setting the viscosity in creeping flow, or Navier–Stokes (inertia)
-- Coefficient functions for front propagation or creeping flow
+- Navier–Stokes (inertia) in creeping flow
+- Coefficient functions for front propagation
 - Mixed element shapes or orders in one mesh
 - Heat flux (Neumann) conditions with a non-zero value
 - `"frontal"` for creeping flow, the beam, or the nonlinear general form