# BC for Blockmatrix in time-loop

**URL:** <https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631>\
**Category:** dolfinx\
**Created:** [November 3, 2023, 8:33am UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631 "2023-11-03T08:33:42Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\
**Post date:** [November 3, 2023, 8:33am UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/1 "2023-11-03T08:33:42Z")

</div>

Dear Community,

i am trying to implement the following coupled, scalar, “directed” (by b\in\mathbb{R}^d) advection/transport system with Dirichlet BC on both unknowns:

\begin{aligned} \partial\_t V +b\cdot \nabla \phi &= 0\\ \partial\_t \phi + b\cdot \nabla V &= 0 \end{aligned}

The weak form for the problem on H^1\times H^1 is

\begin{aligned} \left \langle \partial\_t V ,W \right \rangle- \left \langle \phi b, \nabla W \right \rangle &= - \left ( \phi\_\Gamma b\cdot e\_n , W \right )\_\Gamma \quad \forall W\in H^1\\ \left \langle \partial\_t \phi, \psi \right \rangle - \left \langle Vb, \nabla \psi \right \rangle &= - \left ( V\_\Gamma b\cdot e\_n, \psi \right )\_\Gamma \quad \forall \psi\in H^1 \end{aligned}

Or fully discrete, with implicit Euler:

\begin{aligned} \begin{pmatrix} M & E \\ F & N \end{pmatrix} \begin{pmatrix} \vec V ^{n+1}\\ \vec\phi^{n+1} \end{pmatrix} = \begin{pmatrix} M \vec V^n - \vec r^{n+1} \\ F\vec V^n -\vec s^{n+1} \end{pmatrix} \end{aligned}

with:

\begin{aligned} M\_{ij} &= \left \langle W\_j, W\_i \right \rangle , \qquad &&N\_{ij} = \left \langle \psi\_j,\psi\_i \right \rangle \\ E\_{ij} &= \left \langle \psi\_j b, \nabla W\_i \right \rangle , \qquad &&F\_{ij} = \left \langle W\_j b,\nabla \psi\_i\right \rangle \\ r\_i^{n+1} &= \left ( \phi\_{\Gamma,j}^{n+1} b \cdot e\_n, W\_i \right )\_\Gamma , \qquad && s\_i^{n+1} = \left ( V\_{\Gamma,j}^{n+1} b\cdot e\_n, \psi\_i \right )\_\Gamma \end{aligned}

I am following the [Navier-Stokes example](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_navier-stokes.html), BUT: I use `MixedElement()` in the following!

```auto
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import plot
from dolfinx.fem import (Constant, Function, dirichletbc, form, 
                         functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block,assemble_vector, create_vector_block,assemble_vector_block, create_vector
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.mesh import *
from ufl import div, dx, ds, grad, inner, FacetNormal, dot, 
import pyvista
from mpi4py import MPI
from petsc4py import PETSc

# parameters
dt = 0.1
t = 0
T_end = 1*dt
num_steps = int(T_end/dt)

# mesh, spaces, functions
msh = create_unit_square(MPI.COMM_WORLD, 20,20)
x = ufl.SpatialCoordinate(msh)
n = FacetNormal(msh)
P1_elem = element("Lagrange", "triangle", 1) 
P1_mixed = ufl.MixedElement(P1_elem,P1_elem)
X = functionspace(msh, P1_mixed) # mixed space P1 x P1
XV = X.sub(0).collapse()[0]
Xp = X.sub(1).collapse()[0]
W,q = ufl.TestFunctions(X)
V,p = ufl.TrialFunctions(X)
V_old = Function(XV)
V_new = Function(XV)
p_old = Function(Xp)
p_new = Function(Xp)
bfield = Constant(msh, (+1.0,0.0))

# init cond
def initcond_V(x):
    return (x[0]-0.6)**2 + (x[1]-0.3)**2 < 0.1**2
def initcond_p(x):
    return (x[0]-0.2)**2 + (x[1]-0.7)**2 < 0.1**2
V_old.interpolate(initcond_V)
p_old.interpolate(initcond_p)

# BC
facets = locate_entities_boundary(msh, dim=1,marker=lambda x: np.logical_or.reduce((
                                  np.isclose(x[0], 0.0),
                                  np.isclose(x[0], 1.0),
                                  np.isclose(x[1], 0.0),
                                  np.isclose(x[1], 1.0))))
V_dofs = locate_dofs_topological((X.sub(0), XV), 1, facets)
p_dofs = locate_dofs_topological((X.sub(1), Xp), 1, facets)
def inflow_BC_expression(x):
    return np.zeros(X.dofmap.index_map.size_global)
BC_V = Function(XV)
BC_V.interpolate(inflow_BC_expression)
BC_p = Function(Xp)
BC_p.interpolate(inflow_BC_expression)

BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)] #?

# weak form implicit euler
M = inner(V,W)*dx
E = inner(bfield*p,grad(W))*dx
F = inner(V*bfield,grad(q))*dx
N = inner(p,q)*dx
a = form([[M, -dt/2*E],
          [-dt/2*F, N]])

# assembly
A = assemble_matrix_block(a, bcs=BCs)
A.assemble()
Vp_vec = A.createVecRight()

# solver
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt

    # Update the right hand side
    l = form([inner(grad(V_old), grad(W)) * dx,
            inner(grad(p_old), grad(q)) * dx])
    L = assemble_vector_block(l,a,bcs=BCs)
    
    # Solve linear problem
    solver.solve(L, Vp_vec)

    # extract solution
    offset = XV.dofmap.index_map.size_global
    V_new.x.array[:] = Vp_vec.array_r[:offset]
    V_new.x.scatter_forward()
    p_new.x.array[:(len(Vp_vec.array_r) - offset)] = Vp_vec.array_r[offset:]
    p_new.x.scatter_forward()

    # Update solution at previous time step
    V_old.x.array = V_new.x.array
    p_old.x.array = p_new.x.array

```

However, it outputs the following error for the “apply lifting”-command, (after tracing back to dolfinx/fem/assemble.py):

```bash
Traceback (most recent call last):
  File ".../transport-system-mixed.py", line 157, in <module>
    p_new.x.array[:] = Vp_vec.array_r[offset:]
    ~~~~~~~~~~~~~ ^^^
ValueError: could not broadcast input array from shape (1323,) into shape (441,)

```

I have questions concerning the following:

- Why does `createVecRight()` create an array that is much too large? (also, shouldnt it be “createVecLeft()”?)
- How do I assemble the matrix correctly by taking into account Dirichlet BC on both unknowns? (I think I cant just assemble the submatrices individually. E.g. assembling submatrix M with BC would not enforce the essential BC on the essential DOFs correctly, as there is submatrix E messing up the respective equations in the resulting system.) I assume the command `BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)]` is wrong. Can i somehow concatenate the Dirichlet values before?
- When do i need `apply_lifting(L, [a], [BCs])` or `set_bc(L,[BCs])` on the right hand side vector?

Thank you very much for your answers!  
Best regards

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [November 3, 2023, 3:03pm UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/2 "2023-11-03T15:03:53Z")

</div>

There are lots of oddities here which don’t make sense to me. Salient points:

1. Your weak formulation only makes sense if \nabla \cdot b = 0.
2. Your discretisation may be unstable with the scheme you’ve chosen.
3. Your boundary conditions are likely ilposed where b \cdot e\_n \> 0.
4. You’ve collapsed the space `X` in to `XV` and `Xp` but you’re still trying to treat it as a monolithic blocked formulation

```python
Xp = X.sub(1).collapse()[0]
...
p_new = Function(Xp)
...
p_new.x.array[:(len(Vp_vec.array_r) - offset)] = Vp_vec.array_r[offset:]

```

Perhaps consider the [Stokes demo](https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_stokes.py) for examples of manipulating block systems.

---

<div class="post-metadata">

**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\
**Post date:** [November 3, 2023, 5:09pm UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/3 "2023-11-03T17:09:12Z")

</div>

Thank you, I highly appreciate your quick answer!

I purposefully opted for these BC, IC and b-vector, for the sake of constructing an [MWE](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21).  
You are right, the problem is convection dominated, which should be considered in the choice of parameters and time discretization. The b-vector is indeed divergence free, which I should have indicated.

Your last point is very insightful for me. Is there a rule of thumb to decide whether to go for a Mixed Element implementation, or one with separate Elements and a “monolithic” system matrix? My corrected code for reference:

```auto
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import plot
from dolfinx.fem import (Constant, Function, dirichletbc, form, 
                         functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block,assemble_vector, create_vector_block,assemble_vector_block, create_vector
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from dolfinx.mesh import *
from ufl import div, dx, ds, grad, inner, FacetNormal, dot, 
import pyvista
from mpi4py import MPI
from petsc4py import PETSc

# parameters
dt = 0.1
t = 0
T_end = 20*dt
num_steps = int(T_end/dt)

# mesh, spaces, functions
msh = create_unit_square(MPI.COMM_WORLD, 20,20)
x = ufl.SpatialCoordinate(msh)
n = FacetNormal(msh)
P1 = element("Lagrange", "triangle", 1) 
XV = functionspace(msh, P1) 
Xp = functionspace(msh, P1) 
V = ufl.TrialFunction(XV)
p = ufl.TrialFunction(Xp)
W = ufl.TestFunction(XV)
q = ufl.TestFunction(Xp)
V_old = Function(XV)
V_new = Function(XV)
p_old = Function(Xp)
p_new = Function(Xp)
bfield = Constant(msh, (0.1,0.0))

# init cond
def initcond_V(x):
    return (x[0]-0.5)**2 + (x[1]-0.5)**2 < 0.1**2
def initcond_p(x):
    return (x[0]-0.2)**2 + (x[1]-0.7)**2 < 0.1**2
V_old.interpolate(initcond_V)
p_old.interpolate(initcond_p)

# BC
facets = locate_entities_boundary(msh, dim=1,marker=lambda x: np.logical_or.reduce((
                                  np.isclose(x[0], 0.0),
                                  np.isclose(x[0], 1.0),
                                  np.isclose(x[1], 0.0),
                                  np.isclose(x[1], 1.0))))
V_dofs = locate_dofs_topological(XV, 1, facets)
p_dofs = locate_dofs_topological(Xp, 1, facets)
def BC_expressionV(x): return np.zeros(XV.dofmap.index_map.size_global)
def BC_expressionp(x): return np.zeros(Xp.dofmap.index_map.size_global)
BC_V = Function(XV)
BC_V.interpolate(BC_expressionV)
BC_p = Function(Xp)
BC_p.interpolate(BC_expressionp)

BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)] #?

# weak form implicit euler
M = inner(V,W)*dx
E = inner(bfield*p,grad(W))*dx
F = inner(V*bfield,grad(q))*dx
N = inner(p,q)*dx
a = form([[M, -dt/2*E],
          [-dt/2*F, N]])

# assembly
A = assemble_matrix_block(a, bcs=BCs)
A.assemble()
Vp_vec = A.createVecLeft()

# solver
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt

    # Update the right hand side
    l = form([inner(grad(V_old), grad(W)) * dx,
            inner(grad(p_old), grad(q)) * dx])
    L = assemble_vector_block(l,a,bcs=BCs)

    # Solve linear problem
    solver.solve(L, Vp_vec)

    # extract solution
    offset = XV.dofmap.index_map.size_global
    V_new.x.array[:] = Vp_vec.array_r[:offset]
    V_new.x.scatter_forward()
    p_new.x.array[:] = Vp_vec.array_r[offset:]
    p_new.x.scatter_forward()

    # Update solution at previous time step
    V_old.x.array[:] = V_new.x.array
    p_old.x.array[:] = p_new.x.array

```

In the [Navier-Stokes](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_navier-stokes.html) and the [Stokes](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_stokes.html) demo examples, I could see the assembly of the system matrix - but only when having Dirichlet BC solely on the first unknown. How can I achieve this for the case with BC on both unknowns? The method I am using above, `BCs = [dirichletbc(BC_V, V_dofs), dirichletbc(BC_p, p_dofs)]`, seems to be wrong to me.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [November 3, 2023, 6:50pm UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/4 "2023-11-03T18:50:30Z")

</div>

Typically you’d choose the matrix structure based on the ease with which you can implement the preconditioner you intend to employ, along with any memory and preformance requirements to which you may be subjected.

Since your demo indicates you want to use standard LU factorisation for the whole system, a pure mixed element matrix is likely most efficient to exploit sparsity patterns with minimal bandwidth.

If you wanted to precondition each block individually then a monolithic block matrix structure would allow you to construct the appropriate preconditioning matrices for each block individually, and apply appropriate (inner) solvers if you require. It would also be very straightforward to solve the whole system simply using LU factorisation if required.

If you wanted to minimise memory use even further you could employ the block nested structure; however, this requires more care in the design of an appropriate iterative solver with preconditioners as standard LU for the whole system are not available.

In each case consider the following lines for setting BCs:

- Mixed element assembly:  
[https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo\_stokes.py#L459-L477](https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_stokes.py#L459-L477)
- Blocked assemby:  
[https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo\_stokes.py#L133-L152](https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_stokes.py#L133-L152)

As a side note: You can technically implement a block preconditioner for a mixed element matrix; however keeping tack of the index sets when using PETSc to solve the system can be a real pain.

---

<div class="post-metadata">

**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\
**Post date:** [November 4, 2023, 9:30am UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/5 "2023-11-04T09:30:03Z")

</div>

Dear nate, thank you for your elaborate answer.

I still have troubles understanding how to generalize the code from the Stokes example you have linked to a case where BC have to be set on both unknowns.

In the Stokes example, we only have BC on the velocity field, not on the pressure. We can therefore define the no-slip boundary and the lid boundary, define Dirichlet BC and collect them in a list

```auto
bcs = [bc0, bc1]

```

My question is: Am I allowed to use this syntax to collect Boundary conditions for different unknowns? E.g., could bc0 be a condition for the velocity, and bc1 for the pressure?

---

<div class="post-metadata">

**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\
**Post date:** [November 5, 2023, 10:55am UTC](https://fenicsproject.discourse.group/t/bc-for-blockmatrix-in-time-loop/12631/6 "2023-11-05T10:55:47Z")

</div>

I found what I needed.

As an example, where Dirichlet BC are set on multiple unknowns, the reader may refer to

- [This discourse article](https://fenicsproject.discourse.group/t/navier-stokes-block-iterative-solver-with-pressure-bc/10607)
- [This subsection of Jørgen Dokken’s tutorial](https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html)