Assembly of a block system


78
views
2
20 days ago by
Davide  
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:
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
Hello, I think I also met the same problem as yours. If you check the type of your H and b, it should be "block.block_mat.block_mat". Can this type be input to the solve function? I highly suspect it.
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.
written 16 days ago by SICHENG SUN  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »