FluidNumerics
diff --git a/‎docs/Models/linear-euler-2d-model.md‎
Lines changed: 21 additions & 16 deletions b/‎docs/Models/linear-euler-2d-model.md‎
Lines changed: 21 additions & 16 deletions
diff --git a/‎docs/Models/linear-euler-2d-pml-model.md‎
Lines changed: 2 additions & 1 deletion b/‎docs/Models/linear-euler-2d-pml-model.md‎
Lines changed: 2 additions & 1 deletion
@@ -18,11 +18,12 @@ $$
     u \\ 
     v \\ 
     p \\
-    c
+    c \\
+    \rho_0
     \end{pmatrix}
 $$
 
-where $\rho$ is a density anomaly, referenced to the density $\rho_0$, $u$ and $v$ are the $x$ and $y$ components of the fluid velocity (respectively), $p$ is the pressure, and $c$ is the speed of sound. The sound speed is carried as a per-node solution variable so that it can vary in space (e.g. across material interfaces); its flux and source are identically zero, so $c$ is held fixed in time at each spatial location.
+where $\rho$ is a density anomaly, referenced to the background density $\rho_0$, $u$ and $v$ are the $x$ and $y$ components of the fluid velocity (respectively), $p$ is the pressure, $c$ is the speed of sound, and $\rho_0$ is the background density. Both the sound speed $c$ and the background density $\rho_0$ are carried as per-node solution variables so that they can vary in space (e.g. across material interfaces); their flux and source are identically zero, so $c$ and $\rho_0$ are held fixed in time at each spatial location. This is entropy-stable for piecewise-constant material regions aligned with element boundaries: element interiors have $\nabla \rho_0 = \nabla c = 0$ (so the flux-divergence form is exact) and the impedance-matched Riemann flux handles the jumps at faces.
 
 When we assume an ideal gas, and a motionless background state, the conservative fluxes are
 
@@ -33,6 +34,7 @@ $$
     \frac{p}{\rho_0} \hat{x} \\
     \frac{p}{\rho_0} \hat{y} \\
     \rho_0 c^2 (u \hat{x} + v \hat{y}) \\
+    \vec{0} \\
     \vec{0}
     \end{pmatrix}
 $$
@@ -50,7 +52,7 @@ $$
 $$
 
 ## Implementation
-The Linear Euler 2D model is implemented as a type extension of the [`DGModel2D` class](../ford/type/dgmodel2d_t.html). The [`LinearEuler2D_t` class](../ford/type/lineareuler2d_t.html) adds a parameter for the reference density $\rho_0$ and overrides `SetNumberOfVariables` (to declare `nvar = 5`), `SetMetadata`, `AdditionalInit`, `entropy_func`, `flux2d`, and `riemannflux2d`. The sound speed is no longer a scalar attribute of the model — it lives in `solution(:,:,:,5)` and can be set independently per node when initializing the simulation.
+The Linear Euler 2D model is implemented as a type extension of the [`DGModel2D` class](../ford/type/dgmodel2d_t.html). The [`LinearEuler2D_t` class](../ford/type/lineareuler2d_t.html) keeps a scalar `rho0` attribute (used as the reference value that fills variable 6 in the built-in initial conditions) and overrides `SetNumberOfVariables` (to declare `nvar = 6` with `nstepped = 4`, so the last two variables are static), `SetMetadata`, `AdditionalInit`, `entropy_func`, `flux2d`, and `riemannflux2d`. The sound speed lives in `solution(:,:,:,5)` and the background density in `solution(:,:,:,6)`; both can be set independently per node when initializing the simulation.
 
 ### Riemann Solver
 The `LinearEuler2D` class is defined using the conservative form of the conservation law. The interface flux is the exact upwind (Godunov) solver for the linearized acoustic system obtained by characteristic decomposition. The normal-flux Jacobian has eigenstructure
@@ -67,29 +69,31 @@ $$
     p^*   = \frac{Z_R p_L + Z_L p_R + Z_L Z_R (u_{n,L} - u_{n,R})}{Z_L + Z_R},
 $$
 
-with acoustic impedance $Z = \rho_0 c$. This reduces correctly to Fresnel reflection / transmission across an impedance jump (e.g. $c_L \neq c_R$ at a material interface). The interface flux is then
+with the per-side acoustic impedance $Z = \rho_0 c$, where each side uses its own background density $\rho_0$ and sound speed $c$. This reduces correctly to Fresnel reflection / transmission across an impedance jump (e.g. $Z_L \neq Z_R$ at a material interface, whether from a density jump, a sound-speed jump, or both). The interface flux is then
 
 $$
     \overleftrightarrow{f}^* \cdot \hat{n} =
     \begin{pmatrix}
