Hi, I’m trying to use DG method to solve 2D shallow water equations. A problem is that I don’t know how to calculate the maximum Eigenvalue α in the Lax-Friedrichs flux using FEniCS. For simplification, the α can be expressed as max(inner(u,n)), where u is an already known vector, n is the outward normal vector for the cells’ edges, max means to find the maximum inner(u,n) over all the cells in the domain. Following simple code shows my problem:
from dolfin import *
parameters["ghost_mode"] = "shared_facet"
# define mesh and FunctionSpace
mesh=UnitSquareMesh(2,2)
V=FunctionSpace(mesh,'DG',1)
# the already known vector
x=interpolate(Expression('x[0]',degree=1),V)
y=interpolate(Expression('2*x[1]',degree=1),V)
u=as_vector((x,y))
# FacetNormal for cells
n=FacetNormal(mesh)
# how to get the max value of inner(u,n) over all the cells in the whole domain?
# print(max(inner(u,n)))
Dear nate, thanks for your reply. Indeed I’m not a math major and not very familiar with the FEM theroy. I gusse, according to the literature I read, I want the global maximum α in Lax-Friedrichs flux but not the local. Here is the relavent part of the literature:
In addition, I’m not very clear how to “interpolate your eigenvalue to a DG0 space and compute the ℓ∞ norm of the underlying vector”. I have tried something as following:
from dolfin import *
parameters["ghost_mode"] = "shared_facet"
# define mesh and FunctionSpace
mesh=UnitSquareMesh(2,2)
V=FunctionSpace(mesh,'DG',1)
# the already known vector
x=interpolate(Expression('1',degree=1),V)
y=interpolate(Expression('-1',degree=1),V)
u=as_vector((x,y))
# FacetNormal for cells
n=FacetNormal(mesh)
# how to get the max value of inner(u,n) over all the cells in the whole domain?
# print(max(inner(u,n)))
V_=FunctionSpace(mesh,'DG',0)
ut=TrialFunction(V_)
v=TestFunction(V_)
A=assemble(ut('+')*jump(v)*dS)
b=assemble(dot(u('+'),n('+'))*jump(v)*dS)
ux=Function(V_)
solve(A,ux.vector(),b)
print(ux.vector()[:])
In theory, all the inner(u,n) at the cells’ interior edges should be 1.414 or -1.414 or 1 or -1, and I suppose (or maybe I understand the FEM wrong) the lenth of the result vector should be 16 since the mesh has 8 shared interior edges for the cells and each edge is shared by 2 cells.
Very very few of us are math majors, but we appreciate (suffer through?) its utility for communication.
I don’t have a local installation of dolfin, but maybe this example will still work.
It computes a space-time FE approximation of the Burgers equation solution. This is nonlinear advection but very similar to the system I assume you’re trying to solve. I use a local Lax-Friedrichs flux in this example, but you could probably modify it for your global dissipation parameter formulation. It should also be reasonably straight forward to modify it for linear advection (as I believe that’s the problem you’re attempting to discretise).
Regarding
As an approximation you could project or interpolate your eigenvalue to a DG0 space and compute the \ell_\infty norm of the underlying vector.
This will compute a cell-wise average of the flux and then compute the maximum over all cells. The pseudocode is something like this with dolfin:
...
# prior code specifying a known flux, u.
...
DG0 = FunctionSpace(mesh, "DG", 0)
a_dg0 = project(u, DG0) # This is expensive, computes the solution of a linear system. Dolfinx has much more powerful and far cheaper interpolation of ufl forms
a_max = a_dg0.vector().norm("linf")
print(a_max) # Check this looks reasonable
