Add two-point vector types and split-form divergence

.github/workflows/fprettify-lint.yml

@@ -5,11 +5,13 @@ on:
     branches:
       - main
     paths-ignore:
+      - 'docs/'
       - 'AUTHORS.md'
       - 'LICENSE.md'
       - 'README.md'
   pull_request:
     paths-ignore:
+      - 'docs/'
       - 'AUTHORS.md'
       - 'LICENSE.md'
       - 'README.md'

.github/workflows/main-docs.yml

@@ -3,6 +3,10 @@ on:
   push:
     branches:
       - main
+    paths:
+      - 'docs/**'
+      - 'mkdocs.yml'
+      - 'self.md'
 
 jobs:
   build:

docs/Learning/ProvableStability.md

@@ -98,7 +98,28 @@ where $s$ is the solution and $u$ is a velocity field. A convex entropy function
 
 
 ## Split form equations
+Split-form equations rewrite the flux divergence so that the resulting discretization satisfies a discrete entropy conservation law. For a scalar conservation law $s_t + (u s)_x = 0$, we can use the product rule to write the equivalent split form:
 
+\begin{equation}
+s_t + \frac{1}{2}\left( (u s)_x + u s_x + s u_x \right) = 0
+\end{equation}
+
+Discretizing the split form with the SBP derivative matrix $D$ and a symmetric two-point flux $f^\#(s_i, s_n)$ leads to the split-form DGSEM:
+
+\begin{equation}
+\frac{ds_i}{dt} = -2 \sum_n D_{\text{split},n,i} \, f^\#(s_i, s_n) - \frac{1}{w_i} \left( f_{\text{Riemann}} - f_{\text{local}} \right)\bigg|_{\partial \Omega_e}
+\end{equation}
+
+where $D_\text{split} = D - \frac{1}{2}M^{-1}B$ is the skew-symmetric split-form derivative operator (see [Split-Form DGSEM](SplitFormDGSEM.md) for details). The skew-symmetry of $D_\text{split}$ ensures the volume integral contributes zero to the discrete entropy rate; all entropy change passes through the surface Riemann flux.
 
 ## Two-point flux
+A two-point flux $\mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R)$ is a symmetric, consistent numerical flux that satisfies the Tadmor entropy conservation condition:
+
+\begin{equation}
+(\mathbf{w}_R - \mathbf{w}_L)^T \mathbf{f}^\# = \Psi_R - \Psi_L
+\end{equation}
+
+where $\mathbf{w} = \partial \eta / \partial \mathbf{s}$ are the entropy variables and $\Psi$ is the entropy flux potential. When used in the split-form volume integral, this condition guarantees that the discrete volume contribution to $d\eta/dt$ vanishes identically.
+
+For linear advection with $f = as$ and $\eta = s^2/2$, the arithmetic mean $f^\#(s_L, s_R) = a(s_L + s_R)/2$ satisfies this condition. For nonlinear systems (Euler, shallow water), deriving the entropy-conserving flux is more involved; see the references in [Split-Form DGSEM](SplitFormDGSEM.md).
 

docs/Learning/SoftwareArchitecture.md

