Question about dgjumps in ngsxfem

NGSolve
1

This question is directed at C. Lehrenfeld.

I was modifying the jupyter notebook MoL_TraceFEM. The original notebook had the bilinear and linear forms being redefined within the time-stepping loop (using SymbolicBFI).

I rewrote it to define the forms once outside the loop and using integrators with dCut and dx on a narrow band. There is no facet based ghost penalty.

This seems to work, but I had to modify the finite element space VhG to include: dgjumps=True. Otherwise, I get an error. But the original code did not require this.

Why is that? I’m guessing it is because of the cut cell integrators. But the other explanations for using dgjumps=True is when you have coupled DoFs across facets (which I don’t here).

-Shawn

2

Hi Shawn,

You are right, if you don’t use ghost penalties, you shouldn’t need the dgjumps=True flag. Can you send me an example?

Best,
Christoph

3

MoL_TraceFEM_example.ipynb (17.6 KB)

This is a modification of a code that you made (I assume). You can run it, then remove the dgjumps=True flag and run it again.

4

Sorry for the late reply.

The problem is a different one.
You change the active domain in every time step. That changes the sparsity pattern of the matrix. In order to respect that you need to reset the Matrix-Graph and reallocate the matrix in every time step, i.e. you need to use a.Assemble(reallocate=True) instead of a.Assemble().
That should suffice to make your code work.

When using dgjumps=True together with no facet_restrictions all facet couplings are reserved in the matrix including all the element-couplings in the active domain which you need. That’s why dgjumps=True also works, but increases the computational effort significantly.

Best,
Christoph

5

So you are saying that if I remove the dgjumps=True, but use a.Assemble(reallocate=True), then it should work? I tried that before, but it did not work. I just checked again, and I get this error:

SparseMatrixTM::AddElementMatrix: illegal dnumsin Assemble BilinearForm 'blf'

I understand about the domain changing, but the original demo example also had a domain changing. Do I need to change one of the options in:

a = RestrictedBilinearForm

What should facet_restriction and check_unused be set to?

-Shawn

6

Ok, I figured it out. I just need to replace RestrictedBilinearForm with BilinearForm. I guess because I don’t have any facet integrals to compute. In the while-loop, the active DoFs are passed to the inverse of a, so it has everything it needs.

So, RestrictedBilinearForm should only be used if there are facet integrals.

Best,

-Shawn

7

Dear Shawn,

Well, you can use the BilinearForm of course. This will reserve the entries for the coupling of all background elements without consideration of the current active element/facet dofs.

The idea of the RestrictedBilinearForm is to avoid the setup of the matrix for the whole domain. Instead you can prescribe the elements and facets (in the case of additional coupling through facets) that should only be taken into account when reserving entries for the matrix. That’s what the corresponding flags are for.

The problem with your script is with the UpdateBand-method. It computes the active elements which are considered for the setup of the matrix pattern. However, the resulting BitArray does not even include the current cut elements. There are some elements marked in ba_new_IF-BitArray which are not part of the ba_band_IF-BitArray. Then, during assembly entries for these cut elements are computed and when trying to write to the matrix this fails as no entries have been reserved for that beforehand.

Using the standard BilinearForm circumvents that issue (at the cost of a larger matrix (costs are only increased for storing, inverting is not more expensive as you only invert on the active_dofs anyway)). However, the origin of the problem is with the marker BitArrays and this described problem from above also implies that you are not applying the normal derivative stabilization on some cut elements which is probably not good for robustness…

Best, Christoph

8

Ok, I’ll look at it a bit more. But I based this off of a demo code that is in ngsxfem. I didn’t change anything with UpdateBand. The original code was also a moving surface problem. But maybe I accidentally changed something.

-Shawn

9

The original demo still runs through (no dgjumps, using RestrictedBLF), but it also creates the RestrictedBLF in every time step (which is less elegant). Not sure where exactly the issue emerged, but you certainly changed a few details.

Best,
Christoph

10

I see the problem with your script now. You create a new CutInfo-Objekt in every time step. Each new CutInfo has its own BitArray so that the BitArrays that are used for the integrals and for the matrices are outdated. I will upload a corrected version of your script in a few minutes.

Best,
Christoph

11

Here are the correction for your script:

  • lset_approx defined before use
  • dgjumps=False
  • CutInfo is not re-created within time loop, but only updated (so that we always work on the same BitArrays)
  • a.Assemble() → a.Assemble(reallocate=True)

This yields the attached file:
MoL_TraceFEM_example_corrected.ipynb (15.5 KB)

Best,
CL

12

Thanks for this. I had made some of these corrections already. But it was still not working because I had done this within the time-loop:

ci.Update(lset_approx)
ba_new_IF.Clear()
ba_new_IF |= ci.GetElementsOfType(IF)

whereas you did this:

ci.Update(lset_approx)
ba_new_IF = ci.GetElementsOfType(IF)

I guess when Clear() is called it completely resets the pointer to that array. So, then it gets orphaned.

Thanks!

-Shawn