Question on interface integrals at material boundaries in DG formulation

Guosheng_Fu · July 25, 2019, 6:02pm

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:uhat = IfPos(n, u, u.Other(bnd=1))
I am using the bilinear form element_boundary formulation

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

Attachment: interface.py

joachim · July 25, 2019, 9:36pm

Hi Guosheng,

element-boundary DG terms on region boundaries are skipped, see

github.com

NGSolve/ngsolve/blob/master/comp/bilinearform.cpp

#include <comp.hpp>
#include <multigrid.hpp>

#include <parallelngs.hpp>
#include "../fem/h1lofe.hpp"

namespace ngcomp
{
  
  // dummy function header 
  void CalcEigenSystem (FlatMatrix<Complex> & elmat, 
                        FlatVector<Complex> & lami, 
                        FlatMatrix<Complex> & evecs)
  { ; }

  void MinusMultAB (SliceMatrix<Complex> a, SliceMatrix<Complex> b, SliceMatrix<Complex> c)
  {
    c = -a * b | Lapack;
  }
  void MinusMultABt (SliceMatrix<Complex> a, SliceMatrix<Complex> b, SliceMatrix<Complex> c)

This file has been truncated. show original

comments in line 4024.

As a workaround, you can multiply with the characteristic function:

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 · July 28, 2019, 11:12pm

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

Attachment: testTwoDomain.py

schruste · July 29, 2019, 5:31am

Dear Guosheng,

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

Best,
Christoph

joachim · July 29, 2019, 5:44am

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

Guosheng_Fu · July 29, 2019, 1:02pm

I see…
I thought the following two lines are mathematically equivalent:

inner_func*(uhat*r*n)*dx(element_boundary=True)

and

(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 https://ngsolve.org/forum/ngspy-forum/218-integral-on-interface-between-two-domains#750
I can do

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

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

joachim · July 29, 2019, 1:13pm

The first version produces 0 element-matrices, and tries to assemble them. But the global matrix has no entries allocated, and thus throws an exception.

The translate to dx/ds syntax should be

a +=uhat.Trace()*vhat.Trace()*ds(definedon=mesh.Boundaries("b2"))

We want a loop over boundary elements, no need for volume elements and thus no skeleton flag.

Guosheng_Fu · July 29, 2019, 1:19pm

Thanks again for the clarification.
Previously, I used ds(skeleton=True,…), it gives me an error.
Now, I tried the same thing in 1D, and it seems .Trace() is not properly implemented yet for 1D Facet Dofs. It gives me this error message:

netgen.libngpy._meshing.NgException: FacetFESpace::GetSFE: unsupported element Point

Guosheng_Fu · July 29, 2019, 1:23pm

Here is the code produce the error message

Attachment: facetfes.py