add SubVolumeElement

src/dg1d.jl

@@ -41,6 +41,7 @@ module dg1d
   include("glgl.jl")
   include("lgl.jl")
   include("spectralelement.jl")
+  include("subvolumeelement.jl")
   include("mesh.jl")
 
   # custom evolution

src/subvolumeelement.jl

+export SubVolumeElement
+
+
+struct SubVolumeElement
+  Npts::Int64
+  z::Vector{Float64}      # flux nodes ∈ [-1,1], include boundary points
+  w::Vector{Float64}      # cell widths defined by z
+  Cs_lhs::Vector{Float64}     # correction coefficients lhs
+  Cs_rhs::Vector{Float64}     # correction coefficients rhs
+end
+
+
+function SubVolumeElement(Npts)
+  @assert Npts > 0
+  dz = 2 / (Npts-1)
+  z = [ -1 + dz * (i-1) for i = 1:Npts ] # include boundary points
+  w = dz * ones(Npts-1)                  # equidistant grid is equally weighted
+  Cs_lhs, Cs_rhs = correction_coefficients(z)
+  return SubVolumeElement(Npts, z, w, Cs_lhs, Cs_rhs)
+end
+
+
+function SubVolumeElement_from_nodes(z)
+  Npts = length(z)
+  zmin, zmax = extrema(z)
+  @assert -1 == zmin && zmax == 1
+  @assert all(isapprox.(z .+ z[end:-1:1], 0.0, atol=1e-15)) """
+  Sub-volume nodes must be symmetric across z = 0.
+  """
+  w = diff(z)
+  Cs_lhs, Cs_rhs = correction_coefficients(z)
+  return SubVolumeElement(Npts, z, w, Cs_lhs, Cs_rhs)
+end
+
+
+function SubVolumeElement_weights(node_w)
+  Npts = length(node_w)
+  @assert all(node_w .> 0) """
+  Sub-volume node weights must be positive.
+  """
+  # construct cell widths from node weights
+  w = 1 ./ abs.(diff(w)) # inverse weighting makes higher weighted nodes have smaller cells
+  w = w .* 2 ./ sum(w) # normalize weights such that sum(w) = 2
+  # construct flux nodes
+  z = [ -1, (-1 + cumsum(w[1:end]))..., 1 ]
+  Cs_lhs, Cs_rhs = correction_coefficients(z)
+  return SubVolumeElement(Npts, z, w, Cs_lhs, Cs_rhs)
+end
+
+
+function SubVolumeElement(element::SpectralElement, divider_wrt_quadrature=false)
+  if divider_wrt_quadrature
+    return SubVolumeElement_from_weights(element.w)
+  else
+    return SubVolumeElement_from_nodes(element.z)
+  end
+end
+
+
+function correction_coefficients(z)
+  Npts = length(z)
+  B = [ (-1)^(j+1) * binomial(Npts-1+j,j) * binomial(Npts,j) for j = 1:Npts ]
+  xis = [ (zi + 1) / 2 for zi in z ] # map [-1,1] to [0,1]
+  monomials_xi = [ xis[j]^i for i = 1:Npts, j = 1:Npts ]
+  Cs = 1 .- transpose(monomials_xi) * B
+  return Cs, Cs[end:-1:1]
+end