Is it possible to calculate a flux value on the edges of a cell in Firedrake?
Hello, I appreciate if you help me with this question. I am trying to implement a postprocessing flux-limiting scheme for a DG advection problem to eliminate the under/overshoots that appears in the vicinity of the front. In order to achieve this, I need to iterate over each mesh element and then on each facet of the element calculate the value of numerical flux. Based on the sign of this Flux value, a constant edge-specific alpha is calculated and the new corrected solutions are then reconstructed for the cell. In firedrake, Is there any way to iterate over facets of each cell, and calculate an integral of an expression (e.g, avg(inner(grad(sol,n))) ) on that facet? This is how I achieved it in fenics: ----------------------------------------------------- facet_domains = MeshFunction('size_t',mesh,mesh.topology().dim()-1) dS = Measure('dS', subdomain_data = facet_domains) elnum = 0 for cell in cells(mesh): flux=np.zeros(3) for fct in facets(cell): flux[fct.index] = assemble( (avg(inner(grad(sol) , n) ) )*dS(fct.index()) print(flux) . . do some operations based on flux[fct.index] values . . elnum += 1 -------------------------------------------------------- Thanks, Sarraf
Hi Sarraf, apologies for the low responses. I have some comments inline below.
On 14 Sep 2020, at 08:01, Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> wrote:
Hello,
I appreciate if you help me with this question.
I am trying to implement a postprocessing flux-limiting scheme for a DG advection problem to eliminate the under/overshoots that appears in the vicinity of the front. In order to achieve this, I need to iterate over each mesh element and then on each facet of the element calculate the value of numerical flux. Based on the sign of this Flux value, a constant edge-specific alpha is calculated and the new corrected solutions are then reconstructed for the cell. In firedrake, Is there any way to iterate over facets of each cell, and calculate an integral of an expression (e.g, avg(inner(grad(sol,n))) ) on that facet?
This is how I achieved it in fenics: ----------------------------------------------------- facet_domains = MeshFunction('size_t',mesh,mesh.topology().dim()-1) dS = Measure('dS', subdomain_data = facet_domains) elnum = 0 for cell in cells(mesh): flux=np.zeros(3) for fct in facets(cell): flux[fct.index] = assemble( (avg(inner(grad(sol) , n) ) )*dS(fct.index()) print(flux) . . do some operations based on flux[fct.index] values . .
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. Lawrence
Dear Lawrence, Thanks for your help.
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
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)
solve(F == 0, fluxes)
# Now fluxes contains your integrals.
Lawrence
Also, how this globally defined fluxes are connected to cells? In other words for a triangular mesh, is it possible to know which 3 components of fluxes.vector().array() are related to cell E? Thank you, Sarraf
On 15 Sep 2020, at 09:18, Mohammad Sarraf Joshaghani <m.sarraf.j@gmail.com> wrote:
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.
solve(F == 0, fluxes)
# Now fluxes contains your integrals.
Lawrence
Also, how this globally defined fluxes are connected to cells? In other words for a triangular mesh, is it possible to know which 3 components of fluxes.vector().array() are related to cell E?
Vt.cell_node_map().values contains for each cell the three fluxes.
participants (2)
-
Lawrence Mitchell
-
Mohammad Sarraf Joshaghani