Question: What is the general approach to assemble differnet PDEs over multiple domains?

[willi-z]

Hello everyone,

I want to model a multiphysic problem where I have different PDEs over multiple subdomains.
A simple case to describe this would be the following to equations:
$$\nabla \cdot \left( K_1 \nabla u_1 \right) = f_1 , x \in \Omega_1$$
and
$$\nabla \cdot \left( K_2 \nabla ln(u_2) \right) + K_1 = f_2 , x \in \Omega_2$$
with
$$\frac{\partial u_1}{\partial n} = 0 $$ and $$\frac{\partial u_2}{\partial n} = 0$$ on the interface between $\Omega_1$ and $\Omega_2$.
The equations above are just two simple functions which are representive for some way more complicated functions. Therefor I did not check if the equtions above actually make sense (My actual problem is closer to a concentration diffusion in a solid and a liquid phase of charged particle in a electric field).
After reading the tutorial and trying the examples I am still unsure how to aproach such a problem. I thought about some thing like:

import skfem as fem


basis = fem.Basis(mesh, element)

@fem.BilinearForm
def a_1(u,v,w):
    return K_1 * dot(grad(u), grad(v))

@fem.BilinearForm
def a_2(u,v,w):
    return K_1 * dot(grad(ln(u)), grad(v))

domain1, domain2 = split_domain(basis)

a = a1.assemble(domain1) + a2.assemble(domain2)

But I was not able to find a function which could help me with split_domain (i only found FacetBasis what seems to do something in this regard).
I would be very gratefull if someone could explain how to approach such a problem in scikit-fem.

[kinnala]

In order to solve a PDE on a subdomain (part of the mesh) you would need to eliminate the degrees-of-freedom that are outside of the subdomain and initialize Basis using elements keyword argument to restrict numerical integration to the subdomain. Here is an example which does exactly that: https://github.com/kinnala/scikit-fem/blob/master/docs/examples/ex26.py

If you want to solve different PDEs on different subdomains then you do similar restrictions for all subdomains. If you want to couple the different PDEs at the interfaces of the subdomains then you first need to decide what type of coupling approach you want to use: iterative or monolithic. Perhaps the easiest to start with is iterative coupling where you solve one PDE at a time and then take its solution to be a boundary condition for the other PDEs, and iterate different PDEs until convergence.

Now the type of the interface condition between the subdomains will decide how the boundary data is applied. If your coupling is Dirichlet type, i.e. something like u_1 = u_2, then you'd simply apply Dirichlet conditions at the interface using condense. If your coupling is something like Neumann type, i.e. the normal derivative of u_1 is affected by u_2, then you'd get some additional terms in the weak formulation and in order to assemble terms at the interfaces between subdomains you can use InteriorFacetBasis.

Your current example is not coupled in any way so I cannot comment further what you should do. I also suggest you start with linear PDEs because it's easier to understand the coupling approach.

[willi-z]

Thank you so much!
I will definitly try this and post my results here.
Until then I hope it is ok, if I leave the question marked unsanswered in case further questions or problems arises.