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

fennecs/src/fields_mod.f90
e3b71186ea57geom_testpardiso

Last Change
Raw File

fields_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: beam
!
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for initializing the magnetic field and solve the Poisson equation
!------------------------------------------------------------------------------
MODULE fields
USE constants
USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, splrz, nlperiod, phinorm, nlPhis, nrank, rhs, mpirank, mpisize, step, it2d
USE beam, ONLY: parts
USE bsplines
USE mumps_bsplines
USE mpi
IMPLICIT NONE
DOUBLE PRECISION, allocatable, SAVE :: matcoef(:,:), sol(:), vec1(:), vec2(:)
TYPE(mumps_mat), SAVE :: femat ! Finite Element Method matrix
CONTAINS
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid
!
!---------------------------------------------------------------------------
SUBROUTINE fields_init
USE basic, ONLY: pot, nlperiod, nrank, rhs
USE bsplines
USE mumps_bsplines
USE mpihelper
INTEGER :: nrz(2)
! Set up 2d spline splrz used in the FEM
CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod)
! Calculate dimension of splines
nrz(1)=nz
nrz(2)=nr
CALL get_dim(splrz, nrank, nrz, femorder)
ALLOCATE(matcoef(nrank(1),nrank(2)))
ALLOCATE(pot((nr+1)*(nz+1)))
ALLOCATE(rhs(nrank(1)*nrank(2)))
ALLOCATE(sol(nrank(1)*nrank(2)))
ALLOCATE(Br((nr+1)*(nz+1)),Bz((nr+1)*(nz+1)))
ALLOCATE(Athet((nr+1)*(nz+1)))
ALLOCATE(Er((nr+1)*(nz+1)),Ez((nr+1)*(nz+1)))
! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed)
CALL magnet
! Calculate and factorise FEM matrix (depends only on mesh)
CALL fematrix
CALL factor(femat)
! Initialize the MPI column overlap type for rhs
CALL init_coltypempi(nrank(2), nrank(1), MPI_DOUBLE_PRECISION, rhsoverlap_type)
END SUBROUTINE fields_init
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Construct the right hand side vector used in the FEM Poisson solver
!
!---------------------------------------------------------------------------
SUBROUTINE rhscon
USE bsplines
USE mumps_bsplines
USE mpi
USE basic, ONLY: omegap, omegac, V, potinn, potout, ierr, nrank, rhs, nplasma
USE beam, ONLY: parts
USE mpihelper
INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch
INTEGER,DIMENSION(:),ALLOCATABLE::zleft, rleft
DOUBLE PRECISION :: chargecoeff
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:,:), fun2(:,:,:)
nbunch=100 ! Particle bunch size used when calling basfun
ALLOCATE(zleft(nbunch),rleft(nbunch))
ALLOCATE(fun(1:femorder(1)+1,0:0,nbunch),fun2(1:femorder(2)+1,0:0,nbunch)) ! Arrays keeping values of b-splines at gauss node
rhs=0 ! Reset rhs
chargecoeff=0.5/pi*omegap/omegac*V/nplasma ! Normalized charge density simulated by each macro particle
! Assemble rhs
IF (parts%Nptot .ne. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, i, iend, irow, jcol, mu, k) PRIVATE(fun, fun2) REDUCTION(+:rhs)
DO i=1,parts%Nptot,nbunch
! Avoid segmentation fault by accessing non relevant data
iend=min(i+nbunch-1,parts%Nptot)
! Localize the particle
CALL locintv(splrz%sp1, parts%Z(i:iend), zleft)
CALL locintv(splrz%sp2, parts%R(i:iend), rleft)
! Compute the value of the splines at the particles positions
CALL basfun(parts%Z(i:iend), splrz%sp1, fun, zleft+1)
CALL basfun(parts%R(i:iend), splrz%sp2, fun2, rleft+1)
DO k=1,nbunch
DO jw=1,(femorder(2)+1)
DO it=1,(femorder(1)+1)
irow=zleft(k)+it
jcol=rleft(k)+jw
mu=irow+(jcol-1)*(nrank(1))
! Add contribution of particle nbunch to rhs grid point mu
rhs(mu)=rhs(mu)+fun(it,0,k)*fun2(jw,0,k)*chargecoeff
END DO
END DO
END DO
END DO
!$OMP END PARALLEL DO
END IF
! If we are using MPI parallelism, reduce the rhs on the root process
IF(mpisize .gt. 1) THEN
IF (mpirank.eq.0) THEN
CALL MPI_REDUCE(MPI_IN_PLACE, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, &
& MPI_SUM, 0, MPI_COMM_WORLD, ierr)
ELSE
CALL MPI_REDUCE(rhs, rhs, nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, &
& MPI_SUM, 0, MPI_COMM_WORLD, ierr)
END IF
END IF
! Apply the Dirichlet boundary conditions on the coaxial insert (if present)
IF(rgrid(0).ne.0.0_db) rhs(1:nrank(1))=potinn
! Apply the Dirichlet boundary conditions on the cylinder walls
rhs((nrank(2)-1)*nrank(1)+1:nrank(2)*nrank(1))=potout
DEALLOCATE(fun,fun2, zleft, rleft)
END SUBROUTINE rhscon
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Solves Poisson equation using FEM. Distributes the result on all MPI workers and interpolate the electric forces
!> for each particle.
!
!---------------------------------------------------------------------------
SUBROUTINE poisson
USE bsplines
USE mumps_bsplines
USE futils
INTEGER:: partend, i, ierr, jder(2)
jder(1)=0
jder(2)=0
partend=parts%Nptot
! Only the root process solves Poisson
IF(mpirank.eq.0) CALL bsolve(femat,rhs,sol)
! Broadcast the solution to all MPI workers
CALL MPI_Bcast(sol(1), nrank(1)*nrank(2), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
matcoef= reshape(sol, (/nrank(1),nrank(2)/))
! Evaluate the electric potential on the grid points
CALL gridval(splrz, vec1, vec2, pot, jder, matcoef)
! Evaluate the electric potential and field at the particles position
CALL getgrad(splrz, parts%Z(1:partend), parts%R(1:partend), parts%pot(1:partend), parts%EZ(1:partend), parts%ER(1:partend))
IF (nlphis) THEN
! We consider self-consistent interactions
Do i=1,partend
parts%EZ(i)=-parts%Ez(i)
parts%ER(i)=-parts%ER(i)
END DO
ELSE
! Self-consistent electric field is neglected
Do i=1,partend
parts%EZ(i)=0
parts%ER(i)=0
END DO
END IF
IF( mpirank .eq. 0 .and. modulo(step,it2d) .eq. 0) THEN
! On the root process, compute the electric field for diagnostic purposes
CALL gridval(splrz, vec1, vec2, Ez, (/1,0/))
CALL gridval(splrz, vec1, vec2, Er, (/0,1/))
Ez=-Ez
Er=-Er
END IF
END SUBROUTINE poisson
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Constucts the FEM matrix using bsplines initialized in fields_init
!---------------------------------------------------------------------------
SUBROUTINE fematrix
USE bsplines
USE mumps_bsplines
DOUBLE PRECISION, ALLOCATABLE :: xgauss(:,:), wgauss(:), zg(:), rg(:), wzg(:), wrg(:)
DOUBLE PRECISION, ALLOCATABLE :: f(:,:), aux(:)
DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
DOUBLE PRECISION, ALLOCATABLE :: fun(:,:), fun2(:,:)
DOUBLE PRECISION :: contrib
INTEGER, ALLOCATABLE :: idert(:), iderw(:)
INTEGER :: i, j, k, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms
kterms=2 !Number of terms in weak form
ALLOCATE(fun(1:femorder(1)+1,0:1),fun2(1:femorder(2)+1,0:1))!Arrays keeping values of b-splines at gauss node
ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays
ALLOCATE(f((femorder(1)+1)*(femorder(2)+1),2),aux(femorder(1)+1)) !Auxiliary arrays ordering bsplines
ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms)) !Pointers on the order of derivatives
! Constuction of auxiliary array ordering bsplines in given interval
DO i=1,(femorder(1)+1)
aux(i)=i
END DO
DO i=1,(femorder(2)+1)
f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),1)=aux
f((i-1)*(femorder(1)+1)+1:i*(femorder(1)+1),2)=i
END DO
! Initialisation of FEM matrix
CALL init(nrank(1)*nrank(2),2,femat)
! Assemble FEM matrix
DO j=1,nr ! Loop on r position
DO i=1,nz ! Loop on z position
! Computation of gauss weight and position in r and z direction for gaussian integration
CALL get_gauss(splrz%sp1, ngauss(1), i, zg, wzg)
CALL get_gauss(splrz%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
DO igauss=1,ngauss(1)*ngauss(2) ! Loop on gaussian weights and positions
CALL basfun(xgauss(igauss,1), splrz%sp1, fun, i)
CALL basfun(xgauss(igauss,2), splrz%sp2, fun2, j)
CALL coefeq(xgauss(igauss,:), idert, iderw, coefs)
DO iterm=1,kterms ! Loop on the two integration dimensions
DO jt=1,(1+femorder(1))*(femorder(2)+1)
DO iw=1,(1+femorder(1))*(femorder(2)+1)
irow=i+f(jt,1)-1; jcol=j+f(jt,2)-1
irow2=i+f(iw,1)-1; jcol2=j+f(iw,2)-1
mu=irow+(jcol-1)*nrank(1)
mu2=irow2+(jcol2-1)*nrank(1)
contrib=fun(f(jt,1),idert(iterm))* fun(f(iw,1),idert(iterm))* &
& fun2(f(jt,2),iderw(iterm))*fun2(f(iw,2),iderw(iterm))* &
& wgauss(igauss)*xgauss(igauss,2)
CALL updtmat(femat, mu, mu2, contrib)
END DO
END DO
END DO
END DO
END DO
END DO
! Impose Dirichlet boundary conditions on FEM matrix
CALL fe_dirichlet
DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg)
DEALLOCATE(f, aux)
DEALLOCATE(idert, iderw, coefs, fun,fun2)
END SUBROUTINE fematrix
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Imposes the dirichlet boundary conditions on the FEM matrix.
!---------------------------------------------------------------------------
SUBROUTINE fe_dirichlet
DOUBLE PRECISION, ALLOCATABLE :: arr(:)
INTEGER :: i
ALLOCATE(arr(nrank(1)*nrank(2)))
DO i=1,nrank(1)*nrank(2)
CALL getrow(femat,i,arr)
END DO
DO i=1,nrank(1)
IF(rgrid(0).ne.0.0_db) THEN
arr=0; arr(i)=1;
CALL putrow(femat,i,arr)
END IF
arr=0; arr(nrank(1)*nrank(2)+1-i)=1 ;
CALL putrow(femat,nrank(1)*nrank(2)+1-i,arr)
END DO
DEALLOCATE(arr)
END SUBROUTINE fe_dirichlet
!________________________________________________________________________________
SUBROUTINE coefeq(x, idt, idw, c)
DOUBLE PRECISION, INTENT(in) :: x(2)
INTEGER, INTENT(out) :: idt(:), idw(:)
DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
c(1) = x(2)
idt(1) = 0
idw(1) = 1
c(2) = x(2)
idt(2) = 1
idw(2) = 0
END SUBROUTINE coefeq
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the magnetic field on the grid according to a magnetic mirror.
!---------------------------------------------------------------------------
SUBROUTINE magnet
USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm
USE constants, ONLY: Pi
DOUBLE PRECISION :: rg, zg,halfLz, MirrorRatio
INTEGER :: i, rindex
halfLz=(zgrid(nz)+zgrid(0))/2
MirrorRatio=(Rcurv-1)/(Rcurv+1)
DO i=1,(nr+1)*(nz+1)
rindex=(i-1)/(nz+1)
rg=rgrid(rindex)
zg=zgrid(i-rindex*(nz+1)-1)-halfLz
Br(i)=-B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm
Bz(i)=B0*(1-MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm
Athet(i)=0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm))
END DO
END SUBROUTINE magnet
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the zero order
FUNCTION bessi0(x)
DOUBLE PRECISION :: bessi0,x
DOUBLE PRECISION :: ax
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0,0.2659732d0,0.360768d-1,0.45813d-2/
DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2 ,0.916281d-2, &
& -0.2057706d-1, 0.2635537d-1,-0.1647633d-1, 0.392377d-2 /
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))
else
ax=abs(x)
y=3.75/ax
bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
endif
return
END FUNCTION bessi0
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the first order
FUNCTION bessi1(x)
DOUBLE PRECISION :: bessi1,x
DOUBLE PRECISION :: ax
DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9,y
SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9
DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0,0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/
DATA q1,q2,q3,q4,q5,q6,q7,q8,q9 /0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2,-0.1031555d-1, &
& 0.2282967d-1, -0.2895312d-1, 0.1787654d-1,-0.420059d-2/
if (abs(x).lt.3.75) then
y=(x/3.75)**2
bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))))
else
ax=abs(x)
y=3.75/ax
bessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*(q7+y*(q8+y*q9))))))))
if(x.lt.0.)bessi1=-bessi1
endif
return
END FUNCTION bessi1
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Free the memory used by the fields module
!---------------------------------------------------------------------------
SUBROUTINE clean_fields
Use bsplines
DEALLOCATE(matcoef)
DEALLOCATE(pot)
DEALLOCATE(rhs)
DEALLOCATE(sol)
DEALLOCATE(Br,Bz)
DEALLOCATE(Er,Ez)
DEALLOCATE(vec1,vec2)
Call DESTROY_SP(splrz)
END SUBROUTINE clean_fields
END MODULE fields
c4science · Help