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