python tutorial on multiphysics?

[Guosheng Fu]

Hello all,

I am wondering if there is python tutorials available on multiphysics problems.
Say, Stokes-Darcy coupling.
I am interested to see how to define/use finite elements on subdomains.

Best,
Guosheng

[Joachim Schoeberl]

Hi Guosheng,

I added a chapter on that in the docu:
HOW-TO ---> Spaces and forms on sub-domains, see
https://ngsolve.org/docu/nightly/how_to/howto_definedon.html

The Stokes-Darcy coupling I would do by the H(div)-conforming HDG for Stokes, and just dropping the viscosity-term (and facet term) in the Darcy domain.

Best, Joachim

[Guosheng Fu]

Hi Joachim,

See. This is very good.

Now, I have trouble imposing interface conditions.

Here I am solving a two-Domain poisson equation with a interface. In the spirit of Nitsche-xfem, we defined a piecewise-H1 space and impose interface condition weakly.

The FE space is a compound one:

V1 = H1(mesh, order=order, dirichlet=[1,3,4], definedon=[1])
V2 = H1(mesh, order=order, dirichlet=[5,6,7], definedon=[2])
V = FESpace([V1,V2])

the final numerical solution shall be u = u1+u2 with u1 only defined on subdomain 1, and u2 on subdomain 2.

Then, we implement the Nitsche-interface term

u1,u2 = V.TrialFunction()
v1,v2 = V.TestFunction()
n = specialcf.normal(2)
h = specialcf.mesh_size

a = BilinearForm(V)
a += SymbolicBFI (grad(u1)*grad(v1),definedon=[1])
a += SymbolicBFI (grad(u2)*grad(v2),definedon=[2])
# nitsche term
cf1 = 0.5*(grad(u1)+grad(u1.Other()) + grad(u2)+grad(u2.Other())) \
    *n*(v1-v1.Other()+v2-v2.Other())
cf2 = 0.5*(grad(v1)+grad(v1.Other()) + grad(v2)+grad(v2.Other())) \
    *n*(u1-u1.Other()+u2-u2.Other())
cf3 = 2*order*order/h*(v1-v1.Other()+v2-v2.Other())*(u1-u1.Other()+u2-u2.Other())
a += SymbolicBFI (-cf1-cf2+cf3, BND, skeleton=True, definedon=[2])

Here the average and jump terms take into accound contributions from each subdomain.
But the interface term is not implemented correctly. I don't know what when wrong.
The code is attached.

Best,
Guosheng

[schruste]

Dear Guosheng,

The skeleton-BFI together with the BND flag is only implemented for domain boundaries. It expects that only one aligned volume element exists. In the implementation it gets a list of neighbors (here two) and only takes the first one (it does not expect a second in the list).

The case with two neighbors is a "regular" skeleton-BFI with VOL flag. Your case, where you only want to treat jumps across an interface is not treated straight-forwardly, yet.

Attached is a version with a cheap work-around which works:

We make use of the (only recently provided) SetDefinedOnElements method which marks only some facets where the integrator should be applied (similar to the definedon-flag - difference is that the decision is not done domain-wise, but element/facet-wise).
In the attached file I do the following: In a "hacked way" I determine all facets in the mesh which are on the interface and store this in a BitArray (see output).
Later, I use the standard SymbolicBFI(…, skeleton=True, VOL) and combine it with the SetDefinedOnElements-method where I enter the preprocessed BitArray. Now, the skeleton-integrator does what you expected it to do.

Note also that I changed your CompoundFESpace to
V = FESpace([V1,V2], flags={"dgjumps":True})
because we need to add the Nitsche-couplings to the matrix.

Explanation:
In NGSolve, so far, the sparsity pattern of a matrix is done before the assembly. The sparsity pattern is based on overlapping supports. For the DG/Nitsche-terms we have additional couplings (non-zero entries in the matrix). To also consider possible couplings through facets we need to provide the "dgjumps"-flag.
Note that this introduces much more nonzero entries than necessary in your case. You would only need these additional couplings along the interface. The dgjumps-flag introduces them through every facet.
More clever ways could be:

• Change the way the sparsity pattern of the matrix is set up (requires some coding)

• Use an additional facet unknown for the coupling! (FacetFESpace with definedonbound=...)

I also cleaned up the Nitsche formulation a bit...

Best,
Christoph

[Guosheng Fu]

Hi Christoph,

Thanks for the code!

Now I have two follow-up questions.
1, When we do SymbolicBFI with VOL and skeleton=True, there will be two elements sharing the facet, and two normal directions. What's the convention in ngsolve to choose the normal direction?
Based on your Nitsche term, I guess the convention might be choose the direction so that its first argument is positive? Is it correct?
It seems I asked you the same question before... but I forgot the answer

2, Concerning your suggestion on using FacetFESpace on the interface, I got stuct with defining the element.
I defined a facet fespace V3, with flag definedonbound = [2], where 2 is the number of the interface, but the FreeDofs of the facet fe space is still on all the facets...

For the lagrange multiplier, I used the ip-hdg form with interface unknown as u restricted on it. In the form, the coupling on the interface shall be (V1<-->V3) and (V2<-->V3), but I still got sparsity patten error without using dgjumps flag.

Best,
Guosheng

[schruste]

Hey Guosheng,

There is always one pronounced direction. I am actually not sure what the policy is for that. I think it is always from the element with the lower element number to the higher one. So it essentially depends on the mesh generator. I am not sure if there is a rule, but in your example the elements on the left always have a smaller number then those on the right. But I have to pass that question to Joachim or another Netgen-expert, sorry.

Concerning the Mortaring with the FacetFESpace:
Your version still requires dgjumps because you use skeleton=True integrals. Those will assemble Facet-Matrices which involve all the d.o.f.s of the neighboring elements. If some of them are zero or not is not respected. NGSolve assumes the worst case, i.e. that the Facet matrix is full. In this case the sparsity pattern does not fit unless you have dgjumps activated. So, as it is implemented right now, you cannot use skeleton=True (+VOL) bilinearform integrator if you assemble a matrix without the dgjumps flag. The way around it, is to write the interface integrals as two one-sided integrals (element-to-facet) and (element-to-facet) to circumvent (element-to-element).

You are right that also this version is not a standard case for NGSolve, but all the tools that you need to hack it in are there.
Again, there are two ways:

• You want to have one integrator that integrate only on the domain interface. That requires new integrals (change in bilinearform.cpp). It's not too hard but requires some insight and some C++ coding.

• There is a faster and simpler approach (based only on python), which is however less efficient and less elegant.

I wil explain the second approach that you should be able to try yourself.

1. Use the VOL-integrals with the element_boundary=True flag and assemble the terms in the usual HDG way. However, introduce a facet-indicator function (similar to the dom_ind_average in my last example) so that only the interface facets are really used. Formally this way you require the facetspace everywhere, but we accept this for now. This way (which is not too hard to implement) you should be able to get the right integrals.
   Of course the integrator will run over all elements and assemble many zeros for this part, but that is the price we accept to pay in this approach.

2.
Now, in terms of the linear systems you should mark only the facet d.o.f.s that are on the interface as FreeDofs (all others should have zero entry in the respective columns and rows anyway). To this end you have to find out which d.o.f.s are at the interface. You can do this in different ways. The simplest way is that you count the facets (w.g. number of d.o.f.s in the indicator-facet-space) and use the information from your indicator function which tells you which facet is an interface facet. Then, you can calculate the d.o.f.s that are on the interface (note that you will have the lowest order d.o.f.s first (number of facet) and then k (in 2D) dofs per facet afterwards.

Hope this helps,
Christoph