Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

A bug on boundary condition for i-tutorial unit 3.3.1, and related questions

More
3 years 9 months ago #2973 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
Code:
u.n=0
. The numerical flux on the boundary is
Code:
phat = p uhat.n = 0
(Here maybe make a note that

p.Other()

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 )

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
More
3 years 9 months ago #2977 by Guosheng Fu
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)
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.108 seconds