Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

test_intergrid0.f90
View Options

!>
!> @file test_intergrid0.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>
!>
PROGRAM main
!
! intergrid transfer using *serial* multigrid module:
! - restriction of rhs
! - prolongation of sol
!
USE multigrid, ONLY : grid2d, mg_info, &
& get_lmax, create_grid_fd, ibc_transf, &
& prolong, restrict
IMPLICIT NONE
DOUBLE PRECISION, PARAMETER :: pi=4.0d0*ATAN(1.0d0)
DOUBLE PRECISION :: Lx, Ly, kx, ky, icrosst, beta, miome
INTEGER :: nx, ny, levels
CHARACTER(len=4) :: prb
LOGICAL :: nldebug
!
DOUBLE PRECISION :: dx, dy
DOUBLE PRECISION, ALLOCATABLE :: x(:),y(:)
!
TYPE(mg_info) :: info ! info for MG
TYPE(grid2d), ALLOCATABLE :: grids(:)
!
INTEGER :: i, l
!
NAMELIST /parameters/ prb, nx, ny, levels, Lx, Ly, kx, ky, icrosst, beta, &
& miome, nldebug
!--------------------------------------------------------------------------------
!
! Default inputs
!
nx=32
ny=32
levels = 2
kx=1
ky=1
icrosst=1.0d0
Lx = 1.0D0
Ly = 1.0D0
miome = 200d0
beta = 0d0
prb = 'dddd'
nldebug = .FALSE.
!
READ(*,parameters)
WRITE(*,parameters)
!
! Fine grid
!
dx = lx/REAL(nx,8)
dy = ly/REAL(ny,8)
ALLOCATE(x(0:nx), y(0:ny))
x = dx * [(i,i=0,nx)]
y = dy * [(i,i=0,ny)]
WRITE(*,'(a/10(1pe12.3))') 'x =', x
WRITE(*,'(a/10(1pe12.3))') 'y =', y
!
! Create array of grids
!
levels = MIN(levels, get_lmax(nx), get_lmax(ny))
WRITE(*,'(a,i4)') 'Number of levels', levels
ALLOCATE(grids(levels))
info%nu1 = 1
info%nu2 = 1
info%mu = 1
info%nu0 = 1
info%levels = levels
info%relax = 'jac'
info%omega = 1
CALL create_grid_fd(x, y, grids, info, mat_type='cds', debug=nldebug)
!
! Set BC on grid transfer matrices
!
IF(prb.EQ.'dddd') CALL ibc_transf(grids,1,3) ! Direction X
CALL ibc_transf(grids,2,3) ! Direction Y
!
! Define RHS at l=1, compute RHS at l=2,...,levels by "restriction".
!
grids(1)%f(:,:) = frhs(grids(1)%x,grids(1)%y)
DO l=2,levels
grids(l)%f = restrict(grids(l)%matp, grids(l-1)%f)
grids(l)%f = 0.25d0*grids(l)%f ! Scaling for FD
END DO
!
! Define SOL at l=levels, compute SOL at l=levels-1,..,1 by "prolongation"
!
grids(levels)%v(:,:) = fsol(grids(levels)%x,grids(levels)%y)
DO l=levels-1,1,-1
grids(l)%v = prolong(grids(l+1)%matp, grids(l+1)%v)
END DO
!
IF(nldebug) THEN
DO l=1,levels
WRITE(*,'(a,i3)') '==== Level', l
WRITE(*,'(a)') 'f ='
DO i=0,grids(l)%n(1)
WRITE(*,'(10f8.3)') grids(l)%f(i,:)
END DO
WRITE(*,'(a)') 'v ='
DO i=0,grids(l)%n(1)
WRITE(*,'(10f8.3)') grids(l)%v(i,:)
END DO
END DO
END IF
!
! Epilogue
!
CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
!+++
FUNCTION fdense(x)
!
! Return density
!
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: fdense(SIZE(x))
fdense(:) = 0.5d0*beta*miome*EXP( -(x(:)-Lx/3d0)**2d0/(Lx/2d0)**2d0 );
END FUNCTION fdense
!+++
FUNCTION frhs(x,y)
!
! Return RHS
!
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: frhs(SIZE(x),SIZE(y))
DOUBLE PRECISION :: c, s, d(SIZE(x))
DOUBLE PRECISION :: corr
INTEGER :: j
corr = 1.d0+icrosst**2/4.0d0
d(:) = fdense(x(:))
IF(prb.EQ.'dddd') THEN
DO j=1,SIZE(y)
c = COS(2.0d0*pi*ky*y(j)/Ly)
s = SIN(2.0d0*pi*ky*y(j)/Ly)
frhs(:,j) = -4.0d0*pi*pi*( (kx*kx/Lx/Lx+corr*ky*ky/Ly/Ly) * SIN(2.0d0*pi*kx*x(:)/Lx)*s &
& -icrosst*kx*ky/Lx/Ly * COS(2.0d0*pi*kx*x(:)/Lx)*c ) &
& + d(:) * SIN(2.0d0*pi*kx*x(:)/Lx)*s
END DO
ELSE IF (prb.EQ.'nndd') THEN
DO j=1,SIZE(y)
c = COS(2.0d0*pi*ky*y(j)/Ly)
s = SIN(2.0d0*pi*ky*y(j)/Ly)
frhs(:,j) = -4.0d0*pi*pi*( (kx*kx/Lx/Lx+corr*ky*ky/Ly/Ly) * COS(2.0d0*pi*kx*x(:)/Lx)*s &
& +icrosst*kx*ky/Lx/Ly * SIN(2.0d0*pi*kx*x(:)/Lx)*c ) &
& + d(:) * COS(2.0d0*pi*kx*x(:)/Lx)*s
END DO
END IF
END FUNCTION frhs
!+++
FUNCTION fsol(x,y)
!
! Return analytical solution
!
DOUBLE PRECISION, INTENT(in) :: x(:), y(:)
DOUBLE PRECISION :: fsol(SIZE(x),SIZE(y))
DOUBLE PRECISION :: c
INTEGER :: j
IF(prb.EQ.'dddd') THEN
DO j=1,SIZE(y)
c = SIN(2.0d0*pi*ky*y(j)/Ly)
fsol(:,j) = SIN(2.0d0*pi*kx*x(:)/Lx)*c
END DO
ELSE IF (prb.EQ.'nndd') THEN
DO j=1,SIZE(y)
c = SIN(2.0d0*pi*ky*y(j)/Ly)
fsol(:,j) = COS(2.0d0*pi*kx*x(:)/Lx)*c
END DO
END IF
END FUNCTION fsol
!+++
SUBROUTINE h5file
USE futils
USE csr, ONLY : putmat
CHARACTER(len=128) :: file='test_intergrid0.h5'
INTEGER :: fid
INTEGER :: l
CHARACTER(len=64) :: dsname
CALL creatf(file, fid, real_prec='d')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NY', ny)
CALL attach(fid, '/', 'KX', kx)
CALL attach(fid, '/', 'KY', ky)
CALL attach(fid, '/', 'LX', Lx)
CALL attach(fid, '/', 'LY', Ly)
CALL attach(fid, '/', 'BETA', beta)
CALL attach(fid, '/', 'LEVELS', levels)
CALL attach(fid, '/', 'PRB', prb)
CALL attach(fid, '/', 'NLDEBUG', nldebug)
CALL creatg(fid, '/mglevels')
DO l=1,levels
WRITE(dsname,'("/mglevels/level.",i2.2)') l
CALL creatg(fid, TRIM(dsname))
IF(l.GT.1) THEN
CALL putmat(fid, TRIM(dsname)//'/matpx', grids(l)%matp(1))
CALL putmat(fid, TRIM(dsname)//'/matpy', grids(l)%matp(2))
END IF
CALL putarr(fid, TRIM(dsname)//'/x', grids(l)%x)
CALL putarr(fid, TRIM(dsname)//'/y', grids(l)%y)
CALL putarr(fid, TRIM(dsname)//'/f', grids(l)%f)
CALL putarr(fid, TRIM(dsname)//'/v', grids(l)%v)
END DO
CALL closef(fid)
END SUBROUTINE h5file
!+++
END PROGRAM main
c4science · Help