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

fennecs/src/geometry_mod.f90
e4669700ed0dgeom_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
USE splinebound
use weighttypes
IMPLICIT NONE
type innerspline
Integer, Allocatable:: k(:) ! Index in reduced set
real(kind=db), Allocatable:: weight(:) ! geomtric weight at relevant cell
end type innerspline
type test_params
real(kind=db):: z0
real(kind=db):: r0
real(kind=db):: Lz
real(kind=db):: Lr
end type
Integer, save:: testkr=1, testkz=1
Logical, save:: nlweb=.false.
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
Integer, save, Allocatable:: linkedspline(:,:) ! Array containing the lowerleft linked spline in case of boundary spline
Real(kind=db), save, allocatable:: gridwgeom(:,:) ! Stores the geometric weight at the grid points
Real(kind=db), save, allocatable:: gridwdom(:,:) ! Stores the domain 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
type(test_params), save:: test_pars
PROCEDURE(geom_eval), POINTER:: dirichlet_weight => NULL()
PROCEDURE(geom_eval), POINTER:: domain_weight => NULL()
PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL()
PUBLIC:: geom_weight, dom_weight
ABSTRACT INTERFACE
SUBROUTINE gtilde_eval(z,r,g,w)
USE constants
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: g(:,0:)
Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:)
END SUBROUTINE
SUBROUTINE geom_eval(z,r,w,wupper)
USE constants
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db), INTENT(OUT), OPTIONAL:: wupper(:,0:)
END SUBROUTINE
END INTERFACE
INTERFACE geom_weight
MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2
END INTERFACE geom_weight
INTERFACE dom_weight
MODULE PROCEDURE dom_weight0, dom_weight1, dom_weight2, dom_weight3
END INTERFACE dom_weight
NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, z_a, z_b ,walltype, L_r, L_z, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight, testkr, testkz
contains
SUBROUTINE read_geom(Fileid, rnorm, splrz, Potinn, Potout)
Use mpi
Use bsplines
use basic, ONLY: phinorm
use weighttypes
type(spline2d):: splrz
Real(kind=db):: rnorm, Potinn, Potout
Integer:: Fileid, mpirank, ierr, istat
character(len=1000) :: line
!type(spline_boundary):: border0
!Real(kind=db):: cpoints(2,3),t(1), distance, pos
!Real(kind=db), ALLOCATABLE::result(:,:,:)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
Rewind(Fileid)
READ(Fileid, geomparams, iostat=istat)
if (istat.gt.0) then
if(mpirank .eq. 0) then
backspace(Fileid)
read(Fileid,fmt='(A)') line
write(*,'(A)') &
'Invalid line in geomparams: '//trim(line)
end if
call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
stop
end if
if(mpirank .eq. 0) WRITE(*, geomparams)
!! Normalizations and initialization of geometric variables
r_a=r_a/rnorm
r_b=r_b/rnorm
z_a=z_a/rnorm
z_b=z_b/rnorm
r_bLeft=r_bLeft/rnorm
r_bRight=r_bRight/rnorm
if(r_a .eq. 0 .and. r_b .eq.0) then
!! in case no geom_params have been defined
r_a=splrz%sp2%knots(0)
r_b=splrz%sp2%knots(splrz%sp2%nints)
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
Phidown=Potinn
Phiup=Potout
L_r=L_r/rnorm
L_z=L_z/rnorm
SELECT CASE (abs(walltype))
CASE (2)
! coaxial cylinder and top ellipse with cylinder extensions
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w2
domain_weight=>geom_w2
CASE (3)
! Two ellipses with "parallel" tangents with cylinders
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w3
domain_weight=>geom_w3
CASE (4)
! Two ellipses with same radii with cylinders
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w4
domain_weight=>geom_w4
CASE (5)
! Two ellipses with same radii with cylinders and
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w5
domain_weight=>geom_w5
CASE (6)
! circular coaxial tilted ellipse right Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w6
domain_weight=>geom_w6
CASE (7)
! circular coaxial tilted ellipse right and left Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w7
domain_weight=>geom_w7
CASE (8)
! circular coaxial tilted ellipse right and left Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w8
domain_weight=>geom_w8
CASE (9)
! Geometry defined as a spline curve
total_gtilde=>gspline
Dirichlet_weight=>geom_spline
domain_weight=>geom_splinetot
call read_splinebound(Fileid,the_domain, splrz, rnorm, phinorm, Potinn,Potout)
CASE (10)
! square section disc
total_gtilde=>gUpDown
domain_weight=>geom_w10
Dirichlet_weight=>geom_w10
CASE DEFAULT
! Ellipse as in gt170 standard weight and straight coaxial configs
total_gtilde=>gstd
Dirichlet_weight=>geom_weightstd
domain_weight=>geom_weightstd
END SELECT
! If we are lauching a test case, we load the test_pars variable
if(walltype.lt.0) then
test_pars%Lr=(splrz%sp2%knots(splrz%sp2%nints)-splrz%sp2%knots(0))/testkr
test_pars%Lz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/testkz
test_pars%r0=0.5*(splrz%sp2%knots(splrz%sp2%nints)+splrz%sp2%knots(0))
test_pars%z0=0.5*(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0))
total_gtilde=>gtest
end if
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)
end if
ALLOCATE(gridwgeom((nrz(1)+1)*(nrz(2)+1),0:2))
gridwgeom=0
ALLOCATE(gridwdom((nrz(1)+1)*(nrz(2)+1),0:2))
gridwdom=0
ALLOCATE(gtilde((nrz(1)+1)*(nrz(2)+1),0:2))
gtilde=0
Call geom_weight (vec1, vec2, gridwgeom)
CALL total_gtilde(vec1, vec2, gtilde, gridwgeom)
Call dom_weight (vec1, vec2, gridwdom)
end Subroutine init_geom
Subroutine geom_diag(File_handle, str, rnorm)
Use mpi
Use futils
use weighttypes
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_a", z_a*RNORM)
Call attach(File_handle, trim(grpname), "z_b", z_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), "L_r", test_pars%Lr*RNORM)
Call attach(File_handle, trim(grpname), "L_z", test_pars%Lz*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',gridwdom)
Call putarr(File_handle, trim(grpname)//'/dirichletweight',gridwgeom)
Call putarr(File_handle, trim(grpname)//'/gtilde',gtilde)
Call putarr(File_handle, trim(grpname)//'/ctype',gridcelltype)
Call putarr(File_handle, trim(grpname)//'/linked_s',linkedspline)
Call putarr(File_handle, trim(grpname)//'/bsplinetype',bsplinetype)
END IF
End subroutine geom_diag
Subroutine buildetilde(spl2, bsplinetype, celltype)
USE mumps_bsplines
Use basic, ONLY: rnorm
type(spline2d):: spl2
Integer:: bsplinetype(:)
Integer:: celltype(1:,1:)
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(:)
integer, allocatable:: igridmesh(:), jgridmesh(:)
!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(linkedspline(nrank(1),nrank(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(igridmesh(nrank(1)*nrank(2)), jgridmesh(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
igridmesh(i*nrank(1)+1:(i+1)*nrank(1))=(/ (j,j=0,nrank(1)-1)/)
jgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=i
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
linkedspline(icellz,icellr)=i
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 positive splines in this spline domain 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)! calculate distance between center of Interior splines and spline j
d=(zgridmesh(j)-zgridmesh)**2+(rgridmesh(j)-rgridmesh)**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=real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db)
if(indexdistance .gt. ((norder(1)+1)**2+(norder(2)+1)**2)*2)then
Write(*,'(a)') 'Warning on system conditioning, the number of radial or axial points could be too low!'
Write(*,'(a,1f6.2,a,2(1pe12.4))') 'Distance found: ', sqrt(indexdistance), ' at (z,r): ',zgridmesh(j)*rnorm, rgridmesh(j)*rnorm
!stop
end if
linkedspline(jcellz+1,jcellr+1)=k
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)))
!eij=eij*(jcellr-icellr-m)/(n-m)
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*(jcellz-icellz-l)/(i-l)
end do
linkedi=innersplinelist%k(k+i+n*nrank(1)) ! equivalent to findloc, necessary for ifort 17
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.300)
!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).lt.1e-13) 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), ndegree(2)
Integer:: i, j, mu, imin, imax, jmin, jmax
Integer, Allocatable:: splinespan(:,:)
call get_dim(spl2, nrank, nrz, ndegree)
! by default, all splines have part of their support outside the domain D
bsptype=0
Allocate(splinespan(ndegree(1)+1,ndegree(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 for correct behavior in boundaries
splinespan=-1
! find the axial span of this spline in cell indices
imin=max(1,i-ndegree(1))
imax=min(nrz(1),i)
! find the radial span of this spline in cell indices
jmin=max(1,j-ndegree(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
! if at least one cell is fully in the domain the spline is an inside spline
bsptype(mu)=1
else if (.not. ANY( splinespan==0 )) then
! if all the cells are outside, this is an outside spline
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 Dirichlet_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 Dirichlet_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 Dirichlet_weight(z,r,w)
!
End SUBROUTINE geom_weight2
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> general placeholder function used in fields_mod to compute the domain weight for weighted_bsplines
!---------------------------------------------------------------------------
SUBROUTINE dom_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 domain_weight(ztmp,rtmp,wtmp)
w=wtmp(1,1)
End SUBROUTINE dom_weight0
SUBROUTINE dom_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 domain_weight(ztmp,rtmp,wtmp)
w=wtmp(1,:)
End SUBROUTINE dom_weight1
SUBROUTINE dom_weight2(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
call domain_weight(z,r,w)
!
End SUBROUTINE dom_weight2
SUBROUTINE dom_weight3(z,r,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:)
Real(kind=db):: wtmp(size(w,1),1:1)
call domain_weight(z,r,wtmp)
w=wtmp(:,1)
!
End SUBROUTINE dom_weight3
SUBROUTINE geom_weightstd(z,r,w,wupper)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db), OPTIONAL:: wupper(:,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 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)
! check if cell is fully inside, outside or on the boundary of the domain
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 dom_weight(zeval,reval,weights)
ctype=int(sign(1.0_db,weights(5,1)))
If(weights(1,1)*weights(2,1) .le. 0 ) then
ctype=0
return
End If
If(weights(1,1)*weights(3,1) .le. 0 ) then
ctype=0
return
End If
If(weights(2,1)*weights(4,1) .le. 0 ) then
ctype=0
return
End If
If(weights(3,1)*weights(4,1) .le. 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,wlu,wld, 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
! We check if the boundary goes through the cell upper and lower limit
Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown)
Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup)
call dom_weight(zgrid(i),rgrid(j),wld)
call dom_weight(zgrid(i),rgrid(j+1),wlu)
select case ( directionup+directiondown)
Case (0) ! The intersections are only on the left and right cell boundaries
! or the upper and lower limits are full boundaries
nbzpoints=2
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(i),zgrid(i+1)/)
Case(1)
if(directiondown.eq.1)then
if( wlu .gt. 0) then
! the lower left corner is inside
nbzpoints=3
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(i), xptdown,zgrid(i+1)/)
else
nbzpoints=2
Allocate(zpoints(nbzpoints))
if(wld.gt.0)then
! the upper left corner is inside
zpoints=(/zgrid(i),xptdown/)
else
zpoints=(/xptdown,zgrid(i+1)/)
end if
end if
else
if(wld .gt. 0) then
! the lower left corner is inside
nbzpoints=3
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(i), xptup,zgrid(i+1)/)
else
nbzpoints=2
Allocate(zpoints(nbzpoints))
if(wlu.gt.0)then
! the upper left corner is inside
zpoints=(/zgrid(i),xptup/)
else
zpoints=(/xptup,zgrid(i+1)/)
end if
end if
end if
Case(2)
nbzpoints=4
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(i),min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /)
!If(wld.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /)
!else
! nbzpoints=2
! Allocate(zpoints(nbzpoints))
! If(wld.ge.0) zpoints=(/zgrid(i),min(xptdown,xptup) /)
! If(wld.lt.0) zpoints=(/max(xptdown,xptup), zgrid(i+1) /)
!end if
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 dom_weight(zg(1),rgrid(j),wld)
rmin=rgrid(j)
if (wld .le. 0) rmin = rgrid(j+1)
Do k=1,ngauss(1)
direction=2
Call Find_crosspointdico((/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
Real(kind=db):: w1, w2, w3
Real(kind=db):: x1, x2, x3
Integer, Intent(INOUT):: direction
Integer:: i
! calculates the position of the boundary ( where the weight changes sign )
! between x(1) and x(2) at 2nd coordinate y
xptm=x(1)
xpt=x(2)
i=0
! direction=1 finds cross-point along z
if(direction .eq.1) Call dom_weight(xptm,y,wold)
! direction=2 finds cross-point along r
if(direction .eq.2) Call dom_weight(y,xptm,wold)
if(direction .eq.1) Call dom_weight(xpt,y,w)
if(direction .eq.2) Call dom_weight(y,xpt,w)
! if the weight doesn't change sign there is no cross point
if(w*wold.gt.0) then
direction=0
return
End If
! Find the cross-point
Do while( i .le.100 .and. abs(w).gt.1e-9)
!Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold
xptm=xpt-w*(xpt-xptm)/(w-wold)
wold=w
temp=xptm
xptm=xpt
xpt=temp
i=i+1
if(direction .eq.1) Call dom_weight(xpt,y,w)
if(direction .eq.2) Call dom_weight(y,xpt,w)
End do
if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0
End Subroutine
Subroutine Find_crosspointdico(x,y,xpt, direction)
Real(kind=db):: x(2), y
Real(kind=db):: xpt
Real(kind=db):: w1, w2, w3
Real(kind=db):: x1, x2, x3
Integer, Intent(INOUT):: direction
Integer:: i
! calculates the position of the boundary ( where the weight changes sign )
! between x(1) and x(2) at 2nd coordinate y
x1=x(1)
x2=x(2)
i=0
select case(direction)
case(1)
! direction=1 finds cross-point along z
Call dom_weight(x1,y,w1)
Call dom_weight(x2,y,w2)
case(2)
! direction=2 finds cross-point along r
Call dom_weight(y,x1,w1)
Call dom_weight(y,x2,w2)
end select
! if the weight doesn't change sign there is no cross point
if(w1*w2.gt.0) then
direction=0
xpt=x2
return
End If
! Find the cross-point
Do while( i .le. 500 .and. abs(x2-x1).gt.1e-13)
x3=0.5*(x1+x2)
i=i+1
select case(direction)
case(1)
! direction=1 finds cross-point along z
Call dom_weight(x3,y,w3)
case(2)
! direction=2 finds cross-point along r
Call dom_weight(y,x3,w3)
end select
if(w1*w3.gt.0)then! we are in a region were there is no change of sign
x1=x3
w1=w3
else
x2=x3
w2=w3
end if
!if(abs(w3).lt.1e-14.and. w3 .ge.0)Exit
End do
xpt=x3
if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0
End Subroutine
Subroutine gstd(z,r,gtilde,w)
! 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), INTENT(IN),OPTIONAL::w(:,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))
if (above1.eq.0) then
gtilde=0
gtilde(:,0)=Phiup
RETURN
end if
! 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 gUpDown(z,r,gtilde,w)
! 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), INTENT(IN),OPTIONAL::w(:,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
call Dirichlet_weight(z,r,belowtmp,abovetmp)
! 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 gUpDown
Subroutine gtest(z,r,gtilde,w)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(:,0:)
Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:)
Real(kind=db):: wtmp(size(gtilde,1),0:size(gtilde,2)-1)
!
if(present(w)) then
wtmp=w
else
call Dirichlet_weight(z,r,wtmp)
end if
gtilde(:,0)=(sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))+2)
If(size(gtilde,2) .gt. 1) then ! first derivative
gtilde(:,1)=pi/(test_pars%Lz)*cos(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(:,0))-wtmp(:,1)*gtilde(:,0)
gtilde(:,2)=pi/(test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*cos(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(:,0))-wtmp(:,2)*gtilde(:,0)
End If
gtilde(:,0)=gtilde(:,0)*(1-wtmp(:,0))
End subroutine
Subroutine ftest(z,r,f)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT)::f(:,0:)
!
f(:,0)=(pi/test_pars%Lz)**2*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))&
& + (pi/test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))&
& *( -1/r*(cos(pi*(r-test_pars%r0)/(test_pars%Lr))) + (pi/test_pars%Lr)*sin(pi*(r-test_pars%r0)/(test_pars%Lr)))
!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, tempmat1, tempmat2
Integer:: fullrank, reducedrank
Integer:: i,j, k
call to_mat(full)
fullrank=full%rank
reducedrank=nbreducedspline
call init(reducedrank,2,reduced)
call init(fullrank,2,tempmat1)
call init(reducedrank,2,tempmat2)
tempmat1=mmx_mumps_mat_loc(full,etildet,(/fullrank,fullrank/),(/fullrank,reducedrank/),.false.)
tempmat2=mmx_mumps_mat_loc(etildet,tempmat1,(/fullrank,reducedrank/),(/fullrank,reducedrank/),.true.)
do i=1,reducedrank
do k=tempmat2%irow(i),tempmat2%irow(i+1)-1
j=tempmat2%cols(k)
call putele(reduced, i, j, tempmat2%val(k))
end do
end do
end subroutine
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
FUNCTION mmx_mumps_mat_loc(mata, matb, ranka, rankb, transpose) RESULT(matc)
!
! Return product matc=mata*matb
! All matrices are represented in csr sparse representation
!
Use mumps_bsplines
TYPE(mumps_mat) :: mata,matb,matc
INTEGER:: ranka(2), rankb(2) ! nb of (rows,columns) for each matrix
INTEGER :: n, info, m, request, sort
LOGICAL:: transpose
Character(1):: T
!
m=ranka(1)
n=ranka(2)
T='N'
if (transpose)then
T='T'
m=ranka(2)
n=ranka(1)
end if
!#ifdef MKL
!
if(associated(matc%val)) deallocate(matc%val)
if(associated(matc%cols)) deallocate(matc%cols)
if(associated(matc%irow)) deallocate(matc%irow)
allocate(matc%irow(m+1))
allocate(matc%cols(1))
allocate(matc%val(1))
request=1
sort=7
CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), &
& mata%val, mata%cols, mata%irow(1), &
& matb%val, matb%cols, matb%irow(1), &
& matc%val, matc%cols, matc%irow(1), &
& 1, info)
if(info .gt. 0) WRITE(*,'(a,i)') " Error in mmx_mumps_mat_loc: ", info
if(associated(matc%val)) deallocate(matc%val)
if(associated(matc%cols)) deallocate(matc%cols)
allocate(matc%val(matc%irow(m+1)-1))
allocate(matc%cols(matc%irow(m+1)-1))
request=2
CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), &
& mata%val, mata%cols, mata%irow(1), &
& matb%val, matb%cols, matb%irow(1), &
& matc%val, matc%cols, matc%irow(1), &
& 1, info)
if(info .ne. 0) WRITE(*,'(a,i)') " Error in mmx_mumps_mat_loc: ", info
if (transpose)then
matc%rank=ranka(2)
else
matc%rank=ranka(1)
end if
!
END FUNCTION mmx_mumps_mat_loc
END MODULE geometry
c4science · Help