@@ -78,7 +78,7 @@ and in 3-D
 F(0:N,0:N,0:N,0:N,1:nEl,1:nvar,1:3)
 ```
 
-The first dimension refers to the directional component (physical or computational) of the vector. The second dimension is a computational coordinate direction and variations in this dimensions are equivalent to looping over $n$ in the two-point vector representation shown above. The remaining array dimensions of size `0:N` are computational coordinate directions and are the usual `(i,j,k)` dimensions, in this order. This is followed by the dimension over the solution variables and then the elements. 
+The first dimension is a computational coordinate direction and variations in this dimension are equivalent to looping over $n$ in the two-point vector representation shown above. The next dimensions are `(i,j)` (2-D) or `(i,j,k)` (3-D) computational coordinate directions, followed by the elements, solution variables, and the directional component of the vector `(...,nEl,nvar,idir)`, respectively.
 
 
 ### Tensors

docs/Learning/SplitFormDGSEM.md

@@ -0,0 +1,166 @@
+# Split-Form Entropy-Conserving DGSEM
+
+This page documents the split-form discontinuous Galerkin spectral element method (DGSEM) implemented in SELF's `ECDGModel2D_t` and `ECDGModel3D_t` classes. The formulation follows Gassner, Winters, and Kopriva (2016) and is aligned with the [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) implementation for curved meshes.
+
+## Mathematical Background
+
+### Conservation Law
+We consider a system of conservation laws on a domain $\Omega$:
+
+$$
+\mathbf{s}_t + \vec{\nabla} \cdot \vec{\mathbf{f}}(\mathbf{s}) = \mathbf{q}
+$$
+
+where $\mathbf{s}$ is the state vector, $\vec{\mathbf{f}}$ is the flux, and $\mathbf{q}$ is a source term.
+
+### Entropy Conservation
+A convex entropy function $\eta(\mathbf{s})$ satisfies the entropy conservation law $\eta_t + \vec{\nabla} \cdot \vec{\Psi} = 0$ (see [Provable Stability](ProvableStability.md)). A numerical scheme is **entropy-conserving** if it satisfies a discrete analogue of this conservation law, and **entropy-stable** if entropy can only decrease (mimicking physical dissipation at discontinuities).
+
+### Two-Point Entropy-Conserving Flux
+An entropy-conserving two-point flux $\mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R)$ is a symmetric function of two states satisfying the Tadmor condition:
+
+$$
+(\mathbf{w}_R - \mathbf{w}_L)^T \mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R) = \Psi_R - \Psi_L
+$$
+
+where $\mathbf{w} = \partial \eta / \partial \mathbf{s}$ are the entropy variables and $\Psi$ is the entropy flux potential.
+
+**Example:** For linear advection $f = a s$ with entropy $\eta = s^2/2$, the arithmetic mean $f^\#(s_L, s_R) = a(s_L + s_R)/2$ is entropy-conserving.
+
+## The Split-Form Derivative Operator
+
+### Standard Derivative Matrix $D$
+The standard DGSEM derivative matrix $D$ satisfies $D \cdot \mathbf{1} = 0$ (derivative of a constant is zero) and the SBP property:
+
+$$
+M D + D^T M = B
+$$
+
+where $M = \text{diag}(w_0, \ldots, w_N)$ is the mass matrix (quadrature weights) and $B$ is the boundary operator.
+
+### Split-Form Matrix $D_\text{split}$
+The split-form derivative matrix is defined as:
+
+$$
+D_\text{split} = D - \tfrac{1}{2} M^{-1} B
+$$
+
+**Key property:** $D_\text{split}$ is skew-symmetric under the $M$ inner product:
+
+$$
+M D_\text{split} + D_\text{split}^T M = 0
+$$
+
+This follows directly from the SBP property: $M D_\text{split} + D_\text{split}^T M = (M D + D^T M) - B = B - B = 0$.
+
+**$D_\text{split}$ is NOT an SBP operator.** Its weighted symmetric part is zero (not the boundary term $B$). This means the volume integral formed with $D_\text{split}$ contributes **zero** to the entropy rate. All entropy change passes through the surface term.
+
+For Gauss-Lobatto-Legendre (GLL) nodes, $D_\text{split}$ reduces to $(D - D^T)/2$, and SELF requires GLL nodes for `TwoPointVector` types.
+
+In SELF, the split-form matrix is stored as `interp%dSplitMatrix` in the Lagrange class:
+
+```fortran
+dSplitMatrix(ii,i) = dMatrix(ii,i) &
+  - 0.5*(bMatrix(i,2)*bMatrix(ii,2) - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i)
+```
+
+### Relationship to Trixi.jl
+Trixi.jl defines `derivative_split = 2D - M^{-1}B` and applies it **without** a factor of 2. SELF defines `dSplitMatrix = D - 0.5 M^{-1}B` and multiplies by 2 in the divergence formula. These are algebraically identical:
+
+$$
+2 \cdot D_\text{split}^\text{SELF} = D_\text{split}^\text{Trixi}
+$$
+
+
+## The EC-DGSEM Semidiscretization
+
+### Full Tendency
+The `ECDGModel2D_t` computes the time tendency as:
+
+$$
+\frac{d\mathbf{s}}{dt} = \mathbf{q} - \frac{1}{J} \left[ 2 \sum_n D_{\text{split},n,i} \, \tilde{F}^r(n,i,j) \;+\; M^{-1} B^T \mathbf{f}_\text{Riemann} \right]
+$$
+
+where:
+
+- $\tilde{F}^r(n,i,j)$ is a **scalar contravariant** two-point flux for direction $r$
+- $\mathbf{f}_\text{Riemann}$ is the surface Riemann flux (e.g., Local Lax-Friedrichs)
+- $J$ is the element Jacobian
+
+### Equivalence to the Strong-Form Penalty
+Using $D_\text{split}$ for the volume combined with the plain Riemann flux on the surface is algebraically identical to using the standard $D$ for the volume with a **penalty** $(f_\text{Riemann} - f_\text{local})$ on the surface:
+
+$$
+2 D_\text{split} \tilde{F} + M^{-1} B^T f_\text{Riemann} = 2 D \tilde{F} + M^{-1} B^T (f_\text{Riemann} - \tilde{F}_\text{boundary})
+$$
+
+SELF uses the left-hand form (matching Trixi.jl's `SurfaceIntegralWeakForm`).
+
+## Contravariant Two-Point Fluxes on Curved Meshes
+
+### Convention
+Following Trixi.jl for curved meshes, `interior(n,i,j,iEl,iVar,r)` stores a **scalar** contravariant two-point flux for the $r$-th computational direction. Each direction uses its own node pairing and its own averaged metric:
+
+| Direction $r$ | Node pair | Averaged metric |
+|:---:|---|---|
+| $\xi^1$ | $(i,j) \leftrightarrow (n,j)$ | $\overline{J\mathbf{a}^1} = \tfrac{1}{2}\left(J\mathbf{a}^1_{i,j} + J\mathbf{a}^1_{n,j}\right)$ |
+| $\xi^2$ | $(i,j) \leftrightarrow (i,n)$ | $\overline{J\mathbf{a}^2} = \tfrac{1}{2}\left(J\mathbf{a}^2_{i,j} + J\mathbf{a}^2_{i,n}\right)$ |
+| $\xi^3$ | $(i,j,k) \leftrightarrow (i,j,n)$ | $\overline{J\mathbf{a}^3} = \tfrac{1}{2}\left(J\mathbf{a}^3_{i,j,k} + J\mathbf{a}^3_{i,j,n}\right)$ |
+
+The contravariant flux is the dot product of the averaged metric with the physical EC flux:
+
+$$
+\tilde{F}^r(n,i,j) = \sum_d \overline{J a^r_d} \; f^\#_d(\mathbf{s}_L, \mathbf{s}_R)
+$$
+
+where $f^\#_d$ is the $d$-th physical component of the entropy-conserving two-point flux evaluated between the direction-specific node pair.
+
+### Why Not Physical-Space Storage?
+A single `interior(n,i,j,...,d)` array with $d$ indexing physical directions cannot correctly serve both the $\xi^1$ and $\xi^2$ sums, because the same index $n$ refers to different physical node partners in each direction. Storing pre-projected scalar contravariant fluxes (one per computational direction) avoids this cross-contamination entirely.
+
+### Implementation
+`TwoPointFluxMethod` in `ECDGModel2D_t` computes the contravariant projection:
+
+```fortran
+! xi^1: pair (i,j)-(nn,j), project onto avg(Ja^1)
+sR = solution(nn,j,iel,:)
+Fphys = twopointflux2d(sL, sR)         ! physical EC flux (nvar x 2)
+Fc = sum_d avg(Ja^1_d) * Fphys(:,d)    ! scalar contravariant
+twoPointFlux%interior(nn,i,j,...,1) = Fc
+
+! xi^2: pair (i,j)-(i,nn), project onto avg(Ja^2)
+sR = solution(i,nn,iel,:)
+Fphys = twopointflux2d(sL, sR)
+Fc = sum_d avg(Ja^2_d) * Fphys(:,d)
+twoPointFlux%interior(nn,i,j,...,2) = Fc
+```
+
+`MappedDivergence` then simply applies `Divergence / J` -- no metric operations inside.
+
+## Implementing a Concrete EC Model
+
+To create an entropy-conserving model, extend `ECDGModel2D_t` (or `ECDGModel3D_t`) and override:
+
+| Procedure | Purpose |
+|-----------|---------|
+| `SetNumberOfVariables` | Set `this%nvar` |
+| `twopointflux2d(sL, sR)` | Return the physical EC two-point flux (nvar x 2 array) |
+| `riemannflux2d(sL, sR, dsdx, nhat)` | Return the surface Riemann flux (nvar array). Use LLF for symmetric dissipation. |
+| `entropy_func(s)` | Return the scalar entropy $\eta(\mathbf{s})$ for diagnostics |
+| `SetMetadata` | (Optional) name and units for each variable |
+
+The base class handles the split-form volume integral, metric averaging, surface correction, time integration (inherited from `DGModel2D_t`), and I/O.
+
+**Example:** `ECAdvection2D_t` implements entropy-conserving linear advection with:
+- `twopointflux2d`: arithmetic mean $f^\# = \mathbf{a}(s_L + s_R)/2$
+- `riemannflux2d`: Local Lax-Friedrichs $f_\text{LLF} = \tfrac{1}{2}\left[u_n(s_L+s_R) - |\mathbf{a}|(s_R-s_L)\right]$
+- `entropy_func`: $\eta = s^2/2$
+
+
+## References
+
+1. Gassner, G. J., Winters, A. R., & Kopriva, D. A. (2016). Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations. *Journal of Computational Physics*, 327, 39-66.
+
+2. Winters, A. R., Kopriva, D. A., Gassner, G. J., & Hindenlang, F. (2021). Construction of Modern Robust Nodal Discontinuous Galerkin Spectral Element Methods for the Compressible Navier-Stokes Equations. *Lecture Notes in Computational Science and Engineering*, Springer.
+
+3. Ranocha, H., Schlottke-Lakemper, M., Winters, A. R., Faulhaber, E., Chan, J., & Gassner, G. J. (2022). Adaptive numerical simulations with Trixi.jl: A case study of Julia for scientific computing. *Proceedings of the JuliaCon Conferences*.

src/SELF_DGModel1D_t.f90

@@ -307,7 +307,7 @@ subroutine CalculateEntropy_DGModel1D_t(this)
       do i = 1,this%solution%interp%N+1
         J = this%geometry%dxds%interior(i,iel,1)
         s(1:this%solution%nvar) = this%solution%interior(i,iel,1:this%solution%nvar)
-        e = e+this%entropy_func(s)*J
+        e = e+this%entropy_func(s)*J*this%solution%interp%qWeights(i)
       enddo
     enddo
 

src/SELF_DGModel2D_t.f90

@@ -343,7 +343,9 @@ subroutine CalculateEntropy_DGModel2D_t(this)
         do i = 1,this%solution%interp%N+1
           jac = abs(this%geometry%J%interior(i,j,iel,1))
           s = this%solution%interior(i,j,iel,1:this%nvar)
-          e = e+this%entropy_func(s)*jac
+          e = e+this%entropy_func(s)*jac* &
+              this%solution%interp%qWeights(i)* &
+              this%solution%interp%qWeights(j)
         enddo
       enddo
     enddo

src/SELF_DGModel3D_t.f90

@@ -344,7 +344,10 @@ subroutine CalculateEntropy_DGModel3D_t(this)
           do i = 1,this%solution%interp%N+1
             jac = abs(this%geometry%J%interior(i,j,k,iel,1))
             s = this%solution%interior(i,j,k,iel,1:this%nvar)
-            e = e+this%entropy_func(s)*jac
+            e = e+this%entropy_func(s)*jac* &
+                this%solution%interp%qWeights(i)* &
+                this%solution%interp%qWeights(j)* &
+                this%solution%interp%qWeights(k)
           enddo
         enddo
       enddo

src/SELF_ECAdvection2D_t.f90

@@ -0,0 +1,131 @@
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+!
+! Maintainers : support@fluidnumerics.com
+! Official Repository : https://github.com/FluidNumerics/self/
+!
+! Copyright © 2024 Fluid Numerics LLC
+!
+! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+!
+! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+!
+! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
+!    the documentation and/or other materials provided with the distribution.
+!
+! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
+!    this software without specific prior written permission.
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+!
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+
+module SELF_ECAdvection2D_t
+  !! Entropy-conserving linear scalar advection equation in 2-D.
+  !!
+  !! Solves:
+  !!   du/dt + u_a * du/dx + v_a * du/dy = 0
+  !!
+  !! Entropy function: eta(u) = u^2 / 2  (same as Trixi.jl)
+  !!
+  !! EC two-point flux (arithmetic mean, entropy-conserving for eta = u^2/2):
+  !!   F^EC(uL, uR) = (u_a * (uL+uR)/2,  v_a * (uL+uR)/2)
+  !!
+  !! Surface Riemann flux (Local Lax-Friedrichs / Rusanov):
+  !!   F_Riemann = 0.5 * (a.n) * (uL+uR) - 0.5 * |a| * (uR-uL)
+  !! where |a| = sqrt(u^2 + v^2) is the maximum wave speed.
+  !! This is dissipative (entropy-stable) and reduces to the central flux
+  !! when uL = uR (no dissipation at no-normal-flow or mirror boundaries).
+
+  use SELF_ECDGModel2D_t
+  use SELF_mesh
+
+  implicit none
+
+  type,extends(ECDGModel2D_t) :: ECAdvection2D_t
+
+    real(prec) :: u ! x-component of advection velocity
+    real(prec) :: v ! y-component of advection velocity
+
+  contains
+
+    procedure :: entropy_func => entropy_func_ECAdvection2D_t
+    procedure :: twopointflux2d => twopointflux2d_ECAdvection2D_t
+    procedure :: riemannflux2d => riemannflux2d_ECAdvection2D_t
+
+    ! Mirror BC: sets extBoundary = interior boundary state so that
+    ! sR = sL at all domain faces.  Used to eliminate boundary dissipation
+    ! when testing entropy conservation.
+    procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_ECAdvection2D_t
+
+  endtype ECAdvection2D_t
+
+contains
+
+  pure function entropy_func_ECAdvection2D_t(this,s) result(e)
+    !! Quadratic entropy: eta(u) = u^2 / 2
+    class(ECAdvection2D_t),intent(in) :: this
+    real(prec),intent(in) :: s(1:this%nvar)
+    real(prec) :: e
+
+    e = 0.5_prec*s(1)*s(1)
+    if(.false.) e = e+this%u ! suppress unused-dummy-argument warning
+
+  endfunction entropy_func_ECAdvection2D_t
+
+  pure function twopointflux2d_ECAdvection2D_t(this,sL,sR) result(flux)
+    !! Arithmetic-mean two-point flux for linear advection.
+    !! Entropy-conserving with respect to eta(u) = u^2/2.
+    class(ECAdvection2D_t),intent(in) :: this
+    real(prec),intent(in) :: sL(1:this%nvar)
+    real(prec),intent(in) :: sR(1:this%nvar)
+    real(prec) :: flux(1:this%nvar,1:2)
+    ! Local
+    real(prec) :: savg
+
+    savg = 0.5_prec*(sL(1)+sR(1))
+    flux(1,1) = this%u*savg
+    flux(1,2) = this%v*savg
+
+  endfunction twopointflux2d_ECAdvection2D_t
+
+  pure function riemannflux2d_ECAdvection2D_t(this,sL,sR,dsdx,nhat) result(flux)
+    !! Local Lax-Friedrichs (Rusanov) Riemann flux for linear advection.
+    !! Entropy-stable: provides symmetric dissipation at element interfaces.
+    !! lambda_max = sqrt(u^2 + v^2) is the maximum wave speed.
+    class(ECAdvection2D_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:2)
+    real(prec),intent(in) :: nhat(1:2)
+    real(prec) :: flux(1:this%nvar)
+    ! Local
+    real(prec) :: un,lam
+
+    un = this%u*nhat(1)+this%v*nhat(2)
+    lam = sqrt(this%u*this%u+this%v*this%v)
+    flux(1) = 0.5_prec*(un*(sL(1)+sR(1))-lam*(sR(1)-sL(1)))
+    if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning
+
+  endfunction riemannflux2d_ECAdvection2D_t
+
+  pure function hbc2d_NoNormalFlow_ECAdvection2D_t(this,s,nhat) result(exts)
+    !! Mirror boundary condition: sets extBoundary = interior state.
+    !! With the LLF Riemann flux, this gives sR = sL at the boundary,
+    !! so the Riemann flux reduces to the central flux (a.n)*s — no
+    !! dissipation.  Use this BC when testing entropy conservation.
+    class(ECAdvection2D_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)
+
+    exts = s
+    if(.false.) exts(1) = exts(1)+nhat(1) ! suppress unused-dummy-argument warning
+
+  endfunction hbc2d_NoNormalFlow_ECAdvection2D_t
+
+endmodule SELF_ECAdvection2D_t