Homec4science
Diffusion FENNECS lite (3Dmpimode) Browse

fennecslite/src/geometry_mod.f90
6b37315661933Dmpimode

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 extended b-splines and
!> can load the definition of the geometry from input files
!> This module is based on the theory by K. Hollig and book "Finite element methods with b-splines"
!> SIAM Frontiers in applied mathematics
!------------------------------------------------------------------------------
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 !< parameters defining the manufactured test solutionwhen negative weighttypes is used
real(kind=db):: z0
real(kind=db):: r0
real(kind=db):: Lz
real(kind=db):: Lr
end type
real(kind=db), save:: testkr=1, testkz=1
Logical, save:: nlweb=.true. !< use weighted extended b-splines and not only weighted b-splies
Integer, save :: walltype = 0 !< type of geometric weight to use (see readgeom)
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:: gridwdir(:,:) !< Stores the Dirichlet 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
Real(kind=db), save, allocatable:: fsolution(:,:) ! Stores the manufactured solution
type(mumps_mat):: etilde !> Matrix of extendend web splines definition 4.9 p48 of Hollig's book
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()!< Function evaluating the weight for Dirichelt boundary conditions
PROCEDURE(geomtot_eval), POINTER:: domain_weight => NULL() !< function giving the limits of the simulation domain
PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() !< Computes the parameter gtilde used to impose dirichlet boundary conditions with phi=uh+gtilde
PUBLIC:: geom_weight, dom_weight
ABSTRACT INTERFACE
SUBROUTINE gtilde_eval(r,z,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), OPTIONAL :: wupper(0:,:)
END SUBROUTINE
SUBROUTINE geomtot_eval(z,r,w,idwall)
USE constants
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(0:,:)
INTEGER, optional, INTENT(OUT):: idwall(:)
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, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight, testkr, testkz
contains
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Read the input parameters to initialize the geometry module from the standard input file
!> @param[in] Fileid Text file id of the input file containing namelists
!> @param[in] rnorm distance normalization constant
!> @param[in] splrz bspline structure used by the FEM comming form bspline library
!> @param[in] Potinn Normalized electric potential on the inner boundary
!> @param[in] Potout Normalized electric potential on the outer boundary
!---------------------------------------------------------------------------
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 ! normalisation variable for distances
Real(kind=db):: Potinn ! potential at inner electrode from basic
Real(kind=db):: Potout ! potential at outer electrode from basic
Integer:: Fileid, mpirank, ierr, istat
character(len=1000) :: line
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 definedwe take the defaults from the grid limits
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
SELECT CASE (abs(walltype))
CASE (2)
! coaxial cylinder and top ellipse with cylinder extensions
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w2
domain_weight=>geom_rvaschtot
CASE (3)
! Two ellipses with "parallel" tangents with cylinders
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w3
domain_weight=>geom_rvaschtot
CASE (4)
! Two ellipses with same radii with cylinders
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w4
domain_weight=>geom_rvaschtot
CASE (5)
! Two ellipses with same radii with cylinders and
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w5
domain_weight=>geom_rvaschtot
CASE (6)
! circular coaxial tilted ellipse right Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w6
domain_weight=>geom_rvaschtot
CASE (7)
! circular coaxial tilted ellipse right and left Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w7
domain_weight=>geom_rvaschtot
CASE (8)
! circular coaxial tilted ellipse right and left Dirichlet
total_gtilde=>gUpDown
Dirichlet_weight=>geom_w8
domain_weight=>geom_rvaschtot
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)
CASE (10)
! square section disc
total_gtilde=>gUpDown
domain_weight=>geom_rvaschtot
Dirichlet_weight=>geom_w10
CASE (11)
! square section disc
total_gtilde=>gUpDown
domain_weight=>geom_rvaschtot
Dirichlet_weight=>geom_w11
CASE (12)
! square section disc
total_gtilde=>gUpDown
domain_weight=>geom_rvaschtot
Dirichlet_weight=>geom_w12
CASE DEFAULT
! Ellipse as in gt170 standard weight and straight coaxial configs
total_gtilde=>gstd
Dirichlet_weight=>geom_weightstd
domain_weight=>geom_rvaschtot
END SELECT
! If we are lauching a test case, we load the test_pars variable used for the manufactured solution
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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Initialises the module and precomputes the matrix e used for extended splines
!> Classify every grid cell (inner, outer, boundary)
!> @param[in] splr2 bspline structure used by the FEM comming form bspline library
!> @param[in] vec1 Axial Meshgrid array to precompute the weigths on the crid points
!> @param[in] vec2 Radial Meshgrid array to precompute the weigths on the crid points
!---------------------------------------------------------------------------
SUBROUTINE geom_init(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(2)+1))
Allocate(rgrid(1:nrz(1)+1))
zgrid=spl2%sp2%knots(0:nrz(2))
rgrid=spl2%sp1%knots(0:nrz(1))
! create a table of the classification for each cell
Allocate(gridcelltype(nrz(1),nrz(2)))
Call classifycells(rgrid, zgrid, gridcelltype)
if (nlweb) then
! if we use extended splines, we need to build the e matrix linking inner and outer splines
Allocate(bsplinetype(nrank(1)*nrank(2)))
Call classifyspline(spl2, gridcelltype, bsplinetype)
Call buildetilde(spl2,bsplinetype,gridcelltype)
end if
! Precompute the domain and dirichlet weights at the grid/cell positions
ALLOCATE(gridwdir(0:2,(nrz(1)+1)*(nrz(2)+1)))
gridwdir=0
ALLOCATE(gridwdom(0:2,(nrz(1)+1)*(nrz(2)+1)))
gridwdom=0
ALLOCATE(gtilde(0:2,(nrz(1)+1)*(nrz(2)+1)))
gtilde=0
ALLOCATE(fsolution(0:2,(nrz(1)+1)*(nrz(2)+1)))
fsolution=0
Call geom_weight (vec1, vec2, gridwdir)
Call dom_weight (vec1, vec2, gridwdom)
! Precompute the gtilde at the grid cell position
CALL total_gtilde(vec1, vec2, gtilde, gridwdir)
! Precompute the manufactured solution at the grid cell position
if (walltype .lt. 0) then
CALL manufsolution(vec1, vec2, fsolution)
end if
end Subroutine geom_init
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Save this module run parameters to the h5 result file in the correct group
!> @param[in] File_handle h5 file id of the result file
!> @param[in] parentgroup h5 parent group where to save the simulation parameters
!> @param[in] rnorm distance normalization constant
!---------------------------------------------------------------------------
Subroutine geom_diag(File_handle, parentgroup, rnorm)
use mpi
Use futils
use weighttypes
Integer:: File_handle
Real(kind=db):: rnorm
Character(len=*):: parentgroup
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(parentgroup),"/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',transpose(gridwdom))
Call putarr(File_handle, trim(grpname)//'/dirichletweight',transpose(gridwdir))
Call putarr(File_handle, trim(grpname)//'/gtilde',transpose(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)
Call putarr(File_handle, trim(grpname)//'/fsolution',transpose(fsolution))
END IF
End subroutine geom_diag
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Construct the e matrix used to link inner and boundary b-splines
!> @param[in] spl2 bi-variate b-spline structure
!> @param[in] bsplinetype 1D array classifying the individual b-splines as inner, boundary, outer
!> @param[in] celltype 2D array classifying the individual grid cells as inner, boundary, outer
!---------------------------------------------------------------------------
Subroutine buildetilde(spl2, bsplinetype, celltype)
USE mumps_bsplines
Use basic, ONLY: rnorm, mpirank
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(rgrid(1:nrz(1)+1))
Allocate(zgrid(1:nrz(2)+1))
rgrid=spl2%sp1%knots(0:nrz(1))
zgrid=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.
! Compute the center of the 2D splines for the Lagrange interpolation
call calcsplinecenters(spl2%sp1,cr,hr)
call calcsplinecenters(spl2%sp2,cz,hz)
Do i=0,nrank(2)-1
rgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cr
zgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cz(i+1)
hmesh(i*nrank(1)+1:(i+1)*nrank(1))=hz(i+1)*hr
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 memory for etilde and its transpose
call init(nrank(1)*nrank(2),nbreducedspline,etilde)
call init(nrank(1)*nrank(2),nbreducedspline,etildet)
k=1
! Compute the terms eii for the inner b-splines
Do i=1,nrank(1)*nrank(2)
if(bsplinetype(i) .ne. 1) cycle ! span of this spline is almost completely outside D
! We consider only splines with
! one cell of bspline i completely in D
icellr=mod(i-1,nrank(1))+1
icellz=(i-1)/(nrank(1))+1
wgeomi=1
outer: do l=max(1,icellr-norder(1)),min(nrz(1),icellr)
do m=max(1,icellz-norder(2)),min(nrz(2),icellz)
if (celltype(l,m).eq.1) then
call geom_weight((rgrid(l)+rgrid(l+1))/2,(zgrid(m)+zgrid(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(icellr,icellz)=i
Ibsplinetype(i)=icellr+norder(1) .le. nrank(1) .and. icellz+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
! Compute the terms eij for the outer b-splines
!!$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)
icellr=mod(k-1,nrank(1))
icellz=(k-1)/(nrank(1))
jcellr=mod(j-1,nrank(1))
jcellz=(j-1)/(nrank(1))
indexdistance=real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db)
if(d(k) .gt. ((norder(1)+2)**2+(norder(2)+2)**2)*2 .and. mpirank .eq. 0)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
! Compute the Lagrange polynomia linking spline i and spline j
linkedspline(jcellr+1,jcellz+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*(zgridmesh(j)-zgridmesh(k+m*nrank(1)))/(zgridmesh(k+n*nrank(1))-zgridmesh(k+m*nrank(1)))
!eij=eij*real((jcellr-icellr-m),db)/real((n-m),db)
end do
do l=0,norder(1)
if( i .eq. l ) cycle
eij=eij*(rgridmesh(j)-rgridmesh(k+l))/(rgridmesh(k+i)-rgridmesh(k+l))
!eij=eij*real((jcellz-icellz-l),db)/real((i-l),db)
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))
! add the polynomia to the etilde matrix
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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Routine to compute the center of the 2d b-splines used in Lagrange interpolation
!> @param[in] spl 1D b-spline structure
!> @param[out] ctrs 1D array of spline centers
!> @param[out] heights 1D array of the maximum amplitude of each spline function
!---------------------------------------------------------------------------
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(spl%knots(i-order))
x2=spl%knots(i-1)-1e5*Epsilon(spl%knots(i-1))
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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Classify each grid cell as inner (1), boundary (0), or outer (-1) cell and store it in the
!> gridctype 2D array
!> @param[in] zgrid array of axial limits of the grid-cells
!> @param[in] rgrid array of radial limits of the grid-cells
!> @param[out] gridctype 2D array classifying the individual grid cells as inner, boundary, outer
!---------------------------------------------------------------------------
Subroutine classifycells(rgrid,zgrid,gridctype)
Real(kind=db):: zgrid(1:),rgrid(1:)
Integer:: gridctype(:,:)
Integer::i,j
gridctype=-1
!$OMP parallel do private(i,j)
Do j=1,size(zgrid,1)-1
Do i=1,size(rgrid,1)-1
! Determines the type inner/boundary/outer for each cell
Call classifycell(zgrid(j:j+1),rgrid(i:i+1),&
& gridctype(i,j))
End Do
End do
!$OMP end parallel do
End subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> For each spline i determines if its support is inside(1), partially inside (0) or outside (-1) of the simulation domain
!> @param[in] spl2 structure storing the bi-variate b-splines
!> @param[in] gridctype 2D array classifying the individual grid cells as inner, boundary, outer
!> @param[out] bsptype array of b-spline classification
!---------------------------------------------------------------------------
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 call the correct subroutine pointed by dirichlet_weight
!
!> @param[in] z axial position
!> @param[in] r radial position
!> @param[out] w geometric weight
!---------------------------------------------------------------------------
SUBROUTINE geom_weight0(r,z,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(r,z,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(0:size(w,1)-1,1:1)
ztmp=z
rtmp=r
call Dirichlet_weight(ztmp,rtmp,wtmp)
w=wtmp(:,1)
End SUBROUTINE geom_weight1
SUBROUTINE geom_weight2(r,z,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
!> This functions call the correct subroutine pointed by domain_weight
!> if the weight is negative return the id of the closest wall.
!
!> @param[in] z axial position
!> @param[in] r radial position
!> @param[out] w geometric weight
!> @param[out] idwall unique identifier of the boundary on which the particle is lost
!---------------------------------------------------------------------------
SUBROUTINE dom_weight0(r,z,w,idwall)
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)
INTEGER:: idwalltmp(1:1)
INTEGER, optional, INTENT(OUT):: idwall
ztmp=z
rtmp=r
call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp)
w=wtmp(1,1)
if (present(idwall)) idwall=idwalltmp(1)
End SUBROUTINE dom_weight0
SUBROUTINE dom_weight1(r,z,w,idwall)
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(0:size(w,1)-1,1:1)
INTEGER:: idwalltmp(1:1)
INTEGER, optional, INTENT(OUT):: idwall
ztmp=z
rtmp=r
call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp)
w=wtmp(:,1)
if (present(idwall)) idwall=idwalltmp(1)
End SUBROUTINE dom_weight1
SUBROUTINE dom_weight2(r,z,w,idwall)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(0:,:)
INTEGER, optional, INTENT(OUT):: idwall(:)
call domain_weight(z,r,w,idwall=idwall)
!
End SUBROUTINE dom_weight2
SUBROUTINE dom_weight3(r,z,w,idwall)
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:)
Real(kind=db):: wtmp(1:1,size(w,1))
INTEGER, optional, INTENT(OUT):: idwall(:)
call domain_weight(z,r,wtmp,idwall=idwall)
!Write(*,*) 'wtmp, size', wtmp, size(wtmp)
w=wtmp(1,:)
!
End SUBROUTINE dom_weight3
SUBROUTINE geom_weightstd(z,r,w,wupper)
! return the geometric weight for a coaxial configuration
! or central conductor + ellipse if walltype ==1
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(0:,:)
Real(kind=db), OPTIONAL:: wupper(0:,:)
Real(kind=db):: walltmp(0:size(w,1)-1,size(w,2)), elliptmp(0:size(w,1)-1,size(w,2))
Real(kind=db):: squareroot(size(w,2)), denom(size(w,2))
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 (present(wupper))then
w=walltmp
wupper=elliptmp
return
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,1) .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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Classify each grid cell as inner (1), boundary (0), or outer (-1) cell and store it in the
!> gridctype 2D array
!> @param[in] zgrid array of axial limits of the grid-cells
!> @param[in] rgrid array of radial limits of the grid-cells
!> @param[out] ctype classification of the cell
!---------------------------------------------------------------------------
SUBROUTINE classifycell(z, r, ctype)
! classify if cell is fully inside, outside or on the boundary of the domain
! by calculating the weight on each corner and the cell.
! It is assumed that the cells are sufficiently small such that there is no sharp edge entering the cell
! where an inner portion of only one cell edge is outside of the domain
real(kind=db), INTENT(IN):: r(2), z(2)
INTEGER, INTENT(OUT):: ctype
Real(kind=db)::weights(1:1,5)
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(reval,zeval,weights)
ctype=int(sign(1.0_db,weights(1,5)))
If(weights(1,1)*weights(1,2) .le. 0 ) then
ctype=0
return
End If
If(weights(1,1)*weights(1,3) .le. 0 ) then
ctype=0
return
End If
If(weights(1,2)*weights(1,4) .le. 0 ) then
ctype=0
return
End If
If(weights(1,3)*weights(1,4) .le. 0 ) then
ctype=0
return
End If
If(weights(1,3)*weights(1,2) .le. 0 ) then
ctype=0
return
End If
If(weights(1,1)*weights(1,4) .le. 0 ) then
ctype=0
return
End If
end subroutine
! ##################################################################
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> calculates the gauss quadrature integration points for the FEM method for any cell type
!> takes care of boundary cells as well and limit the integration boundaries accordingly
!> returns the gauss quadrature points and weights for cell (i,j)
!> @param[in] spl2 structure storing the bi-variate b-spline
!> @param[in] ngauss 2D array of number of gauss ponts in z and r
!> @param[in] i axial cell index
!> @param[in] j radial cell index
!> @param[out] xgauss 2D array of evaluation points z=xgauss(:,1), r=xgauss(:,2)
!> @param[out] xgauss 1D array of evaluation weights
!> @param[out] gausssize total number of gauss points for this cell
!> @param[out] celltype return the type of the cell inner (1), boundary (0), outer(-1)
!---------------------------------------------------------------------------
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, intent(out)::xgauss(:,:), wgauss(:)
Integer:: i,j, ngauss(2)
Real(kind=db),Allocatable:: zpoints(:)
Real(kind=db):: rg(ngauss(1)),zg(ngauss(2)), wrg(ngauss(1)), wzg(ngauss(2))
Integer:: k, l, directiondown, directionup, nbzpoints, direction, ctype
Logical:: hasmaxpoint
Real(kind=db):: xptup, xptdown,wlu,wld, rmin, rmax, zmin, zmax
type(spline1d):: splz, splr
splr=spl2%sp1
splz=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(PRESENT(celltype)) celltype=ctype
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, rg, wrg)
CALL get_gauss(spl2%sp2, ngauss(2), j, zg, wzg)
! 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)=rg
xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=zg(k)
wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wzg(k)*wrg
END DO
Else If(ctype.eq.0) then ! we have a boundary cell
directiondown=1
directionup=1
! We check if the boundary goes through the cell upper and lower limit
Call Find_crosspointdico(zgrid(j:j+1),rgrid(i),xptdown,directiondown)
Call Find_crosspointdico(zgrid(j:j+1),rgrid(i+1),xptup,directionup)
call dom_weight(rgrid(i),zgrid(j),wld)
call dom_weight(rgrid(i+1),zgrid(j),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(j:j+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(j), xptdown,zgrid(j+1)/)
else
nbzpoints=2
Allocate(zpoints(nbzpoints))
if(wld.gt.0)then
! the upper left corner is inside
zpoints=(/zgrid(j),xptdown/)
else
zpoints=(/xptdown,zgrid(j+1)/)
end if
end if
else
if(wld .gt. 0) then
! the lower left corner is inside
nbzpoints=3
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(j), xptup,zgrid(j+1)/)
else
nbzpoints=2
Allocate(zpoints(nbzpoints))
if(wlu.gt.0)then
! the upper left corner is inside
zpoints=(/zgrid(j),xptup/)
else
zpoints=(/xptup,zgrid(j+1)/)
end if
end if
end if
Case(2)
nbzpoints=2
Allocate(zpoints(nbzpoints))
zpoints=(/zgrid(j), zgrid(j+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
Allocate(xgauss(ngauss(1)*ngauss(2),2))
Allocate(wgauss(ngauss(1)*ngauss(2)))
gausssize=ngauss(1)*ngauss(2)
! Compute gauss points
!CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg)
!Call gauleg(zpoints(l),zpoints(l+1),zg,wzg,ngauss(1))
call gauleg(rgrid(i),rgrid(i+1),rg,wrg,ngauss(1))
Do k=1,ngauss(1)
! We test if the lower or upper side is in the domain
call dom_weight(rg(k),zgrid(j),wld)
zmin=zgrid(j)
if (wld .lt. 0) zmin = zgrid(j+1)
direction=1
Call Find_crosspointdico(zgrid(j:j+1),rg(k),zmax,direction) ! We compute the radial limits at each z position
Call gauleg(min(zmin,zmax),max(zmin,zmax),zg,wzg,ngauss(2)) ! We obtain the gauss w and pos for these boundaries
if(direction .eq. 0.and. wld .lt. 0) then
wzg=0
end if
xgauss(k : ngauss(2)*ngauss(1) : ngauss(2),2)= zg
xgauss(k : ngauss(2)*ngauss(1) : ngauss(2),1)= rg(k)
wgauss(k : ngauss(2)*ngauss(1) : ngauss(2)) = wzg*wrg(k)
End do
return
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(2))
Do k=1,ngauss(2)
! We test if the lower or upper side is in the domain
call dom_weight(rgrid(i),zg(k),wld)
rmin=rgrid(i)
if (wld .lt. 0) rmin = rgrid(i+1)
direction=2
Call Find_crosspointdico(rgrid(i:i+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(1)) ! We obtain the gauss w and pos for these boundaries
if(direction .eq. 0.and. wld .lt. 0) then
wrg=0
end if
xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(2),1) = rg
xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(2),2) = zg(k)
wgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(2)) = wrg*wzg(k)
End do
End Do
Else
gausssize=0
Allocate(xgauss(1:1,2))
Allocate(wgauss(1:1))
! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration
xgauss(1,2)=(zgrid(j)+zgrid(j+1))*0.5
xgauss(1,1)=(rgrid(i)+rgrid(i+1))*0.5
wgauss(1)=0
End If
End Subroutine calc_gauss
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> computes quickly if a particle is inside the domain D and returns the id of the closest boundary
!> if the particle is outside the domain
!> @param[in] z axial particle position
!> @param[in] r radial particle position
!> @param[in] i axial cell index
!> @param[in] j radial cell index
!> @param[out] idwall id of the closest boundary if the particle is outside the domain
!> @param[out] inside .true. if the particle is inside D
!---------------------------------------------------------------------------
subroutine is_insidegeom(r,z,i,j,idwall,inside)
REAL(kind=db), intent(in):: z,r
INTEGER,INTENT(in):: i,j
INTEGER:: idwall
logical:: inside
Real(kind=db):: weight
if(i.gt.size(gridcelltype,1)-1.or. i.lt.0 .or. j.gt.size(gridcelltype,2)-1.or. j.lt.0 )then
inside=.false.
idwall=0
return
end if
select case(gridcelltype(i+1,j+1))
case(-1)
inside=.false.
idwall=0
case(1)
inside=.true.
idwall=0
case(0)
call dom_weight(r, z, weight, idwall=idwall)
inside=weight.gt.0
end select
end subroutine is_insidegeom
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> calculates the boundary limit between x(1) and x(2) using Newton's method
!> @param[in] x axial/radial initial guesses for the boundary position
!> @param[in] y radial/axial position at which the boundary is searched
!> @param[out] xpt position of the boundary found
!> @param[inout] direction (1) searches boundary along z, (2) searches boundary along r (0) no boundary found between x(1) and x(2)
!---------------------------------------------------------------------------
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
! 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(y,xptm,wold)
! direction=2 finds cross-point along r
if(direction .eq.2) Call dom_weight(xptm,y,wold)
if(direction .eq.1) Call dom_weight(y,xpt,w)
if(direction .eq.2) Call dom_weight(xpt,y,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.2) Call dom_weight(xpt,y,w)
if(direction .eq.1) Call dom_weight(y,xpt,w)
End do
if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0
End Subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> calculates the boundary limit between x(1) and x(2) using dichotomy method
!> @param[in] x axial/radial initial guesses for the boundary position
!> @param[in] y radial/axial position at which the boundary is searched
!> @param[out] xpt position of the boundary found
!> @param[inout] direction (1) searches boundary along z, (2) searches boundary along r (0) no boundary found between x(1) and x(2)
!---------------------------------------------------------------------------
Subroutine Find_crosspointdico(x,y,xpt, direction)
! calculates the boundary limit between x(1) and x(2) using dichotomy method
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(2)
! direction=1 finds cross-point along r
Call dom_weight(x1,y,w1)
Call dom_weight(x2,y,w2)
case(1)
! direction=2 finds cross-point along z
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 .or. (w1.eq.0 .and.w2.eq.0)) then
direction=0
xpt=x2
return
End If
! Find the cross-point
Do while( i .le. 1000 .and. abs((x2-x1)/(x(2)-x(1))).gt.1e-14)
x3=0.5*(x1+x2)
i=i+1
select case(direction)
case(2)
! direction=1 finds cross-point along r
Call dom_weight(x3,y,w3)
case(1)
! direction=2 finds cross-point along z
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 .gt. x(2) .or. xpt .lt. x(1) ) direction=0
End Subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> returns the total weight which is the same as the Dirichlet weight for a boundary
!> defined with Rvaschev functions
!> @param[in] z axial position
!> @param[in] r radial position
!> @param[out] w dirichlet weight
!> @param[out] idwall id of the closest boundary if the particle is outside the domain
!---------------------------------------------------------------------------
SUBROUTINE geom_rvaschtot(z,r,w,idwall)
! returns the total weight which is the same as the Dirichlet weight for a boundary
! defined with Rvaschev functions
Use splinebound, ONLY: spline_w
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(0:,:)
Real(kind=db), allocatable:: w2(:,:)
Real(kind=db), allocatable:: w3(:,:)
INTEGER, optional, INTENT(OUT):: idwall(:)
INTEGER:: sw2
sw2=size(w,2)
if(present(idwall)) then
idwall=0
allocate(w2(1:size(w,1),1:size(w,2)))
allocate(w3(1:size(w,1),1:size(w,2)))
call Dirichlet_weight(z,r,w2,w3)
! gives total weight
where(w2(1,:).le.0) idwall=1
where(w3(1,:).le.0) idwall=2
call Combine(w2, w3, w, -1)
else
call Dirichlet_weight(z,r,w)
end if
End SUBROUTINE geom_rvaschtot
Subroutine gstd(r,z,gtilde,w)
! g tilde function added to rhs of poisson solver for the standard coaxial configuration
! for the default weight function
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(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,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,1) .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(r,z,gtilde,w)
! g tilde function added to rhs of poisson solver by combining two boundaries set at different potentials
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(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,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,1) .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(r,z,gtilde,w)
! calculates the Poisson gtilde term for testing the solver on a new geometry
! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2
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(0:size(gtilde,1)-1,size(gtilde,2))
!
!gtilde=0
!return
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)
!gtilde(0,:)=1+r**6+2*z**6
!gtilde(0,:)=exp(r*z)
If(size(gtilde,1) .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,:)
! gtilde(2,:)=6*r**5*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:)
! gtilde(1,:)=12*z**5*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:)
! !gtilde(2,:)=z*exp(r*z)*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:)
! !gtilde(1,:)=r*exp(r*z)*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:)
End If
gtilde(0,:)=gtilde(0,:)*(1-wtmp(0,:))
End subroutine
Subroutine manufsolution(r,z,gtilde)
! calculates the Poisson gtilde term for testing the solver on a new geometry
! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: gtilde(0:,:)
!
!gtilde=0
!return
gtilde(0,:)=(sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))+2)
!gtilde(0,:)=1+r**6+2*z**6
!gtilde(0,:)=exp(r*z)
If(size(gtilde,1) .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))
gtilde(2,:)=-pi/(test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*cos(pi*(r-test_pars%r0)/(test_pars%Lr))
! gtilde(2,:)=6*r**5*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:)
! gtilde(1,:)=12*z**5*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:)
! !gtilde(2,:)=z*exp(r*z)*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:)
! !gtilde(1,:)=r*exp(r*z)*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:)
End If
gtilde(0,:)=gtilde(0,:)
End subroutine
Subroutine ftest(r,z,f)
! calculates the Poisson source term for testing the solver on a new geometry
! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT)::f(0:,:)
!f(0,:)=1
!return
!
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,:)=-36*r**4-60*z**4
!f(0,:)=-(z*exp(r*z)/r+z**2*exp(r*z)+r**2*exp(r*z))
End Subroutine ftest
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Reduce the Finite element matrix in the LHS from full spline space to reduced web-spline space
!> If A is the standard FEM matrix on the full grid space then the web-spline FEM matrix is etilde A etildet
!> with etildet the transposed etilde matrix
!> @param[in] full full FEM matrix on the weighted b-spline domain
!> @param[out] reduced reduced matrix on the web-spline space
!---------------------------------------------------------------------------
subroutine Reducematrix(full,reduced)
!
use mumps_bsplines
type(mumps_mat):: full, reduced, tempmat1, tempmat2
Integer:: fullrank, reducedrank
Integer:: i,j,k, dimsmult1(2), dimsmult2(2)
Real(kind=db):: val
logical:: error
call to_mat(full)
fullrank=full%rank
reducedrank=nbreducedspline
!WRITE(*,*) "# web-splines", nbreducedspline
!WRITE(*,*) "# splines", fullrank
dimsmult1=(/fullrank,fullrank/)
dimsmult2=(/fullrank,reducedrank/)
call init(reducedrank,2,reduced,nlpos=.false.)
call init(fullrank,2,tempmat1)
call init(reducedrank,2,tempmat2)
tempmat1=mmx_mumps_mat_loc(full,etildet,dimsmult1,dimsmult2,.false.)
tempmat2=mmx_mumps_mat_loc(etildet,tempmat1,dimsmult2,dimsmult2,.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
call to_mat(reduced)
end subroutine
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Com putes the sparse matrix matrix multiplication using MKL A B=C
!> @param[in] mata sparse matrix A in csr representation
!> @param[in] matb sparse matrix B in csr representation
!> @param[in] ranka Rank of matrix A
!> @param[in] rankb Rank of matrix B
!> @param[in] transpose Is A transposed
!> @param[out] matc sparse matrix C in csr representation
!---------------------------------------------------------------------------
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))
! Compute the size of the resulting matrix and the indices of the non-zero elements
request=1
sort=7
CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), &
& mata%val, mata%cols, mata%irow, &
& matb%val, matb%cols, matb%irow, &
& matc%val, matc%cols, matc%irow, &
& 1, info)
if(info .gt. 0) WRITE(*,'(a,i4)') " 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))
matc%val=0
request=2
! Do the actual matrix matrix multiplication
CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), &
& mata%val, mata%cols, mata%irow, &
& matb%val, matb%cols, matb%irow, &
& matc%val, matc%cols, matc%irow, &
& 1, info)
if(info .ne. 0) WRITE(*,'(a,i4)') " 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