Homec4science
Diffusion SPClibs (master) Browse

spclibs/examples/moments.f90
253e659ccb68master

Last Change
Raw File

moments.f90
View Options

!>
!> @file moments.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>
!>
!
! Compute moments of f(x), using its Spline representation
!
MODULE globals
USE bsplines
USE matrix
IMPLICIT NONE
DOUBLE PRECISION, PARAMETER :: pi = 3.14159265359d0
DOUBLE PRECISION, ALLOCATABLE :: xgrid(:), rhs(:), sol(:)
DOUBLE PRECISION, ALLOCATABLE :: finteg(:,:), moms(:)
TYPE(spline1d), SAVE :: splx
TYPE(gbmat), SAVE :: mat
CONTAINS
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE dismat(spl, mat)
!
! Assembly FE matrix mat using spline spl
!
TYPE(spline1d), INTENT(in) :: spl
TYPE(gbmat), INTENT(inout) :: mat
INTEGER :: nrank, nx, nidbas, ngauss
INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION :: contrib
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER, ALLOCATABLE :: idert(:), iderw(:) ! Derivative order
DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, nrank, nx, nidbas)
ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
CALL coefeq(xgauss(igauss), idert, iderw, coefs)
DO iterm=1,kterms
DO jt=0,nidbas
DO iw=0,nidbas
contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
& fun(iw,iderw(iterm)) * wgauss(igauss)
irow=i+iw; jcol=i+jt
CALL updtmat(mat, irow, jcol, contrib)
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
DEALLOCATE(iderw, idert, coefs)
!
CONTAINS
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Mass matrix
!
c(1) = 1.d0
idt(1) = 0
idw(1) = 0
END SUBROUTINE coefeq
END SUBROUTINE dismat
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE disrhs(spl, rhs)
!
! Assenbly the RHS using spline spl
!
TYPE(spline1d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
INTEGER :: nrank, nx, nidbas, ngauss
INTEGER :: i, igauss
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION:: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl, nrank, nx, nidbas)
!
ALLOCATE(fun(nidbas+1,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
rhs(1:nrank) = 0.0d0
!
DO i=1,nx
CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
DO igauss=1,ngauss
CALL basfun(xgauss(igauss), spl, fun, i)
contrib = wgauss(igauss)*rhseq(xgauss(igauss))
rhs(i:i+nidbas) = rhs(i:i+nidbas) + contrib*fun(:,1)
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
!
CONTAINS
DOUBLE PRECISION FUNCTION rhseq(v)
DOUBLE PRECISION, INTENT(in) :: v
rhseq = SQRT(1.0d0/(2.0d0*pi)) * EXP(-0.5d0*v**2)
END FUNCTION rhseq
END SUBROUTINE disrhs
END MODULE globals
!
!================================================================================
PROGRAM main
USE bsplines
USE globals
IMPLICIT NONE
INTEGER :: nidbas, nx, nmoms, ngauss, rank, kl, ku
INTEGER :: i
DOUBLE PRECISION :: a, b, dx
!
! Input
!
WRITE(*,'(a)') 'Enter, nidbas, a, b, nx, nmoms'
READ(*,*) nidbas, a, b, nx, nmoms
!
! Equidistant mesh
!
ALLOCATE(xgrid(0:nx))
dx = (b-a)/REAL(nx)
DO i=0, nx
xgrid(i) = a + i*dx
END DO
WRITE(*,'(a/(8(1pe12.4)))') 'XGRID', xgrid
!
! Set up spline
!
ngauss = nidbas+1
CALL set_spline(nidbas, ngauss, xgrid, splx)
!
! Mass matrix
!
CALL get_dim(splx, rank) ! Rank of the FE Mass matrix
kl = nidbas
ku = kl
CALL init(kl, ku, rank, 1, mat)
WRITE(*,'(a,3i6)') 'kl, ku, rank', kl, ku, rank
CALL dismat(splx, mat)
!
! RHS
!
ALLOCATE(rhs(rank), sol(rank))
CALL disrhs(splx, rhs)
!
! Solve for Spline coefs
!
CALL factor(mat)
CALL bsolve(mat, rhs, sol)
WRITE(*,'(a/(8(1pe12.4)))') 'SOL', sol
WRITE(*,'(a,1pe20.12)') ' Integral of sol using FINTG', fintg(splx, sol)
!
! Moments
!
ALLOCATE(finteg(rank,0:nmoms), moms(0:nmoms))
CALL calc_integ(splx, finteg)
DO i=0,nmoms
moms(i) = DOT_PRODUCT(sol(:), finteg(:,i))
END DO
WRITE(*,'(a,i3)') 'Moments of orders from 0 to', nmoms
WRITE(*,'(8(1pe20.12))') moms
!
DEALLOCATE(finteg, moms)
DEALLOCATE(rhs, sol)
DEALLOCATE(xgrid)
CALL destroy(mat)
CALL destroy_sp(splx)
END PROGRAM main
c4science · Help