Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/ibcmat.f90
03540012fedefeature/fft-split

Last Change
Raw File

ibcmat.f90
View Options

!>
!> @file ibcmat.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
SUBROUTINE ibcmat(mat, ny)
!
! Apply BC on matrix
!
USE matrix
IMPLICIT NONE
TYPE(gbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: ny
INTEGER :: kl, ku, nrank, i, j
DOUBLE PRECISION, ALLOCATABLE :: zsum(:), arr(:)
INTEGER :: i0, ii
INTEGER :: i0_arr(ny), col(ny)
!===========================================================================
! 1.0 Prologue
!
kl = mat%kl
ku = mat%ku
nrank = mat%rank
ALLOCATE(zsum(nrank), arr(nrank))
!
i0 = nrank - ku
WRITE(*,'(a,i6)') 'Estimated i0', i0
DO i=1,ny
arr = 0.0d0
col(i) = nrank-ny+i
CALL getcol(mat, nrank-ny+i, arr)
DO ii=1,nrank
i0_arr(i)=ii
IF(arr(ii) .NE. 0.0d0) EXIT
END DO
END DO
!!$ WRITE(*,'(a/(10i6))') 'col', col
!!$ WRITE(*,'(a/(10i6))') 'i0_arr', i0_arr
!
!===========================================================================
! 2.0 Unicity at the axis
!
! The vertical sum on the NY-th row
!
zsum = 0.0d0
DO i=1,ny
arr = 0.0d0
CALL getrow(mat, i, arr)
DO j=1,ny+ku
zsum(j) = zsum(j) + arr(j)
END DO
END DO
CALL putrow(mat, ny, zsum)
!
! The horizontal sum on the NY-th column
!
zsum = 0.0d0
DO j=1,ny
arr = 0.0d0
CALL getcol(mat, j, arr)
DO i=ny,ny+kl
zsum(i) = zsum(i) + arr(i)
END DO
END DO
CALL putcol(mat, ny, zsum)
!
! The away operator
!
DO j = 1,ny-1
arr = 0.0d0; arr(j) = 1.0d0
CALL putcol(mat, j, arr)
END DO
!
DO i = 1,ny-1
arr = 0.0d0; arr(i) = 1.0d0
CALL putrow(mat, i, arr)
END DO
!
!===========================================================================
! 3.0 Dirichlet on right boundary
!
DO j = nrank, nrank-ny+1, -1
arr = 0.0d0; arr(j) = 1.0d0
CALL putcol(mat, j, arr)
END DO
!
DO i = nrank, nrank-ny+1, -1
arr = 0.0d0; arr(i) = 1.0d0
CALL putrow(mat, i, arr)
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(zsum, arr)
!
END SUBROUTINE ibcmat
!+++
SUBROUTINE ibcrhs(rhs, ny)
!
! Apply BC on RHS
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(inout) :: rhs(:)
INTEGER, INTENT(in) :: ny
INTEGER :: nrank
DOUBLE PRECISION :: zsum
!===========================================================================
! 1.0 Prologue
!
nrank = SIZE(rhs,1)
!===========================================================================
! 2.0 Unicity at the axis
!
! The vertical sum
!
zsum = SUM(rhs(1:ny))
rhs(ny) = zsum
rhs(1:ny-1) = 0.0d0
!===========================================================================
! 3.0 Dirichlet on right boundary
!
rhs(nrank-ny+1:nrank) = 0.0d0
END SUBROUTINE ibcrhs
!+++
SUBROUTINE ibcrhs3(rhs, ny)
!
! Apply BC on RHS
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(inout) :: rhs(:,:)
INTEGER, INTENT(in) :: ny
INTEGER :: nrank, nz, k
DOUBLE PRECISION :: zsum
!===========================================================================
! 1.0 Prologue
!
nrank = SIZE(rhs,1)
nz = SIZE(rhs,2)
!===========================================================================
! 2.0 Unicity at the axis
!
! The vertical sum
!
DO k=1,nz
zsum = SUM(rhs(1:ny,k))
rhs(ny,k) = zsum
rhs(1:ny-1,k) = 0.0d0
END DO
!===========================================================================
! 3.0 Dirichlet on right boundary
!
DO k=1,nz
rhs(nrank-ny+1:nrank,k) = 0.0d0
END DO
END SUBROUTINE ibcrhs3
c4science · Help