This repository is no longer maintained; please see StartUpDG.jl or Trixi.jl for my current Julia-based DG codes.
This package provides utilities for flux differencing and computation of flux differencing Jacobians in terms of derivatives of flux functions. The code based in part on the preprint "Efficient computation of Jacobian matrices for entropy stable summation-by-parts schemes".
The routines are meant to be fairly general, but specialize depending on whether the operators are general arrays or SparseMatrixCSC (to capitalize on sparsity).
using LinearAlgebra
using FluxDiffUtils
using Test
# make 3-field solution
u = collect(LinRange(-1,1,4))
U = (u,u,u)
# define a flux fS
avg(x,y) = @. .5*(x+y)
function fS(UL,UR)
uL,vL,wL = UL
uR,vR,wR = UR
Fx = avg(uL,uR),avg(vL,vR),avg(wL,wR)
Fy = avg(uL,uR),avg(vL,vR),avg(wL,wR)
return SVector{3}(Fx),SVector{3}(Fy)
end
# jacobians w.r.t. (uR,vR)
dfS(UL,UR) = ([.5 0 0; 0 .5 0; 0 0 .5], [.5 0 0; 0 .5 0; 0 0 .5])
A_list = (A->A+A').(ntuple(x->randn(4,4),2)) # make symmetric to check formula
rhs = hadamard_sum(A_list,fS,U)
jac = hadamard_jacobian(A_list, dfS, U)
# jac = hadamard_jacobian(A_list, :sym, df, U) # optimized version
jac11_exact = sum((A->.5*(A + diagm(vec(sum(A,dims=1))))).(A_list))
@test norm(jac11_exact-jac[1,1]) < 1e-12
# converts tuple-block storage of jac to a global matrix
jac_global = blockcat(size(jac,2),jac)