Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fluidnumericsJoe merged 1 commit into from
Jul 13, 2026
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 21 additions & 28 deletions docs/Models/linear-euler-3d-model.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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} =
Expand All @@ -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}
Expand All @@ -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
Expand Down
160 changes: 102 additions & 58 deletions src/SELF_LinearEuler3D_t.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
!!
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
Loading
Loading