Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_testmumps) Browse

fennecs/src/geometry_mod.f90
865582e1b728geom_testmumps

Last Change
Raw File

geometry_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: geometry
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for handling geometries with non constant radius using b-splines interpolation
!> This module defines ways to comupte the weight function needed for weighted b-splines and
!> can load the definition of the geometry from input files
!------------------------------------------------------------------------------
MODULE geometry
USE constants
USE bsplines
USE mumps_bsplines
IMPLICIT NONE
type innerspline
Integer, Allocatable:: k(:) ! Index in reduced set
real(kind=db), Allocatable:: weight(:) ! geomtric weight at relevant cell
end type innerspline
Real(kind=db),ALLOCATABLE, SAVE ::IVand (:,:)
Real(kind=db),save:: z_0=0, r_0=0, invr_r=0, invr_z=0, r_a=0, r_b=0, z_r=1, r_r=1
Real(kind=db), save:: Phidown=0,Phiup=0
Real(kind=db), save:: L_r=0.16, L_z=0.24
Logical, save:: nlweb=.false.
Integer, save::Interior=-1
Integer, save:: above1=1, above2=-1
Integer, save :: walltype = 0
Integer, save, Allocatable:: bsplinetype(:) ! Array containing the inner/outer type for each bspline
Integer, save, Allocatable:: gridcelltype(:,:) ! Array containing the inner/outer type for each gridcell
Real(kind=db), save, allocatable:: gridwgeom(:,:) ! Stores the geometric weight at the grid points
Real(kind=db), save, allocatable:: gtilde(:,:) ! Stores the extension to the domain of the boundary conditions
type(mumps_mat):: etilde ! Matrix of extendend web splines
type(mumps_mat):: etildet ! Transpose of Matrix of extendend web splines
integer,save :: nbreducedspline ! Number of splines in the reduced set
PROCEDURE(geom_eval), POINTER:: total_weight => NULL()
PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL()
PUBLIC:: geom_weight
ABSTRACT INTERFACE
SUBROUTINE gtilde_eval(z,r,g)
USE constants
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: g(:,0:)
END SUBROUTINE
SUBROUTINE geom_eval(z,r,g)
USE constants
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: g(:,0:)
END SUBROUTINE
END INTERFACE
!
! INTERFACE
! SUBROUTINE gtest(z,r,w)
! USE constants
! Real(kind=db), INTENT(IN):: r(:),z(:)
! Real(kind=db), INTENT(OUT):: w(:,0:)
! END SUBROUTINE
! END INTERFACE
INTERFACE geom_weight
MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2
END INTERFACE geom_weight
NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb, above1, above2
contains
SUBROUTINE read_geom(Fileid, rnorm, rgrid, Potinn, Potout)
Use mpi
Real(kind=db):: rnorm, rgrid(:), Potinn, Potout
Integer:: Fileid, mpirank, ierr
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
Rewind(Fileid)
READ(Fileid, geomparams)
if(mpirank .eq. 0) WRITE(*, geomparams)
!! Normalizations and initialization of geometric variables
r_a=r_a/rnorm
r_b=r_b/rnorm
if(r_a .eq. 0 .and. r_b .eq.0) then
!! in case no geom_params have been defined
r_a=rgrid(1)
r_b=rgrid(size(rgrid,1))
end if
z_0=z_0/rnorm
r_0=r_0/rnorm
z_r=z_r/rnorm
r_r=r_r/rnorm
invr_r=1/r_r
invr_z=1/z_r
Interior=-1
above1=1
above2=-1
Phidown=Potinn
Phiup=Potout
L_r=L_r/rnorm
L_z=L_z/rnorm
SELECT CASE (walltype)
CASE(-1)
! test case
total_gtilde=>gtest
total_weight=>geom_weightstd
CASE (2)
! coaxial cylinder and top ellipse with cylinder extensions
total_gtilde=>g2
total_weight=>geom_w2
CASE (3)
! Two ellipses with cylinders
total_gtilde=>g3
total_weight=>geom_w3
CASE DEFAULT
! Ellipse as in gt170 standard weight
total_weight=>geom_weightstd
total_gtilde=>gstd
END SELECT
end subroutine read_geom
SUBROUTINE init_geom(spl2, vec1, vec2)
type(spline2d):: spl2
real(kind=db):: vec1(:),vec2(:)
Real(kind=db), Allocatable:: zgrid(:),rgrid(:)
Integer:: nrank(2), nrz(2)
Call get_dim(spl2, nrank, nrz)
! Obtain grid data drom the spline structure
Allocate(zgrid(1:nrz(1)+1))
Allocate(rgrid(1:nrz(2)+1))
zgrid=spl2%sp1%knots(0:nrz(1))
rgrid=spl2%sp2%knots(0:nrz(2))
! create a table of the calssification for each cell
Allocate(gridcelltype(nrz(1),nrz(2)))
Call classifycell(zgrid, rgrid, gridcelltype)
if (nlweb) then
Allocate(bsplinetype(nrank(1)*nrank(2)))
Call classifyspline(spl2, gridcelltype, bsplinetype)
Call buildetilde(spl2,bsplinetype,gridcelltype,etilde)
end if
ALLOCATE(gridwgeom((nrz(1)+1)*(nrz(2)+1),0:2))
gridwgeom=0
ALLOCATE(gtilde((nrz(1)+1)*(nrz(2)+1),0:2))
gtilde=0
CALL total_gtilde(vec1, vec2, gtilde)
Call geom_weight (vec1, vec2, gridwgeom)
end Subroutine init_geom
Subroutine geom_diag(File_handle, str, rnorm)
Use mpi
Use futils
Integer:: File_handle
Real(kind=db):: rnorm
Character(len=*):: str
CHARACTER(len=128):: grpname
Integer:: ierr, mpirank
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
IF(mpirank .eq. 0) THEN
Write(grpname,'(a,a)') trim(str),"/geometry"
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "r_a", r_a*RNORM)
Call attach(File_handle, trim(grpname), "r_b", r_b*RNORM)
Call attach(File_handle, trim(grpname), "z_0", z_0*RNORM)
Call attach(File_handle, trim(grpname), "r_0", r_0*RNORM)
Call attach(File_handle, trim(grpname), "r_r", r_r*RNORM)
Call attach(File_handle, trim(grpname), "z_r", z_r*RNORM)
Call attach(File_handle, trim(grpname), "interior", interior)
Call attach(File_handle, trim(grpname), "above1", above1)
Call attach(File_handle, trim(grpname), "above2", above2)
Call attach(File_handle, trim(grpname), "walltype", walltype)
Call putarr(File_handle, trim(grpname)//'/geomweight',gridwgeom)
END IF
End subroutine geom_diag
Subroutine buildetilde(spl2, bsplinetype, celltype, etilde)
USE mumps_bsplines
type(spline2d):: spl2
Integer:: bsplinetype(:)
Integer:: celltype(:,:)
type(mumps_mat):: etilde
Logical, Allocatable:: Ibsplinetype(:)
Integer:: i,j,k, icellz, icellr, jcellz, jcellr, nrank(2), nrz(2), norder(2)
real(kind=db), allocatable:: zgrid(:), rgrid(:)
real(kind=db):: wgeomi, indexdistance, eij
integer:: l,m, n, linkedi
type(innerspline):: innersplinelist
real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:)
!if(allocated(etilde)) deallocate(etilde)
nbreducedspline=count(bsplinetype .eq. 1)
Call get_dim(spl2,nrank,nrz,norder)
! Obtain grid data drom the spline structure
Allocate(zgrid(1:nrz(1)+1))
Allocate(rgrid(1:nrz(2)+1))
zgrid=spl2%sp1%knots(0:nrz(1))
rgrid=spl2%sp2%knots(0:nrz(2))
allocate(innersplinelist%k(nrank(1)*nrank(2)),innersplinelist%weight(nrank(1)*nrank(2)))
allocate(Ibsplinetype(nrank(1)*nrank(2)))
allocate(rgridmesh(nrank(1)*nrank(2)), zgridmesh(nrank(1)*nrank(2)))
allocate(hmesh(nrank(1)*nrank(2)))
allocate(d(size(bsplinetype)))
Ibsplinetype=.False.
call calcsplinecenters(spl2%sp1,cz,hz)
call calcsplinecenters(spl2%sp2,cr,hr)
Do i=0,nrank(2)-1
zgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cz
rgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cr(i+1)
hmesh(i*nrank(1)+1:(i+1)*nrank(1))=hr(i+1)*hz
End do
! allocate etilde and its transpose
call init(nrank(1)*nrank(2),nbreducedspline,etilde)
call init(nrank(1)*nrank(2),nbreducedspline,etildet)
k=1
Do i=1,nrank(1)*nrank(2)
if(bsplinetype(i) .ne. 1) cycle ! span of this spline is completely outside D
! one cell of bspline i is completely in D
icellz=mod(i-1,nrank(1))+1
icellr=(i-1)/(nrank(1))+1
outer: do l=max(1,icellz-norder(1)),min(nrz(1),icellz)
do m=max(1,icellr-norder(2)),min(nrz(2),icellr)
if (celltype(l,m).eq.1) then
call geom_weight((zgrid(l)+zgrid(l+1))/2,(rgrid(m)+rgrid(m+1))/2,wgeomi)
EXIT outer
end if
end do
end do outer
call putele(etilde,k,i,1/wgeomi)
call putele(etildet,i,k,1/wgeomi)
innersplinelist%k(i)=k
innersplinelist%weight(i)=1/wgeomi
k=k+1
!Ibsplinetype(i)=icellz-norder(1) .gt. 0 .and. icellr-norder(2) .gt. 0 .and. icellz .le. nrz(1) .and. icellr .le. nrz(2)
Ibsplinetype(i)=icellz+norder(1) .le. nrank(1) .and. icellr+norder(2) .le. nrank(2)
if(.not. Ibsplinetype(i)) cycle
do m=0,norder(2)
do l=0,norder(1)
! Check if all linked splines are inner splines
Ibsplinetype(i)=Ibsplinetype(i) .and. (bsplinetype(i+l+m*nrank(1)) .eq. 1)
end do
end do
end do
!!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(d,k,icellz,icellr,jcellr,jcellz,indexdistance,n,i,m,l,eij,linkedi)
Do j=1,nrank(1)*nrank(2)
! find the closest b-spline fully in D for the web method
if(bsplinetype(j) .ne. 0) cycle
d=0
where(Ibsplinetype)
d=(rgridmesh(j)-rgridmesh)**2+(zgridmesh(j)-zgridmesh)**2
end where
k=minloc(d,1,MASK=Ibsplinetype)
icellz=mod(k-1,nrank(1))
icellr=(k-1)/(nrank(1))
jcellz=mod(j-1,nrank(1))
jcellr=(j-1)/(nrank(1))
indexdistance=sqrt(real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db))
if(indexdistance .gt. 10)then
Write(*,*) 'Error on system conditioning, the number of radial or axial points is too low!'
Write(*,*) 'Distance found: ', indexdistance
!stop
end if
do n=0,norder(2)
do i=0,norder(1)
eij=1
!eij=1
do m=0,norder(2)
if(n.eq.m) cycle
eij=eij*(rgridmesh(j)-rgridmesh(k+m*nrank(1)))/(rgridmesh(k+n*nrank(1))-rgridmesh(k+m*nrank(1)))
end do
do l=0,norder(1)
if( i .eq. l ) cycle
eij=eij*(zgridmesh(j)-zgridmesh(k+l))/(zgridmesh(k+i)-zgridmesh(k+l))
!eij=eij*(jcellr-icellr-m)/(n-m)*(jcellz-icellz-l)/(i-l)
end do
linkedi=innersplinelist%k(k+i+n*nrank(1)) ! equivalent to findloc, necessary for ifort 17
!findloc(innersplinelist%mu,k,1)
eij=eij*innersplinelist%weight(k+i+n*nrank(1))
call putele(etilde,linkedi,j,eij)
call putele(etildet,j,linkedi,eij)
end do
end do
end do
!!$OMP End parallel do
call to_mat(etilde)
call to_mat(etildet)
end subroutine
Subroutine calcsplinecenters(spl,ctrs,heights)
use bsplines
type(spline1d):: spl
real(kind=db), allocatable:: ctrs(:), heights(:)
integer:: nrank, nx, order, i, left1, j, left2, left3
real(kind=db):: x1, x2, x3
real(kind=db), allocatable:: fun1(:,:), fun2(:,:), fun3(:,:)
real(kind=db),allocatable:: init_heights(:)
call get_dim(spl, nrank, nx, order)
if (allocated(ctrs)) deallocate(ctrs)
if (allocated(heights)) deallocate(heights)
allocate(heights(nrank))
allocate(ctrs(nrank))
allocate(fun1(1:order+1,0:1), fun2(1:order+1,0:1), fun3(1:order+1,0:1))
allocate(init_heights(nrank))
init_heights=1.0
ctrs(:)=spl%knots(0:nrank-1)
call gridval(spl, ctrs, heights, 0, init_heights)
if (order .gt.1) then
do i=2,nrank-1
x1=spl%knots(i-order)+1e5*Epsilon(x1)
x2=spl%knots(i-1)-1e5*Epsilon(x2)
left1=max(0,i-order)
left2=min(nx-2,i-2)
call basfun(x1,spl,fun1, left1+1)
!Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold
call basfun(x2,spl,fun2, left2+1)
fun3=1
j=0
Do while( j .le.100)
!Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold
x3=(x1+x2)/2
call locintv(spl, x3, left3)
call basfun(x3,spl,fun3, left3+1)
if( (x2-x1)/2.lt.1e-12) exit
if(abs(fun3(i-left3,1)).lt.1e-15) exit
if(fun3(i-left3,1)*fun1(i-left1,1).le.0) then
fun2=fun3
x2=x3
left2=left3
else
fun1=fun3
x1=x3
left1=left3
end if
j=j+1
End do
ctrs(i)=x3
heights(i)=fun3(i-left3,0)
end do
end if
end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Subroutine classifycell(zgrid,rgrid,gridctype)
Real(kind=db):: zgrid(:),rgrid(:)
Integer:: gridctype(:,:)
Integer::i,j
gridctype=-1
!$OMP parallel do private(i)
Do j=1,size(rgrid,1)-1
Do i=1,size(zgrid,1)-1
! Determines the type inner/boundary/outer for each cell
Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),&
& gridctype(i,j))
End Do
End do
!$OMP end parallel do
End subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine classifyspline(spl2,gridctype,bsptype)
type(spline2d):: spl2
Integer:: gridctype(:,:)
Integer:: bsptype(:)
Integer:: nrank(2), nrz(2), norder(2)
Integer:: i, j, mu, imin, imax, jmin, jmax
Integer, Allocatable:: splinespan(:,:)
call get_dim(spl2, nrank, nrz, norder)
! by default, all splines have part of their support outside the domain D
bsptype=0
Allocate(splinespan(norder(1)+1,norder(2)+1))
Do mu=1,nrank(1)*nrank(2)
! scan the spline space
! obtain the axial spline index
i=mod(mu-1,nrank(1))+1
! obtain the radial spline index
j=(mu-1)/nrank(1)+1
! by default all cells are outside
splinespan=-1
! find the axial span of this spline in cell indices
imin=max(1,i-norder(1))
imax=min(nrz(1),i)
! find the radial span of this spline in cell indices
jmin=max(1,j-norder(2))
jmax=min(nrz(2),j)
! obtain the cell type on which the spline is defined
splinespan(1:(imax-imin+1),1:(jmax-jmin+1))=gridctype(imin:imax,jmin:jmax)
if ( ANY( splinespan==1 ) ) then
bsptype(mu)=1
else if (sum(splinespan).eq. -(norder(1)+1)*(1+norder(2))) then
bsptype(mu)=-1
end if
end do
End subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> general placeholder function used in fields_mod to compute the weight for weighted_bsplines
!> This functions combine an ellipse function and an horizontal wall
!
!---------------------------------------------------------------------------
SUBROUTINE geom_weight0(z,r,w)
Real(kind=db), INTENT(IN):: r,z
Real(kind=db), INTENT(OUT):: w
Real(kind=db):: rtmp(1:1), ztmp(1:1), wtmp(1:1,1:1)
ztmp=z
rtmp=r
call total_weight(ztmp,rtmp,wtmp)
w=wtmp(1,1)
End SUBROUTINE geom_weight0
SUBROUTINE geom_weight1(z,r,w)
Real(kind=db), INTENT(IN):: r,z
Real(kind=db), INTENT(OUT):: w(0:)
Real(kind=db):: rtmp(1:1), ztmp(1:1)
Real(kind=db):: wtmp(1:1,0:size(w,1)-1)
ztmp=z
rtmp=r
call total_weight(ztmp,rtmp,wtmp)
w=wtmp(1,:)
End SUBROUTINE geom_weight1
SUBROUTINE geom_weight2(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
call total_weight(z,r,w)
!
End SUBROUTINE geom_weight2
SUBROUTINE geom_weightstd(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1)
Real(kind=db):: squareroot(size(w,1)), denom(size(w,1))
call cyllweight(z,r,walltmp, r_a, above1)
SELECT CASE (walltype)
CASE (0)
call cyllweight(z,r,elliptmp, r_b, above2)
CASE DEFAULT
call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior)
END SELECT
if(interior.eq.0 .and. above2.eq.0) then
w=walltmp(:,:)
else if (above1 .eq. 0) then
w=elliptmp(:,:)
else
denom=walltmp(:,0)**2+elliptmp(:,0)**2
squareroot=sqrt(denom)
w(:,0)=walltmp(:,0)+elliptmp(:,0)-squareroot ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=walltmp(:,1)+elliptmp(:,1)-(elliptmp(:,1)*elliptmp(:,0)+walltmp(:,1)*walltmp(:,0))/&
& squareroot ! z derivative of w
w(:,2)=walltmp(:,2)+elliptmp(:,2)-(elliptmp(:,2)*elliptmp(:,0)+walltmp(:,2)*walltmp(:,0))/&
& squareroot ! r derivative of w
End If
End if
End SUBROUTINE geom_weightstd
SUBROUTINE geom_w2(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1)
Real(kind=db):: squareroot(size(w,1)), denom(size(w,1))
call cyllweight(z,r,walltmp, r_a, above1)
call ellipsewall(z,r,elliptmp,r_b,above2,r_0, z_0, invr_z, invr_r, Interior)
denom=walltmp(:,0)**2+elliptmp(:,0)**2
squareroot=sqrt(denom)
w(:,0)=walltmp(:,0)+elliptmp(:,0)-squareroot ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=walltmp(:,1)+elliptmp(:,1)-(elliptmp(:,1)*elliptmp(:,0)+walltmp(:,1)*walltmp(:,0))/&
& squareroot ! z derivative of w
w(:,2)=walltmp(:,2)+elliptmp(:,2)-(elliptmp(:,2)*elliptmp(:,0)+walltmp(:,2)*walltmp(:,0))/&
& squareroot ! r derivative of w
End If
End SUBROUTINE geom_w2
SUBROUTINE geom_w3(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db):: belowtmp(size(w,1),0:size(w,2)-1), abovetmp(size(w,1),0:size(w,2)-1)
Real(kind=db):: squareroot(size(w,1)), denom(size(w,1))
! Weight functions necessary for calculation of g
! Ellipse and bottom cylinder combination
call ellipsewall(z,r,belowtmp,r_a,1,r_b, z_0, 1/(z_r+r_b-r_a), 1/(z_r+r_b-r_a), 1)
! Ellipse and top cylinder combination
call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1)
denom=abovetmp(:,0)**2+belowtmp(:,0)**2
squareroot=sqrt(denom)
w(:,0)=belowtmp(:,0)+abovetmp(:,0)-squareroot ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=belowtmp(:,1)+abovetmp(:,1)-(abovetmp(:,1)*abovetmp(:,0)+belowtmp(:,1)*belowtmp(:,0))/&
& squareroot ! z derivative of w
w(:,2)=belowtmp(:,2)+abovetmp(:,2)-(abovetmp(:,2)*abovetmp(:,0)+belowtmp(:,2)*belowtmp(:,0))/&
& squareroot ! r derivative of w
End If
End SUBROUTINE geom_w3
Subroutine Boundary_limits(pts)
Real(kind=db), Allocatable, INTENT(out) :: pts(:,:)
If(allocated(pts))Deallocate(pts)
allocate(pts(4,2))
pts(1,1:2)=(/z_0-invr_z,r_0/)
pts(2,1:2)=(/z_0,r_0+invr_r/)
pts(3,1:2)=(/z_0+invr_z,r_0/)
pts(4,1:2)=(/z_0,r_0-invr_r/)
End subroutine Boundary_limits
SUBROUTINE classification(z, r, ctype)
real(kind=db), INTENT(IN):: r(2), z(2)
INTEGER, INTENT(OUT):: ctype
Real(kind=db)::weights(5,1:1)
Real(kind=db):: zeval(5),reval(5)
zeval=(/ z(1),z(2),z(1),z(2), (z(2)+z(1))/2 /)
reval=(/ r(1),r(1),r(2),r(2), (r(2)+r(1))/2 /)
CAll geom_weight(zeval,reval,weights)
ctype=int(sign(1.0_db,weights(5,1)))
If(weights(1,1)*weights(2,1) .lt. 0 ) then
ctype=0
return
End If
If(weights(1,1)*weights(3,1) .lt. 0 ) then
ctype=0
return
End If
If(weights(2,1)*weights(4,1) .lt. 0 ) then
ctype=0
return
End If
If(weights(3,1)*weights(4,1) .lt. 0 ) then
ctype=0
return
End If
end subroutine
! ##################################################################
Subroutine calc_gauss(spl2, ngauss, i, j, xgauss, wgauss, gausssize, celltype)
type(spline2d), INTENT(IN):: spl2
Integer, Intent(out):: gausssize
INTEGER, Intent(out), Optional :: celltype
Real(kind=db), Allocatable::rgrid(:),zgrid(:)
Real(kind=db), ALLOCATABLE::xgauss(:,:), wgauss(:)
Integer:: i,j, ngauss(2)
Real(kind=db),Allocatable:: pts(:,:), zpoints(:)
Real(kind=db):: zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))
Integer:: k, l, directiondown, directionup, nbzpoints, direction, ctype
Logical:: hasmaxpoint
Real(kind=db):: xptup, xptdown, w, rmin, rmax
type(spline1d):: splz, splr
splz=spl2%sp1
splr=spl2%sp2
Allocate(zgrid(1:splz%nints+1))
Allocate(rgrid(1:splr%nints+1))
zgrid=splz%knots(0:splz%nints)
rgrid=splr%knots(0:splr%nints)
hasmaxpoint=.false.
If(allocated(xgauss)) deallocate(xgauss)
if(allocated(wgauss)) deallocate(wgauss)
!Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),ctype)
ctype=gridcelltype(i,j)
If (ctype .ge. 1) then ! we have a normal internal cell
Allocate(xgauss(ngauss(1)*ngauss(2),2))
Allocate(wgauss(ngauss(1)*ngauss(2)))
gausssize=ngauss(1)*ngauss(2)
!Computation of gauss weight and position in r and z direction for gaussian integration
CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg)
CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg)
! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration
DO k=1,ngauss(2)
xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg
xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k)
wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg
END DO
Else If(ctype.eq.0) then ! we have a boundary cell
!Call Boundary_Limits(pts)
!Do k=1,size(pts,1)
! hasmaxpoint=hasmaxpoint .or. (pts(k,1) .ge. zgrid(i) .and. pts(k,1) .le. zgrid(i+1) &
! & .and. pts(k,2) .ge. rgrid(j) .and. pts(k,2) .le. rgrid(j+1))
!End do
If(.not. hasmaxpoint) Then
directiondown=1
directionup=1
Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown)
Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup)
select case ( directionup+directiondown)
Case (0) ! The intersections are only on the left and right cell boundaries
nbzpoints=2
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(i),zgrid(i+1)/)
Case(1)
nbzpoints=3
Allocate(zpoints(nbzpoints))
if(directiondown.eq.1) zpoints=(/zgrid(i), xptdown,zgrid(i+1)/)
if(directionup .eq. 1) zpoints=(/zgrid(i), xptup,zgrid(i+1)/)
Case(2)
nbzpoints=3
Allocate(zpoints(nbzpoints))
call geom_weight(zgrid(i),rgrid(j),w)
If(w.ge.0) zpoints=(/zgrid(1),min(xptdown,xptup),max(xptdown,xptup) /)
If(w.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /)
End select
Allocate(xgauss(ngauss(1)*ngauss(2)*(nbzpoints-1),2))
Allocate(wgauss(ngauss(1)*ngauss(2)*(nbzpoints-1)))
gausssize=ngauss(1)*ngauss(2)*(nbzpoints-1)
! Compute gauss points
Do l=1,nbzpoints-1
!CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg)
Call gauleg(zpoints(l),zpoints(l+1),zg,wzg,ngauss(1))
! We test if the lower or upper side is in the domain
call geom_weight(zg(1),rgrid(j),w)
rmin=rgrid(j)
if (w .lt. 0) rmin = rgrid(j+1)
Do k=1,ngauss(1)
direction=2
Call Find_crosspoint((/rgrid(j),rgrid(j+1)/),zg(k),rmax,direction) ! We compute the radial limits at each z position
Call gauleg(min(rmin,rmax),max(rmin,rmax),rg,wrg,ngauss(2)) ! We obtain the gauss w and pos for these boundaries
xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),1) = zg(k)
xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),2) = rg
wgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1)) = wrg*wzg(k)
End do
End Do
End If
Else
gausssize=1
Allocate(xgauss(1:1,2))
Allocate(wgauss(1:1))
!Computation of gauss weight and position in r and z direction for gaussian integration
CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg)
CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg)
! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration
xgauss(1,1)=(zg(1)+zg(ngauss(1)))*0.5
xgauss(1,2)=(rg(1)+rg(ngauss(2)))*0.5
wgauss(1)=0
End If
If(PRESENT(celltype)) celltype=ctype
End Subroutine calc_gauss
Subroutine Find_crosspoint(x,y,xpt, direction)
Real(kind=db):: x(2), y
Real(kind=db):: xptm, xpt, temp
Real(kind=db):: w, wold
Integer, Intent(INOUT):: direction
Integer:: i
xptm=x(1)
xpt=x(2)
i=0
if(direction .eq.1) Call geom_weight(xptm,y,wold)
if(direction .eq.2) Call geom_weight(y,xptm,wold)
if(direction .eq.1) Call geom_weight(xpt,y,w)
if(direction .eq.2) Call geom_weight(y,xpt,w)
if(w*wold.gt.0) then
direction=0
return
End If
Do while( i .le.100 .and. w*wold .lt. 0)
!Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold
if(direction .eq.1) Call geom_weight(xpt,y,w)
if(direction .eq.2) Call geom_weight(y,xpt,w)
xptm=xpt-w*(xpt-xptm)/(w-wold)
wold=w
temp=xptm
xptm=xpt
xpt=temp
i=i+1
End do
if(xpt .ge. x(2) .or. xpt .le. x(1)) direction=0
End Subroutine
SUBROUTINE cyllweight(z,r,w, r_lim, above)
Real(kind=db):: z(:),r(:),w(:,0:), r_lim
Integer:: above
w(:,0)=above*(r-r_lim) ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=0 ! z derivative of w
w(:,2)=above ! r derivative of w
End If
End subroutine cyllweight
SUBROUTINE ellipsewall(z,r,w,r_lim,above,r_0, z_0, invr_z, invr_r, Interior)
Real(kind=db):: z(:),r(:),w(:,0:)
Integer:: Interior, above
Real(kind=db):: r_lim, r_0,z_0,invr_z,invr_r
Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1)
Real(kind=db):: squareroot(size(w,1)), denom(size(w,1))
call cyllweight(z,r,walltmp, r_lim, above)
call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior)
denom=walltmp(:,0)**2+elliptmp(:,0)**2
squareroot=sqrt(denom)
w(:,0)=walltmp(:,0)+elliptmp(:,0)+Interior*squareroot ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=walltmp(:,1)+elliptmp(:,1)+Interior*(elliptmp(:,1)*elliptmp(:,0)+walltmp(:,1)*walltmp(:,0))/&
& squareroot ! z derivative of w
w(:,2)=walltmp(:,2)+elliptmp(:,2)+Interior*(elliptmp(:,2)*elliptmp(:,0)+walltmp(:,2)*walltmp(:,0))/&
& squareroot ! r derivative of w
End If
END SUBROUTINE ellipsewall
Subroutine ellipseweight(z,r,w, r_0, z_0, invr_z, invr_r, Interior)
Real(kind=db):: z(:),r(:),w(:,0:)
Real(kind=db):: r_0,z_0,invr_z,invr_r
Integer:: Interior
w(:,0)=Interior*(1-(r-r_0)**2*invr_r**2-(z-z_0)**2*invr_z**2) ! weight at position r,z
If(size(w,2) .gt. 1) then ! first derivative
w(:,1)=-2*(z-z_0)*invr_z**2*Interior ! z derivative of w
w(:,2)=-2*(r-r_0)*invr_r**2*Interior ! r derivative of w
End If
End subroutine ellipseweight
Subroutine gstd(z,r,gtilde)
! standard g function added to rhs of poisson solver
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1)
Real(kind=db):: denom(size(r,1))
! Weight functions necessary for calculation of g
! coaxial insert
call cyllweight(z,r,belowtmp, r_a, above1)
SELECT CASE (walltype)
CASE (0)
! top cylinder
call cyllweight(z,r,abovetmp, r_b, above2)
CASE DEFAULT
! Ellipse as in gt170
call ellipseweight(z,r,abovetmp,r_0,z_0,invr_z,invr_r,Interior)
END SELECT
! Extension to the whole domain of the boundary conditions
! constructed by weight multiplied with boundary value
gtilde(:,0)=(Phidown*abovetmp(:,0) + Phiup*belowtmp(:,0) ) / &
& (abovetmp(:,0)+belowtmp(:,0))
If(size(gtilde,2) .gt. 2) then ! first derivative
denom=(abovetmp(:,0)+belowtmp(:,0))**2
gtilde(:,1)=(Phiup-Phidown)*(belowtmp(:,1)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,1)) / &
& denom ! Axial derivative
gtilde(:,2)=(Phiup-Phidown)*(belowtmp(:,2)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,2)) / &
& denom ! Radial derivative
End If
End subroutine
Subroutine g2(z,r,gtilde)
! standard g function added to rhs of poisson solver
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1)
Real(kind=db):: denom(size(r,1))
! Weight functions necessary for calculation of g
! coaxial insert
call cyllweight(z,r,belowtmp, r_a, above1)
! Ellipse and top cylinder combination
call ellipsewall(z,r,abovetmp,r_b,above2,r_0,z_0,invr_z,invr_r,Interior)
! Extension to the whole domain of the boundary conditions
! constructed by weight multiplied with boundary value
gtilde(:,0)=(Phidown*abovetmp(:,0) + Phiup*belowtmp(:,0) ) / &
& (abovetmp(:,0)+belowtmp(:,0))
If(size(gtilde,2) .gt. 2) then ! first derivative
denom=(abovetmp(:,0)+belowtmp(:,0))**2
gtilde(:,1)=(Phiup-Phidown)*(belowtmp(:,1)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,1)) / &
& denom ! Axial derivative
gtilde(:,2)=(Phiup-Phidown)*(belowtmp(:,2)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,2)) / &
& denom ! Radial derivative
End If
End subroutine
SUBROUTINE g3(z,r,gtilde)
! standard g function added to rhs of poisson solver
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1)
Real(kind=db):: denom(size(r,1))
! Weight functions necessary for calculation of g
! Ellipse and bottom cylinder combination
call ellipsewall(z,r,belowtmp,r_a,1,r_b, z_0, 1/(z_r+r_b-r_a), 1/(z_r+r_b-r_a), 1)
! Ellipse and top cylinder combination
call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1)
! Extension to the whole domain of the boundary conditions
! constructed by weight multiplied with boundary value
gtilde(:,0)=(Phidown*abovetmp(:,0) + Phiup*belowtmp(:,0) ) / &
& (abovetmp(:,0)+belowtmp(:,0))
If(size(gtilde,2) .gt. 2) then ! first derivative
denom=(abovetmp(:,0)+belowtmp(:,0))**2
gtilde(:,1)=(Phiup-Phidown)*(belowtmp(:,1)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,1)) / &
& denom ! Axial derivative
gtilde(:,2)=(Phiup-Phidown)*(belowtmp(:,2)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,2)) / &
& denom ! Radial derivative
End If
END SUBROUTINE g3
Subroutine gtest(z,r,gtilde)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
!
gtilde(:,0)=sin(2*pi*z/L_z)*cos(2*pi*r/L_r) + 2
If(size(gtilde,2) .gt. 1) then ! first derivative
gtilde(:,1)=2*pi/L_z*cos(2*pi*z/L_z)*cos(2*pi*r/L_r)
gtilde(:,2)=-2*pi/L_r*sin(2*pi*z/L_z)*sin(2*pi*r/L_r)
End If
End subroutine
Subroutine ftest(z,r,f)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT)::f(:,0:)
!
f(:,0)=2*pi*sin(2*pi*z/L_z)*( &
& 2*pi*cos(2*pi*r/L_r)*(1/L_z**2+1/L_r**2) &
& +1/(r*L_r)*sin(2*pi*r/L_r) )
End Subroutine ftest
subroutine Reducematrix(full,reduced)
use mumps_bsplines
type(mumps_mat):: full, reduced, tempmat
Real(kind=db), allocatable:: temparray(:), rsltarray(:)
Integer:: fullrank, reducedrank
Integer:: i,j
call to_mat(full)
fullrank=full%rank
reducedrank=nbreducedspline
call init(reducedrank,2,reduced)
call init(fullrank,2,tempmat)
!$OMP parallel private(j), firstprivate(temparray,rsltarray), default(shared)
allocate(temparray(fullrank),rsltarray(fullrank))
!$OMP do
do i=1,fullrank
call getcol(etildet,i,temparray)
rsltarray=vmx_mumps_mat_loc(full,temparray)
DO j=1,fullrank
CALL putele_sploc(tempmat%mat%row(i), j, rsltarray(j), .false.)
END DO
end do
!$OMP end do
!$OMP barrier
!$OMP single
call to_mat(tempmat)
!$OMP end single
!$OMP do
do i=1,reducedrank
call getcol(etildet,i,temparray)
rsltarray=vmx_mumps_mat_loc(tempmat,temparray)
DO j=1,reducedrank
CALL putele_sploc(reduced%mat%row(i), j, rsltarray(j), .false.)
END DO
end do
!$OMP end do
deallocate(temparray,rsltarray)
!$OMP end parallel
end subroutine
!!###############################################################################################
! subroutines imported from mumps_bsplines and sparse lib to have thread-private routines
SUBROUTINE putele_sploc(arow, j, val, nlforce_zero)
!
! Put (overwrite) element j of row arow or insert it in an increasing "index"
!
USE sparse
TYPE(sprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: val
LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero
!
TYPE(elt), TARGET :: pre_root
TYPE(elt), POINTER :: t, p
LOGICAL :: rmnode
!
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
!
! Remove node which has zero val or not?
! But never create new node with zero val
!
rmnode = .TRUE.
IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero
!
DO WHILE( ASSOCIATED(t%next) )
p => t%next
IF( p%index .EQ. j ) THEN
IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val!
t%next => p%next
arow%nnz = arow%nnz-1
arow%row0 => pre_root%next ! In case the head is altered
DEALLOCATE(p)
ELSE
p%val = val
END IF
RETURN
END IF
IF( p%index .GT. j ) EXIT
t => t%next
END DO
!
! Never create new node with zero val
!
IF(ABS(val).GT.EPSILON(0.0d0)) THEN
ALLOCATE(p)
p = elt(j, val, t%next)
t%next => p
arow%nnz = arow%nnz+1
arow%row0 => pre_root%next
END IF
END SUBROUTINE putele_sploc
!
SUBROUTINE getcol_loc(amat, j, arr)
!
! Get a column from sparse matrix
!
USE mumps_bsplines
TYPE(mumps_mat), INTENT(in) :: amat
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(out) :: arr(:)
INTEGER :: i
!
DO i=amat%istart,amat%iend
CALL getele_loc(amat, i, j, arr(i))
END DO
END SUBROUTINE getcol_loc
SUBROUTINE getele_loc(mat, i, j, val)
!
! Get element (i,j) of sparse matrix
!
USE mumps_bsplines
TYPE(mumps_mat) :: mat
INTEGER, INTENT(in) :: i, j
DOUBLE PRECISION, INTENT(out) :: val
INTEGER :: iget, jget
INTEGER :: k, s, e
!
iget = i
jget = j
IF(mat%nlsym) THEN
IF( i.GT.j ) THEN ! Lower triangular part
iget = j
jget = i
END IF
END IF
!
val = 0.0d0 ! Assume zero val if outside
IF( iget.LT.mat%istart .OR. iget.GT.mat%iend ) RETURN
!
IF(ASSOCIATED(mat%mat)) THEN
CALL getele(mat%mat, iget, jget, val)
ELSE
s = mat%irow(iget) - mat%nnz_start + 1
e = mat%irow(iget+1) - mat%nnz_start
k = isearch(mat%cols(s:e), jget)
IF( k.GE.0 ) THEN
val =mat%val(s+k)
ELSE
val = 0.0d0 ! Assume zero val if not found
END IF
END IF
END SUBROUTINE getele_loc
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION vmx_mumps_mat_loc(mat, xarr) RESULT(yarr)
!
! Return product mat*x
!
Use mumps_bsplines
TYPE(mumps_mat) :: mat
DOUBLE PRECISION, INTENT(in) :: xarr(:)
DOUBLE PRECISION :: yarr(SIZE(xarr))
DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0
CHARACTER(len=6) :: matdescra
INTEGER :: n
!
n = mat%rank
!
!#ifdef MKL
IF(mat%nlsym) THEN
matdescra = 'sun'
ELSE
matdescra = 'g'
END IF
!
CALL mkl_dcsrmv('N', n, n, alpha, matdescra, mat%val, &
& mat%cols, mat%irow(1), mat%irow(2), xarr, &
& beta, yarr)
!#else
! yarr = 0.0d0
! DO i=1,n
! DO j=mat%irow(i), mat%irow(i+1)-1
! yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j))
! END DO
! IF(mat%nlsym) THEN
! DO j=mat%irow(i)+1, mat%irow(i+1)-1
! yarr(mat%cols(j)) = yarr(mat%cols(j)) &
! & + mat%val(j)*xarr(i)
! END DO
! END IF
! END DO
! Write(*,*) "linear mult")
!#endif
!
END FUNCTION vmx_mumps_mat_loc
END MODULE geometry
c4science · Help