DG and the L2 space

Main Forum Old-Forum - do not post here
1

Dear NGSolve community,

I am currently dealing with a function that solves an advection-dominated equation,
which I let NGSolve approximate in L2 space with DG.
Analytically, this function is smooth enough to have a gradient which I want to calculate.
In this context I played a little with the L2 space and came up with behaviour I do not understand,
so I kindly ask for elaboration on this:

  1. Gradients of functions from L2 with order=1 seem to evaluate to zero, but increasing the order
    gives non-zero gradients. I would understand that the gradient isn’t implemented since it isn’t expected to exists analytically - although in a discrete sense one could always calculate it element-wise, right? -, but am puzzled with the behaviour for higher orders. Is this related with different basis functions?

  2. I can evaluate L2 functions on vertices and also the gradient of higher order L2 functions although
    it isn’t continous across elements (I found discontinuous gradient - Kunena). How can I understand the value at the vertex (which is shared by multiple elements)?

  3. I found that when doing mesh refinement (code below), the gradient of higher order L2 functions
    vanishes. Is this intended and if yes, is there a proper way to work around this?

from ngsolve import *       
from netgen.geom2d import unit_square
                            
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))                  
fes = L2(mesh, order=4)     
phi = GridFunction(fes, nested = True)
phi.Set(atan(sqrt((x-0.5)**2 + (y-0.5)**2)-0.2))
               
print("phi: ", max(phi.vec), ", ", min(phi.vec))
for v in mesh.vertices:
    print("grad ", cf_gphi(mesh(v.point[0], v.point[1])))
               
for i in range(3):
    mesh.Refine()                     
    phi.Update()                                      
    fes.Update()

cf_gphi = grad(phi)
for v in mesh.vertices:
    print("grad ", cf_gphi(mesh(v.point[0], v.point[1])))
                         
print("phi: ", max(phi.vec), ", ", min(phi.vec))

I would be very grateful if you could help me out here: DG is kind of new to me.

All the best,
Philipp

2
  1. This code
from ngsolve import *       
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.5))                  
fes = L2(mesh, order=1)     
phi = GridFunction(fes, nested = True)
phi.Set(x)
cf_gphi = grad(phi)
for v in mesh.vertices:
    print("grad ", cf_gphi(mesh(v.point[0], v.point[1])))

gives this output

grad  (1.0000000000000002, 1.1102230246251384e-16)
grad  (0.9999999999999997, -2.220446049250313e-16)
grad  (1.0000000000000004, 1.1102230246251383e-16)
grad  (1.0, 0.0)
grad  (1.0000000000000002, 1.1102230246251384e-16)
grad  (0.9999999999999997, -2.220446049250313e-16)
grad  (1.0000000000000004, 1.1102230246251383e-16)
grad  (1.0000000000000002, 1.1102230246251384e-16)

If you set the order to 0, the base functions are constant on each element, so the local gradient is 0.

  1. When you evaluate a gridfunction at a point, what actually happens is this: First, some element that contains the point is found, the basis functions that belong to that element are evaluated at this point and the values are added according to the gridfunction’s coefficients for that element. Vertices are contained in multiple elements, so evaluation will give you the value the restriction of the gridfunction to one of these elements has in that point. This is not mathematically well defined for arbitrary L2 functions, but will give you the right value if what you evaluate is continuous.

  2. This behavior is not intuitive. Take a look not at the gradient of the gridfunction but at the gridfunction itself. You will see that after refining, the gridfunction is piecewise constant, so it’s “gradient” is zero. What happens here is that

   phi.Update()

only updates the lowest order DOF for each element, so the constant base function. The higher order ones are just not implemented.

A workaround would be to create a copy of the unrefined mesh and interpolate to the refined level manually (see attachment).

Best,
Lukas

Attachment: discg.py

3

Hello Lukas,

  1. I somehow expected the default order of L2 to be one and just wrote

fes = L2(mesh), which is kinda stupid. I should’ve checked. Sorry!

  1. and 3) Thank you very much. This clarifies a lot and is very helpful.

All the best,
Philipp