Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

test_gvec1d.f90
View Options

!>
!> @file test_gvec1d.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>
!>
!
! Test implementation of 1D vectors with arbirary
! vector bounds and ghost cell width.
!
! T.M. Tran (09/2013)
!
MODULE gvector
USE iso_fortran_env, ONLY : rkind => real64
IMPLICIT NONE
PRIVATE
PUBLIC :: rkind, gvector_1d, disp, norm2, &
& OPERATOR(+), OPERATOR(-), OPERATOR(*), &
& ASSIGNMENT(=)
TYPE gvector_1d
INTEGER :: s, e, g
REAL(rkind), ALLOCATABLE :: val(:)
END TYPE gvector_1d
INTERFACE gvector_1d
MODULE PROCEDURE constructor
END INTERFACE gvector_1d
INTERFACE OPERATOR(+)
MODULE PROCEDURE add_scal
MODULE PROCEDURE add_vec
END INTERFACE OPERATOR(+)
INTERFACE OPERATOR(-)
MODULE PROCEDURE minus_vec
MODULE PROCEDURE substract_vec
END INTERFACE OPERATOR(-)
INTERFACE OPERATOR(*)
MODULE PROCEDURE scale_left
MODULE PROCEDURE scale_right
END INTERFACE OPERATOR(*)
INTERFACE ASSIGNMENT(=)
MODULE PROCEDURE from_scal
MODULE PROCEDURE from_vec
END INTERFACE ASSIGNMENT(=)
INTERFACE norm2
module procedure norm2_gvector_1d
END INTERFACE norm2
CONTAINS
FUNCTION constructor(s, e, g) RESULT(res)
INTEGER, INTENT(in) :: s, e
INTEGER, OPTIONAL, INTENT(in) :: g
TYPE(gvector_1d) :: res
INTEGER :: lb, ub
res%g= 0
IF(PRESENT(g)) res%g=g
res%s=s
res%e=e
lb = res%s-res%g
ub = res%e+res%g
ALLOCATE(res%val(lb:ub))
res%val = -9999.0
END FUNCTION constructor
FUNCTION add_vec(lhs, rhs) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: lhs, rhs
TYPE(gvector_1d) :: res
res = gvector_1d(lhs%s, lhs%e, lhs%g)
res%val(res%s:res%e) = lhs%val(res%s:res%e) + rhs%val(res%s:res%e)
END FUNCTION add_vec
FUNCTION add_scal(lhs, rhs) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: lhs
REAL(rkind), INTENT(in) :: rhs
TYPE(gvector_1d) :: res
res = gvector_1d(lhs%s, lhs%e, lhs%g)
res%val(res%s:res%e) = lhs%val(res%s:res%e) + rhs
END FUNCTION add_scal
FUNCTION minus_vec(this) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: this
TYPE(gvector_1d) :: res
res = gvector_1d(this%s, this%e, this%g)
res%val(res%s:res%e) = -this%val(res%s:res%e)
END FUNCTION minus_vec
FUNCTION substract_vec(lhs, rhs) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: lhs, rhs
TYPE(gvector_1d) :: res
res = gvector_1d(lhs%s, lhs%e, lhs%g)
res = lhs + (-rhs)
END FUNCTION substract_vec
FUNCTION scale_left(lhs, rhs) RESULT(res)
REAL(rkind), INTENT(in) :: lhs
TYPE(gvector_1d), INTENT(in) :: rhs
TYPE(gvector_1d) :: res
res = gvector_1d(rhs%s, rhs%e, rhs%g)
res%val(res%s:res%e) = lhs * rhs%val(res%s:res%e)
END FUNCTION scale_left
FUNCTION scale_right(lhs, rhs) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: lhs
REAL(rkind), INTENT(in) :: rhs
TYPE(gvector_1d) :: res
res = gvector_1d(lhs%s, lhs%e, lhs%g)
res%val(res%s:res%e) = rhs * lhs%val(res%s:res%e)
END FUNCTION scale_right
SUBROUTINE from_vec(lhs, rhs)
TYPE(gvector_1d), INTENT(inout) :: lhs
REAL(rkind), INTENT(in) :: rhs(:)
INTEGER :: n
n = lhs%e - lhs%s + 1
IF(SIZE(rhs) .NE. n) THEN
PRINT*, 'from_vec: sizes of rhs and lhs not equal!'
STOP
END IF
lhs%val(lhs%s:lhs%e) = rhs(1:n)
END SUBROUTINE from_vec
SUBROUTINE from_scal(lhs, rhs)
TYPE(gvector_1d), INTENT(inout) :: lhs
REAL(rkind), INTENT(in) :: rhs
lhs%val(lhs%s:lhs%e) = rhs
END SUBROUTINE from_scal
SUBROUTINE disp(str,this)
CHARACTER(len=*), INTENT(in) :: str
TYPE(gvector_1d), INTENT(in) :: this
WRITE(*,'(/a,3i6)') str//': s, e, g =', this%s, this%e, this%g
WRITE(*,'(10(1pe12.3))') this%val
END SUBROUTINE disp
FUNCTION norm2_gvector_1d(this) RESULT(res)
TYPE(gvector_1d), INTENT(in) :: this
REAL(rkind) :: res
res = NORM2(this%val(this%s:this%e))
END FUNCTION norm2_gvector_1d
END MODULE gvector
PROGRAM main
USE gvector
IMPLICIT NONE
INTEGER :: s=0, e=5, g=1
INTEGER :: i, lb, ub
REAL(rkind) :: a=0.1
TYPE(gvector_1d) :: v1, v2, v3
!
lb = s-g
ub = e+g
v1 = gvector_1d(s, e, g)
v1%val(s:e) = [ (i, i=s,e) ]
CALL disp('v1', v1)
!
v2 = v1 + a*v1
CALL disp('v1+a*v1', v2)
!
v3 = v1 - v1*a
CALL disp('v1-v1*a', v3)
!
WRITE(*,'(a,1pe12.3)') 'norm of v1 =', NORM2(v1)
WRITE(*,'(a,1pe12.3)') 'norm of v1-a*v1 =', NORM2(v1-a*v1)
!
v1 = 0.0d0
CALL disp('Should be all zero', v1)
v2 = [ 1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0 ]
CALL disp('Should be (1. 2. 3. 4. 5. 6.)', v2)
END PROGRAM main
c4science · Help