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