Re: [firedrake] Is it possible to calculate a flux value on the edges of a cell in Firedrake?
Reinstating the firedrake mailing list cc (please do keep it cc'd since others may have more useful comments, or wish to know the answer).
On 15 Sep 2020, at 10:57, Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> wrote:
Thanks a lot for your prompt response and help.
The inner(grad(sol),n)*v in last two terms give me an error: UFL:ERROR Shapes do not match: <Grad id=4986788992> and <Zero id=4998074176>.
But the projection was okay when changed them to sol*v. Do you know where is this error coming from? (FYI: my sol Function is defined on DG1 functionspace)
I always mix up the definitions of jump(foo, n). I actually prefer just to write out what you want explicitly in terms of the restricted values.
I think I was not clear here. My code works when define F as: F = (fluxes*w)('+')*dS + fluxes*w*ds - avg(sol)*w('+')*dS - sol*w*ds‘
But when change it to a desired flux that contains grad(sol) such as below: F = (fluxes*w)('+')*dS + fluxes*w*ds - avg(inner(grad(sol), n))*w('+')*dS - inner(grad(sol), n)*w*ds
It returns error even when changing to grad(sol) only on ds term. This is the error I get: UFL:ERROR Shapes do not match: <Grad id=4986788992> and <Zero id=4998074176>.
Somehow one of the two terms was simplified to zero (perhaps n is zero)? This from firedrake import * mesh = UnitSquareMesh(1, 1) V = FunctionSpace(mesh, "DG", 1) Vt = FunctionSpace(mesh, "HDiv Trace", 0) v = TestFunction(Vt) sol = Function(V) n = FacetNormal(mesh) inner(grad(sol), n)*v*ds assemble(inner(grad(sol), n)*v*ds) works for me, but perhaps I am doing something different. Thanks, Lawrence
On Sep 15, 2020, at 5:03 AM, Lawrence Mitchell <wence@gmx.li> wrote: .
I think I was not clear here. My code works when define F as: F = (fluxes*w)('+')*dS + fluxes*w*ds - avg(sol)*w('+')*dS - sol*w*ds‘
But when change it to a desired flux that contains grad(sol) such as below: F = (fluxes*w)('+')*dS + fluxes*w*ds - avg(inner(grad(sol), n))*w('+')*dS - inner(grad(sol), n)*w*ds
It returns error even when changing to grad(sol) only on ds term. This is the error I get: UFL:ERROR Shapes do not match: <Grad id=4986788992> and <Zero id=4998074176>.
Somehow one of the two terms was simplified to zero (perhaps n is zero)?
This
from firedrake import * mesh = UnitSquareMesh(1, 1) V = FunctionSpace(mesh, "DG", 1) Vt = FunctionSpace(mesh, "HDiv Trace", 0) v = TestFunction(Vt) sol = Function(V) n = FacetNormal(mesh) inner(grad(sol), n)*v*ds assemble(inner(grad(sol), n)*v*ds)
works for me, but perhaps I am doing something different. That is right my n had been simplified to zero by mistake. It works now.
Thank you Lawrence. May I expand on my previous question? You earlier showed me a workaround for integration on the interface facets.
We can't directly integrate like this. You can do a whole-field operation to get the fluxes into a DG0 trace space:
Vt = FunctionSpace(mesh, "HDiv Trace", 0)
fluxes = Function(Vt)
# you have to code this trace projection by hand:
v = TestFunction(Vt)
F = (fluxes*v)('+')*dS + fluxes*v*ds - avg(inner(grad(sol), n))*v('+')*dS - inner(grad(sol), n))*v*ds solve(F == 0, fluxes)
# Now fluxes contains your integrals.
Explanation: The goal is to implement a flux-limiting algorithm which is cell-based and: i) For cell E, It requires to know the flux integral H_{E,E’} on the facets of cell E and the adjacent cell E’ (e.g., H_{E,E’} = int_dS avg(inner(grad(sol),n)) + jump(sol) ) . ii) for cell E, it assumes that n(‘+’) is always unit outward normal from the cell. The consequence of (ii) is that for two adjacent cells E and E’, H_{E,E’} = -H_{E’,E} (flux goes from E to E’ is equal opposite sign to flux goes from E’ to E). Projecting fluxes integration to DG0 trace space gave us the value H at the interface and using Vt. cell_node_map().values we know that this H belongs to which two adjacent cells. However, we don’t know if H=H_{E,E’} or H=-H_{E,E’}. This is because (I assume) in firedrake it is arbitrary that which side of the facet is chosen positive and which side is chosen negative. To solve this sign issue at the interface of E and E’, I need to know if the side of cell E is take (+), if yes then H=H_{E_E’}. Otherwise I should reverse the sign of H. Question: Is it possible to somehow extract that at the interface of two cells, which of E or E’ side is taken as (+)? Is there a way of defining the FacetNormal on internal facets such that n(+) for all facets point from left to right (i.e., n(+)[0] be positive )? In other words, change +,- signs and make them regular in domain. (FYI: my mesh is triangular) Thanks, Sarraf
On 17 Sep 2020, at 16:55, Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> wrote:
[...]
We can't directly integrate like this. You can do a whole-field operation to get the fluxes into a DG0 trace space:
Vt = FunctionSpace(mesh, "HDiv Trace", 0)
fluxes = Function(Vt)
# you have to code this trace projection by hand:
v = TestFunction(Vt)
F = (fluxes*v)('+')*dS + fluxes*v*ds - avg(inner(grad(sol), n))*v('+')*dS - inner(grad(sol), n))*v*ds solve(F == 0, fluxes)
# Now fluxes contains your integrals.
Explanation: The goal is to implement a flux-limiting algorithm which is cell-based and:
i) For cell E, It requires to know the flux integral H_{E,E’} on the facets of cell E and the adjacent cell E’ (e.g., H_{E,E’} = int_dS avg(inner(grad(sol),n)) + jump(sol) ) .
ii) for cell E, it assumes that n(‘+’) is always unit outward normal from the cell.
The consequence of (ii) is that for two adjacent cells E and E’, H_{E,E’} = -H_{E’,E} (flux goes from E to E’ is equal opposite sign to flux goes from E’ to E).
Projecting fluxes integration to DG0 trace space gave us the value H at the interface and using Vt.
The flux you get with this procedure is single-valued on the edge.
cell_node_map().values we know that this H belongs to which two adjacent cells. However, we don’t know if H=H_{E,E’} or H=-H_{E,E’}. This is because (I assume) in firedrake it is arbitrary that which side of the facet is chosen positive and which side is chosen negative.
To solve this sign issue at the interface of E and E’, I need to know if the side of cell E is take (+), if yes then H=H_{E_E’}. Otherwise I should reverse the sign of H.
Question: Is it possible to somehow extract that at the interface of two cells, which of E or E’ side is taken as (+)?
Is there a way of defining the FacetNormal on internal facets such that n(+) for all facets point from left to right (i.e., n(+)[0] be positive )? In other words, change +,- signs and make them regular in domain. (FYI: my mesh is triangular)
You can't do this. If you need to do some kind of upwinding the usual trick is to pick based on 0.5*(abs(u.n) + u.n) Which is u.n on the upwind side and zero on the downwind side. I guess therefore you want to include a sign-flip based on the direction of n? I don't quite see how to do that. Lawrence
Thank you for your prompt reply.
On Sep 17, 2020, at 11:02 AM, Lawrence Mitchell <wence@gmx.li> wrote: …
Explanation: The goal is to implement a flux-limiting algorithm which is cell-based and:
i) For cell E, It requires to know the flux integral H_{E,E’} on the facets of cell E and the adjacent cell E’ (e.g., H_{E,E’} = int_dS avg(inner(grad(sol),n)) + jump(sol) ) .
ii) for cell E, it assumes that n(‘+’) is always unit outward normal from the cell.
The consequence of (ii) is that for two adjacent cells E and E’, H_{E,E’} = -H_{E’,E} (flux goes from E to E’ is equal opposite sign to flux goes from E’ to E).
Projecting fluxes integration to DG0 trace space gave us the value H at the interface and using Vt.
The flux you get with this procedure is single-valued on the edge.
cell_node_map().values we know that this H belongs to which two adjacent cells. However, we don’t know if H=H_{E,E’} or H=-H_{E,E’}. This is because (I assume) in firedrake it is arbitrary that which side of the facet is chosen positive and which side is chosen negative.
To solve this sign issue at the interface of E and E’, I need to know if the side of cell E is take (+), if yes then H=H_{E_E’}. Otherwise I should reverse the sign of H.
Question: Is it possible to somehow extract that at the interface of two cells, which of E or E’ side is taken as (+)?
Is there a way of defining the FacetNormal on internal facets such that n(+) for all facets point from left to right (i.e., n(+)[0] be positive )? In other words, change +,- signs and make them regular in domain. (FYI: my mesh is triangular)
You can't do this.
If you need to do some kind of upwinding the usual trick is to pick based on
0.5*(abs(u.n) + u.n)
Which is u.n on the upwind side and zero on the downwind side.
I guess therefore you want to include a sign-flip based on the direction of n? I don't quite see how to do that.
Yes, exactly. I want to do a sign-flip operation on H (the single-valued fluxes at the edge) based on the direction of n(+). I think it is achievable by first accessing the layout of + and - cells at edges in domain. I wonder if there is any direct function to access this layout mapping? If it is not directly possible, I think this layout can be calculated in two steps: I) project inner(Constant((1,0)) * n(+) ) into DG0 trace space. Get direction of n(+) at the interface. For example, If projection is positive, we know that n(+)[0] is pointing from left to right. II) compare the x-component of mid-points of two cells E and E' adjacent to the interface. If cell E midpoint is less than cell E’ midpoint then we know that cell E is positive and cell E’ is negative. Is there a way to extract cell midpoints? Thanks, Sarraf
participants (2)
-
Lawrence Mitchell
-
Mohammad Sarraf Joshaghani