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

fennecs/src/geometry_mod.f90
a2847b9464fageom_testpardiso

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 pardiso_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(pardiso_mat):: etilde ! Matrix of extendend web splines
type(pardiso_mat):: etildet ! Transpose of Matrix of extendend web splines
integer,save :: nbreducedspline ! Number of splines in the reduced set
PUBLIC:: geom_weight
INTERFACE geom_weight
MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2
END INTERFACE geom_weight
NAMELIST /geom_params/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb
contains
SUBROUTINE init_geom(Fileid, spl2, Potinn, Potout, rnorm, vec1, vec2)
type(spline2d):: spl2
Real(kind=db):: Potinn, Potout, rnorm
real(kind=db):: vec1(:),vec2(:)
Real(kind=db), Allocatable:: zgrid(:),rgrid(:)
Integer:: Fileid, nrank(2), nrz(2)
Rewind(Fileid)
READ(Fileid, geom_params)
WRITE(*, geom_params)
!! Normalizations and initialization of geometric variables
r_a=r_a/rnorm
r_b=r_b/rnorm
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
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 CompgBoundary(vec1,vec2,gtilde)
Call geom_weight(vec1, vec2, gridwgeom)
end Subroutine init_geom
Subroutine diag_geom(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
Subroutine buildetilde(spl2, bsplinetype, celltype, etilde)
USE pardiso_bsplines
type(spline2d):: spl2
Integer:: bsplinetype(:)
Integer:: celltype(:,:)
type(pardiso_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, eij, indexdistance, eij1
integer:: l,m, n, linkedi
type(innerspline):: innersplinelist
real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:)
real(kind=db),allocatable :: etildemat(:,:)
!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(nbinspline,nrank(1)*nrank(2)))
!etilde=0
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 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)
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. 20)then
Write(*,*) 'Error on system conditioning, the number of radial or axial points is too low'
stop
end if
do n=0,norder(2)
do i=0,norder(1)
eij1=1
!eij=1
do m=0,norder(2)
if(n.eq.m) cycle
eij1=eij1*(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
eij1=eij1*(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)
eij1=eij1*innersplinelist%weight(k+i+n*nrank(1))
!eij=eij*innersplinelist%weight(k+i+n*nrank(1))
call putele(etilde,linkedi,j,eij1)
call putele(etildet,j,linkedi,eij1)
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,j)
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)
bsptype=0
Allocate(splinespan(norder(1)+1,norder(2)+1))
Do mu=1,nrank(1)*nrank(2)
i=mod(mu-1,nrank(1))+1
j=(mu-1)/nrank(1)+1
splinespan=-1
imin=max(1,i-norder(1))
imax=min(nrz(1),i)
jmin=max(1,j-norder(2))
jmax=min(nrz(2),j)
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), walltmp(1:1,1:1), elliptmp(1:1,1:1)
ztmp=z
rtmp=r
call wallweight(ztmp,rtmp,walltmp, r_a, above1)
if (walltype .ne. 0) then
call ellipseweight(ztmp,rtmp,elliptmp,Interior)
else
call wallweight(ztmp,rtmp,elliptmp, r_b, above2)
end if
if(interior.eq.0 .and. above2.eq.0) then
w=walltmp(1,1)
else if (above1 .eq. 0) then
w=elliptmp(1,1)
else
w=walltmp(1,1)+elliptmp(1,1)-sqrt(walltmp(1,1)**2+elliptmp(1,1)**2) ! weight at position r,z
end if
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), walltmp(1:1,0:size(w,1)-1), elliptmp(1:1,0:size(w,1)-1)
Real(kind=db):: squareroot, denom
ztmp=z
rtmp=r
call wallweight(ztmp,rtmp,walltmp, r_a, above1)
if (walltype .ne. 0) then
call ellipseweight(ztmp,rtmp,elliptmp,Interior)
else
call wallweight(ztmp,rtmp,elliptmp, r_b, above2)
end if
if(interior.eq.0 .and. above2.eq.0) then
w=walltmp(1,:)
else if (above1 .eq. 0) then
w=elliptmp(1,:)
else
denom=walltmp(1,0)**2+elliptmp(1,0)**2
squareroot=sqrt(denom)
w(0)=walltmp(1,0)+elliptmp(1,0)-squareroot ! weight at position r,z
If(size(w,1) .gt. 1) then ! first derivative
w(1)=walltmp(1,1)+elliptmp(1,1)-(elliptmp(1,1)*elliptmp(1,0)+walltmp(1,1)*walltmp(1,0))/&
& squareroot ! z derivative of w
w(2)=walltmp(1,2)+elliptmp(1,2)-(elliptmp(1,2)*elliptmp(1,0)+walltmp(1,2)*walltmp(1,0))/&
& squareroot ! r derivative of w
End If
If (size(w,1).gt.3) then ! second derivative
w(3)=walltmp(1,3)+elliptmp(1,3)-((elliptmp(1,1)*elliptmp(1,0)-walltmp(1,1)*walltmp(1,0))**2+&
& (elliptmp(1,3)*elliptmp(1,0)+walltmp(1,3)*walltmp(1,0))*denom)/&
& (squareroot*denom)! zz derivative of w
w(4)=walltmp(1,4)+elliptmp(1,4)-((elliptmp(1,4)*elliptmp(1,0)-walltmp(1,4)*walltmp(1,0))*denom+&
& elliptmp(1,1)*elliptmp(1,2)*walltmp(1,0)**2 + walltmp(1,1)*walltmp(1,2)*elliptmp(1,0)**2- &
& elliptmp(1,1)*elliptmp(1,0)*walltmp(1,2)*walltmp(1,0) - &
& elliptmp(1,2)*elliptmp(1,0)*walltmp(1,1)*walltmp(1,0) )/&
& (squareroot*denom)! zr derivative of w
w(5)=walltmp(1,5)+elliptmp(1,5)-((elliptmp(1,2)*elliptmp(1,0)-walltmp(1,2)*walltmp(1,0))**2+&
& (elliptmp(1,5)*elliptmp(1,0)+walltmp(1,5)*walltmp(1,0))*denom)/&
& (squareroot*denom)! rr derivative of w
End If
End if
End SUBROUTINE geom_weight1
SUBROUTINE geom_weight2(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 wallweight(z,r,walltmp, r_a, above1)
if (walltype .ne. 0) then
call ellipseweight(z,r,elliptmp,Interior)
else
call wallweight(z,r,elliptmp, r_b, above2)
end if
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
If (size(w,2).gt.3) then ! second derivative
w(:,3)=walltmp(:,3)+elliptmp(:,3)-((elliptmp(:,1)*elliptmp(:,0)-walltmp(:,1)*walltmp(:,0))**2+&
& (elliptmp(:,3)*elliptmp(:,0)+walltmp(:,3)*walltmp(:,0))*denom)/&
& (squareroot*denom)! zz derivative of w
w(:,4)=walltmp(:,4)+elliptmp(:,4)-((elliptmp(:,4)*elliptmp(:,0)-walltmp(:,4)*walltmp(:,0))*denom+&
& elliptmp(:,1)*elliptmp(:,2)*walltmp(:,0)**2 + walltmp(:,1)*walltmp(:,2)*elliptmp(:,0)**2- &
& elliptmp(:,1)*elliptmp(:,0)*walltmp(:,2)*walltmp(:,0) - &
& elliptmp(:,2)*elliptmp(:,0)*walltmp(:,1)*walltmp(:,0) )/&
& (squareroot*denom)! zr derivative of w
w(:,5)=walltmp(:,5)+elliptmp(:,5)-((elliptmp(:,2)*elliptmp(:,0)-walltmp(:,2)*walltmp(:,0))**2+&
& (elliptmp(:,5)*elliptmp(:,0)+walltmp(:,5)*walltmp(:,0))*denom)/&
& (squareroot*denom)! rr derivative of w
End If
End if
End SUBROUTINE geom_weight2
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=.false.
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
hasmaxpoint=.False.
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 wallweight(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
If (size(w,2).gt.3) then ! second derivative
w(:,3)= 0 ! zz derivative
w(:,4)= 0 ! zr derivative
w(:,5)= 0 ! rr derivative
End If
End subroutine wallweight
Subroutine ellipseweight(z,r,w, Interior)
Real(kind=db):: z(:),r(:),w(:,0:)
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
If (size(w,2).gt.3) then ! second derivative
w(:,3)=-2*invr_z**2*Interior ! zz derivative
w(:,4)=0 ! zr derivative
w(:,5)=-2*invr_r**2*Interior ! rr derivative
End If
End subroutine ellipseweight
Subroutine CompgBoundary(z,r,gtilde)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
if (walltype .eq. -1) then
call gtest(z,r,gtilde)
else
call gstd(z,r,gtilde)
end if
End subroutine
Subroutine gstd(z,r,gtilde)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
Real(kind=db):: walltmp(size(r,1),0:size(gtilde,2)-1), elliptmp(size(r,1),0:size(gtilde,2)-1)
Real(kind=db):: denom(size(r,1))
call wallweight(z,r,walltmp, r_a, above1)
if (walltype .eq. 1) then
call ellipseweight(z,r,elliptmp,Interior)
else
call wallweight(z,r,elliptmp, r_b, above2)
end if
! Extension to the whole domain of the boundary conditions
gtilde(:,0)=(Phidown*elliptmp(:,0) + Phiup*walltmp(:,0) ) / &
& (elliptmp(:,0)+walltmp(:,0))
If(size(gtilde,2) .gt. 1) then ! first derivative
denom=(elliptmp(:,0)+walltmp(:,0))**2
gtilde(:,1)=(Phiup-Phidown)*(walltmp(:,1)*elliptmp(:,0)-walltmp(:,0)*elliptmp(:,1)) / &
& denom ! Axial derivative
gtilde(:,2)=(Phiup-Phidown)*(walltmp(:,2)*elliptmp(:,0)-walltmp(:,0)*elliptmp(:,2)) / &
& denom ! Radial derivative
End If
End subroutine
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
END MODULE geometry
c4science · Help