Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/examples/pde1dp_mod.f90
03540012fedefeature/refactoring

Last Change
Raw File

pde1dp_mod.f90
View Options

!>
!> @file pde1dp_mod.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>
!>
MODULE pde1dp_mod
!
USE bsplines
USE matrix
IMPLICIT NONE
DOUBLE PRECISION, ALLOCATABLE :: bcoef(:)
TYPE(spline1d), SAVE :: splx
!
CONTAINS
SUBROUTINE meshdist(c, x, nx)
!
! Construct an 1d non-equidistant mesh given a
! mesh distribution function defined in FDIST
!
DOUBLE PRECISION, INTENT(in) :: c(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
INTEGER :: nintg
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
DOUBLE PRECISION :: a, b, dx, f0, f1, scal
INTEGER :: i, k
!
a=x(0)
b=x(nx)
nintg = 10*nx
ALLOCATE(xint(0:nintg), fint(0:nintg))
!
! Mesh distribution
!
dx = (b-a)/REAL(nintg)
xint(0) = a
fint(0) = 0.0d0
f1 = fdist(xint(0))
DO i=1,nintg
f0 = f1
xint(i) = xint(i-1) + dx
f1 = fdist(xint(i))
fint(i) = fint(i-1) + 0.5*(f0+f1)
END DO
!
! Normalization
!
scal = REAL(nx) / fint(nintg)
fint(0:nintg) = fint(0:nintg) * scal
!
! Obtain mesh point by (inverse) interpolation
!
k = 1
DO i=1,nintg-1
IF( fint(i) .GE. REAL(k) ) THEN
x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
& (k-fint(i))
k = k+1
END IF
END DO
!
DEALLOCATE(xint, fint)
CONTAINS
DOUBLE PRECISION FUNCTION fdist(x)
DOUBLE PRECISION, INTENT(in) :: x
fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
END FUNCTION fdist
END SUBROUTINE meshdist
!+++
SUBROUTINE dismat(spl, mat)
!
! Assembly FE matrix (with periodic BC) mat using spline spl
!
TYPE(spline1d), INTENT(in) :: spl
TYPE(periodic_mat), INTENT(inout) :: mat
INTEGER :: dim, 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, dim, nx, nidbas)
ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
! Weak form
!
kterms = mat%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=MODULO(i+iw-1,nx) + 1 ! Periodic BC
jcol=MODULO(i+jt-1,nx) + 1
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)
END SUBROUTINE dismat
!
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.0d0
idt(1) = 0
idw(1) = 0
END SUBROUTINE coefeq
!+++
SUBROUTINE disrhs(spl, rhs)
!
! Assenbly the RHS using spline spl
!
TYPE(spline1d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
INTEGER :: dim, nrank, nx, nidbas, ngauss
INTEGER :: i, igauss, it, irow
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
DOUBLE PRECISION:: contrib
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
nrank = SIZE(rhs)
CALL get_dim(spl, dim, nx, nidbas)
!
ALLOCATE(fun(0:nidbas,1)) ! needs only basis functions (no derivatives)
!
! Gauss quadature
!
CALL get_gauss(spl, ngauss)
ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
! 2.0 Assembly loop
!
rhs(:) = 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))
DO it=0,nidbas
irow=MODULO(i+it-1,nx) + 1 ! Periodic BC
rhs(irow) = rhs(irow) + contrib*fun(it,1)
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(fun)
DEALLOCATE(xgauss, wgauss)
END SUBROUTINE disrhs
!
DOUBLE PRECISION FUNCTION rhseq(x)
DOUBLE PRECISION, INTENT(in) :: x
DOUBLE PRECISION :: xarr(1), farr(1)
INTEGER, SAVE :: icall =0
xarr(1) = x
IF( icall.EQ.0 ) THEN
icall = icall+1
CALL gridval(splx, xarr, farr, 0, bcoef)
ELSE
CALL gridval(splx, xarr, farr, 0)
END IF
rhseq = farr(1)
END FUNCTION rhseq
END MODULE pde1dp_mod
c4science · Help