# How to get mass matrix of composite FE spaces?

Hello, I need to study problems related to fluid free surfaces based on the level set method.

My ideas are as follows:

we consider the following two PDE problems at the same time:

Find (u,p):[0,T] \to (H_{0,D}^1)^d \times L^2, s.t.

\begin{align} \int_{\Omega} \partial_t u \cdot v + \int_{\Omega} \nu(c) \nabla u \nabla v + u \cdot \nabla u v - \int_{\Omega} \operatorname{div}(v) p &= \int f v && \forall v \in (H_{0,D}^1)^d, \\ - \int_{\Omega} \operatorname{div}(u) q &= 0 && \forall q \in L^2, \\ \quad u(t=0) & = u_0 \end{align}

with \nu(c) depends on the solution c of the unsteady scalar linear transport problem, which is the level set equation of free surface

Find c: [0,T] \to V_D := \{ c \in L^2(\Omega), u \cdot \nabla c \in L^2(\Omega), c|_{\Gamma_{in}} = c_D\}, s.t.

\int_{\Omega} \partial_t c d + u \cdot \nabla c d = \int_{\Omega} f d \qquad \forall d \in V_0 = \{ c \in L^2(\Omega), u \cdot \nabla c \in L^2(\Omega), c|_{\Gamma_{in}} = 0\}

We use a Taylor-Hood discretization of degree k. The finite element space for the (vectorial) velocity \mathbf{u} and the pressure p are:

\begin{align} \mathbf{u} \in \mathbf{V} &= \{ v \in H^1(\Omega) | v|_T \in \mathcal{P}^k(T) \}^2 \\ p \in Q &= \{ q \in H^1(\Omega) | q|_T \in \mathcal{P}^{k-1}(T) \} \end{align}

Stokes problem for initial values

Find \mathbf{u} \in \mathbf{V}, p \in Q so that

\begin{align} \int_{\Omega} \nu(c) \nabla \mathbf{u} : \nabla \mathbf{v} - \int_{\Omega} \operatorname{div}(\mathbf{v}) p &= \int \mathbf{f} \cdot \mathbf{v} && \forall \mathbf{v} \in \mathbf{V}, \\- \int_{\Omega} \operatorname{div}(\mathbf{u}) q &= 0 && \forall q \in Q, \end{align}

with boundary conditions:

• the inflow boundary data (‘inlet’) with mean value 1
• “do-nothing” boundary conditions on the ouflow (‘outlet’) and
• homogeneous Dirichlet conditions on all other boundaries (‘wall’).

To solve the Stokes (and later the Navier-Stokes) problem, we introduce the bilinear form:

a((u,p),(v,q)) := \int_{\Omega} \nu(c) \nabla u : \nabla v - \operatorname{div}(v) p - \operatorname{div}(u) q

IMEX time discretization for the unsteady Navier Stokes equations

For the time integration we consider an semi-implicit Euler method (IMEX) where the convection is treated only explicitly and the remaining part implicitly:

Find (\mathbf{u}^{n+1},p^{n+1}) \in X = \mathbf{V} \times Q, s.t. for all (\mathbf{v},q) \in X = \mathbf{V} \times Q

\begin{align*} \underbrace{m(\mathbf{u}^{n+1},\mathbf{v}) + \Delta t ~\cdot~a((\mathbf{u}^{n+1},p^{n+1}),(\mathbf{v},q))}_{ \to M^\ast} ~=~m(\mathbf{u}^{n},\mathbf{v}) - \Delta t ~\cdot~c(\mathbf{u}^{n}; \mathbf{u}^{n},\mathbf{v}) \end{align*}

with

\begin{align*} m(\mathbf{u},\mathbf{v}) = \int_{\Omega} \mathbf{u} \cdot \mathbf{v} \end{align*}

and

\begin{align*} c(\mathbf{w},\mathbf{u},\mathbf{v}) = \int_{\Omega} \mathbf{w} \cdot \nabla \mathbf{u} \cdot \mathbf{v} \end{align*}

We prefer the incremental form (as it homogenizes the boundary conditions):

\begin{align*} m(\delta \mathbf{u}^{n+1},\mathbf{v}) + \Delta t ~\cdot~a((\delta \mathbf{u}^{n+1},\delta p^{n+1}),(\mathbf{v},q)) ~=~ \Delta t (-a((\mathbf{u}^{n},p^n),(\mathbf{v},q)) -c(\mathbf{u}^{n}; \mathbf{u}^{n},\mathbf{v})) \end{align*}

