BC for Blockmatrix in time-loop

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, BUT: I use MixedElement() in the following!

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):

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

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
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 for examples of manipulating block systems.

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.
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:

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 and the Stokes 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.

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:

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.

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

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?

I found what I needed.

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