- Thank you received: 6
A bug on boundary condition for i-tutorial unit 3.3.1, and related questions
- Guosheng Fu
- Topic Author
- Offline
- Elite Member
Less
More
4 years 4 months ago #2973
by Guosheng Fu
A bug on boundary condition for i-tutorial unit 3.3.1, and related questions was created by Guosheng Fu
Hello,
I am looking into traceoperator/geom_free implementation of DG.
My starting point is this i-tutorial ngsolve.org/docu/latest/i-tutorials/unit...hlight=traceoperator
In version 1, you are using the wall boundary condition
. The numerical flux on the boundary is
(Here maybe make a note that
By increasing final time tend = 0.3, we find the boundary condition is indeed correctly implemented.
Then in Version 2, phat is the average on interior edges, but phat=p on boundary edges. Hence, to have a correct implementation of the boundary term Btr in block 11, I added the following hack to have the correct boundary condition:
With this correction, the version 2 results are identical to version 1 up to time tend=0.3. (Otherwise, we observe something wrong when wave hit the boundary)
Unfortunately, the same trick does not work for the geom_free version.
Do you have any suggestions to fix the boundary condition for Version 3?
Another question: does traceoperator and/or geom_free work on a periodic mesh?
I have some trouble getting consistent results on periodic mesh for the traceoperator.
Best,
Guosheng
I am looking into traceoperator/geom_free implementation of DG.
My starting point is this i-tutorial ngsolve.org/docu/latest/i-tutorials/unit...hlight=traceoperator
In version 1, you are using the wall boundary condition
Code:
u.n=0
Code:
phat = p
uhat.n = 0
returns p, not zero on the boundary edge, as it is sometimes confusing what the magic Other() operator is doing on domain boundary. See this discussion ngsolve.org/forum/ngspy-forum/28-inconsi...-in-bnd-integral#146 )p.Other()
By increasing final time tend = 0.3, we find the boundary condition is indeed correctly implemented.
Then in Version 2, phat is the average on interior edges, but phat=p on boundary edges. Hence, to have a correct implementation of the boundary term Btr in block 11, I added the following hack to have the correct boundary condition:
Code:
# A HACK
gf = GridFunction(FacetFESpace(mesh, dirichlet=".*"))
gf.Set(1, BND)
Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u)
Btr += 0.5 * (1+gf)*phat * (v*n) * dx(element_boundary=True)
With this correction, the version 2 results are identical to version 1 up to time tend=0.3. (Otherwise, we observe something wrong when wave hit the boundary)
Unfortunately, the same trick does not work for the geom_free version.
Do you have any suggestions to fix the boundary condition for Version 3?
Another question: does traceoperator and/or geom_free work on a periodic mesh?
I have some trouble getting consistent results on periodic mesh for the traceoperator.
Best,
Guosheng
- Guosheng Fu
- Topic Author
- Offline
- Elite Member
Less
More
- Thank you received: 6
4 years 4 months ago #2977
by Guosheng Fu
Replied by Guosheng Fu on topic A bug on boundary condition for i-tutorial unit 3.3.1, and related questions
I like the HDG-like implementation of an explicit DG scheme.
The traceop can be explicitly realized via a matrix-vector form (although the performance is a bit slower than traceop, but it handles the boundary condition correctly)
But the geom_free flag still does not work for a periodic facet fespace: the following operator shall allow for geom_free implementation, but it fails for a periodic facet fespace VF.
The traceop can be explicitly realized via a matrix-vector form (although the performance is a bit slower than traceop, but it handles the boundary condition correctly)
Code:
# test/trial functions
u,v = V.TnT()
uhat, vhat = VF.TnT()
### mass matrix for uhat (diagonal)
Mtr = BilinearForm(VF)
Mtr += uhat*vhat*dx(element_boundary=True)
Mtr.Assemble()
invuhat = Mtr.mat.CreateSmoother(VF.FreeDofs())
### u-->uhat (uhat=average of u)
Bhat = BilinearForm(trialspace=V, testspace=VF)
Bhat += u*vhat*dx(element_boundary=True)
Bhat.Assemble()
traceop = invuhat @ Bhat.mat
But the geom_free flag still does not work for a periodic facet fespace: the following operator shall allow for geom_free implementation, but it fails for a periodic facet fespace VF.
Code:
Btr = BilinearForm(trialspace=VF, testspace=W, geom_free=True)
Btr += uhat * (r*n) * dx(element_boundary=True)
Btr.Assemble()
Time to create page: 0.100 seconds