Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/examples/disrhs.f90
03540012fedefeature/fft-split

Last Change
Raw File

disrhs.f90
View Options

!>
!> @file disrhs.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>
!>
SUBROUTINE disrhs(mbess, spl, rhs)
!
! Assembly the RHS using 2d spline spl
!
USE bsplines
IMPLICIT NONE
INTEGER, INTENT(in) :: mbess
TYPE(spline2d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
INTEGER :: n1, nidbas1, ndim1, ng1
INTEGER :: n2, nidbas2, ndim2, ng2
INTEGER :: i, j, ig1, ig2, k1, k2, i1, j2, ij, nrank
DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:)
DOUBLE PRECISION:: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
CALL get_dim(spl%sp2, ndim2, n2, nidbas2)
!
ALLOCATE(fun1(0:nidbas1,1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun2(0:nidbas2,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(spl%sp1, ng1)
CALL get_gauss(spl%sp2, ng2)
WRITE(*,'(/a, 2i3)') 'Gauss points and weights, ngauss =', ng1, ng2
ALLOCATE(xg1(ng1), wg1(ng1))
ALLOCATE(xg2(ng1), wg2(ng1))
!===========================================================================
! 2.0 Assembly loop
!
nrank = SIZE(rhs)
rhs(1:nrank) = 0.0d0
!
DO i=1,n1
CALL get_gauss(spl%sp1, ng1, i, xg1, wg1)
DO ig1=1,ng1
CALL basfun(xg1(ig1), spl%sp1, fun1, i)
DO j=1,n2
CALL get_gauss(spl%sp2, ng2, j, xg2, wg2)
DO ig2=1,ng2
CALL basfun(xg2(ig2), spl%sp2, fun2, j)
contrib = wg1(ig1)*wg2(ig2) * &
& rhseq(xg1(ig1),xg2(ig2), mbess)
DO k1=0,nidbas1
i1 = i+k1
DO k2=0,nidbas2
j2 = MODULO(j+k2-1,n2) + 1
ij = j2 + (i1-1)*n2
rhs(ij) = rhs(ij) + contrib*fun1(k1,1)*fun2(k2,1)
END DO
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(xg1, wg1, fun1)
DEALLOCATE(xg2, wg2, fun2)
!
CONTAINS
DOUBLE PRECISION FUNCTION rhseq(x1, x2, m)
DOUBLE PRECISION, INTENT(in) :: x1, x2
INTEGER, INTENT(in) :: m
rhseq = REAL(4*(m+1),8)*x1**(m+1)*COS(REAL(m,8)*x2)
END FUNCTION rhseq
END SUBROUTINE disrhs
!
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE disrhs3(mbess, npow, spl, rhs)
!
! Assembly the RHS using 3d spline spl
!
USE bsplines
IMPLICIT NONE
INTEGER, INTENT(in) :: mbess, npow
TYPE(spline2d1d), TARGET :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:,:)
!
TYPE(spline1d), POINTER :: sp1, sp2, sp3
INTEGER :: n1, nidbas1, ndim1, ng1
INTEGER :: n2, nidbas2, ndim2, ng2
INTEGER :: n3, nidbas3, ndim3, ng3
INTEGER :: i, j, k, ig1, ig2, ig3, k1, k2, k3, i1, j2, ij, kk, nrank
DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xg3(:), wg3(:), fun3(:,:)
DOUBLE PRECISION:: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
sp1 => spl%sp12%sp1
sp2 => spl%sp12%sp2
sp3 => spl%sp3
!
CALL get_dim(sp1, ndim1, n1, nidbas1)
CALL get_dim(sp2, ndim2, n2, nidbas2)
CALL get_dim(sp3, ndim3, n3, nidbas3)
!
ALLOCATE(fun1(0:nidbas1,1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun2(0:nidbas2,1)) ! needs only basis functions (no derivatives)
ALLOCATE(fun3(0:nidbas3,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(sp1, ng1)
CALL get_gauss(sp2, ng2)
CALL get_gauss(sp3, ng3)
WRITE(*,'(/a, 3i3)') 'Gauss points and weights, ngauss =', ng1, ng2, ng3
ALLOCATE(xg1(ng1), wg1(ng1))
ALLOCATE(xg2(ng2), wg2(ng2))
ALLOCATE(xg3(ng3), wg3(ng3))
!===========================================================================
! 2.0 Assembly loop
!
nrank = SIZE(rhs,1)
rhs(1:nrank,1:n3) = 0.0d0
!
DO i=1,n1
CALL get_gauss(sp1, ng1, i, xg1, wg1)
DO ig1=1,ng1
CALL basfun(xg1(ig1), sp1, fun1, i)
DO j=1,n2
CALL get_gauss(sp2, ng2, j, xg2, wg2)
DO ig2=1,ng2
CALL basfun(xg2(ig2), sp2, fun2, j)
DO k=1,n3
CALL get_gauss(sp3, ng3, k, xg3, wg3)
DO ig3=1,ng3
CALL basfun(xg3(ig3), sp3, fun3, k)
contrib = wg1(ig1)*wg2(ig2)*wg3(ig3) * &
& rhseq(xg1(ig1), xg2(ig2), xg3(ig3), mbess, npow)
DO k1=0,nidbas1
i1 = i+k1
DO k2=0,nidbas2
j2 = MODULO(j+k2-1,n2) + 1
ij = j2 + (i1-1)*n2
DO k3=0,nidbas3
kk = MODULO(k+k3-1,n3) + 1
rhs(ij,kk) = rhs(ij, kk) + &
& contrib*fun1(k1,1)*fun2(k2,1)*fun3(k3,1)
END DO
END DO
END DO
END DO
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(xg1, wg1, fun1)
DEALLOCATE(xg2, wg2, fun2)
DEALLOCATE(xg3, wg3, fun3)
!
CONTAINS
DOUBLE PRECISION FUNCTION rhseq(x1, x2, x3, m, n)
DOUBLE PRECISION, INTENT(in) :: x1, x2, x3
INTEGER, INTENT(in) :: m, n
rhseq = REAL(4*(m+1),8)*x1**(m+1)*COS(REAL(m,8)*x2)*COS(x3)**n
END FUNCTION rhseq
END SUBROUTINE disrhs3
c4science · Help