-    \rho_0\,u_n^* \\
-    p^*\, n_x / \rho_0 \\
-    p^*\, n_y / \rho_0 \\
-    \rho_0\,\overline{c^2}\,u_n^* \\
+    \overline{\rho_0}\,u_n^* \\
+    p^*\, n_x / \overline{\rho_0} \\
+    p^*\, n_y / \overline{\rho_0} \\
+    \overline{\rho_0}\,\overline{c^2}\,u_n^* \\
+    0 \\
     0
     \end{pmatrix}, \qquad
+    \overline{\rho_0} = \tfrac{1}{2}(\rho_{0,L} + \rho_{0,R}), \quad
     \overline{c^2} = \tfrac{1}{2}(c_L^2 + c_R^2).
 $$
 
-The pressure-flux uses an averaged $c^2$ as a pragmatic treatment of the non-conservative product $\rho_0 c^2 \nabla \cdot \vec{v}$ at a face where $c$ jumps. The previously used local Lax-Friedrichs solver was found to over-dissipate the tangential and entropy modes and to fail to stably handle impedance mismatch at high polynomial order (aliasing instability at material interfaces), and has been replaced by the characteristic flux above. Details are in [`self_lineareuler2d_t.f90`](../ford/sourcefile/self_lineareuler2d_t.f90.html).
+The reconstructed density/momentum/pressure fluxes use the face-averaged $\overline{\rho_0}$ and $\overline{c^2}$ as a pragmatic treatment of the non-conservative products $p/\rho_0$ and $\rho_0 c^2 \nabla \cdot \vec{v}$ at a face where the coefficients jump; the physically important reflection/transmission is carried exactly by the per-side impedances above. The previously used local Lax-Friedrichs solver was found to over-dissipate the tangential and entropy modes and to fail to stably handle impedance mismatch at high polynomial order (aliasing instability at material interfaces), and has been replaced by the characteristic flux above. Details are in [`self_lineareuler2d_t.f90`](../ford/sourcefile/self_lineareuler2d_t.f90.html).
 
 ### Boundary conditions
-Boundary conditions are managed through the [extensible boundary-condition system](../Learning/BoundaryConditions.md). Each model registers its hyperbolic boundary conditions inside `AdditionalInit` by calling `hyperbolicBCs%RegisterBoundaryCondition(id, name, fn)`. `LinearEuler2D_t` registers a no-normal-flow handler out of the box; the GPU build additionally registers a radiation handler that runs on the device. Prescribed boundary conditions are registered by user code (or by an example subclass) so that the external state can be set as a function of position and time.
+Boundary conditions are managed through the [extensible boundary-condition system](../Learning/BoundaryConditions.md). Each model registers its hyperbolic boundary conditions inside `AdditionalInit` by calling `hyperbolicBCs%RegisterBoundaryCondition(id, name, fn)`. `LinearEuler2D_t` registers both a no-normal-flow and a radiation handler out of the box (CPU); the GPU build overwrites both registrations with equivalent device kernels. Prescribed boundary conditions are registered by user code (or by an example subclass) so that the external state can be set as a function of position and time.
 
 The built-in boundary identifiers used with the mesh generators are
 
 * `SELF_BC_RADIATION` — set the external state on the boundary to zero in the Riemann solver (open/non-reflecting).
-* `SELF_BC_NONORMALFLOW` — reflect the velocity vector about the boundary normal and prolong $\rho$, $p$, and $c$. This produces a reflecting (free-slip) wall and works for arbitrarily oriented normals.
+* `SELF_BC_NONORMALFLOW` — reflect the velocity vector about the boundary normal and prolong $\rho$, $p$, $c$, and $\rho_0$. This produces a reflecting (free-slip) wall and works for arbitrarily oriented normals.
 * `SELF_BC_PRESCRIBED` — use a user-registered handler to fill the external state.
 
 As an example, when using the built-in structured mesh generator,
@@ -115,17 +119,18 @@ integer :: bcids(1:4)
     See the [Structured Mesh documentation](../MeshGeneration/StructuredMesh.md) for details on using the `structuredmesh` procedure, and the [Boundary Condition System](../Learning/BoundaryConditions.md) for how to register new BC handlers.
 
 !!! note
-    To set a prescribed state as a function of position and time, create a type-extension of `LinearEuler2D` and register a custom BC method against `SELF_BC_PRESCRIBED` from `AdditionalInit`. Remember that your handler must also fill `solution%extBoundary(:,:,:,5)` with the appropriate sound-speed value at the boundary — the planewave examples show this pattern.
+    To set a prescribed state as a function of position and time, create a type-extension of `LinearEuler2D` and register a custom BC method against `SELF_BC_PRESCRIBED` from `AdditionalInit`. Remember that your handler must also fill `solution%extBoundary(:,:,:,5)` (sound speed) and `solution%extBoundary(:,:,:,6)` (background density) with the appropriate values at the boundary — the planewave examples show this pattern.
 
-### Setting the sound speed
+### Setting the sound speed and background density
 
