Homec4science
Diffusion SPClibs (master) Browse

spclibs/multigrid/src/two_grid.f90
03540012fedemaster

Last Change
Raw File

two_grid.f90
View Options

!>
!> @file two_grid.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
!
! Check some properties of grid transfer
!
USE multigrid
USE math_util, ONLY : root_bessj
IMPLICIT NONE
!
INTEGER :: nx=8, nidbas=1, ngauss=2, alpha=0
INTEGER :: modem=22, modep=10
INTEGER :: levels=2
INTEGER :: l, nrank
DOUBLE PRECISION :: sigma=1.0d0, kmode=10.0, pi=4.0d0*ATAN(1.0d0)
DOUBLE PRECISION, ALLOCATABLE :: v_prolong(:)
!
TYPE(grid1d) :: gridx(2)
!
NAMELIST /newrun/ nx, nidbas, ngauss, sigma, kmode, modem, modep, alpha
!--------------------------------------------------------------------------------
! 1. Prologue
! Inputs
!
READ(*,newrun)
WRITE(*,newrun)
!
! Create grids
!
CALL create_grid(nx, nidbas, ngauss, alpha, gridx)
WRITE(*,'(a/(20i6))') 'Number of intervals in grids', (gridx(l)%n, l=1,levels)
!
! Create FE matrice and set BC u(0)=u(1)=0
!
DO l=1,levels
CALL femat(gridx(l)%spl, gridx(l)%mata, coefeq)
!
! Left Dirichlet BC (only for Cartesian geometry)
IF(alpha .EQ. 0) THEN
CALL ibcmat(1, gridx(l)%mata)
END IF
!
! Right Dirichlet BC
CALL ibcmat(gridx(l)%mata%rank, gridx(l)%mata)
!
! BC on grid transfer operator
IF(l.GT.1) THEN
WHERE( ABS(gridx(l)%transf%val) < 1.d-8) gridx(l)%transf%val=0.0d0
IF(alpha .EQ. 0) gridx(l)%transf%val(2:,1)=0.0d0
gridx(l)%transf%val(1:gridx(l-1)%rank-1,gridx(l)%rank)=0.0d0
END IF
END DO
!
! Construct RHS and set BC only on the finest grid
!
nrank = gridx(1)%rank
CALL disrhs(gridx(1)%spl, gridx(1)%f, rhs)
!
! Left Dirichlet BC (only for Cartesian geometry)
IF(alpha .EQ. 0) THEN
gridx(1)%f(1) = 0.0d0
END IF
!
! Right Dirichlet BC
gridx(1)%f(nrank) = 0.0d0
!
! RHS on coarse grid by restriction
!
gridx(2)%f = restrict(gridx(2)%transf,gridx(1)%f)
!--------------------------------------------------------------------------------
! 2. Direct solutions
!
DO l=1,levels
CALL direct_solve(gridx(l), gridx(l)%v)
WRITE(*,'(a,i3/(10(1pe12.3)))') 'Sol at level', l, gridx(l)%v
END DO
!
! Prolongation of coarse solution
!
ALLOCATE(v_prolong(SIZE(gridx(1)%v)))
!
v_prolong = prolong(gridx(2)%transf, gridx(2)%v)
WRITE(*,'(a,i3/(10(1pe12.3)))') 'Prolong. sol.', l, v_prolong
WRITE(*,'(a,1pe12.3)') 'Error ||V_prolong-V_fine||', normf(gridx(1)%matm, v_prolong-gridx(1)%v)
!--------------------------------------------------------------------------------
! 9. Epilogue
!
! Creata HDF5 file
!
CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
SUBROUTINE h5file
USE futils
CHARACTER(len=128) :: file='two_grid.h5'
INTEGER :: fid
INTEGER :: l
CHARACTER(len=64) :: dsname
CALL creatf(file, fid, real_prec='d')
CALL attach(fid, '/', 'NX', nx)
CALL attach(fid, '/', 'NIDBAS', nidbas)
CALL attach(fid, '/', 'SIGMA', sigma)
CALL attach(fid, '/', 'KMODE', kmode)
CALL attach(fid, '/', 'ALPHA', alpha)
CALL attach(fid, '/', 'LEVELS', levels)
CALL creatg(fid, '/mglevels')
DO l=1,levels
WRITE(dsname,'("/mglevels/level.",i2.2)') l
CALL creatg(fid, TRIM(dsname))
CALL putmat(fid, TRIM(dsname)//'/mata', gridx(l)%mata)
IF(l.GT.1) THEN
CALL putarr(fid, TRIM(dsname)//'/matp', gridx(l)%transf%val)
CALL attach(fid, TRIM(dsname)//'/matp', 'M', gridx(l)%transf%mrows)
CALL attach(fid, TRIM(dsname)//'/matp', 'N', gridx(l)%transf%ncols)
END IF
CALL putarr(fid, TRIM(dsname)//'/f', gridx(l)%f)
CALL putarr(fid, TRIM(dsname)//'/v', gridx(l)%v)
END DO
CALL putarr(fid, '/v_prolong', v_prolong)
END SUBROUTINE h5file
FUNCTION rhs(x)
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION :: rhs
DOUBLE PRECISION :: nump
SELECT CASE (alpha)
CASE(0) ! Cartesian geometry
rhs = SIN(pi*kmode*x)
CASE(1) ! Cylindrical
nump = root_bessj(modem, modep)
rhs = x * nump**2 * bessel_jn(modem, nump*x)
END SELECT
END FUNCTION rhs
FUNCTION sol(x)
DOUBLE PRECISION, INTENT(in) :: x(:)
DOUBLE PRECISION :: sol(SIZE(x))
DOUBLE PRECISION :: nump
SELECT CASE (alpha)
CASE(0) ! Cartesian geometry
sol(:) = 1.0d0/((pi*kmode)**2+sigma)*SIN(pi*kmode*x(:))
CASE(1) ! Cylindrical
nump = root_bessj(modem, modep)
sol(:) = bessel_jn(modem, nump*x(:))
END SELECT
END FUNCTION sol
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
SELECT CASE (alpha)
CASE(0) ! Cartesian geometry
c(1) = 1.0d0
idt(1) = 1
idw(1) = 1
c(2) = sigma
idt(2) = 0
idw(2) = 0
CASE(1) ! Cylindrical
c(1) = x
idt(1) = 1
idw(1) = 1
c(2) = REAL(modem,8)**2/x
idt(2) = 0
idw(2) = 0
END SELECT
END SUBROUTINE coefeq
END PROGRAM main
c4science · Help