- Thank you received: 6
Question on interface integrals at material boundaries in DG formulation
- Guosheng Fu
- Topic Author
- Offline
- Elite Member
Less
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:
I am using the bilinear form element_boundary formulation
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
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))
Code:
a = -uhat*v*dx(element_boundary=True, definedon=mesh.Materials("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:
5 years 3 months ago #1735
by joachim
Replied by joachim on topic Question on interface integrals at material boundaries in DG formulation
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:
Joachim
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
- Guosheng Fu
- Topic Author
- Offline
- Elite Member
Less
More
- Thank you received: 6
5 years 3 months ago #1741
by Guosheng Fu
Replied by Guosheng Fu on topic Question on interface integrals at material boundaries in DG formulation
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
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
Attachments:
5 years 3 months ago #1742
by schruste
Replied by schruste on topic Question on interface integrals at material boundaries in DG formulation
Dear Guosheng,
You simply missed to restrict your element_boundary integrals to the interior domain.
Best,
Christoph
You simply missed to restrict your element_boundary integrals to the interior domain.
Best,
Christoph
The following user(s) said Thank You: Guosheng Fu
5 years 3 months ago #1743
by joachim
Replied by joachim on topic Question on interface integrals at material boundaries in DG formulation
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
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
- Guosheng Fu
- Topic Author
- Offline
- Elite Member
Less
More
- Thank you received: 6
5 years 3 months ago #1745
by Guosheng Fu
Replied by Guosheng Fu on topic Question on interface integrals at material boundaries in DG formulation
I see...
I thought the following two lines are mathematically equivalent:
and
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
This works great. But I would like to change SymbolicBFI to dx/ds language, so, I use
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
I thought the following two lines are mathematically equivalent:
Code:
inner_func*(uhat*r*n)*dx(element_boundary=True)
Code:
(uhat*r*n)*dx(element_boundary=True, definedon=mesh.Materials("inner"))
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"))
Code:
a +=uhat*vhat*dx(skeleton=True, definedon=mesh.Boundaries("b2"))
Thanks in advance. Your timely response really helped me a lot.
Best,
Guosheng
Time to create page: 0.109 seconds