-Because $c$ is a solution variable, you initialize it the same way you initialize $\rho$, $u$, $v$, and $p$:
+Because $c$ and $\rho_0$ are solution variables, you initialize them the same way you initialize $\rho$, $u$, $v$, and $p$:
 
 ```fortran
 this%solution%interior(i,j,iel,5) = c_value_at_this_node
+this%solution%interior(i,j,iel,6) = rho0_value_at_this_node
 ```
 
-For a uniform background, set every node to the same constant. For a piecewise-constant medium (e.g. bone embedded in marrow), assign the local material's $c$. The `SphericalSoundWave` initializer takes the (uniform) sound speed as an explicit argument:
+For a uniform background, set every node to the same constants. For a piecewise-constant medium (e.g. bone embedded in marrow), assign the local material's $c$ and $\rho_0$. The `SphericalSoundWave` initializer takes the (uniform) sound speed as an explicit argument and fills the background density from the scalar `this%rho0`:
 
 ```fortran
 call model%SphericalSoundWave(rhoprime=1.0e-2_prec, Lr=0.1_prec, &
@@ -157,4 +162,4 @@ For examples, see any of the following
 * [`examples/lineareuler2d_spherical_soundwave_closeddomain.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/linear_euler2d_spherical_soundwave_closeddomain.f90) - Simulation with a gaussian pressure and density anomaly as an initial condition in a domain with no-normal-flow boundary conditions on all sides. Demonstrates uniform sound speed via the `SphericalSoundWave` initializer.
 * [`examples/linear_euler2d_planewave_propagation.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/linear_euler2d_planewave_propagation.f90) - Gaussian plane wave that propagates at a $45^\circ$ angle through a square domain. The initial and boundary conditions are an exact plane-wave solution to the linear Euler equations. The example subclass carries its own `c` attribute and writes it into `solution(...,5)` for both the initial condition and the prescribed boundary state.
 * [`examples/linear_euler2d_planewave_reflection.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/linear_euler2d_planewave_reflection.f90) - Gaussian plane wave reflected off a wall at $x=1$ via the method of images. Combines prescribed boundary conditions with no-normal-flow on the reflecting side.
-* [`examples/linear_euler2d_boneandmarrow.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/linear_euler2d_boneandmarrow.f90) - Heterogeneous-medium test on a HOHQMesh ISM-MM mesh tagged with three materials (muscle/bone/marrow). Each material is mapped to a representative sound speed written into `solution(...,5)`, and a Gaussian acoustic pulse refracts and reflects at the material interfaces. Exercises the impedance-matched Riemann solver across $c$ discontinuities.
+* [`examples/linear_euler2d_boneandmarrow.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/linear_euler2d_boneandmarrow.f90) - Heterogeneous-medium test on a HOHQMesh ISM-MM mesh tagged with three materials (muscle/bone/marrow). Each material is mapped to a representative sound speed and background density, written into `solution(...,5)` and `solution(...,6)`, and a Gaussian acoustic pulse refracts and reflects at the material interfaces. Exercises the impedance-matched Riemann solver across $\rho_0 c$ (impedance) discontinuities.
@@ -76,7 +76,8 @@ Both paths are covered step-by-step in the [PML tutorial](../Tutorials/LinearEul
 * `SetMetadata` — registers names and units for the four PML auxiliaries (plus `sigma_x`, `sigma_y`).
 * `AdditionalInit` — allocates `sigma_x` and `sigma_y` and registers PML-aware no-normal-flow and radiation boundary handlers.
 * `AdditionalFree` — releases the damping fields.
-* `flux2d`, `riemannflux2d` — reuse the parent linear-Euler flux for variables 1–5 and return zero for variables 6–9 (auxiliaries carry no flux).
+* `flux2d`, `riemannflux2d` — use the parent linear-Euler acoustic flux form for variables 1–5 and return zero for variables 6–9 (auxiliaries carry no flux). Note the PML variant uses the **scalar** background density `rho0`: unlike the parent model (which carries `rho0` as a per-node solution variable in slot 6), the PML model repurposes slot 6 for the auxiliary `phi_rho`.
+* `entropy_func` — acoustic energy using the scalar `rho0` and the per-node sound speed `s(5)`. This override is required because the parent's `entropy_func` reads `rho0` from `s(6)`, which the PML model uses for `phi_rho`.
 * `sourcemethod` — implements the ADE source term node-by-node. We override `sourcemethod` rather than `source2d` because the pure `source2d(s, dsdx)` signature has no access to the per-node $\sigma_x(i,j,iel)$, $\sigma_y(i,j,iel)$ lookups.
 
 A new procedure `SetPMLProfile(x_interior_min, x_interior_max, y_interior_min, y_interior_max, pml_width, sigma_max, ramp_exponent)` populates `sigma_x` and `sigma_y` from the per-element material tags and the geometric ramp.