Linear Euler (3D)#
Definition#
The SELF_LinearEuler3D_t module defines the LinearEuler3D_t class. In SELF, models are posed in the form of a generic conservation law
where \(\vec{s}\) is a vector of solution variables, \(\overleftrightarrow{f}\) is the conservative flux, and \(\vec{q}\) are non-conservative source terms.
For the Linear Euler equations in 3-D
where \(u\), \(v\), and \(w\) are the \(x\), \(y\), and \(z\) components of the fluid velocity (respectively), \(p\) is the pressure, \(c\) is the speed of sound, \(\rho_0\) is the background density, and \(\sigma\) is the relaxation rate (the inverse of a local damping time scale) used to build sponge/damping layers. The density anomaly is not carried as a solution variable: for a motionless background state it is slaved to the pressure through the acoustic relation \(\rho = p/c^2\) and never feeds back into the velocity or pressure dynamics, so only the velocity components and the pressure are forward-stepped. If the density anomaly is needed as a diagnostic, it can be recovered pointwise as \(\rho = p/c^2\). The sound speed \(c\), the background density \(\rho_0\), and the relaxation rate \(\sigma\) are carried as solution variables so that heterogeneous (spatially varying) media and sponge/damping layers are supported; they have zero flux and zero source, so they are static in time and set by the initial condition (mirroring the 2-D model). This is entropy-stable for piecewise-constant material regions aligned with element boundaries. When we assume an ideal gas, and a motionless background state, the conservative fluxes are
The source term is the linear relaxation ("sponge"/damping) term
With \(\sigma = 0\) everywhere - the default, since the solution array is zero-initialized - the source vanishes identically and the model reduces to the undamped linear Euler system. See Sponge layers and damping below.
To track stability of the Euler equation, the total entropy function is
Implementation#
The Linear Euler 3D model is implemented as a type extension of the DGModel3D class. The LinearEuler3D_t class keeps scalar rho0 and c attributes (used as the reference values that fill variables 6 and 5 in the built-in initial conditions) and overrides SetNumberOfVariables (to declare nvar = 7 with nstepped = 4, so the last three variables are static), SetMetadata, AdditionalInit, entropy_func, flux3d, riemannflux3d, and SourceMethod. The sound speed lives in solution(:,:,:,:,5), the background density in solution(:,:,:,:,6), and the relaxation rate in solution(:,:,:,:,7); all three can be set independently per node when initializing the simulation.
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 per-side acoustic impedances \(Z_L = \rho_{0,L} c_L\) and \(Z_R = \rho_{0,R} c_R\) evaluated from the per-node background density and sound speed on either side of the face, the interface normal velocity and pressure are
and the normal flux is
Because the interface states are resolved with the per-side impedances, material interfaces in a heterogeneous field (density and/or sound-speed jumps) produce the physically correct transmission and reflection. The face-averaged \(\overline{\rho_0}\) and \(\overline{c^2}\) are used to reconstruct the momentum/pressure fluxes. The details for this implementation can be found in self_lineareuler3d_t.f90
Boundary conditions#
When initializing the mesh for your Euler 3D equation solver, you can change the boundary conditions to
SELF_BC_RADIATIONto set the external acoustic state on model boundaries to 0 in the Riemann solver, prolonging \(c\), \(\rho_0\) and \(\sigma\) from the interior side (open/non-reflecting).SELF_BC_PRESCRIBEDto set a prescribed external state, via a handler you register yourself.
Warning
Unlike the 2-D model, LinearEuler3D_t does not register a SELF_BC_NONORMALFLOW handler — AdditionalInit registers radiation only, on both the CPU and the GPU path. A face tagged SELF_BC_NONORMALFLOW therefore receives no exterior-state update from this model, and in particular no \(\sigma\) prolongation. Use SELF_BC_PRESCRIBED with your own reflecting handler until a 3-D no-normal-flow handler is implemented. This no longer fails quietly: the first ForwardStep warns that those faces carry a bcid with no registered boundary condition (see Diagnosing unhandled boundaries).
As an example, when using the built-in structured mesh generator, you can do the following
type(Mesh3D),target :: mesh
integer :: bcids(1:6)
bcids(1:6) = (/&
SELF_BC_RADIATION,& ! Bottom boundary condition
SELF_BC_RADIATION,& ! South boundary condition
SELF_BC_RADIATION,& ! East boundary condition
SELF_BC_PRESCRIBED,& ! North boundary condition
SELF_BC_RADIATION,& ! West boundary condition
SELF_BC_RADIATION & ! Top boundary condition
/)
call mesh%StructuredMesh(nxPerTile=5,nyPerTile=5,nzPerTile=5,&
nTileX=2,nTileY=2,nTileZ=2,&
dx=0.1_prec,dy=0.1_prec,dz=0.1_prec,bcids)
Note
See the Structured Mesh documentation for details on using the structuredmesh procedure
Note
To set a prescribed state as a function of position and time, you can create a type-extension of the LinearEuler3D class and override the hbc3d_Prescribed
The no-normal-flow boundary condition#
Note
This section describes the intended formulation. As noted above, LinearEuler3D_t does not yet register a handler for it.
To set the no-normal-flow boundary condition in SELF, we set the external state that is used as input to a Riemann solver. To determine the three components of the velocity field, we use the following conditions
- \(\vec{u}_{ext}\cdot \hat{n} = -\vec{u}_{in}\cdot \hat{n}\)
- \(\vec{u}_{ext}\cdot \hat{t}_1 = \vec{u}_{in}\cdot \hat{t}_1\)
- \(\vec{u}_{ext}\cdot \hat{t}_2 = \vec{u}_{in}\cdot \hat{t}_2\)
where \(\hat{n}\) is the outward pointing unit normal vector and \(\hat{t}_1\) and \(\hat{t}_2\) are mutually orthogonal vectors that are tangent to the boundary surface.
Sponge layers and damping#
The relaxation rate \(\sigma\) (solution variable 7, units s\(^{-1}\)) is the inverse of a local damping time scale. Wherever \(\sigma > 0\), the source term
pulls the acoustic state back toward the motionless, unperturbed background. Because the same rate is applied to all three velocity components and to the pressure, the local acoustic energy density decays as \(e^{-2\sigma t}\) and the local characteristic structure is untouched: the damping introduces no impedance mismatch of its own between a damped and an undamped region. That is what makes it usable as a sponge layer — a shell of elements next to an open boundary in which outgoing waves are absorbed before they can reflect back into the region of interest.
\(\sigma\) is set exactly the way \(c\) and \(\rho_0\) are set — it is a static, per-node solution variable with zero flux and zero source of its own, and it is excluded from the time-stepped variables, so once written it is preserved bitwise for the whole run:
\(\sigma = 0\) — the default, since the solution array is zero-initialized — recovers the undamped equations exactly, so an existing setup produces the same results as before without being changed.
The variable count did change, though, and that is a public-surface change: nvar goes from 6 to 7. Two things follow for existing code.
- A fixed-length array sized to
nvarneeds one more entry. The in-tree case isRefinementIndicator3D'sSetEnergyWeights, which validatessize(w) == nVarwhenEstimateruns, so a 6-element weight vector now fails the check — give \(\sigma\) a zero weight, as the background fields already have. The same applies to user code that declares a state vector asreal(prec) :: s(1:6). - Pickup files written before \(\sigma\) existed still load.
ReadModelskips a solution variable whose dataset the file does not hold, leaves it at its initialized value, and reports which variable it skipped — so a legacy restart comes back with \(\sigma = 0\), the configuration it was written under.
A prescribed-boundary handler does not need updating: \(\sigma\) carries no flux, so its exterior state never enters the tendency.
Choosing a profile#
Two rules of thumb:
- Ramp \(\sigma\) smoothly from zero. A jump in \(\sigma\) at the inner edge of the layer partially reflects the incoming wave. Start the layer at \(\sigma = 0\) and grow it toward the boundary, e.g. quadratically or cubically in the normalized depth into the layer.
- Make the layer a few elements thick and pick \(\sigma_{max}\) from the transit time. A wave crossing a layer of thickness \(L\) at speed \(c\) is attenuated in amplitude by \(\exp\left(-\int \sigma\, \mathrm{d}t\right)\); for a quadratic ramp that integral is \(\sigma_{max} L / (3c)\). Aim for a value of a few (an amplitude reduction of \(10^{-2}\) to \(10^{-3}\)) rather than for the largest \(\sigma\) that is stable — an over-strong layer reflects more than it absorbs. Keep \(\sigma \Delta t \lesssim 1\) so that the explicit time integrator resolves the relaxation.
A layer of thickness layer next to every domain boundary, with a quadratic ramp up to sigma_max, is written as
do iel = 1,mesh%nElem
do k = 1,modelobj%solution%N+1
do j = 1,modelobj%solution%N+1
do i = 1,modelobj%solution%N+1
x = modelobj%geometry%x%interior(i,j,k,iel,1,1)
y = modelobj%geometry%x%interior(i,j,k,iel,1,2)
z = modelobj%geometry%x%interior(i,j,k,iel,1,3)
! distance to the nearest domain boundary
d = min(x-xmin,xmax-x,y-ymin,ymax-y,z-zmin,zmax-z)
if(d < layer) then
! zero at the inner edge of the layer, sigma_max at the boundary
modelobj%solution%interior(i,j,k,iel,7) = sigma_max*((layer-d)/layer)**2
else
modelobj%solution%interior(i,j,k,iel,7) = 0.0_prec
endif
enddo
enddo
enddo
enddo
call modelobj%solution%UpdateDevice()
Pair the layer with SELF_BC_RADIATION on the outer boundaries: the sponge removes most of the outgoing energy and the radiation condition handles the remainder. SELF_BC_NONORMALFLOW also works — the sponge then absorbs the wave before it can reflect off the wall — which is useful when the mesh has no open boundary.
The relaxation is not restricted to boundary layers: a spatially uniform \(\sigma\) acts as a bulk absorber, and a \(\sigma\) supported on an arbitrary subregion damps only there.
Warning
\(\sigma\) is expected to be non-negative. A negative \(\sigma\) amplifies the solution and is not a supported configuration.
GPU Acceleration#
When building SELF with GPU acceleration enabled, the Linear Euler (3-D) model overrides the following DGModel3D type-bound procedures
BoundaryFluxFluxMethodSourceMethodSetBoundaryConditionSetGradientBoundaryCondition
These methods are one-level above the usual pure function type-bound procedures used to define the riemann solver, flux, source terms, and boundary conditions. These procedures need to be overridden with calls to GPU accelerated kernels to make the solver fully resident on the GPU.
The relaxation source term is fully device-resident: SourceMethod launches the sourcemethod_LinearEuler3D_gpu kernel, which reads \(\sigma\) from solution variable 7 on the device, so enabling a sponge layer adds no host-device traffic.
Out-of-the-box, the no-normal-flow and radiation boundary conditions are GPU accelerated. However, prescribed boundary conditions are CPU-only. We have opted to keep the prescribed boundary conditions CPU-only so that their implementation remains easy-to-use. This implies that some data is copied between host and device every iteration when prescribed boundary conditions are enabled.
Note
In simulations where no prescribed boundaries are used, or your prescribed boundaries are time independent, you can disable prescribed boundary conditions by explicitly setting modelobj % prescribed_bcs_enabled = .false.. This can improve the time-to-solution for your simulation by avoiding unnecessary host-device memory movement. An example of this feature is shown in examples/lineareuler3d_spherical_soundwave_radiation.f90
Example usage#
For examples, see any of the following
examples/linear_euler3d_spherical_soundwave_radiation.f90- Implements a simulation with a gaussian pressure anomaly as an initial condition in a domain with radiation boundary conditions on all sides. Uses theSphericalSoundWaveinitializer, which fills the uniform sound speed and background density fields.test/lineareuler3d_sponge_damping.f90- Sponge-layer regression test. A Gaussian pulse is released at the center of a cubic domain with radiation boundaries on all sides; the outermost shell of elements carries a quadratically ramped \(\sigma\). The same problem is run twice, with \(\sigma = 0\) and with the sponge, and the two are compared. Shows the profile-setting pattern described in Sponge layers and damping.