How to apply boundary conditions to left-hand matrix

Errors
1

Hi, I want to replace the “dolfin” with the latest “dolfinx” in my old project, but I met a problem: the result from “fenics.assemble_system” is different compared with “dolfinx.fem.assemble_matrix” in dolfinx. I find that the result from “dolfinx.fem.assemble_matrix” is equal to the result of “fenics.assemble_system” without boundary conditions. So how can I apply the boundary conditions to the left-hand side matrix?

looking forward to anyone’s help.
Best regard.

2

See: Application of Dirichlet boundary conditions — FEniCS Tutorial @ Sorbonne
as lifting is equivalent to assemble_system in legacy dolfin (allthough it is not properly documented in the all code).

1 Like
3

I am wondering what is difference between A.zeroRowsLocal and cpp.fem.petsc.insert_diagonal? It seem to result in a more correct result with A.zeroRowsLocal.

I see you do the following to manually apply boundary conditions for the LHS matrix A, where zero the selected row in local order of matrix and leave diagnoal 1:

for bc in bcs:
    dofs, _ = bc._cpp_object.dof_indices()
    A.zeroRowsLocal(dofs, diag=1)

I was directed to see the following insert_diagonal when calling dolfinx.fem.petsc.assemble_matrix. I guess it is also used to apply bc with changing entities and diagonal:

if a.function_spaces[0] is a.function_spaces[1]:
    A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
    A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
    _cpp.fem.petsc.insert_diagonal(A, a.function_spaces[0], _bcs, diagonal)
4

These operations do different things.
ZeroRowsLocal zeros a set of rows given by indices local to the process, and optionally sets a diagonal value.
Insert_diagonal changes the diagonal, But doesn’t Zero of diagonal entries.
Im not dure what you mean by one giving a more correct result than the other.

That does assembly with BCs as input, so rows and columns are already zeroed, as illustrated by the lifting/symmetric application of bcs

5

Thanks for the explanation.

My problem is that I did not know how to manually apply bc for A correctly. In my applications, I need to first get the matrix A and then modify A with boundary conditions applied. So my code looks like

# 1. Get A and b
A = assemble_matrix(a, bcs)
A.assemble()
b = assemble_vector(L)
apply_lifting(b, [a], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) 
set_bc(b, self.bcs)

# 2. Code to modify A and b, e.g: 
A += delta_A
b += delta_b

# 3. Apply bc for A and b again
_bcs = [bc._cpp_object for bc in bcs]
if a.function_spaces[0] is a.function_spaces[1]:
	_cpp.fem.petsc.insert_diagonal(A, a.function_spaces[0], _bcs, diagonal=1.0)
A.assemble()
set_bc(b, bcs)

It does not give the desired values on the boundary, for instance, a non-zero velocity boundary value with no-slip bc. However, if I replace the above Part 3 by the following code, the result shows a possibly correct value on boundary:

# 3. Apply bc for A and b again
for bc in self.bcs:
	dofs, _ = bc._cpp_object.dof_indices()
	A.zeroRowsLocal(dofs, diag=1)
A.assemble()
set_bc(b, self.bcs)
6

Is it only accomplished by

cpp.fem.petsc.insert_diagonal(A, a.function_spaces[0], _bcs, diagonal)

or some other functions?

7

The zeroing happens when you add bcs to assemble_matrix.

The tutorial I have referenced above:

shows how the bcs are applied correctly to both system.

As you have not provided a minimal reproducible example, I can be of no further help.