Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/multigrid/src/test_relax.f90
c723281d03b6feature/refactoring

Last Change
Raw File

test_relax.f90
View Options

!>
!> @file test_relax.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/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
! Test different relaxations
!
USE multigrid
USE math_util, ONLY : root_bessj
IMPLICIT NONE
!
INTEGER :: nx=8, nidbas=1, alpha=0, nits=40
INTEGER :: modem=22, modep=10
CHARACTER(len=4) :: relax='jac '
LOGICAL :: nlfixed = .FALSE.
DOUBLE PRECISION :: sigma=1.0d0, kmode=10.0, pi=4.0d0*ATAN(1.0d0)
DOUBLE PRECISION :: omega=2.0d0/3.0d0
INTEGER :: ngauss, i, nrank, its
DOUBLE PRECISION, ALLOCATABLE :: u_exact(:), u_direct(:), u_calc(:)
DOUBLE PRECISION, ALLOCATABLE :: sol_direct(:)
DOUBLE PRECISION, ALLOCATABLE :: sol_calc(:), err(:), resid(:), errdisc(:)
DOUBLE PRECISION :: errdisc_dir
!
TYPE(grid1d) :: gridx(1)
!
NAMELIST /newrun/ nx, nidbas, ngauss, sigma, kmode, modem, modep, &
& alpha, relax, omega, nits, nlfixed
!--------------------------------------------------------------------------------
! 1. Prologue: read input, construct matrix and RHS
READ(*,newrun)
WRITE(*,newrun)
!
! Set grid
!
CALL create_grid(nx, nidbas, ngauss, alpha, gridx)
!
! Create FE matrice and set BC u(0)=u(1)=0
!
CALL femat(gridx(1)%spl, gridx(1)%mata, coefeq)
nrank = gridx(1)%rank
!
! Left Dirichlet BC (only for Cartesian geometry)
IF(alpha .EQ. 0) THEN
CALL ibcmat(1, gridx(1)%mata)
END IF
!
! Right Dirichlet BC
CALL ibcmat(nrank, gridx(1)%mata)
!
! Construct RHS
!
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
!--------------------------------------------------------------------------------
! 2. Direct solution
!
! Direct solutions
!
ALLOCATE(sol_direct(nrank))
CALL direct_solve(gridx(1), sol_direct)
!
! Grid values
!
ALLOCATE(u_exact(0:nx))
ALLOCATE(u_direct(0:nx))
ALLOCATE(u_calc(0:nx))
!
u_exact = sol(gridx(1)%x)
CALL gridval(gridx(1)%spl, gridx(1)%x, u_direct, 0, sol_direct)
errdisc_dir = disc_err(gridx(1)%spl, sol_direct, sol)
WRITE(*,'(a,1pe12.3)') 'Discretization error', errdisc_dir
!--------------------------------------------------------------------------------
! 3. Relaxation
!
ALLOCATE(sol_calc(nrank))
ALLOCATE(err(0:nits))
ALLOCATE(errdisc(0:nits))
ALLOCATE(resid(0:nits))
!
! Initial guess
sol_calc(:) = 0.0d0
IF(nlfixed) THEN
sol_calc(:) = sol_direct(:)
END IF
err(0) = normf(gridx(1)%matm, sol_calc-sol_direct)
errdisc(0) = disc_err(gridx(1)%spl, sol_calc, sol)
resid(0) = residue(gridx(1)%mata, gridx(1)%f, sol_calc)
!
! Iterations
DO its=1,nits
SELECT CASE (TRIM(relax))
CASE('jac')
CALL jacobi(gridx(1)%mata, omega, 1, sol_calc, gridx(1)%f)
CASE('gs')
CALL gs(gridx(1)%mata, 1, sol_calc, gridx(1)%f)
END SELECT
err(its) = normf(gridx(1)%matm, sol_calc-sol_direct)
errdisc(its) = disc_err(gridx(1)%spl, sol_calc, sol)
resid(its) = residue(gridx(1)%mata, gridx(1)%f, sol_calc)
END DO
CALL gridval(gridx(1)%spl, gridx(1)%x, u_calc, 0, sol_calc)
!
WRITE(*,'(/a4,3a12)') 'its', 'error', 'residue', 'disc. err'
WRITE(*,'(i4,3(1pe12.3))') 0, err(0), resid(0), errdisc(0)
WRITE(*,'((i4,6(1pe12.3)))') (its, err(its), resid(its), errdisc(its), &
& err(its)/err(its-1), resid(its)/resid(its-1), &
& errdisc(its)/errdisc(its-1), its=1,nits,MAX(1,nits/10))
!
CALL h5file
!--------------------------------------------------------------------------------
CONTAINS
!+++
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
!+++
SUBROUTINE h5file
USE futils
CHARACTER(len=128) :: file='test_relax.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, '/', 'KX', kmode)
CALL attach(fid, '/', 'SIGMA', sigma)
CALL attach(fid, '/', 'ALPHA', alpha)
CALL attach(fid, '/', 'RELAX', relax)
CALL attach(fid, '/', 'OMEGA', omega)
CALL attach(fid, '/', 'NITS', nits)
CALL attach(fid, '/', 'MODEM', modem)
CALL attach(fid, '/', 'MODEP', modep)
!
! Solutions at finest grid
!
CALL creatg(fid, '/solutions')
CALL putarr(fid, '/solutions/xg', gridx(1)%x)
CALL putarr(fid, '/solutions/direct', u_direct)
CALL putarr(fid, '/solutions/anal', u_exact)
CALL putarr(fid, '/solutions/calc', u_calc)
!
CALL creatg(fid, '/relaxation')
CALL putarr(fid, '/relaxation/errdisc', errdisc)
CALL putarr(fid, '/relaxation/resid', resid)
!
! Store FE matrix
!
CALL putmat(fid, '/MATA', gridx(1)%mata)
!
CALL closef(fid)
END SUBROUTINE h5file
!+++
END PROGRAM main
c4science · Help