Homec4science
Diffusion SPClibs (feature/refactoring) Browse

spclibs/multigrid/src/transfer1d.f90
c723281d03b6feature/refactoring

Last Change
Raw File

transfer1d.f90
View Options

!>
!> @file transfer1d.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>
!>
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
!--------------------------------------------------------------------------------
PROGRAM main
USE multigrid
IMPLICIT NONE
!
INTEGER :: nx=8, nidbas=1, ngauss=4, alpha=0, modem=0
DOUBLE PRECISION :: sigma=1.0d0
LOGICAL :: nlper=.FALSE.
INTEGER :: j
!
TYPE(grid1d) :: gridx(2)
TYPE(gemat) :: prolong_mat, restrict_mat, coarse_mat
DOUBLE PRECISION, ALLOCATABLE :: arow(:), temp(:,:)
!
NAMELIST /newrun/ nx, nidbas, ngauss, sigma, alpha, modem, nlper
!--------------------------------------------------------------------------------
READ(*,newrun)
WRITE(*,newrun)
!
! Set up fine and coarse grids
!
CALL create_grid(nx, nidbas, ngauss, alpha, gridx, period=nlper)
CALL printmat('** Prolongation matrix **', gridx(2)%transf)
!
! Restriction matrix = transpose of prolongation matrix
!
CALL mcopy(gridx(2)%transf, prolong_mat)
CALL init(prolong_mat%mrows, 1, restrict_mat, mrows=prolong_mat%ncols)
restrict_mat%val = TRANSPOSE(prolong_mat%val)
!
! Compute femat on fine and coarse grids
!
IF(nlper) THEN
CALL femat(gridx(1)%spl, gridx(1)%matap, coefeq)
CALL printmat('** FE matrix on fine mesh **', gridx(1)%matap)
CALL femat(gridx(2)%spl, gridx(2)%matap, coefeq)
CALL printmat('** FE matrix on coarse mesh **', gridx(2)%matap)
ELSE
CALL femat(gridx(1)%spl, gridx(1)%mata, coefeq)
CALL printmat('** FE matrix on fine mesh **', gridx(1)%mata)
CALL femat(gridx(2)%spl, gridx(2)%mata, coefeq)
CALL printmat('** FE matrix on coarse mesh **', gridx(2)%mata)
END IF
!
! Compute coarse FE matrix using transfer matrix
!
IF(nlper) THEN
CALL init(gridx(2)%matap%rank, 1, coarse_mat)
ALLOCATE(temp(gridx(1)%matap%rank,gridx(2)%matap%rank))
DO j=1,gridx(2)%matap%rank
temp(:,j) = vmx(gridx(1)%matap,prolong_mat%val(:,j))
END DO
coarse_mat%val = vmx(restrict_mat,temp)
DEALLOCATE(temp)
ELSE
CALL init(gridx(2)%mata%rank, 1, coarse_mat)
coarse_mat%val = vmx(restrict_mat,vmx(gridx(1)%mata,prolong_mat%val))
END IF
CALL printmat('** Coarse FE matrix using transfer operators **', coarse_mat)
!
! Compute the diff of Ac - R*Af*P
!
IF(nlper) THEN
coarse_mat%val = coarse_mat%val - gridx(2)%matap%val
ELSE
ALLOCATE(arow(gridx(2)%mata%rank))
DO j=1,gridx(2)%mata%rank
CALL getcol(gridx(2)%mata, j, arow)
coarse_mat%val(:,j) = coarse_mat%val(:,j)-arow(:)
END DO
DEALLOCATE(arow)
END IF
WRITE(*,'(a,1pe12.3)') 'Diff =', MAXVAL(ABS(coarse_mat%val))
CONTAINS
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(:)
SELECT CASE (alpha)
CASE(0) ! Cartesian geometry
c(1) = 1.0d0
idt(1) = 1
idw(1) = 1
c(2) = sigma
idt(2) = 0
idw(2) = 0
CASE(1)
c(1) = x
idt(1) = 1
idw(1) = 1
c(2) = modem**2/x
idt(2) = 0
idw(2) = 0
CASE default
WRITE(*,'(a,i0,a)') 'COEFEQ: alpha ', alpha, ' not defined!'
END SELECT
END SUBROUTINE coefeq
!--------------------------------------------------------------------------------
END PROGRAM main
c4science · Help