How to use DiagonalPolicy to set essential BC in mixed Blocksystem

[markusrenoldner]

Dear community,

we are currently implementing a Dirichlet problem that involves multiple function spaces and results in a blockmatrixsystem. We did not understand how FormLinearSystem() worked, so we implemented a small example and found this:

Suppose i want to solve the following LSE:

$$a\cdot x = b$$

$$ \begin{pmatrix} 1 & 1 &0 &0 &0 \\\ 1 & 1 & 0& 0& 0\\\ 0& 0 &1&0 &0 \\\ 0& 0 & 0& 0&1 \\\ 1& 0& 0& 1&0 \end{pmatrix} \cdot \begin{pmatrix} x_0\\\ x_1\\\ x_2\\\ 5\\\ x_4 \end{pmatrix} = \begin{pmatrix} 11\\\ 11\\\ 12\\\ 13\\\ 14 \end{pmatrix} $$

so the matrix $a$ has four independent row/column vectors, and the unknown vector $x$ has one value that is fixed. The solution of this system is $x_0=9$, $x_1=2$, $x_2=12$, and $x_4=13$.
Now I call

a.FormLinearSysem(essdof, x, b, A, X, B)

where essdof = {3}. Then I print the resulting matrix $A$, the vector $X$ and right hand side $B$ and get

$$ A\cdot X = B$$

$$ \begin{pmatrix} 1 & 1 &0 &0 &0 \\\ 1 & 1 & 0& 0& 0\\\ 0& 0 &1&0 &0 \\\ 0& 0 & 0& 1&0 \\\ 1& 0& 0& 0&0 \end{pmatrix} \cdot \begin{pmatrix} x_0\\\ x_1\\\ x_2\\\ 5\\\ x_4 \end{pmatrix} = \begin{pmatrix} 11 \\\ 11\\\ 12\\\ 5\\\ 9 \end{pmatrix} $$

The solution is now $x_0=9$, $x_1=2$, and $x_2=12$. But the value of $x_4$ is undetermined. Apparently, there is something called

mfem::Operator::DiagonalPolicy

that influences how the command FormLinearSysem() treats certain entries of the system matrix. We dont understand this behavior and similar we dont know whether and how to use this DiagonalPolicy. Our hypothesis is that the zeros diagonal elements that frequently appear in blocksystems are not handled like we expect.
Any help is much appreciated.

[mlstowell]

Hello, @markusrenoldner,

Unfortunately, I think this may be a poorly chosen test matrix. If you make it symmetric by adding a '1' in the upper right entry your example should work fine.

When applying Dirichlet constraints in a FEM context we remove the rows and columns of the matrix which correspond to Dirichlet entries and we modify the right hand side vector so that the solution has the desired Dirichlet values in the proper locations. The DiagonalPolicy controls how the diagonal entries of the matrix are replaced on these Dirichlet rows and columns. The standard way to do this is to place a '1' on the diagonal and place the desired Dirichlet value in the right hand side (this is the default which corresponds to a DiagonalPolicy of DIAG_ONE). If the matrix entries are very different from 1 the this has the potential to create problems with some solvers so we offer the choice to keep the original diagonal entry and set b_i = a_{ii}x_i (this is DiagonalPolicy of DIAG_KEEP). The third option is DIAG_ZERO. This is a special purpose choice which zeros the rows and columns and leaves a zero on the diagonal. This third choice is useful for preparing off-diagonal matrices in block systems.

[markusrenoldner]

Ah I see thank you I will check that out.

Our FEM problem is something like the following: Given
$u|_{\partial\Omega}$ and $w|_{\partial\Omega}$ as well as
$u^0$ and $w^0$
find $u\in$ X and $w \in$ Y st.
$$m \cdot x = b$$

$$\iff \begin{bmatrix}A & B \\ C& 0\end{bmatrix} \cdot \begin{bmatrix} u^{k+1} \\ w^{k+1} \end{bmatrix} = \begin{bmatrix} A u^{k} + B w^{k} +f\\ 0 \end{bmatrix}$$

