DG and the L2 space

More
4 years 4 months ago - 4 years 4 months ago #2979 by philipp
DG and the L2 space was created by philipp
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 ngsolve.org/forum/ngspy-forum/1090-discontinuous-gradient#2660 ). 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?
Code:
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
Last edit: 4 years 4 months ago by philipp.
More
4 years 4 months ago #2982 by lkogler
Replied by lkogler on topic DG and the L2 space
1) This code
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
Code:
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.

2) 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.

3) 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
Code:
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
Attachments:
The following user(s) said Thank You: philipp
More
4 years 4 months ago - 4 years 4 months ago #2986 by philipp
Replied by philipp on topic DG and the L2 space
Hello Lukas,

1) I somehow expected the default order of L2 to be one and just wrote
Code:
fes = L2(mesh)
, which is kinda stupid. I should've checked. Sorry!

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

All the best,
Philipp
Last edit: 4 years 4 months ago by philipp.
Time to create page: 0.123 seconds