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