$A,B,C$ are SparseMatrix objects, $m$ is a BlockOperator, $u$ and $w$ are Gridfunctions, the unknown vector $x$ is a Vector created from concatenating the gridfunctions. $f$ is a LinearForm and the rest of the right hand side is created from computing the matrixproduct with the gridfunctions. The ess_dof array is created by calling GetBoundaryTrueDofs() on both FEM spaces and then concatenating both resulting arrays like this:

mfem::Array<int> essdof;
mfem::Array<int> essdof2;
X.GetBoundaryTrueDofs(essdof); 
Y.GetBoundaryTrueDofs(essdof2);
for (int i=0; i<essdof2.Size(); i++) {
    essdof2[i] += X.GetNDofs() ;
}
ess_dof.Append(essdof2);

Then FormLinearSystem is called on the Blockoperator $m$ using the essdof array. This results in the new system $M\cdot Z=B$. The following steps are:

mfem::MINRES(*M, B, Z, 0, iter, tol*tol, tol*tol);

m.RecoverFEMSolution(Z, b, x);
mfem::Array<int> u_dofs (u.Size());
std::iota(&u_dofs[0], &u_dofs[u.Size()], 0);
x.GetSubVector(u_dofs, u);

Is it necessary to produce a symmetric matrix by multiplying the system by the transposed of $m$ before calling MINRES?

[v-dobrev]

MINRES requires a symmetric matrix. You can try using GMRES instead, unless $C = B^T$ and A is symmetric.

[markusrenoldner]

Thank you.

Just to document an alternative solution: Enforcing boundary conditions by adding equations:

$$ \begin{bmatrix} A\\ D \end{bmatrix} \cdot \begin{bmatrix} u\\ \tilde{u} \end{bmatrix}= \begin{bmatrix} f\\ \tilde{g} \end{bmatrix} $$

where $D$ is a sparse $n\times m$-matrix where $n$ (number of rows) is the number of essential DOFs, and $m$ is the the number of cols of $A$. $D$ has a $1$ in each row, in the column of the index of the essential DOFs of $u$. E.g.:

$$ \begin{bmatrix} 1 & 0&0 &0 &0&0\\ 1 & 1&0 &0 &0&0\\ 1 & 0&0 &0 &1&0\\ 1 & 0&0 &0 &0&1\\ \end{bmatrix} $$

this choice of $D$ fixes $u_0$, $u_1$, $u_4$, and $u_5$.

The vector $\tilde{g}$ contains the dirichlet values and $\tilde{u}$ simply contains the essential DOFs.

[v-dobrev]

If I understand you correctly, I don't think the above makes sense in the standard case, e.g. mfem's example 1. The reason is that the added equations simply make the system over-determined, in general. What mfem does is to replace equations (rows in $A$) with new equations of the form $x_j = g_j$ and at the same time eliminate the unknown $x_j$ from all other equations using the value $g_j$. This corresponds to the following modification of the matrix and r.h.s. performed by FormLinearSystem (in the case of DIAG_ONE policy):

$$ \begin{pmatrix} A_i & A_{ib} \\ A_{bi} & A_b \end{pmatrix} \rightarrow \begin{pmatrix} A_i & 0 \\ 0 & I \end{pmatrix} \qquad \begin{pmatrix} b_i \\ b_b \end{pmatrix} \rightarrow \begin{pmatrix} b_i - A_{ib} x_b \\ x_b \end{pmatrix} $$

where the block index $i$ denotes the "interior" (non-essential) dofs, and the block index $b$ denotes the "boundary" (essential) dofs. Clearly, this modified system is equivalent to the reduced FEM system with imposed essential b.c.: $A_i x_i = b_i - A_{ib} x_b$.

[markusrenoldner]

Yes that is correct indeed!

I should have mentioned, that the problem I am solving was inherently underdetermined and only by adding equations, the problem was posed in a meaningful way.

Thank you for pointing out the premises of our approach.