Question on interface integrals at material boundaries in DG formulation

More
5 years 3 months ago #1732 by Guosheng Fu
Hello,

I would like to perform DG calculation on a subdomain of a two-material mesh, but I have big trouble understanding the numerical flux, in particular the .Other() operator.

Below is a minimal test to illustrate my confusion.
I have a 1D domain with 4 elements, left half of them are labelled mat1, and right half are labelled mat2.
So, I want to solve u_t + u_x = 0 only on mat1 using DG-P0 formulation.
the numerical flux is taken to be upwinding:
Code:
uhat = IfPos(n, u, u.Other(bnd=1))
I am using the bilinear form element_boundary formulation
Code:
a = -uhat*v*dx(element_boundary=True, definedon=mesh.Materials("mat1"))
But the code fails to produce a constant 1 solution on mat1.
Actually, when I apply the vector u=[1,1,0,0] to the bilinear form, I would assume to get a zero-vector, but instead, I get w=[0,1,0,0].
It is clear that the bilinear form forgot to integrate the boundary term at the interface shared by mat1 and mat2: for test function v=1 at the second cell, the bilinear form returns [u0] instead of [u0-u1] as my numerical flux suggest. Do you know what's going on?

P.S. I also tried to use skeleton-based formulation, but I don't know how to treat the interface properly...

Best,
Guosheng
Attachments:
More
5 years 3 months ago #1735 by joachim
Hi Guosheng,

element-boundary DG terms on region boundaries are skipped, see
github.com/NGSolve/ngsolve/blob/master/comp/bilinearform.cpp
comments in line 4024.

As a workaround, you can multiply with the characteristic function:
Code:
mat1func = CoefficientFunction( [1 if mat=="mat1" else 0 for mat in mesh.GetMaterials() ] ) a += -mat1func*uhat*v*n*dx(element_boundary=True) # , definedon=mesh.Materials("mat1"))

Joachim
More
5 years 3 months ago #1741 by Guosheng Fu
Hi Joachim,

Thank you for the clarification. It would be better if we can include the material interface...

Here is another issue I encounter:
I want to do HDG on a subdomain, with static condensation on.
Attached is the code, and here is the set-up:
I have two subdomains, with label "inner" and "outer".
Then, I have three finite element spaces, the vector DG space fes_q, and the scalar DG space fes_u, both live on the whole domain, and another hybrid FESpace live only on the skeleton of subdomain "inner". (I don't want to use HDG on "outer" subdomain)
Then, I solve the eqn
-lap u + u = 1 on domain "inner" with zero flux bc on the material interface.
The solution shall be u = 1 on iinner" and zero on "outer".
To solve the problem, I have to manually set those DG dofs on subdomain "outer" as not FreeDOfs. (line 40-45)
Iin the attached code, when condensation is set off, the results are correst, but when condensation is set on, the solver produce nan values for those dofs in subdomain "outer". I noticed that line 59 on a.inner_solve create the issue. Do you know what's going on here?

Best,
Guosheng
More
5 years 3 months ago #1742 by schruste
Dear Guosheng,

You simply missed to restrict your element_boundary integrals to the interior domain.

Best,
Christoph
The following user(s) said Thank You: Guosheng Fu
More
5 years 3 months ago #1743 by joachim
my answer above applies only to DG, not to HDG.
HDG is much nicer since we can assemble element by element, and don't need the neighbouring elements at the same time as in DG. The neighbouring elements outside the Region were the issue in your first post.

Joachim
More
5 years 3 months ago #1745 by Guosheng Fu
I see...
I thought the following two lines are mathematically equivalent:
Code:
inner_func*(uhat*r*n)*dx(element_boundary=True)
and
Code:
(uhat*r*n)*dx(element_boundary=True, definedon=mesh.Materials("inner"))
But the first gives a buggy result, while the second produce the correct result. Can you explain the hidden difference?

Now, I have another (minor) question on the subject:
I want to add the bilinear form <uhat, vhat> on interface edges, following the post ngsolve.org/forum/ngspy-forum/218-integr...ween-two-domains#750
I can do
Code:
a += SymbolicBFI(uhat.Trace()*vhat.Trace(), definedon=mesh.Boundaries("b2"))
This works great. But I would like to change SymbolicBFI to dx/ds language, so, I use
Code:
a +=uhat*vhat*dx(skeleton=True, definedon=mesh.Boundaries("b2"))
The code runs, but the above two lines produce different results. Again, I thought these two lines are identical... So, what's the correct equivalent of the above SymbolicBFI in dx/ds language?

Thanks in advance. Your timely response really helped me a lot.

Best,
Guosheng
Time to create page: 0.109 seconds