### Assembly of a block system

78

views

2

Hello, I am having some problems in assembling a system with a block matrix and block vector. To make things simple let us say I want to solve the heat equations

\begin{align}

v_t = \Delta v, \quad x\in\Omega\\

w_t = d\Delta w, \quad x\in\Omega\\

\end{align}

and boundary conditions

\begin{align}

- n\cdot\nabla v = w, \quad x\in\partial\Omega\\

- n\cdot\nabla w = v, \quad x\in\partial\Omega\\

\end{align}

This leads to the weak formulation

\begin{align}

\int_{\Omega} v_t\phi dx + \int_{\Omega} \nabla v\cdot\nabla\phi dx + \int_{\partial\Omega} w \phi ds = 0\\

\int_{\Omega} w_t\psi dx + d\int_{\Omega} \nabla w\cdot\nabla\psi dx + d\int_{\partial\Omega} v \psi ds = 0\\

\end{align}

Since I want to use an implicit method for the time discretisation, for example

\begin{align}

\frac{1}{\Delta t}\int_{\Omega} v^n\phi dx + \int_{\Omega} \nabla v^n\cdot\nabla\phi dx + \int_{\partial\Omega} w^n \phi ds = \frac{1}{\Delta t}\int_{\Omega} v^{n-1}\phi dx \\

\frac{1}{\Delta t}\int_{\Omega} w^n\psi dx + d\int_{\Omega} \nabla w^n\cdot\nabla\psi dx + d\int_{\partial\Omega} v^n \psi ds = \frac{1}{\Delta t}\int_{\Omega} w^{n-1}\psi dx

\end{align}

I would like to assembly a single system. I am trying to use cbc.block, however I did not understand very well how to do it. I copy the code below:

This returns the error:

What am I doing wrong?

Another question is how to include Dirichlet conditions, for example, if boundary conditions would be as:

\begin{align*}

v = v_0, \quad x\in\Gamma_D\\

w = w_0, \quad x\in\Gamma_D\\

- n\cdot\nabla v = w, \quad x\in\Gamma_N\\

- n\cdot\nabla w = v, \quad x\in\Gamma_N\\

\end{align*}

I would insert the following line

Hope some of you could help!

Thank you

Davide

\begin{align}

v_t = \Delta v, \quad x\in\Omega\\

w_t = d\Delta w, \quad x\in\Omega\\

\end{align}

and boundary conditions

\begin{align}

- n\cdot\nabla v = w, \quad x\in\partial\Omega\\

- n\cdot\nabla w = v, \quad x\in\partial\Omega\\

\end{align}

This leads to the weak formulation

\begin{align}

\int_{\Omega} v_t\phi dx + \int_{\Omega} \nabla v\cdot\nabla\phi dx + \int_{\partial\Omega} w \phi ds = 0\\

\int_{\Omega} w_t\psi dx + d\int_{\Omega} \nabla w\cdot\nabla\psi dx + d\int_{\partial\Omega} v \psi ds = 0\\

\end{align}

Since I want to use an implicit method for the time discretisation, for example

\begin{align}

\frac{1}{\Delta t}\int_{\Omega} v^n\phi dx + \int_{\Omega} \nabla v^n\cdot\nabla\phi dx + \int_{\partial\Omega} w^n \phi ds = \frac{1}{\Delta t}\int_{\Omega} v^{n-1}\phi dx \\

\frac{1}{\Delta t}\int_{\Omega} w^n\psi dx + d\int_{\Omega} \nabla w^n\cdot\nabla\psi dx + d\int_{\partial\Omega} v^n \psi ds = \frac{1}{\Delta t}\int_{\Omega} w^{n-1}\psi dx

\end{align}

I would like to assembly a single system. I am trying to use cbc.block, however I did not understand very well how to do it. I copy the code below:

```
from fenics import *
import numpy as np
from dolfin import *
mesh = RectangleMesh(Point(0.0,0.0), Point(5,7.5), 10, 15)
tol = 1E-14
n = FacetNormal(mesh)
V = FunctionSpace(mesh, 'P', 1)
W = FunctionSpace(mesh, 'P', 1)
v = TrialFunction(V)
w = TrialFunction(W)
phi = TestFunction(V)
psi = TestFunction(W)
Mv = assemble(v*phi*dx)
Mw = assemble(w*psi*dx)
Kv = assemble(inner(grad(v),grad(phi))*dx)
Kw = assemble(inner(grad(w),grad(psi))*dx)
MGammav = assemble(w*phi*ds)
MGammaw = assemble(v*psi*ds)
v_n = interpolate(Constant(2.0), V)
w_n = interpolate(Constant(1.0), W)
t = 0
Dt = 0.1
from block import *
H = block_mat([[(1/Dt)*Mv + Kv, MGammav], [MGammaw, (1/Dt)*Mw + Kw]])
v = Function(V)
w = Function(W)
while t<10 :
t += Dt
sol = block_vec([v, w])
B = block_vec([(1/Dt)*Mv*v_n.vector(), (1/Dt)*Mw*w_n.vector()])
solve(H, sol, B)
v = sol[0]
w = sol[1]
v_n.assign(v)
w_n.assign(w)
```

This returns the error:

*Traceback (most recent call last):**File "test.py", line 51, in <module>**solve(H, sol, B)**File "/usr/local/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 310, in solve**return cpp.la_solve(*args)**File "/usr/local/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 4898, in la_solve**return _la.la_solve(*args)**TypeError: in method 'la_solve', argument 1 of type 'dolfin::GenericLinearOperator const &'**Aborted*What am I doing wrong?

Another question is how to include Dirichlet conditions, for example, if boundary conditions would be as:

\begin{align*}

v = v_0, \quad x\in\Gamma_D\\

w = w_0, \quad x\in\Gamma_D\\

- n\cdot\nabla v = w, \quad x\in\Gamma_N\\

- n\cdot\nabla w = v, \quad x\in\Gamma_N\\

\end{align*}

I would insert the following line

```
bcV = DirichletBC(V, v0, boundary_markers, 1)
bcW = DirichletBC(W, w0, boundary_markers, 1)
```

but how to apply these to the system? Would this work?```
bc_V.apply(H,B)
bc_W.apply(H,B)
```

Hope some of you could help!

Thank you

Davide

Community: FEniCS Project

Please login to add an answer/comment or follow this question.

Thus I think maybe we need to transform block matrix to petsc matrix.

You can find it at here:

https://fenicsproject.org/qa/10752/cbc-block-matrices-to-petsc-matrices/

and

https://www.allanswered.com/post/bevkr/

However, I am not very clear about what's the next step.