Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/multigrid/src/fdmat_mod.f90
03540012fedefeature/fft-split

Last Change
Raw File

fdmat_mod.f90
View Options

!>
!> @file fdmat_mod.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>
!>
MODULE fdmat_mod
!
USE multigrid
IMPLICIT NONE
!
INTERFACE fdmat
MODULE PROCEDURE fdmat_stencil
MODULE PROCEDURE fdmat_gen, fdmat_csr, fdmat_cds
END INTERFACE fdmat
INTERFACE ibc_fdmat
MODULE PROCEDURE ibc_fdmat_stencil
MODULE PROCEDURE ibc_fdmat_gen, ibc_fdmat_csr, ibc_fdmat_cds
END INTERFACE ibc_fdmat
INTERFACE ibc_rhs
MODULE PROCEDURE ibc_rhs_g
END INTERFACE ibc_rhs
!
CONTAINS
!--------------------------------------------------------------------------------
SUBROUTINE fdmat_stencil(grid, fdense, icrosst, mat)
!
! Construct model GBS FD partitioned matrix
!
USE iso_fortran_env, ONLY : rkind => real64
USE stencil, ONLY : stencil_2d
USE parmg, ONLY : grid2_type
!
TYPE(grid2_type) :: grid
REAL(rkind), INTENT(in) :: icrosst
TYPE(stencil_2d) :: mat
INTERFACE
FUNCTION fdense(x)
USE iso_fortran_env, ONLY : rkind => real64
REAL(rkind), INTENT(in) :: x(:)
REAL(rkind) :: fdense(SIZE(x))
END FUNCTION fdense
END INTERFACE
!
REAL(rkind) :: dx, dy
REAL(rkind),ALLOCATABLE :: dense(:)
REAL(rkind) :: stencil_arr(-1:1,-1:1), zdiag(-1:1,-1:1), corr
INTEGER :: nx, ny, i, j, k, d(2)
!
! Grid properties
!
nx = grid%npt(1)-1
ny = grid%npt(2)-1
dx = grid%x(1)-grid%x(0) ! Assume equidistant grid
dy = grid%y(1)-grid%y(0)
ALLOCATE(dense(0:nx)) ! electron densoty vary only along x
!
! Stencil array
!
stencil_arr = 0.0d0
zdiag = 0.0d0
corr = 1.d0+icrosst**2/4.0d0
stencil_arr(0,0) = -2.0d0/dx/dx-2.0d0/dy/dy*corr
stencil_arr(-1,0) = 1.0/dx/dx
stencil_arr(1,0) = 1.0/dx/dx
stencil_arr(0,-1) = 1.0/dy/dy*corr
stencil_arr(0,1) = 1.0/dy/dy*corr
stencil_arr(-1,-1) = icrosst*1.0/4.0/dx/dy
stencil_arr(1,-1) = icrosst*(-1.0/4.0/dx/dy)
stencil_arr(-1,1) = icrosst*(-1.0/4.0/dx/dy)
stencil_arr(1,1) = icrosst*(1.0/4.0/dx/dy)
zdiag(0,0) = 1.0d0
dense(:) = fdense(grid%x(:))
!
! Assemble the stencil by scanning local grid points
!
DO k=0,mat%npoints-1
d(:) = mat%id(k,:)
DO j=mat%s(2),mat%e(2)
DO i=mat%s(1),mat%e(1)
mat%val(i,j,k) = stencil_arr(d(1),d(2)) + &
& zdiag(d(1),d(2))*dense(i)
END DO
END DO
END DO
!
DEALLOCATE(dense)
END SUBROUTINE fdmat_stencil
!--------------------------------------------------------------------------------
SUBROUTINE ibc_rhs_g(f, s0, e0, prb)
!
! Impose BC on rhs
!
USE iso_fortran_env, ONLY : rkind => real64
USE gvector, ONLY : gvector_2d
!
TYPE(gvector_2d), INTENT(inout) :: f
CHARACTER(len=*), INTENT(in) :: prb
INTEGER, INTENT(in) :: s0(2), e0(2)
INTEGER :: s(2), e(2)
!
s = f%s
e = f%e
!
IF(s(1).EQ.s0(1)) THEN
IF(prb(1:1).EQ.'d') THEN ! West face
f%val(s(1),s(2):e(2)) = 0.0_rkind
ELSE
f%val(s(1),s(2):e(2)) = 0.5_rkind*f%val(s(1),s(2):e(2))
END IF
END IF
IF(e(1).EQ.e0(1)) THEN
IF(prb(2:2).EQ.'d') THEN ! East face
f%val(e(1),s(2):e(2)) = 0.0_rkind
ELSE
f%val(e(1),s(2):e(2)) = 0.5_rkind*f%val(e(1),s(2):e(2))
END IF
END IF
IF(s(2).EQ.s0(2)) THEN
IF(prb(3:3).EQ.'d') THEN ! South face
f%val(s(1):e(1),s(2)) = 0.0_rkind
ELSE
f%val(s(1):e(1),s(2)) = 0.5_rkind*f%val(s(1):e(1),s(2))
END IF
END IF
IF(e(2).EQ.e0(2)) THEN
IF(prb(4:4).EQ.'d') THEN ! North face
f%val(s(1):e(1),e(2)) = 0.0_rkind
ELSE
f%val(s(1):e(1),e(2)) = 0.5_rkind*f%val(s(1):e(1),e(2))
END IF
END IF
END SUBROUTINE ibc_rhs_g
!--------------------------------------------------------------------------------
SUBROUTINE ibc_fdmat_stencil(mat, prb)
!
! Impose BC on matrix
!
USE iso_fortran_env, ONLY : rkind => real64
USE stencil, ONLY : stencil_2d
!
TYPE(stencil_2d), INTENT(inout) :: mat
CHARACTER(len=*), INTENT(in) :: prb
!
INTEGER :: s0(2), e0(2), s(2), e(2)
!
s0 = mat%s0
e0 = mat%e0
s = mat%s
e = mat%e
!
! Neumann BC
! WARNING: Divide the stencil by 2 => should do the same for RHS!
!
! N
! 6---7---8
! | | |
! W 4---0---5 E Numbering of stencil
! | | |
! 1---2---3
! S
!
IF(s(1).EQ.s0(1) .AND. prb(1:1).EQ.'n') THEN ! West face
mat%val(s(1),s(2):e(2),1) = 0.0_rkind
mat%val(s(1),s(2):e(2),3) = 0.0_rkind
mat%val(s(1),s(2):e(2),4) = 0.0_rkind
mat%val(s(1),s(2):e(2),5) = 2.0d0*mat%val(s(1),s(2):e(2),5)
mat%val(s(1),s(2):e(2),6) = 0.0_rkind
mat%val(s(1),s(2):e(2),8) = 0.0_rkind
mat%val(s(1),s(2):e(2),:) = 0.5_rkind*mat%val(s(1),s(2):e(2),:)
END IF
IF(e(1).EQ.e0(1) .AND. prb(2:2).EQ.'n') THEN ! East face
mat%val(e(1),s(2):e(2),1) = 0.0_rkind
mat%val(e(1),s(2):e(2),3) = 0.0_rkind
mat%val(e(1),s(2):e(2),4) = 2.0d0*mat%val(e(1),s(2):e(2),4)
mat%val(e(1),s(2):e(2),5) = 0.0_rkind
mat%val(e(1),s(2):e(2),6) = 0.0_rkind
mat%val(e(1),s(2):e(2),8) = 0.0_rkind
mat%val(e(1),s(2):e(2),:) = 0.5_rkind*mat%val(e(1),s(2):e(2),:)
END IF
IF(s(2).EQ.s0(2) .AND. prb(3:3).EQ.'n') THEN ! South face
mat%val(s(1):e(1),s(2),1) = 0.0_rkind
mat%val(s(1):e(1),s(2),2) = 0.0_rkind
mat%val(s(1):e(1),s(2),3) = 0.0_rkind
mat%val(s(1):e(1),s(2),6) = 0.0_rkind
mat%val(s(1):e(1),s(2),7) = 2.0d0*mat%val(s(1):e(1),s(2),7)
mat%val(s(1):e(1),s(2),8) = 0.0_rkind
mat%val(s(1):e(1),s(2),:) = 0.5_rkind*mat%val(s(1):e(1),s(2),:)
END IF
IF(e(2).EQ.e0(2) .AND. prb(4:4).EQ.'n') THEN ! North face
mat%val(s(1):e(1),e(2),1) = 0.0_rkind
mat%val(s(1):e(1),e(2),2) = 2.0d0*mat%val(s(1):e(1),e(2),2)
mat%val(s(1):e(1),e(2),3) = 0.0_rkind
mat%val(s(1):e(1),e(2),6) = 0.0_rkind
mat%val(s(1):e(1),e(2),7) = 0.0_rkind
mat%val(s(1):e(1),e(2),8) = 0.0_rkind
mat%val(s(1):e(1),e(2),:) = 0.5_rkind*mat%val(s(1):e(1),e(2),:)
END IF
!
! Dirichlet BC
!
IF(s(1).EQ.s0(1) .AND. prb(1:1).EQ.'d') THEN ! West face
mat%val(s(1),s(2):e(2),:) = 0.0_rkind
mat%val(s(1),s(2):e(2),0) = 1.0_rkind
END IF
IF(e(1).EQ.e0(1) .AND. prb(2:2).EQ.'d') THEN ! East face
mat%val(e(1),s(2):e(2),:) = 0.0_rkind
mat%val(e(1),s(2):e(2),0) = 1.0_rkind
END IF
IF(s(2).EQ.s0(2) .AND. prb(3:3).EQ.'d') THEN ! South face
mat%val(s(1):e(1),s(2),:) = 0.0_rkind
mat%val(s(1):e(1),s(2),0) = 1.0_rkind
END IF
IF(e(2).EQ.e0(2) .AND. prb(4:4).EQ.'d') THEN ! North face
mat%val(s(1):e(1),e(2),:) = 0.0_rkind
mat%val(s(1):e(1),e(2),0) = 1.0_rkind
END IF
END SUBROUTINE ibc_fdmat_stencil
!--------------------------------------------------------------------------------
SUBROUTINE fdmat_gen(grid, fdense, icrosst, noinit)
!
! Generic version
!
TYPE(grid2d), INTENT(inout) :: grid
DOUBLE PRECISION, INTENT(in) :: icrosst
LOGICAL, INTENT(in), OPTIONAL :: noinit
INTERFACE
FUNCTION fdense(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: fdense(SIZE(x))
END FUNCTION fdense
END INTERFACE
!
IF(ALLOCATED(grid%mata)) THEN
CALL fdmat_csr(grid, fdense, icrosst, grid%mata, noinit)
ELSE
CALL fdmat_cds(grid, fdense, icrosst, grid%mata_cds, noinit)
END IF
END SUBROUTINE fdmat_gen
!--------------------------------------------------------------------------------
SUBROUTINE fdmat_cds(grid, fdense, icrosst, mat, noinit)
!
! Construct FD matrix
!
TYPE(grid2d), INTENT(in) :: grid
DOUBLE PRECISION, INTENT(in) :: icrosst
TYPE(cds_mat), INTENT(inout) :: mat
LOGICAL, INTENT(in), OPTIONAL :: noinit
INTERFACE
FUNCTION fdense(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: fdense(SIZE(x))
END FUNCTION fdense
END INTERFACE
!
INTEGER :: n, nx, ny
INTEGER :: kl, ku, k
INTEGER :: ix, ix2, jx, iy, jy, iy2, irow
DOUBLE PRECISION :: lx, ly, dx, dy, mele
DOUBLE PRECISION :: dense(0:grid%n(1))
DOUBLE PRECISION :: stencil(-1:1,-1:1), zdiag(-1:1,-1:1)
DOUBLE PRECISION :: corr
LOGICAL :: run_init
INTEGER, ALLOCATABLE :: dists(:)
!--------------------------------------------------------------------------------
run_init = .TRUE.
IF(PRESENT(noinit)) run_init = .NOT.noinit
!
! Grid properties
!
nx = grid%n(1)
ny = grid%n(2)
dx = grid%x(1) - grid%x(0)
dy = grid%y(1) - grid%y(0)
lx = grid%x(nx)
ly = grid%y(ny)
n = PRODUCT(grid%rank) ! Rank of matrix
!
! Stencil
!
stencil = 0.0d0
zdiag = 0.0d0
!
corr = 1.d0+icrosst**2/4.0d0
stencil(0,0) = -2.0d0/dx/dx-2.0d0/dy/dy*corr
stencil(-1,0) = 1.0/dx/dx
stencil(1,0) = 1.0/dx/dx
stencil(0,-1) = 1.0/dy/dy*corr
stencil(0,1) = 1.0/dy/dy*corr
stencil(-1,-1) = icrosst*1.0/4.0/dx/dy
stencil(1,-1) = icrosst*(-1.0/4.0/dx/dy)
stencil(-1,1) = icrosst*(-1.0/4.0/dx/dy)
stencil(1,1) = icrosst*(1.0/4.0/dx/dy)
zdiag(0,0) = 1.0d0
!
! 9-point stencil "diagonal storage"
!
kl=4
ku=4
ALLOCATE(dists(-kl:ku))
DO iy2=-1,1
DO ix2=-1,1
k=3*iy2+ix2
dists(k) = iy2*(nx+1) + ix2
END DO
END DO
!
IF(run_init) THEN
CALL init(n, dists, 1, mat)
END IF
!
! Assemble matrix by scanning all grid points
!
dense(:) = fdense(grid%x(:))
DO iy=0,ny
DO ix=0,nx
irow = iy*(nx+1)+ix+1
DO iy2=-1,1
jy=iy+iy2
IF(jy.GE.0 .AND. jy.LE.ny) THEN
DO ix2=-1,1
jx=ix+ix2
IF(jx.GE.0 .AND.jx.LE.nx) THEN
mele = stencil(ix2,iy2) + zdiag(ix2,iy2)*dense(ix)
k=3*iy2+ix2
mat%val(irow,k) = mele
END IF
END DO
END IF
END DO
END DO
END DO
!
DEALLOCATE(dists)
END SUBROUTINE fdmat_cds
!--------------------------------------------------------------------------------
SUBROUTINE fdmat_csr(grid, fdense, icrosst, mat, noinit)
!
! Construct FD matrix
!
TYPE(grid2d), INTENT(in) :: grid
DOUBLE PRECISION, INTENT(in) :: icrosst
TYPE(csr_mat), INTENT(inout) :: mat
LOGICAL, INTENT(in), OPTIONAL :: noinit
INTERFACE
FUNCTION fdense(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: fdense(SIZE(x))
END FUNCTION fdense
END INTERFACE
!
INTEGER :: n, nx, ny
INTEGER :: ix, ix2, jx, iy, jy, iy2, irow, icol
DOUBLE PRECISION :: lx, ly, dx, dy, mele
DOUBLE PRECISION :: dense(0:grid%n(1))
DOUBLE PRECISION :: stencil(-1:1,-1:1), zdiag(-1:1,-1:1)
DOUBLE PRECISION :: corr
LOGICAL :: run_init
!--------------------------------------------------------------------------------
run_init = .TRUE.
IF(PRESENT(noinit)) run_init = .NOT.noinit
!
! Grid properties
!
nx = grid%n(1)
ny = grid%n(2)
dx = grid%x(1) - grid%x(0)
dy = grid%y(1) - grid%y(0)
lx = grid%x(nx)
ly = grid%y(ny)
n = PRODUCT(grid%rank) ! Rank of matrix
!
! Stencil
!
stencil = 0.0d0
zdiag = 0.0d0
!
corr = 1.d0+icrosst**2/4.0d0
stencil(0,0) = -2.0d0/dx/dx-2.0d0/dy/dy*corr
stencil(-1,0) = 1.0/dx/dx
stencil(1,0) = 1.0/dx/dx
stencil(0,-1) = 1.0/dy/dy*corr
stencil(0,1) = 1.0/dy/dy*corr
stencil(-1,-1) = icrosst*1.0/4.0/dx/dy
stencil(1,-1) = icrosst*(-1.0/4.0/dx/dy)
stencil(-1,1) = icrosst*(-1.0/4.0/dx/dy)
stencil(1,1) = icrosst*(1.0/4.0/dx/dy)
zdiag(0,0) = 1.0d0
!
! Create CSR matrix
IF(run_init) THEN
CALL init(n, 1, mat)
END IF
!
! Assemble matrix by scanning all grid points
!
dense(:) = fdense(grid%x(:))
DO iy=0,ny
DO ix=0,nx
irow = numb(ix,iy)
DO iy2=-1,1
jy=iy+iy2
IF(jy.GE.0 .AND. jy.LE.ny) THEN
DO ix2=-1,1
jx=ix+ix2
IF(jx.GE.0 .AND.jx.LE.nx) THEN
icol=numb(jx,jy)
mele = stencil(ix2,iy2) + zdiag(ix2,iy2)*dense(ix)
CALL putele(mat, irow, icol,mele)
END IF
END DO
END IF
END DO
END DO
END DO
!--------------------------------------------------------------------------------
CONTAINS
INTEGER FUNCTION numb(ix,iy)
INTEGER, INTENT(in) :: ix, iy
INTEGER :: stride
stride = grid%rank(1)
numb = iy*stride + (ix+1)
END FUNCTION numb
!--------------------------------------------------------------------------------
END SUBROUTINE fdmat_csr
!++
SUBROUTINE ibc_fdmat_gen(grid, prb)
!
! Generic version
!
TYPE(grid2d), INTENT(inout) :: grid
CHARACTER(len=*), INTENT(in) :: prb
!
IF(ALLOCATED(grid%mata)) THEN
CALL ibc_fdmat_csr(grid, grid%mata, prb)
ELSE
CALL ibc_fdmat_cds(grid, grid%mata_cds, prb)
END IF
!
END SUBROUTINE ibc_fdmat_gen
!++
SUBROUTINE ibc_fdmat_csr(grid, mat, prb)
!
! Impose BC
!
TYPE(grid2d), INTENT(in) :: grid
TYPE(csr_mat), INTENT(inout) :: mat
CHARACTER(len=*), INTENT(in) :: prb
!
DOUBLE PRECISION :: arow(mat%rank)
INTEGER :: nx, ny, nx1, ny1, n, iy, irow, irow1
!--------------------------------------------------------------------------------
nx = grid%n(1)
ny = grid%n(2)
nx1=nx+1
ny1=ny+1
n = nx1*ny1
!
! Dirichelt BC on West/East
!
IF(prb.EQ.'dddd') THEN
DO irow=1,ny*nx1+1,nx1
arow=0.0d0; arow(irow)=1.0d0
CALL putrow(mat, irow, arow)
irow1=irow+nx
arow=0.0d0; arow(irow1)=1.0d0
CALL putrow(mat, irow1, arow)
END DO
!
! Neumann on West/East
! WARNING: Divide the stencil by 2 => should do the same for RHS!
!
ELSE IF(prb.EQ.'nndd') THEN
DO irow=1,ny*nx1+1,nx1
iy = irow/nx1
CALL getrow(mat, irow, arow)
arow(irow+1) = 2.0d0*arow(irow+1)
IF(iy.GT.0) arow(irow-nx) = 0.0d0
IF(iy.LT.ny) arow(irow+nx+2) = 0.0d0
arow(:) = 0.5d0*arow(:)
CALL putrow(mat, irow, arow)
END DO
DO irow=nx1,n,nx1
iy = irow/nx1
CALL getrow(mat, irow, arow)
arow(irow-1) = 2.0d0*arow(irow-1)
IF(iy.GT.0) arow(irow-nx-2) = 0.0d0
IF(iy.LT.ny) arow(irow+nx) = 0.0d0
arow(:) = 0.5d0*arow(:)
CALL putrow(mat, irow, arow)
END DO
!
ELSE
WRITE(*,'(a,a4,a)') 'ibc_mat: prb = ', prb, ' NOT IMPLEMENTED!'
STOP
END IF
!
! Dirichlet BC on South/North sides
!
DO irow=1,nx1
arow=0.0d0; arow(irow)=1.0d0
CALL putrow(mat, irow, arow)
END DO
DO irow1=ny*nx1+1,n
arow=0.0d0; arow(irow1)=1.0d0
CALL putrow(mat, irow1, arow)
END DO
!--------------------------------------------------------------------------------
END SUBROUTINE ibc_fdmat_csr
!++
SUBROUTINE ibc_fdmat_cds(grid, mat, prb)
!
! Impose BC
!
TYPE(grid2d), INTENT(in) :: grid
TYPE(cds_mat), INTENT(inout) :: mat
CHARACTER(len=*), INTENT(in) :: prb
!
INTEGER :: nx, ny, iy, irow
INTEGER :: n, nx1, ny1
!--------------------------------------------------------------------------------
nx = grid%n(1)
ny = grid%n(2)
nx1=nx+1
ny1=ny+1
n = nx1*ny1
!
! 2 == 3 == 4
! | | |
! -1 == 0 == 1
! | | |
! -4 == -3 == -2
!
! Dirichelt BC on West/East
!
IF(prb.EQ.'dddd') THEN
DO irow=1,ny*nx1+1,nx1
mat%val(irow,:) = 0.0d0
mat%val(irow,0) = 1.0d0
mat%val(irow+nx,:) = 0.0d0
mat%val(irow+nx,0) = 1.0d0
END DO
!
! Neumann on West/East
! WARNING: Divide the stencil by 2 => should do the same for RHS!
!
ELSE IF(prb.EQ.'nndd') THEN
DO irow=1,ny*nx1+1,nx1
iy = irow/nx1
IF(iy.GT.0) mat%val(irow,-2)=0.0d0
IF(iy.LT.ny) mat%val(irow,+4)=0.0d0
mat%val(irow,+1)=2.0d0*mat%val(irow,+1)
mat%val(irow,:)=0.5d0*mat%val(irow,:)
END DO
DO irow=nx1,n,nx1
iy = irow/nx1
IF(iy.GT.0) mat%val(irow,-4)=0.0d0
IF(iy.LT.ny) mat%val(irow,+2)=0.0d0
mat%val(irow,-1)=2.0d0*mat%val(irow,-1)
mat%val(irow,:)=0.5d0*mat%val(irow,:)
END DO
!
ELSE
WRITE(*,'(a,a4,a)') 'ibc_mat: prb = ', prb, ' NOT IMPLEMENTED!'
STOP
END IF
!
! Dirichlet BC on South/North sides
!
DO irow=1,nx1
mat%val(irow,:) = 0.0d0
mat%val(irow,0) = 1.0d0
END DO
DO irow=ny*nx1+1,n
mat%val(irow,:) = 0.0d0
mat%val(irow,0) = 1.0d0
END DO
!--------------------------------------------------------------------------------
END SUBROUTINE ibc_fdmat_cds
!++
END MODULE fdmat_mod
c4science · Help