I am solving the poisson equation and I constructed the global stiffness matrix in compressed row storage format. Then I wrote the preconditioned conjugate gradient solver for solving the system of equations. Now my problem is how I can apply the boundary condition when I have the global stiffness matrix in csr format. I can assemble the global stiffness matrix as a sparse matrix and then easily apply the boundary condition by removing some column and row but it's not an efficient way. I would be grateful if you help me with detail.
1
-
$\begingroup$ No,it doesn't, I stored the global stiffness matrix in csr format,it means I build three row vectors. The first row vector store the non_zero element, second one stores the column index of non_zero element and the last one is a pointer to the first non_zero element in each row. I can remove the column an row in the big sparse matrix and apply the boundary condition. But it's not just about solving the poisson equation,it's about memory optimization and solving the system of linear equations parallelized with mpi. So I need to apply the boundary condition in csr format. $\endgroup$– ResaCommented Mar 11, 2021 at 13:12
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
2
The process is not entirely trivial, but you might be interested in watching video lectures 21.6 and 21.65 here: https://www.math.colostate.edu/~bangerth/videos.html
answered Mar 11, 2021 at 10:29
-
$\begingroup$ Thanks, I watched these lectures before but they don't answer my questions. By the way thanks for your help. $\endgroup$– ResaCommented Mar 11, 2021 at 13:41
-
$\begingroup$ @Resa Then I don't know what your question is. That may also be because your question does not, in fact, have a question (or question mark). $\endgroup$ Commented Mar 12, 2021 at 17:11
$\begingroup$
$\endgroup$
5
I do something similar for a structural mechanics problem but I don't remove any rows or columns. I use the penalty method by modifying the diagonal term of the degree of freedom where I want to impose some BC (Dirichlet).
I don't know exactly how you're assembling your global sparse stiffness matrix but I keep a track of where in the global stiffness matrix are the entries belonging of the ith degree of freedom. So basically a lookup table of where the diagonal entry of the ith dof falls and I can modifying exactly that value without disturbing the rest of the matrix. Hope that helps.
-
$\begingroup$ Thanks for your reply, are using csr format for storing non-zero element? $\endgroup$– ResaCommented Mar 11, 2021 at 14:43
-
$\begingroup$ Yes, storing in CSR format and solving using Intel MKL direct sparse solver to solve. I modify only the non-zero values vector, not the column_index and row_ptr vectors. $\endgroup$ Commented Mar 11, 2021 at 15:02
-
$\begingroup$ I assembled the global stiffness matrix such a way that I don't construct the big sparse matrix. So the only thing I have is the non-zero element as a row vector, col_index of each non-zero element and the row_ptr vector as a pointer to the first non-zero element. And I'm using preconditioned conjugate gradient method. I can use some software to apply boundary conditions but I need to know what is going on internally. Do you have any source about this? $\endgroup$– ResaCommented Mar 11, 2021 at 16:39
-
$\begingroup$ Yes I know what you're saying. I'm do the same thing, I never construct the full matrix but only the CSR representation. But with information about how many non-zeros you have for each degree of freedom, you can locate where you're diagonal entries are within the CSR representation. For example if you know the each degree of freedom has 10 entries, you know that the the diagonal entry for the 10th element will be at the 100th position. Highly simplified but hope you get the point I'm trying to make. Check this: link.springer.com/article/10.1007/s42102-020-00041-y $\endgroup$ Commented Mar 11, 2021 at 17:02
-
$\begingroup$ It not trivial, but can be done depending on how much time you can invest in it $\endgroup$ Commented Mar 11, 2021 at 17:11
$\begingroup$
$\endgroup$
2
It's not that hard and there is a fortran code example that illustrates that. It is from John Burkardt collection of code examples. Also C and C++ version is available with a simple search within same collection.
This is coordinate format, or as he calls it sparse triplet. With a simple modification to store pointers to elements starting the row in ia instead of row numbers you get CSR version.
The procedure is to store position of diagonal elements in the array of sparse matrix elements, of size nnz. You access it and set it to one if it corresponds to node/DOF at Dirichlet boundary, while the corresponding element in the right-hand side vector gets value of the Dirichlet BC.
subroutine dirichlet_apply_dsp ( node_num, node_xy, node_condition, &
nz_num, ia, ja, a, f )
!*****************************************************************************80
!
!! DIRICHLET_APPLY_DSP accounts for Dirichlet boundary conditions.
!
! Discussion:
!
! It is assumed that the matrix A and right hand side F have already been
! set up as though there were no boundary conditions. This routine
! then modifies A and F, essentially replacing the finite element equation
! at a boundary node NODE by a trivial equation of the form
!
! A(NODE,NODE) * U(NODE) = NODE_BC(NODE)
!
! where A(NODE,NODE) = 1.
!
! This routine assumes that the coefficient matrix is stored in a
! sparse triplet format.
!
! This routine implicitly assumes that the sparse matrix has a storage
! location for every diagonal element...or at least for those diagonal
! elements corresponding to boundary nodes.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 12 July 2007
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of nodes.
!
! Input, integer ( kind = 4 ) NODE_CONDITION(NODE_NUM), reports the
! condition used to set the unknown associated with the node.
! 0, unknown.
! 1, finite element equation.
! 2, Dirichlet condition;
! 3, Neumann condition.
!
! Input, integer ( kind = 4 ) NZ_NUM, the number of nonzero entries.
!
! Input, integer ( kind = 4 ) IA(NZ_NUM), JA(NZ_NUM), the row and column
! indices of the nonzero entries.
!
! Input/output, real ( kind = 8 ) A(NZ_NUM), the coefficient matrix,
! stored in sparse triplet format; on output, the matrix has been adjusted
! for Dirichlet boundary conditions.
!
! Input/output, real ( kind = 8 ) F(NODE_NUM), the right hand side.
! On output, the right hand side has been adjusted for Dirichlet
! boundary conditions.
!
implicit none
integer ( kind = 4 ) node_num
integer ( kind = 4 ) nz_num
real ( kind = 8 ), dimension(nz_num) :: a
integer ( kind = 4 ) column
integer ( kind = 4 ), parameter :: DIRICHLET = 2
real ( kind = 8 ), dimension(node_num) :: f
integer ( kind = 4 ) ia(nz_num)
integer ( kind = 4 ) ja(nz_num)
integer ( kind = 4 ) node
real ( kind = 8 ), dimension ( node_num ) :: node_bc
integer ( kind = 4 ) node_condition(node_num)
real ( kind = 8 ), dimension(2,node_num) :: node_xy
integer ( kind = 4 ) nz
!
! Retrieve the Dirichlet boundary condition value at every node.
!
call dirichlet_condition ( node_num, node_xy, node_bc ) ! SUPPLY THIS FUNCTION
!
! Consider every matrix entry, NZ.
!
! If the row I corresponds to a boundary node, then
! zero out all off diagonal matrix entries, set the diagonal to 1,
! and the right hand side to the Dirichlet boundary condition value.
!
do nz = 1, nz_num
node = ia(nz)
if ( node_condition(node) == DIRICHLET ) then
column = ja(nz)
if ( column == node ) then
a(nz) = 1.0D+00
f(node) = 0.0D0 ! node_bc(node) FIXME
else
a(nz) = 0.0D+00
end if
end if
end do
return
end
```
-
$\begingroup$ Thanks for your answer but I don't understand why we need to put one instead of diagonal an zero out of diagonal elements $\endgroup$– ResaCommented Mar 12, 2021 at 11:45
-
1$\begingroup$ That's basic linear algebra, that row will multiply solution vector so in that case... $\endgroup$ Commented Mar 12, 2021 at 21:41