LinearEuler3D: heterogeneous sound speed + impedance-matched Riemann flux (#138)

Bring the 3-D linear Euler model up to parity with LinearEuler2D:

- Carry the sound speed as solution variable 6 (nvar 5 -> 6), static in
  time (zero flux/source), so heterogeneous media are supported. The
  scalar attribute `c` remains as the reference value used by initial
  conditions (SphericalSoundWave fills variable 6 with it).
- Replace the local Lax-Friedrichs Riemann solver with the
  impedance-matched (characteristic/Godunov) flux used by the 2-D model:
  interface states resolved with the one-sided impedances Z = rho0*c, so
  material interfaces in a heterogeneous c field get the physically
  correct transmission/reflection.
- GPU backend: fluxmethod reads the per-node c (variable 6) and the
  boundaryflux kernel implements the impedance-matched flux; both drop
  the scalar `c` argument (signatures now mirror the 2-D kernels). The
  radiation BC kernel now takes the interior boundary state and copies c
  through to extBoundary so face Riemann solves see a consistent sound
  speed, with the acoustic perturbation still zeroed.
- Register a CPU radiation BC in AdditionalInit_LinearEuler3D_t
  (previously the radiation BC existed only in the GPU backend); the GPU
  AdditionalInit calls the parent and overwrites it with the kernel
  wrapper, mirroring the 2-D module structure.
- Fix the entropy function: s(3)*(3) typo -> s(3)*s(3), include the w
  contribution correctly, and use the per-node sound speed.
- Update the model documentation accordingly.

Motivation: 3-D ultrasound computed tomography forward modeling
(SELF-UCT AcousticTomography3D) needs spatially varying sound speed with
impedance-matched interface handling, exactly as the 2-D class has.

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: fluidnumericsJoe

docs/Models/linear-euler-3d-model.md

@@ -20,11 +20,12 @@ $$
     u \\ 
     v \\ 
     w \\
-    p
+    p \\
+    c
     \end{pmatrix}
 $$
 
-where $\rho$ is a density anomaly, referenced to the density $\rho_0$, $u$, $v$, and $w$ are the $x$, $y$, and $z$ components of the fluid velocity (respectively), and $p$ is the pressure. When we assume an ideal gas, and a motionless background state, the conservative fluxes are
+where $\rho$ is a density anomaly, referenced to the density $\rho_0$, $u$, $v$, and $w$ are the $x$, $y$, and $z$ components of the fluid velocity (respectively), and $p$ is the pressure. The sound speed $c$ is carried as a solution variable so that heterogeneous (spatially varying) media are supported; it has zero flux and zero source, so it is static in time and is set by the initial condition (mirroring the 2-D model). When we assume an ideal gas, and a motionless background state, the conservative fluxes are
 
 $$
     \overleftrightarrow{f} = 
@@ -33,11 +34,12 @@ $$
     p \hat{x} \\
     p \hat{y} \\
     p \hat{z} \\
-    \rho_0c^2(u \hat{x} + v \hat{y} + w \hat{z})
+    \rho_0c^2(u \hat{x} + v \hat{y} + w \hat{z}) \\
+    0
     \end{pmatrix}
 $$
 
-where $c$ is the (constant) speed of sound. The source term is set to zero.
+The source term is set to zero.
 
 $$
     \vec{q} = \vec{0}
@@ -46,44 +48,35 @@ $$
 To track stability of the Euler equation, the total entropy function is
 
 $$
-    e = \frac{1}{2} \int_V u^2 + v^2 + \frac{p}{\rho_0 c^2} \hspace{1mm} dV
+    e = \frac{1}{2} \int_V \rho_0 (u^2 + v^2 + w^2) + \frac{p^2}{\rho_0 c^2} \hspace{1mm} dV
 $$
 
 ## Implementation
-The Linear Euler 3D model is implemented as a type extension of the [`DGModel3D` class](../ford/type/dgmodel3d_t.html). The [`LinearEuler3D_t` class](../ford/type/lineareuler3d_t.html) adds parameters for the reference density and the speed speed of sound and overrides the `SetMetadata`, `entropy_func`, `flux3d`, and `riemannflux3d` type-bound procedures.
+The Linear Euler 3D model is implemented as a type extension of the [`DGModel3D` class](../ford/type/dgmodel3d_t.html). The [`LinearEuler3D_t` class](../ford/type/lineareuler3d_t.html) adds parameters for the reference density and the reference speed of sound and overrides the `SetMetadata`, `entropy_func`, `flux3d`, and `riemannflux3d` type-bound procedures.
 
 ### Riemann Solver
-The `LinearEuler3D` class is defined using the conservative form of the conservation law. The Riemman solver for the hyperbolic part of Euler equation is the local Lax Friedrichs upwind riemann solver
+The `LinearEuler3D` class is defined using the conservative form of the conservation law. The Riemann solver for the hyperbolic part of the Euler equation is the impedance-matched (characteristic/Godunov) upwind flux, identical in form to the 2-D model's. With the acoustic impedances $Z_L = \rho_0 c_L$ and $Z_R = \rho_0 c_R$ evaluated from the per-node sound speed on either side of the face, the interface normal velocity and pressure are
 
 $$
-    \overleftrightarrow{f}_h^* \cdot \hat{n} = \frac{1}{2}( \overleftrightarrow{f}_L \cdot \hat{n}  + \overleftrightarrow{f}_R \cdot \hat{n}  + c (\vec{s}_L - \vec{s}_R))
+    u_n^* = \frac{Z_L u_{n,L} + Z_R u_{n,R} + (p_L - p_R)}{Z_L + Z_R}, \qquad
+    p^* = \frac{Z_R p_L + Z_L p_R + Z_L Z_R (u_{n,L} - u_{n,R})}{Z_L + Z_R}
 $$
 
-where 
+and the normal flux is
 
 $$
-    \overleftrightarrow{f}_L \cdot \hat{n} = 
+    \overleftrightarrow{f}_h^* \cdot \hat{n} = 
     \begin{pmatrix}
-    \rho_0(u_L n_x + v_L n_y + w_L n_z) \\
-    p_L n_x \\
-    p_L n_y \\
-    p_L n_z \\
-    \rho_0c^2(u_L n_x + v_L n_y + w_L n_z)
-    \end{pmatrix}
-$$
-
-$$
-    \overleftrightarrow{f}_R \cdot \hat{n} = 
-    \begin{pmatrix}
-    \rho_0(u_R n_x + v_R n_y + w_R n_z) \\
-    p_R n_x \\
-    p_R n_y \\
-    p_R n_z \\
-    \rho_0c^2(u_R n_x + v_R n_y + w_R n_z)
-    \end{pmatrix}
+    \rho_0 u_n^* \\
+    p^* n_x / \rho_0 \\
+    p^* n_y / \rho_0 \\
+    p^* n_z / \rho_0 \\
+    \rho_0 \overline{c^2} u_n^* \\
+    0
+    \end{pmatrix}, \qquad \overline{c^2} = \frac{1}{2}(c_L^2 + c_R^2)
 $$
 
-The details for this implementation can be found in [`self_lineareuler3d_t.f90`](../ford/sourcefile/self_lineareuler3d_t.f90.html)
+Because the interface states are resolved with the one-sided impedances, material interfaces in a heterogeneous sound-speed field produce the physically correct transmission and reflection. The details for this implementation can be found in [`self_lineareuler3d_t.f90`](../ford/sourcefile/self_lineareuler3d_t.f90.html)
 
 ### Boundary conditions
 When initializing the mesh for your Euler 3D equation solver, you can change the boundary conditions to 

src/SELF_LinearEuler3D_t.f90

@@ -27,17 +27,20 @@
 module self_LinearEuler3D_t
 !! This module defines a class that can be used to solve the Linear Euler
 !! equations in 3-D. The Linear Euler Equations, here, are the Euler equations
-!! linearized about a motionless background state.
+!! linearized about a motionless background state. The sound speed is carried
+!! as a solution variable (static in time, possibly spatially varying) so that
+!! heterogeneous media are supported, mirroring the 2-D model.
+!!
+!! The solution variables are
 !!
-!! The conserved variables are
-
 !! \begin{equation}
 !! \vec{s} = \begin{pmatrix}
 !!     \rho \\
 !!      u \\
 !!      v \\
 !!      w \\
-!!      p
+!!      p \\
+!!      c
 !!  \end{pmatrix}
 !! \end{equation}
 !!
@@ -49,31 +52,36 @@ module self_LinearEuler3D_t
 !!      \frac{p}{\rho_0} \hat{x} \\
 !!      \frac{p}{\rho_0} \hat{y} \\
 !!      \frac{p}{\rho_0} \hat{z} \\
-!!      c^2 \rho_0 ( u \hat{x} + v \hat{y} + w \hat{z})
+!!      c^2 \rho_0 ( u \hat{x} + v \hat{y} + w \hat{z}) \\
+!!      0
 !!  \end{pmatrix}
 !! \end{equation}
 !!
-!! and the source terms are null.
+!! and the source terms are null. The sound speed variable has zero flux and
+!! zero source; it is set by the initial condition and preserved in time.
+!! The Riemann flux is the impedance-matched (characteristic) flux, identical
+!! in form to the 2-D model's; see riemannflux3d_LinearEuler3D_t.
 !!
 
   use self_model
   use self_dgmodel3D
   use self_mesh
+  use SELF_BoundaryConditions
 
   implicit none
 
   type,extends(dgmodel3D) :: LinearEuler3D_t
     ! Add any additional attributes here that are specific to your model
     real(prec) :: rho0 = 1.0_prec ! Reference density
-    real(prec) :: c = 1.0_prec ! Sound speed
+    real(prec) :: c = 1.0_prec ! Reference sound speed (used to fill variable 6 in initial conditions)
     real(prec) :: g = 0.0_prec ! gravitational acceleration (y-direction only)
 
   contains
     procedure :: SourceMethod => sourcemethod_LinearEuler3D_t
     procedure :: SetNumberOfVariables => SetNumberOfVariables_LinearEuler3D_t
     procedure :: SetMetadata => SetMetadata_LinearEuler3D_t
+    procedure :: AdditionalInit => AdditionalInit_LinearEuler3D_t
     procedure :: entropy_func => entropy_func_LinearEuler3D_t
-    !procedure :: hbc3D_NoNormalFlow => hbc3D_NoNormalFlow_LinearEuler3D_t
     procedure :: flux3D => flux3D_LinearEuler3D_t
     procedure :: riemannflux3D => riemannflux3D_LinearEuler3D_t
     procedure :: SphericalSoundWave => SphericalSoundWave_LinearEuler3D_t
@@ -86,7 +94,7 @@ subroutine SetNumberOfVariables_LinearEuler3D_t(this)
     implicit none
     class(LinearEuler3D_t),intent(inout) :: this
 
-    this%nvar = 5
+    this%nvar = 6
 
   endsubroutine SetNumberOfVariables_LinearEuler3D_t
 
@@ -109,38 +117,65 @@ subroutine SetMetadata_LinearEuler3D_t(this)
     call this%solution%SetName(5,"P") ! Pressure
     call this%solution%SetUnits(5,"kg⋅m⁻¹⋅s⁻²")
 
+    call this%solution%SetName(6,"c") ! Sound speed (static; possibly heterogeneous)
+    call this%solution%SetUnits(6,"m⋅s⁻¹")
+
   endsubroutine SetMetadata_LinearEuler3D_t
 
+  subroutine AdditionalInit_LinearEuler3D_t(this)
+    !! Register the (CPU) radiation boundary condition. GPU builds overwrite
+    !! this registration with the device kernel in AdditionalInit_LinearEuler3D.
+    implicit none
+    class(LinearEuler3D_t),intent(inout) :: this
+    ! Local
+    procedure(SELF_bcMethod),pointer :: bcfunc
+
+    bcfunc => hbc3d_Radiation_LinearEuler3D
+    call this%hyperbolicBCs%RegisterBoundaryCondition( &
+      SELF_BC_RADIATION,"radiation",bcfunc)
+
+  endsubroutine AdditionalInit_LinearEuler3D_t
+
+  subroutine hbc3d_Radiation_LinearEuler3D(bc,mymodel)
+    !! Radiation BC: zero acoustic perturbation in the exterior state; the
+    !! sound speed (variable 6) is copied from the interior side so the
+    !! Riemann solver sees a consistent c.
+    class(BoundaryCondition),intent(in) :: bc
+    class(Model),intent(inout) :: mymodel
+    ! Local
+    integer :: n,i,j,iEl,s
+
+    select type(m => mymodel)
+    class is(LinearEuler3D_t)
+      do n = 1,bc%nBoundaries
+        iEl = bc%elements(n)
+        s = bc%sides(n)
+        do j = 1,m%solution%interp%N+1
+          do i = 1,m%solution%interp%N+1
+            m%solution%extBoundary(i,j,s,iEl,1:5) = 0.0_prec
+            m%solution%extBoundary(i,j,s,iEl,6) = m%solution%boundary(i,j,s,iEl,6) ! c preserved
+          enddo
+        enddo
+      enddo
+    endselect
+
+  endsubroutine hbc3d_Radiation_LinearEuler3D
+
   pure function entropy_func_LinearEuler3D_t(this,s) result(e)
     !! The entropy function is the sum of kinetic and internal energy
     !! For the linear model, this is
     !!
     !! \begin{equation}
-    !!   e = \frac{1}{2} \left( \rho_0*( u^2 + v^2 ) + \frac{P^2}{\rho_0 c^2} \right)
+    !!   e = \frac{1}{2} \left( \rho_0*( u^2 + v^2 + w^2 ) + \frac{P^2}{\rho_0 c^2} \right)
     class(LinearEuler3D_t),intent(in) :: this
     real(prec),intent(in) :: s(1:this%nvar)
     real(prec) :: e
 
-    e = 0.5_prec*this%rho0*(s(2)*s(2)+s(3)*(3)+s(4)*s(4))+ &
-        0.5_prec*(s(5)*s(5)/(this%rho0*this%c*this%c))
+    e = 0.5_prec*this%rho0*(s(2)*s(2)+s(3)*s(3)+s(4)*s(4))+ &
+        0.5_prec*(s(5)*s(5)/(this%rho0*s(6)*s(6)))
 
   endfunction entropy_func_LinearEuler3D_t
 
-  ! pure function hbc3D_NoNormalFlow_LinearEuler3D_t(this,s,nhat) result(exts)
-  !   class(LinearEuler3D_t),intent(in) :: this
-  !   real(prec),intent(in) :: s(1:this%nvar)
-  !   real(prec),intent(in) :: nhat(1:2)
-  !   real(prec) :: exts(1:this%nvar)
-  !   ! Local
-  !   integer :: ivar
-
-  !   exts(1) = s(1) ! density
-  !   exts(2) = (nhat(2)**2-nhat(1)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(3) ! u
-  !   exts(3) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! v
-  !   exts(4) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! w
-  !   exts(5) = s(4) ! p
-
-  ! endfunction hbc3D_NoNormalFlow_LinearEuler3D_t
   subroutine sourcemethod_LinearEuler3D_t(this)
     implicit none
     class(LinearEuler3D_t),intent(inout) :: this
@@ -171,50 +206,58 @@ pure function flux3D_LinearEuler3D_t(this,s,dsdx) result(flux)
     flux(4,2) = 0.0_prec ! z-velocity, y flux; 0
     flux(4,3) = s(5)/this%rho0 ! z-velocity, z flux; p/rho0
 
-    flux(5,1) = this%c*this%c*this%rho0*s(2) ! pressure, x flux : rho0*c^2*u
-    flux(5,2) = this%c*this%c*this%rho0*s(3) ! pressure, y flux : rho0*c^2*v
-    flux(5,3) = this%c*this%c*this%rho0*s(4) ! pressure, y flux : rho0*c^2*w
+    flux(5,1) = s(6)*s(6)*this%rho0*s(2) ! pressure, x flux : rho0*c^2*u
+    flux(5,2) = s(6)*s(6)*this%rho0*s(3) ! pressure, y flux : rho0*c^2*v
+    flux(5,3) = s(6)*s(6)*this%rho0*s(4) ! pressure, z flux : rho0*c^2*w
+
+    flux(6,1) = 0.0_prec ! sound speed, x flux; 0 (c held fixed in time)
+    flux(6,2) = 0.0_prec ! sound speed, y flux; 0
+    flux(6,3) = 0.0_prec ! sound speed, z flux; 0
     if(.false.) flux(1,1) = flux(1,1)+dsdx(1,1) ! suppress unused-dummy-argument warning
 
   endfunction flux3D_LinearEuler3D_t
 
   pure function riemannflux3D_LinearEuler3D_t(this,sL,sR,dsdx,nhat) result(flux)
-    !! Uses a local lax-friedrich's upwind flux
-    !! The max eigenvalue is taken as the sound speed
+    !! Impedance-matched (characteristic/Godunov) Riemann flux, identical in
+    !! form to the 2-D model's. The interface states are resolved with the
+    !! acoustic impedances Z = rho0*c on either side, so material interfaces
+    !! in a heterogeneous sound-speed field are handled with the physically
+    !! correct transmission/reflection:
+    !!
+    !!   un* = (ZL*unL + ZR*unR + (pL - pR)) / (ZL + ZR)
+    !!   p*  = (ZR*pL + ZL*pR + ZL*ZR*(unL - unR)) / (ZL + ZR)
+    !!
+    !! The sound speed variable carries zero flux.
     class(LinearEuler3D_t),intent(in) :: this
     real(prec),intent(in) :: sL(1:this%nvar)
     real(prec),intent(in) :: sR(1:this%nvar)
     real(prec),intent(in) :: dsdx(1:this%nvar,1:3)
     real(prec),intent(in) :: nhat(1:3)
     real(prec) :: flux(1:this%nvar)
     ! Local
-    real(prec) :: fL(1:this%nvar)
-    real(prec) :: fR(1:this%nvar)
-    real(prec) :: u,v,w,p,c,rho0
-
-    u = sL(2)
-    v = sL(3)
-    w = sL(4)
-    p = sL(5)
+    real(prec) :: rho0,cL,cR,ZL,ZR,unL,unR,pL,pR,un_star,p_star,c2_avg
+
     rho0 = this%rho0
-    c = this%c
-    fL(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
-    fL(2) = p*nhat(1)/rho0 ! u
-    fL(3) = p*nhat(2)/rho0 ! v
-    fL(4) = p*nhat(3)/rho0 ! w
-    fL(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure
-
-    u = sR(2)
-    v = sR(3)
-    w = sR(4)
-    p = sR(5)
-    fR(1) = rho0*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! density
-    fR(2) = p*nhat(1)/rho0 ! u
-    fR(3) = p*nhat(2)/rho0 ! v'
-    fR(4) = p*nhat(3)/rho0 ! w
-    fR(5) = rho0*c*c*(u*nhat(1)+v*nhat(2)+w*nhat(3)) ! pressure
-
-    flux(1:5) = 0.5_prec*(fL(1:5)+fR(1:5))+c*(sL(1:5)-sR(1:5))
+    cL = sL(6)
+    cR = sR(6)
+    ZL = rho0*cL
+    ZR = rho0*cR
+
+    unL = sL(2)*nhat(1)+sL(3)*nhat(2)+sL(4)*nhat(3)
+    unR = sR(2)*nhat(1)+sR(3)*nhat(2)+sR(4)*nhat(3)
+    pL = sL(5)
+    pR = sR(5)
+
+    un_star = (ZL*unL+ZR*unR+(pL-pR))/(ZL+ZR)
+    p_star = (ZR*pL+ZL*pR+ZL*ZR*(unL-unR))/(ZL+ZR)
+    c2_avg = 0.5_prec*(cL*cL+cR*cR)
+
+    flux(1) = rho0*un_star
+    flux(2) = p_star*nhat(1)/rho0
+    flux(3) = p_star*nhat(2)/rho0
+    flux(4) = p_star*nhat(3)/rho0
+    flux(5) = rho0*c2_avg*un_star
+    flux(6) = 0.0_prec
     if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning
 
   endfunction riemannflux3D_LinearEuler3D_t
@@ -260,6 +303,7 @@ subroutine SphericalSoundWave_LinearEuler3D_t(this,rhoprime,Lr,x0,y0,z0)
       this%solution%interior(i,j,k,iEl,3) = 0.0_prec
       this%solution%interior(i,j,k,iEl,4) = 0.0_prec
       this%solution%interior(i,j,k,iEl,5) = rho*this%c*this%c
+      this%solution%interior(i,j,k,iEl,6) = this%c ! uniform background sound speed
 
     enddo