Explicit time stepping with a DG formulation for the unsteady scalar linear transport problem

We consider an Upwind DG discretization (in space):

Find c: [0,T] \to V_h := \bigoplus_{T\in\mathcal{T}_h} \mathcal{P}^k(T) so that

\sum_{T} \int_T \partial_t c d - u \cdot \nabla d c + \int_{\partial T} u_n \hat{c} d = \int_{\Omega} f d, \quad \forall d \in V_h.

Here \hat{c} is the Upwind flux, i.e. \hat{c} = c on the outflow boundary part \partial T_{out} = \{ x\in \partial T \mid u(x) \cdot n_T(x) \geq 0 \} of T and \hat{c} = c^{other} else, with c^{other} the value from the neighboring element.

Upwind-in-space combined with an explicit Euler scheme in time yields:

\sum_{T} \int_T c^{n+1} d = \sum_{T} \int_T c^{n} d - \Delta t \sum_{T} \left\{ - \int_T u \cdot \nabla d c + \int_{\partial T} u_n \hat{c} d \right\} + \Delta t \int_{\Omega} f d, \quad \forall d \in V_h,
M c^{n+1} = M c^{n} - \Delta t C(c^n) + \Delta t f

In our first example we set f = 0.

\begin{align*} {m({c}^{n+1},{d})} ~=~m({c}^{n},{d}) - \Delta t ~\cdot~c(\mathbf{u}^{n}; {c}^{n},{d}) \end{align*}

with

\begin{align*} m({c},{d}) = \int_{T} {c} \cdot {d} \end{align*}

and

\begin{align*} c(\mathbf{e}; {c},{d}) = - \int_{T} \mathbf{e} \cdot \nabla {d} \cdot {c} + \int_{\partial T} e_n \hat{c} d \end{align*}

We prefer the incremental form (as it homogenizes the boundary conditions):

\begin{align*} {m({\delta c}^{n+1},{d})} ~=~ - \Delta t ~\cdot~c({c}^{n}; {c}^{n},{d}) \end{align*}

combined form

m((\delta \mathbf{u}^{n+1},\delta c^{n+1}),(\mathbf{v},d)) + \Delta t ~\cdot~a((\delta \mathbf{u}^{n+1},\delta p^{n+1}),(\mathbf{v},q)) ~=~ \Delta t (-a((\mathbf{u}^{n},p^n),(\mathbf{v},q)) -c(\mathbf{u}^{n}; (\mathbf{u}^{n},c^{n}),(\mathbf{v},d)))

with :

m((\mathbf{u},c),(\mathbf{v},d)) = \int_{\Omega} \mathbf{u} \cdot \mathbf{v} + \int_{T} {c} \cdot {d}
a((u,p),(v,q)) := \int_{\Omega} \nu \nabla u : \nabla v - \operatorname{div}(v) p - \operatorname{div}(u) q
c(\mathbf{u}; (\mathbf{u},c),(\mathbf{v}.d))= \int_{\Omega} \mathbf{u} \cdot \nabla \mathbf{u} \cdot \mathbf{v} - \int_{T} \mathbf{u} \cdot \nabla {d} \cdot {c} + \int_{\partial T} u_n \hat{c} d

but when i use the following code, i find that the mass matrix can not be calculated in the above way shown in the formula, can somebody give me some advice?

k = 2
# V = VectorH1(mesh,order=k, dirichlet="Walls|In")
V = VectorH1(mesh,order=k, dirichlet="Walls|In",dirichlety="Bottoms")
Q = H1(mesh,order=k-1)
D = L2(mesh,order=k)
X = V*Q*D

(u,p,c), (v,q,d) = X.TnT()

a = BilinearForm(X, nonassemble = True)
nu = IfPos(c,500,0.02)
a += stokes

dt = 1e-5
mstar = BilinearForm(X)

# D.Mass(1) since use the DG
mstar += InnerProduct(u,v)*dx + dt*stokes + D.Mass(1)


you get the mass matrix for the last component as you get it for the other terms:

c*d*dx


What you tried to do,

D.Mass(1)


gives you a linear operator, not an integral - you cannot add these different types

Joachim

Thanks, but when I tried the above method, I found that I couldn’t get the correct solution. I decided to decouple the equation first, and the problem I encountered now is that the NS equation was calculated incorrectly using the Crank-Nicolson scheme
N_S-demo cn.ipynb (135.2 KB)