Homec4science
Diffusion SPClibs (feature/fft-split) Browse

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

Last Change
Raw File

pde2d_pb.f90
View Options

!>
!> @file pde2d_pb.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>
!>
PROGRAM main
!
! Solving the following 2d PDE using splines:
!
! -d/dx[x d/dx]f - 1/x[d/dy]^2 f = 4(m+1)x^(m+1)cos(my), with f(x=1,y) = 0
! exact solution: f(x,y) = (1-x^2) x^m cos(my)
!
USE bsplines
USE matrix
USE conmat_mod
!
IMPLICIT NONE
INTEGER :: nx, ny, nidbas(2), ngauss(2), mbess, nterms
LOGICAL :: nlppform, nlconmat
INTEGER :: i, j, ij, dimx, dimy, nrank, kl, ku, jder(2), it
DOUBLE PRECISION :: pi, coefx(5), coefy(5)
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, ygrid, rhs, sol
TYPE(spline2d) :: splxy
TYPE(pbmat) :: mat
!
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: arr
DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: bcoef, solcal, solana, errsol
DOUBLE PRECISION :: seconds, mem, dopla
DOUBLE PRECISION :: t0, tmat, tfact, tsolv, tgrid, gflops1, gflops2
INTEGER :: nits=500
!
INTERFACE
SUBROUTINE dismat(spl, mat)
USE bsplines
USE matrix
TYPE(spline2d), INTENT(in) :: spl
TYPE(pbmat), INTENT(inout) :: mat
END SUBROUTINE dismat
SUBROUTINE disrhs(mbess, spl, rhs)
USE bsplines
INTEGER, INTENT(in) :: mbess
TYPE(spline2d), INTENT(in) :: spl
DOUBLE PRECISION, INTENT(out) :: rhs(:)
END SUBROUTINE disrhs
SUBROUTINE meshdist(coefs, x, nx)
DOUBLE PRECISION, INTENT(in) :: coefs(5)
INTEGER, INTENT(iN) :: nx
DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
END SUBROUTINE meshdist
SUBROUTINE ibcmat(mat, ny)
USE matrix
TYPE(pbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: ny
END SUBROUTINE ibcmat
SUBROUTINE ibcrhs(rhs, ny)
DOUBLE PRECISION, INTENT(inout) :: rhs(:)
INTEGER, INTENT(in) :: ny
END SUBROUTINE ibcrhs
SUBROUTINE coefeq(x, y, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:)
END SUBROUTINE coefeq
END INTERFACE
!
NAMELIST /newrun/ nx, ny, nidbas, ngauss, mbess, nlppform, nlconmat, &
& coefx, coefy
!===========================================================================
! 1.0 Prologue
!
! Read in data specific to run
!
nx = 8 ! Number of intervals in x
ny = 8 ! Number of intervals in y
nidbas = (/3,3/) ! Degree of splines
ngauss = (/4,4/) ! Number of Gauss points/interval
mbess = 2 ! Exponent of differential problem
nterms = 2 ! Number of terms in weak form
nlppform = .TRUE. ! Use PPFORM for gridval or not
nlconmat = .TRUE. ! Use CONMAT instead of DISMAT
coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
coefy(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
!
READ(*,newrun)
WRITE(*,newrun)
!
! Define grid on x (=radial) & y (=poloidal) axis
!
pi = 4.0d0*ATAN(1.0d0)
ALLOCATE(xgrid(0:nx), ygrid(0:ny))
xgrid(0) = 0.0d0; xgrid(nx) = 1.0d0
CALL meshdist(coefx, xgrid, nx)
ygrid(0) = 0.0d0; ygrid(ny) = 2.d0*pi
CALL meshdist(coefy, ygrid, ny)
!
WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
WRITE(*,'(a/(10f8.3))') 'YGRID', ygrid(0:ny)
!===========================================================================
! 2.0 Discretize the PDE
!
! Set up spline
!
CALL set_spline(nidbas, ngauss, &
& xgrid, ygrid, splxy, (/.FALSE., .TRUE./), nlppform=nlppform)
!!$ WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in X', splxy%sp1%knots
!!$ WRITE(*,'(/a,/(12(f8.3)))') 'KNOTS in Y', splxy%sp2%knots
!
! FE matrix assembly
!
nrank = (nx+nidbas(1))*ny ! Rank of the FE matrix
kl = (nidbas(1)+1)*ny -1 ! Number of sub-diagnonals
ku = kl ! Number of super-diagnonals
WRITE(*,'(a,5i6)') 'nrank, kl, ku', nrank, kl, ku
!
CALL init(ku, nrank, nterms, mat)
t0 = seconds()
IF(nlconmat) THEN
CALL conmat(splxy, mat, coefeq)
ELSE
CALL dismat(splxy, mat)
END IF
tmat = seconds() - t0
ALLOCATE(arr(nrank))
!
! BC on Matrix
!
IF(nrank.LT.100) &
& WRITE(*,'(a/(10(1pe12.3)))') 'Diag. of matrix before BC', mat%val(ku+1,:)
CALL ibcmat(mat, ny)
IF(nrank.LT.100) &
& WRITE(*,'(a/(10(1pe12.3)))') 'Diag. of matrix after BC', mat%val(ku+1,:)
!
! RHS assembly
!
ALLOCATE(rhs(nrank), sol(nrank))
CALL disrhs(mbess, splxy, rhs)
!
! BC on RHS
!
CALL ibcrhs(rhs, ny)
WRITE(*,'(a,1pe12.3)') 'Mem used so far (MB)', mem()
!===========================================================================
! 3.0 Solve the dicretized system
!
t0 = seconds()
CALL factor(mat)
tfact = seconds() - t0
gflops1 = dopla('DPBTRF',nrank,nrank,kl,ku,0) / tfact / 1.d9
t0 = seconds()
CALL bsolve(mat, rhs, sol)
!
! Backtransform of solution
!
sol(1:ny-1) = sol(ny)
!
! Spline coefficients, taking into account of periodicity in y
! Note: in SOL, y was numbered first.
!
dimx = splxy%sp1%dim
dimy = splxy%sp2%dim
ALLOCATE(bcoef(0:dimx-1, 0:dimy-1))
DO j=0,dimy-1
DO i=0,dimx-1
ij = MODULO(j,ny) + i*ny + 1
bcoef(i,j) = sol(ij)
END DO
END DO
WRITE(*,'(a,2(1pe12.3))') ' Integral of sol', fintg(splxy, bcoef)
!
tsolv = seconds() - t0
gflops2 = dopla('DPBTRS',nrank,1,kl,ku,0) / tsolv / 1.d9
!===========================================================================
! 4.0 Check the solution
!
! Check function values computed with various method
!
ALLOCATE(solcal(0:nx,0:ny), solana(0:nx,0:ny), errsol(0:nx,0:ny))
DO i=0,nx
DO j=0,ny
solana(i,j) = (1-xgrid(i)**2) * xgrid(i)**mbess * COS(mbess*ygrid(j))
END DO
END DO
jder = (/0,0/)
!
! Compute PPFORM at first call to gridval
IF(nlppform) THEN
CALL gridval(splxy, xgrid, ygrid, solcal, jder, bcoef)
END IF
!
WRITE(*,'(/a)') '*** Checking solutions'
t0 = seconds()
DO it=1,nits ! nits iterations for timing
CALL gridval(splxy, xgrid, ygrid, solcal, jder)
END DO
tgrid = (seconds() - t0)/REAL(nits)
errsol = solana - solcal
IF( SIZE(bcoef,2) .LE. 10 ) THEN
CALL prnmat('BCOEF', bcoef)
CALL prnmat('SOLANA', solana)
CALL prnmat('SOLCAL', solcal)
CALL prnmat('ERRSOL', errsol)
END IF
WRITE(*,'(a,2(1pe12.3))') 'Relative discretization errors', &
& norm2(errsol) / norm2(solana)
WRITE(*,'(a,1pe12.3)') 'GRIDVAL2 time (s) ', tgrid
!
! Check derivatives d/dx and d/dy
!
WRITE(*,'(/a)') '*** Checking gradient'
DO i=0,nx
DO j=0,ny
IF( mbess .EQ. 0 ) THEN
solana(i,j) = -2.0d0 * xgrid(i)
ELSE
solana(i,j) = (-(mbess+2)*xgrid(i)**2 + mbess) * &
& xgrid(i)**(mbess-1) * COS(mbess*ygrid(j))
END IF
END DO
END DO
!
jder = (/1,0/)
CALL gridval(splxy, xgrid, ygrid, solcal, jder)
errsol = solana - solcal
WRITE(*,'(a,2(1pe12.3))') 'Error in d/dx', norm2(errsol)
!
DO i=0,nx
DO j=0,ny
solana(i,j) = -mbess * (1-xgrid(i)**2) * xgrid(i)**mbess * SIN(mbess*ygrid(j))
END DO
END DO
!
jder = (/0,1/)
CALL gridval(splxy, xgrid, ygrid, solcal, jder)
errsol = solana - solcal
WRITE(*,'(a,2(1pe12.3))') 'Error in d/dy', norm2(errsol)
!===========================================================================
! 9.0 Epilogue
!
WRITE(*,'(/a)') '---'
WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s) ', tmat
WRITE(*,'(a,1pe12.3)') 'Matrice factorisation time (s)', tfact
WRITE(*,'(a,1pe12.3)') 'Backsolve time (s) ', tsolv
WRITE(*,'(a,2f10.3)') 'Factor/solve Gflop/s', gflops1, gflops2
!
DEALLOCATE(xgrid, rhs, sol)
DEALLOCATE(solcal, solana, errsol)
DEALLOCATE(bcoef)
DEALLOCATE(arr)
CALL destroy_sp(splxy)
CALL destroy(mat)
!
!===========================================================================
!
CONTAINS
SUBROUTINE prntmat(str, a)
DOUBLE PRECISION, DIMENSION(:,:) :: a
CHARACTER(len=*) :: str
INTEGER :: i
WRITE(*,'(a)') TRIM(str)
DO i=1,SIZE(a,1)
WRITE(*,'(10f8.1)') a(i,:)
END DO
END SUBROUTINE prntmat
FUNCTION norm2(x)
!
! Compute the 2-norm of array x
!
IMPLICIT NONE
DOUBLE PRECISION :: norm2
DOUBLE PRECISION, INTENT(in) :: x(:,:)
DOUBLE PRECISION :: sum2
INTEGER :: i, j
!
sum2 = 0.0d0
DO i=1,SIZE(x,1)
DO j=1,SIZE(x,2)
sum2 = sum2 + x(i,j)**2
END DO
END DO
norm2 = SQRT(sum2)
END FUNCTION norm2
SUBROUTINE prnmat(label, mat)
CHARACTER(len=*) :: label
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: mat
INTEGER :: i
WRITE(*,'(/a)') TRIM(label)
DO i=1,SIZE(mat,1)
WRITE(*,'(10(1pe12.3))') mat(i,:)
END DO
END SUBROUTINE prnmat
END PROGRAM main
!
!+++
SUBROUTINE meshdist(c, x, nx)
!
! Construct an 1d non-equidistant mesh given a
! mesh distribution function.
!
IMPLICIT NONE
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
!!$ WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
! 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 of FE matrix mat using spline spl
!
USE bsplines
USE matrix
IMPLICIT NONE
TYPE(spline2d), INTENT(in) :: spl
TYPE(pbmat), INTENT(inout) :: mat
!
INTEGER :: n1, nidbas1, ndim1, ng1
INTEGER :: n2, nidbas2, ndim2, ng2
INTEGER :: i, j, ig1, ig2
INTEGER :: iterm, iw1, iw2, igw1, igw2, it1, it2, igt1, igt2, irow, jcol
DOUBLE PRECISION, ALLOCATABLE :: xg1(:), wg1(:), fun1(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: xg2(:), wg2(:), fun2(:,:,:)
DOUBLE PRECISION:: contrib
!
INTEGER :: kterms ! Number of terms in weak form
INTEGER, ALLOCATABLE :: idert(:,:,:,:), iderw(:,:,:,:) ! Derivative order
DOUBLE PRECISION, ALLOCATABLE :: coefs(:,:,:) ! Terms in weak form
INTEGER, ALLOCATABLE :: left1(:), left2(:)
!===========================================================================
! 1.0 Prologue
!
! Properties of spline space
!
CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
CALL get_dim(spl%sp2, ndim2, n2, nidbas2)
WRITE(*,'(/a, 5i6)') 'n1, nidbas1, ndim1 =', n1, nidbas1, ndim1
WRITE(*,'(a, 5i6)') 'n2, nidbas2, ndim2 =', n2, nidbas2, ndim2
!
! Gauss quadature
!
CALL get_gauss(spl%sp1, ng1)
CALL get_gauss(spl%sp2, ng2)
ALLOCATE(xg1(ng1), wg1(ng1))
ALLOCATE(xg2(ng1), wg2(ng1))
!
! Weak form
!
kterms = mat%nterms
ALLOCATE(idert(kterms,2,ng1,ng2), iderw(kterms,2,ng1,ng2))
ALLOCATE(coefs(kterms,ng1,ng2))
!
ALLOCATE(fun1(0:nidbas1,0:1,ng1)) ! Spline and 1st derivative
ALLOCATE(fun2(0:nidbas2,0:1,ng2)) !
!===========================================================================
! 2.0 Assembly loop
!
ALLOCATE(left1(ng1))
ALLOCATE(left2(ng2))
DO i=1,n1
CALL get_gauss(spl%sp1, ng1, i, xg1, wg1)
left1 = i
CALL basfun(xg1, spl%sp1, fun1, left1)
DO j=1,n2
CALL get_gauss(spl%sp2, ng2, j, xg2, wg2)
left2 = j
CALL basfun(xg2, spl%sp2, fun2, left2)
!
DO ig1=1,ng1
DO ig2=1,ng2
CALL coefeq(xg1(ig1), xg2(ig2), &
& idert(:,:,ig1,ig2), &
& iderw(:,:,ig1,ig2), &
& coefs(:,ig1,ig2))
END DO
END DO
!
DO iw1=0,nidbas1 ! Weight function in dir 1
igw1 = i+iw1
DO iw2=0,nidbas2 ! Weight function in dir 2
igw2 = MODULO(j+iw2-1, n2) + 1
irow = igw2 + (igw1-1)*n2
DO it1=0,nidbas1 ! Test function in dir 1
igt1 = i+it1
DO it2=0,nidbas2 ! Test function in dir 2
igt2 = MODULO(j+it2-1, n2) + 1
jcol = igt2 + (igt1-1)*n2
!-------------
contrib = 0.0d0
DO ig1=1,ng1
DO ig2=1,ng2
DO iterm=1,kterms
contrib = contrib + &
& fun1(iw1,iderw(iterm,1,ig1,ig2),ig1) * &
& fun2(iw2,iderw(iterm,2,ig1,ig2),ig2) * &
& coefs(iterm,ig1,ig2) * &
& fun2(it2,idert(iterm,2,ig1,ig2),ig2) * &
& fun1(it1,idert(iterm,1,ig1,ig2),ig1) * &
& wg1(ig1) * wg2(ig2)
END DO
END DO
END DO
CALL updtmat(mat, irow, jcol, contrib)
!-------------
END DO
END DO
END DO
END DO
END DO
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(xg1, wg1, fun1)
DEALLOCATE(xg2, wg2, fun2)
DEALLOCATE(idert, iderw, coefs)
DEALLOCATE(left1,left2)
!
CONTAINS
SUBROUTINE coefeq(x, y, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Weak form = Int(x*dw/dx*dt/dx + 1/x*dw/dy*dt/dy)dxdy
!
c(1) = x !
idt(1,1) = 1
idt(1,2) = 0
idw(1,1) = 1
idw(1,2) = 0
!
c(2) = 1.d0/x
idt(2,1) = 0
idt(2,2) = 1
idw(2,1) = 0
idw(2,2) = 1
END SUBROUTINE coefeq
END SUBROUTINE dismat
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 ibcmat(mat, ny)
!
! Apply BC on matrix
!
USE matrix
IMPLICIT NONE
TYPE(pbmat), INTENT(inout) :: mat
INTEGER, INTENT(in) :: ny
INTEGER :: kl, ku, nrank, i, j
DOUBLE PRECISION, ALLOCATABLE :: zsum(:), arr(:)
INTEGER :: i0, ii
INTEGER :: i0_arr(ny)
!===========================================================================
! 1.0 Prologue
!
ku = mat%ku
kl = ku
nrank = mat%rank
ALLOCATE(zsum(nrank), arr(nrank))
!
i0 = nrank - ku
WRITE(*,'(a,i6)') 'Estimated i0', i0
DO i=1,ny
CALL getcol(mat, nrank-ny+i, arr)
DO ii=1,nrank
i0_arr(i)=ii
IF(arr(ii) .NE. 0.0d0) EXIT
END DO
END DO
!!$ WRITE(*,'(a/(10i6))') 'i0_arr', i0_arr
!
!===========================================================================
! 2.0 Unicity at the axis
!
! The vertical sum on the NY-th row
!
zsum = 0.0d0
DO i=1,ny
arr = 0.0d0
CALL getrow(mat, i, arr)
DO j=1,ny+ku
zsum(j) = zsum(j) + arr(j)
END DO
END DO
!
zsum(ny) = SUM(zsum(1:ny)) ! using symmetry
CALL putrow(mat, ny, zsum)
!
! The away operator
!
DO i = 1,ny-1
arr = 0.0d0; arr(i) = 1.0d0
CALL putrow(mat, i, arr)
END DO
!===========================================================================
! 3.0 Dirichlet on right boundary
!
DO i = nrank, nrank-ny+1, -1
CALL getcol(mat, i, arr)
arr = 0.0d0; arr(i) = 1.0d0
CALL putrow(mat, i, arr)
END DO
!===========================================================================
! 9.0 Epilogue
!
DEALLOCATE(zsum, arr)
!
END SUBROUTINE ibcmat
!+++
SUBROUTINE ibcrhs(rhs, ny)
!
! Apply BC on RHS
!
IMPLICIT NONE
DOUBLE PRECISION, INTENT(inout) :: rhs(:)
INTEGER, INTENT(in) :: ny
INTEGER :: nrank
DOUBLE PRECISION :: zsum
!===========================================================================
! 1.0 Prologue
!
nrank = SIZE(rhs,1)
!===========================================================================
! 2.0 Unicity at the axis
!
! The vertical sum
!
zsum = SUM(rhs(1:ny))
rhs(ny) = zsum
rhs(1:ny-1) = 0.0d0
!===========================================================================
! 3.0 Dirichlet on right boundary
!
rhs(nrank-ny+1:nrank) = 0.0d0
END SUBROUTINE ibcrhs
!++++
SUBROUTINE coefeq(x, y, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x, y
INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Weak form = Int(x*dw/dx*dt/dx + 1/x*dw/dy*dt/dy)dxdy
!
c(1) = x !
idt(1,1) = 1
idt(1,2) = 0
idw(1,1) = 1
idw(1,2) = 0
!
c(2) = 1.d0/x
idt(2,1) = 0
idt(2,2) = 1
idw(2,1) = 0
idw(2,2) = 1
END SUBROUTINE coefeq
c4science · Help