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