diff --git a/src/fields_mod.f90 b/src/fields_mod.f90
index fc32d3a..d533d19 100644
--- a/src/fields_mod.f90
+++ b/src/fields_mod.f90
@@ -1,1560 +1,1519 @@
!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: fields
!
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for initializing the magnetic field, solving the Poisson equation and computing the moments of the particles distribution function
!------------------------------------------------------------------------------
MODULE fields
USE constants
USE basic, ONLY: nr, nthet, nz, zgrid, rgrid, Br, Bz, Er, Ethet, Ez, femorder, ngauss, nlppform, pot, Athet, &
& splrz, splrz_ext, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d, timera, potxt, Erxt, Ezxt
USE beam, ONLY: partslist
USE bsplines
USE mumps_bsplines
use mpi
Use omp_lib
Use mpihelper, ONLY: db_type
USE particletypes
USE, INTRINSIC :: iso_c_binding
IMPLICIT NONE
INCLUDE 'fftw3.f03'
- REAL(kind=db), allocatable, SAVE :: matcoef(:,:), phi(:,:), dphidthet(:,:), vec1(:), vec2(:)
+ REAL(kind=db), allocatable, SAVE :: matcoef(:,:), phiext(:), phi(:,:), dphidthet(:,:), vec1(:), vec2(:)
REAL(kind=db), ALLOCATABLE, SAVE :: dftrhsreal(:,:), dftrhsimag(:,:), dftphireal(:,:), dftphiimag(:,:)
REAL(kind=db), allocatable, SAVE :: loc_moments(:,:,:), loc_rhs(:,:), omp_loc_rhs(:,:,:), gradgtilde(:), fverif(:), ppformwork(:,:,:)
INTEGER, SAVE:: loc_zspan
TYPE(mumps_mat), SAVE, ALLOCATABLE :: femat(:) !< Finite Element Method matrix for the full domain
TYPE(mumps_mat), SAVE, ALLOCATABLE :: reduccedmat(:) !< Finite Element Method matrix in the redduced web-spline sub-space
INTEGER :: nbmoments = 4 !< Number of moments to be calculated and stored
INTEGER(kind=omp_lock_kind), Allocatable:: mu_lock(:) !< Stores the lock for fields parallelism
CONTAINS
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Calculate magnetic field mirror components in grid points (Davidson analytical formula employed)
!> or load it from magnetfile if present
!---------------------------------------------------------------------------
SUBROUTINE mag_init
USE basic, ONLY: nr, nz,lu_in,cstep
USE magnet
ALLOCATE (Br((nr + 1)*(nz + 1)), Bz((nr + 1)*(nz + 1)))
ALLOCATE (Athet((nr + 1)*(nz + 1)))
call timera(0, "mag_init")
CALL magnet_init(lu_in,cstep)
call timera(1, "mag_init")
end subroutine mag_init
!---------------------------------------------------------------------------
!> @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
USE bsplines
USE geometry
USE mumps_bsplines
USE mpihelper
INTEGER :: nrz(2), i, d2, k1, n1, m
! Auxiliary vectors
ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1)))
DO i=0,nz
vec1(i*(nr+1)+1:(i+1)*(nr+1))=rgrid
vec2(i*(nr+1)+1:(i+1)*(nr+1))=zgrid(i)
END DO
! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential
CALL set_spline(femorder, ngauss, rgrid, zgrid, splrz_ext, nlppform=nlppform, period=nlperiod)
! Set up 2d spline splrz used in the FEM
ALLOCATE(splrz(nthet))
DO m=1,nthet
CALL set_spline(femorder, ngauss, rgrid, zgrid, splrz(m), nlppform=nlppform, period=nlperiod)
END DO
! Allocate the work buffer to calculate the ppform
d2 = splrz(1)%sp2%dim
k1 = splrz(1)%sp1%order
n1 = splrz(1)%sp1%nints
ALLOCATE(ppformwork(d2,k1,n1))
! Calculate dimension of splines
nrz(1) = nz
nrz(2) = nr
CALL get_dim(splrz(1), nrank, nrz, femorder)
! Allocate necessary variables
ALLOCATE(femat(nthet/2+1))
ALLOCATE(reduccedmat(nthet/2+1))
ALLOCATE(matcoef(nrank(1), nrank(2)))
ALLOCATE(pot((nr + 1)*(nz + 1),nthet))
ALLOCATE(Er((nr + 1)*(nz + 1),nthet))
ALLOCATE(Ethet((nr + 1)*(nz + 1),nthet))
ALLOCATE(Ez((nr + 1)*(nz + 1),nthet))
ALLOCATE(potxt((nr + 1)*(nz + 1)))
ALLOCATE(Erxt((nr + 1)*(nz + 1)))
ALLOCATE(Ezxt((nr + 1)*(nz + 1)))
ALLOCATE(rhs(nrank(1)*nrank(2),nthet))
rhs=0
ALLOCATE(dftrhsreal(nrank(1)*nrank(2),nthet/2+1))
ALLOCATE(dftrhsimag(nrank(1)*nrank(2),nthet/2+1))
ALLOCATE(dftphireal(nrank(1)*nrank(2),nthet/2+1))
ALLOCATE(dftphiimag(nrank(1)*nrank(2),nthet/2+1))
ALLOCATE(gradgtilde(nrank(1)*nrank(2)))
gradgtilde = 0
+ ALLOCATE(phiext(nrank(1)*nrank(2)))
ALLOCATE(phi(nrank(1)*nrank(2),nthet))
ALLOCATE(dphidthet(nrank(1)*nrank(2),nthet))
ALLOCATE (volume(nrank(1)*nrank(2)))
volume = 0
ALLOCATE (mu_lock(nrank(1)*nrank(2)))
DO i = 1, nrank(1)*nrank(2)
CALL omp_init_lock(mu_lock(i))
END DO
END SUBROUTINE fields_init
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Set-up the geometry definition and read it from the standard input
!> Precomputes the LHS matrix to solve Poisson abd the RHS effect of the dirichlet boundaries
!---------------------------------------------------------------------------
SUBROUTINE fields_start
USE geometry
USE basic, ONLY: nrank, nthet
implicit none
INTEGER:: i,j, k, m, ierr
DOUBLE PRECISION:: val
! set up the geometry module for setting up non-conforming boundary conditions
call timera(0, "geom_init")
call geom_init(splrz(1), vec1, vec2)
call timera(1, "geom_init")
CALL timera(0, "compute femat")
DO m=1,nthet/2+1
! Initialisation of FEM matrix
CALL init(nrank(1)*nrank(2), 2, femat(m))
! Calculate and factorise FEM matrix (depends only on mesh)
CALL fematrix(femat(m),m)
IF(nlweb) THEN
! Calculate reduced matrix for use of web splines
CALL reducematrix(femat(m), reduccedmat(m))
CALL factor(reduccedmat(m))
ELSE
CALL factor(femat(m))
END IF
END DO
CALL timera(1, "compute femat")
IF(walltype .lt. 0) THEN
ALLOCATE(fverif(nrank(1)*nrank(2)))
fverif = 0
END IF
! Compute the volume of the splines and gtilde for solving E using web-splines
CALL comp_volume
!$OMP PARALLEL
CALL comp_gradgtilde
!$OMP END PARALLEL
CALL vacuum_field
END SUBROUTINE fields_start
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Computes the externally imposed electric field
!---------------------------------------------------------------------------
SUBROUTINE vacuum_field
USE geometry
USE basic, ONLY: pot, rhs
implicit none
INTEGER:: i,m,iend
!$OMP PARALLEL private(i)
!$OMP DO SIMD
DO i=1,size(rhs,1)
rhs(i,1)=-gradgtilde(i)
END DO
!$OMP END DO SIMD
!$OMP END PARALLEL
CALL poisson2D
!$OMP PARALLEL private(i,iend)
CALL vacuum_field_comp(splrz_ext)
!$OMP BARRIER
!$OMP END PARALLEL
END SUBROUTINE vacuum_field
!---------------------------------------------------------------------------
!> @author
!> Pierrick Giroud EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Updates the splinevar variable with the new phi coefficients and calculates
!> Phi Er and Ez on the grid
!---------------------------------------------------------------------------
SUBROUTINE vacuum_field_comp(splinevar)
USE basic, ONLY: nrank, pot, nlend, gradthet, rgrid, nr, nz
USE bsplines, ONLY: spline2d, gridval
USE mumps_bsplines, ONLY: bsolve, vmx
USE futils
USE geometry
type(spline2d) :: splinevar
INTEGER:: ierr, i, j, iend, rindex
INTEGER:: Evalpot(2), EvalEr(2), EvalEz(2)
Evalpot=(/0,0/)
EvalEr=(/1,0/)
EvalEz=(/0,1/)
!$OMP DO SIMD collapse(2)
DO j=1,nrank(2)
DO i=1,nrank(1)
- matcoef(i,j) = phi((j-1)*nrank(1)+i,1)
+ matcoef(i,j) = phiext((j-1)*nrank(1)+i)
END DO
END DO
!$OMP END DO SIMD
! Update the ppform coefficients
CALL updt_ppform2d(splinevar, matcoef)
!$OMP BARRIER
! On the root process, compute the electric field for diagnostic purposes
IF (mpirank .eq. 0) THEN
!$OMP DO
DO i=1,size(potxt),16
iend=min(size(pot,1),i+15)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), potxt(i:iend), Evalpot)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Erxt(i:iend), EvalEr)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Ezxt(i:iend), EvalEz)
Ezxt(i:iend) = -potxt(i:iend)*gridwdir(1,i:iend) - Ezxt(i:iend)*gridwdir(0,i:iend) - gtilde(1,i:iend)
Erxt(i:iend) = -potxt(i:iend)*gridwdir(2,i:iend) - Erxt(i:iend)*gridwdir(0,i:iend) - gtilde(2,i:iend)
potxt(i:iend) = potxt(i:iend)*gridwdir(0,i:iend) + gtilde(0,i:iend)
END DO
!$OMP END DO NOWAIT
END IF
END SUBROUTINE vacuum_field_comp
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Set-up the necessary variables for the communication of moments and rhs grid
!---------------------------------------------------------------------------
SUBROUTINE fields_comm_init(Zbounds)
USE basic, ONLY: nrank,nthet
USE mpihelper
USE omp_lib
INTEGER:: Zbounds(0:), nbthreads
nbthreads = omp_get_max_threads()
loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(2)
if (allocated(loc_moments)) deallocate (loc_moments)
ALLOCATE (loc_moments(nbmoments,loc_zspan*nrank(1),nthet))
if (allocated(loc_rhs)) deallocate (loc_rhs)
ALLOCATE (loc_rhs(loc_zspan*nrank(1),nthet))
if (allocated(omp_loc_rhs)) deallocate (omp_loc_rhs)
ALLOCATE (omp_loc_rhs(nbthreads,loc_zspan*nrank(1),nthet))
IF (mpisize .gt. 1) THEN
CALL init_overlaps(nrank, femorder, Zbounds(mpirank), Zbounds(mpirank + 1), nbmoments)
END IF
END SUBROUTINE fields_comm_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
!
!> @param[in] plist list of the particles type storing the desired specie parameters
!
!---------------------------------------------------------------------------
SUBROUTINE rhscon(plist)
USE bsplines
use mpi
USE basic!, ONLY: rhs, Zbounds
USE beam, ONLY: particles
USE mpihelper
Use geometry
Use omp_lib
type(particles), INTENT(INOUT):: plist(:)
INTEGER:: i,j,k,m
IF(nlphis) THEN
!$OMP DO SIMD
DO i=1,size(loc_rhs,1)
DO m=1,size(loc_rhs,2)
loc_rhs(i,m)=0
END DO
END DO
!$OMP END DO SIMD
! Assemble rhs for each specie
DO i = 1, size(plist, 1)
IF(plist(i)%is_field) CALL deposit_charge(plist(i))
END DO
! Communicate the overlaps
IF(mpisize .gt. 1) CALL rhs_overlap
! ! Add gradgtilde
! !$OMP DO SIMD
! DO i=1,size(loc_rhs,1)
! j=(i-1)/nrank(1)+Zbounds(mpirank)
! k=mod(i-1,nrank(1))+1
! DO m=1,size(loc_rhs,2)
! loc_rhs(i,m)=loc_rhs(i,m)-gradgtilde(k+j*nrank(1))
! END DO
! END DO
! !$OMP END DO SIMD
! If we are using MPI parallelism, reduce the rhs on the root process
IF(mpisize .gt. 1) THEN
CALL rhs_gather(rhs)
ELSE ! in serial copy loc_rhs to rhs
!$OMP DO
DO i=1,size(loc_rhs,1)
DO m=1,size(loc_rhs,2)
rhs(i,m)=loc_rhs(i,m)
END DO
END DO
!$OMP END DO
END IF
ELSE ! We only consider the externally imposed field
!$OMP DO
DO i=1,size(rhs,1)
DO m=1,size(rhs,2)
rhs(i,m)=0 !-gradgtilde(i)
END DO
END DO
!$OMP END DO
END IF
END SUBROUTINE rhscon
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Deposit the particles charges (q) from p on the grid
!
!> @param[in] p the particles type storing the desired specie parameters
!---------------------------------------------------------------------------
SUBROUTINE deposit_charge(p)
USE bsplines
USE mpi
USE constants
USE basic, ONLY: Zbounds, rnorm, phinorm, thetgrid, invdthet
USE beam, ONLY: particles
USE mpihelper
USE geometry
USE omp_lib
TYPE(particles), INTENT(IN) :: p
INTEGER :: irow, jcol, it, jw, mu, i, k, iend, nbunch
INTEGER, DIMENSION(:), ALLOCATABLE :: zindex, rindex, thetindex
REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :)
INTEGER:: num_threads, curr_thread
real(kind=db):: contrib, chargecoeff, wthet
num_threads = omp_get_max_threads()
nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun
nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun
nbunch = min(nbunch, 16) ! Particle bunch size used when calling basfun
chargecoeff = p%weight*p%q/(2*pi*eps_0*phinorm*rnorm) ! Normalized charge density simulated by each macro particle
! Assemble rhs
IF (p%Nploc .gt. 0) THEN
ALLOCATE (zindex(nbunch), rindex(nbunch), thetindex(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
zindex=0
rindex=0
thetindex=0
curr_thread=omp_get_thread_num()+1
!$OMP DO SIMD
DO i=1,size(loc_rhs,1)
omp_loc_rhs(:,i,:)=0
END DO
!$OMP END DO SIMD
!$OMP DO
DO i = 1, p%Nploc, nbunch
! Avoid segmentation fault by accessing non relevant data
iend = min(i + nbunch - 1, p%Nploc)
k = iend - i + 1
! Localize the particle
rindex(1:k) = p%cellindex(1,i:iend)
thetindex(1:k) = p%cellindex(2,i:iend)
zindex(1:k) = p%cellindex(3,i:iend)
! Compute the value of the splines at the particles positions
CALL basfun(p%pos(1,i:iend), splrz(1)%sp1, fun, rindex(1:k) + 1)
CALL basfun(p%pos(3,i:iend), splrz(1)%sp2, fun2, zindex(1:k) + 1)
DO k = 1, (iend - i + 1)
DO jw = 1, (femorder(2) + 1)
jcol = zindex(k) + jw - Zbounds(mpirank)
DO it = 1, (femorder(1) + 1)
irow = rindex(k) + it
mu = irow + (jcol - 1)*(nrank(1))
! Add contribution of particle k to rhs grid point mu
contrib = fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(0,i + k - 1)*chargecoeff
wthet = (p%pos(2,i+k-1)-thetgrid(thetindex(k)))*invdthet
!omp_loc_rhs(curr_thread,mu) = omp_loc_rhs(curr_thread,mu) + contrib
IF(thetindex(k) .eq. nthet-1) THEN
omp_loc_rhs(curr_thread,mu,thetindex(k)+1) = omp_loc_rhs(curr_thread,mu,thetindex(k)+1) + contrib*(1-wthet)
omp_loc_rhs(curr_thread,mu,1) = omp_loc_rhs(curr_thread,mu,1) + contrib*wthet
ELSE
omp_loc_rhs(curr_thread,mu,thetindex(k)+1) = omp_loc_rhs(curr_thread,mu,thetindex(k)+1) + contrib*(1-wthet)
omp_loc_rhs(curr_thread,mu,thetindex(k)+2) = omp_loc_rhs(curr_thread,mu,thetindex(k)+2) + contrib*wthet
END IF
END DO
END DO
END DO
END DO
!$OMP END DO
!$OMP DO SIMD
DO i=1,size(loc_rhs,1)
DO k=1,num_threads
loc_rhs(i,:)=loc_rhs(i,:)+omp_loc_rhs(k,i,:)
END DO
END DO
!$OMP END DO SIMD
DEALLOCATE (fun, fun2, zindex, rindex, thetindex)
END IF
END SUBROUTINE deposit_charge
+!---------------------------------------------------------------------------
+!> @author
+!> Guillaume Le Bars EPFL/SPC
+!
+! DESCRIPTION:
+!> @brief
+!> Do the communication of the local moment matrices between mpi workers
+!> for the overlap grid points
+!---------------------------------------------------------------------------
+SUBROUTINE rhs_overlap
+ USE mpihelper
+ USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc, nthet
+ INTEGER:: ierr, i, m
+
+ !$OMP MASTER
+ DO m=1,nthet
+ CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs(:,m), nrank, femorder, loc_zspan - femorder(2))
+ IF (mpirank .gt. 0) THEN
+ DO i = 1,femorder(2)*nrank(1)
+ loc_rhs(i,m) = loc_rhs(i,m) + rhsoverlap_buffer(i)
+ END DO
+ END IF
+ END DO
+ !$OMP END MASTER
+ !$OMP BARRIER
+END SUBROUTINE rhs_overlap
+
+
+!---------------------------------------------------------------------------
+!> @author
+!> Guillaume Le Bars EPFL/SPC
+!
+! DESCRIPTION:
+!> @brief
+!> Do the communication of the local moment matrices between mpi workers
+!> to reduce the result on the host
+!---------------------------------------------------------------------------
+SUBROUTINE rhs_gather(rhs)
+ USE mpihelper
+ USE basic, ONLY: Zbounds, mpirank, leftproc, rightproc, nthet
+ REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: rhs
+ INTEGER:: ierr, m
+ INTEGER:: displs(mpisize), counts(mpisize)
+ INTEGER:: overlap_type
+ INTEGER:: rcvoverlap_type
+
+ displs = Zbounds(0:mpisize - 1)*nrank(1)
+ counts = (Zbounds(1:mpisize) - Zbounds(0:mpisize - 1))*nrank(1)
+ counts(mpisize) = counts(mpisize) + femorder(1)*nrank(1)
+
+ !$OMP MASTER
+ IF (mpirank .eq. 0) THEN
+ rhs = 0
+ END IF
+ DO m=1,nthet
+ CALL MPI_GATHERV(loc_rhs(:,m), counts(mpirank + 1), db_type, rhs(:,m), counts, displs, db_type, 0, MPI_COMM_WORLD, ierr)
+ END DO
+ !$OMP END MASTER
+ !$OMP BARRIER
+END SUBROUTINE rhs_gather
+
+
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Calculate the 0th 1st and 2nd order moments of the particle p and stores it in moment
!
!> @param[in] p the particles type storing the desired specie parameters
!> @param[out] moment the 2d array storing the calculated moments
!---------------------------------------------------------------------------
SUBROUTINE momentsdiag(p)
USE bsplines
USE mpi
USE beam, ONLY: particles
USE basic
USE mpihelper
USE geometry
type(particles), INTENT(INOUT):: p
INTEGER :: i,istart,iend
!$OMP SINGLE
loc_moments=0 ! Reset the moments matrix
!$OMP END SINGLE
CALL deposit_moments(p)
!$OMP SINGLE
IF(.not. allocated(p%moments)) THEN
IF(mpirank.eq.0) THEN
ALLOCATE(p%moments(4,nrank(1)*nrank(2),nthet))
ELSE
ALLOCATE(p%moments(0,0,0))
END IF
END IF
!$OMP END SINGLE
! If we are using MPI parallelism, reduce the moments on the root process
IF (mpisize .gt. 1) THEN
CALL moments_gather(p%moments)
ELSE
!$OMP SINGLE
p%moments = loc_moments
!$OMP END SINGLE NOWAIT
END IF
END SUBROUTINE momentsdiag
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Deposit the particles moments (n,v,v^2) from p on the grid
!
!> @param[in] p the particles type storing the desired specie parameters
!---------------------------------------------------------------------------
SUBROUTINE deposit_moments(p)
USE bsplines
USE mpi
USE basic, ONLY: Zbounds, thetgrid, invdthet
USE beam, ONLY: particles
USE mpihelper
USE geometry
USE omp_lib
TYPE(particles), INTENT(IN):: p
INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch
INTEGER, DIMENSION(:), ALLOCATABLE :: zindex, rindex, thetindex
REAL(kind=db) :: vr, vthet, vz, coeff, wthet
REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :)
INTEGER:: num_threads
num_threads = omp_get_max_threads()
nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun
nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun
nbunch = min(nbunch, 64) ! Particle bunch size used when calling basfun
! Assemble moments
IF (p%Nploc .gt. 0) THEN
ALLOCATE(zindex(nbunch), rindex(nbunch), thetindex(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
!$OMP DO
DO i = 1, p%Nploc, nbunch
! Avoid segmentation fault by accessing non relevant data
iend = min(i + nbunch - 1, p%Nploc)
k = iend - i + 1
! Localize the particle
rindex(1:k) = p%cellindex(1,i:iend)
thetindex(1:k) = p%cellindex(2,i:iend)
zindex(1:k) = p%cellindex(3,i:iend)
! Compute the value of the splines at the particles positions
CALL basfun(p%pos(1,i:iend), splrz(1)%sp1, fun(:, :, 1:k), rindex(1:k) + 1)
CALL basfun(p%pos(3,i:iend), splrz(1)%sp2, fun2(:, :, 1:k), zindex(1:k) + 1)
DO k = 1, (iend - i + 1)
DO jw = 1, (femorder(2) + 1)
jcol = zindex(k) + jw - Zbounds(mpirank)
DO it = 1, (femorder(1) + 1)
irow = rindex(k) + it
mu = irow + (jcol - 1)*(nrank(1))
coeff = p%weight*fun(it, 0, k)*fun2(jw, 0, k)
wthet = (p%pos(2,i+k-1)-thetgrid(thetindex(k)))*invdthet
! Add contribution of particle nbunch to moments grid point mu
vr = 0.5*(p%U(1,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(1,i + k - 1)/p%Gammaold(i + k - 1))
vz = 0.5*(p%U(3,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(3,i + k - 1)/p%Gammaold(i + k - 1))
vthet = 0.5*(p%U(2,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(2,i + k - 1)/p%Gammaold(i + k - 1))
CALL omp_set_lock(mu_lock(mu))
!loc_moments(1:4,mu)=loc_moments(1:4,mu)+coeff*(/1.0_db,vr,vthet,vz/)
IF(thetindex(k) .eq. nthet-1) THEN
loc_moments(1:4,mu,thetindex(k)+1)=loc_moments(1:4,mu,thetindex(k)+1) + (1-wthet)*coeff*(/1.0_db,vr,vthet,vz/)
loc_moments(1:4,mu,1)=loc_moments(1:4,mu,1) + wthet*coeff*(/1.0_db,vr,vthet,vz/)
ELSE
loc_moments(1:4,mu,thetindex(k)+1)=loc_moments(1:4,mu,thetindex(k)+1) + (1-wthet)*coeff*(/1.0_db,vr,vthet,vz/)
loc_moments(1:4,mu,thetindex(k)+2)=loc_moments(1:4,mu,thetindex(k)+2) + wthet*coeff*(/1.0_db,vr,vthet,vz/)
END IF
CALL omp_unset_lock(mu_lock(mu))
END DO
END DO
END DO
END DO
!$OMP END DO NOWAIT
DEALLOCATE (fun, fun2, zindex, rindex, thetindex)
END IF
END subroutine deposit_moments
-!---------------------------------------------------------------------------
-!> @author
-!> Guillaume Le Bars EPFL/SPC
-!
-! DESCRIPTION:
-!> @brief
-!> Do the communication of the local moment matrices between mpi workers
-!> for the overlap grid points
-!---------------------------------------------------------------------------
-SUBROUTINE rhs_overlap
- USE mpihelper
- USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc, nthet
- INTEGER:: ierr, i, m
-
- DO m=1,nthet
- !$OMP MASTER
- CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs(:,m), nrank, femorder, loc_zspan - femorder(2))
- !$OMP END MASTER
- IF (mpirank .gt. 0) THEN
- !$OMP DO SIMD
- DO i = 1,femorder(2)*nrank(1)
- loc_rhs(i,m) = loc_rhs(i,m) + rhsoverlap_buffer(i)
- END DO
- !$OMP END DO SIMD
- END IF
- !$OMP BARRIER
- END DO
-END SUBROUTINE rhs_overlap
-
-
-!---------------------------------------------------------------------------
-!> @author
-!> Guillaume Le Bars EPFL/SPC
-!
-! DESCRIPTION:
-!> @brief
-!> Do the communication of the local moment matrices between mpi workers
-!> to reduce the result on the host
-!---------------------------------------------------------------------------
-SUBROUTINE rhs_gather(rhs)
- USE mpihelper
- USE basic, ONLY: Zbounds, mpirank, leftproc, rightproc, nthet
- REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: rhs
- INTEGER:: ierr, m
- INTEGER:: displs(mpisize), counts(mpisize)
- INTEGER:: overlap_type
- INTEGER:: rcvoverlap_type
-
- displs = Zbounds(0:mpisize - 1)*nrank(1)
- counts = (Zbounds(1:mpisize) - Zbounds(0:mpisize - 1))*nrank(1)
- counts(mpisize) = counts(mpisize) + femorder(1)*nrank(1)
-
- ! Set communication vector type
- overlap_type = rhsoverlap_type
- rcvoverlap_type = rcvrhsoverlap_type
-
- !$OMP MASTER
- IF (mpirank .eq. 0) THEN
- rhs = 0
- END IF
- DO m=1,nthet
- CALL MPI_GATHERV(loc_rhs(:,m), counts(mpirank + 1), db_type, rhs(:,m), counts, displs, db_type, 0, MPI_COMM_WORLD, ierr)
- END DO
- !$OMP END MASTER
- !$OMP BARRIER
-END SUBROUTINE rhs_gather
-
-
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Do the communication of the local moment matrices between mpi workers
!> for the overlap grid points and reduce the result on the host
!---------------------------------------------------------------------------
SUBROUTINE moments_gather(moment)
USE mpihelper
USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc, nthet
REAL(kind=db), DIMENSION(:,:,:), INTENT(INOUT):: moment
INTEGER:: ierr, i, j, m
INTEGER:: displs(mpisize), counts(mpisize)
displs = Zbounds(0:mpisize - 1)*nrank(1)*4
counts = (Zbounds(1:mpisize) - Zbounds(0:mpisize - 1))*4*nrank(1)
counts(mpisize) = counts(mpisize) + femorder(2)*4*nrank(1)
+ !$OMP MASTER
! Communicate overlaps
DO m=1,nthet
- !$OMP MASTER
CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments(:,:,m), nrank, femorder, loc_zspan - femorder(2))
- !$OMP END MASTER
- !$OMP BARRIER
-
IF (mpirank .gt. 0) THEN
- !!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i)
- !$OMP DO SIMD
DO i = 1, nrank(1)*femorder(2)
loc_moments(1:nbmoments,i,m) = loc_moments(1:nbmoments,i,m) + momentsoverlap_buffer((i-1)*nbmoments+1:i*nbmoments)
END DO
- !$OMP END DO SIMD
END IF
- !$OMP BARRIER
END DO
-
- !$OMP MASTER
- ! Set communication vector type
IF (mpirank .eq. 0) THEN
moment = 0
END IF
DO m=1,nthet
CALL MPI_GATHERV(loc_moments(:,:,m), counts(mpirank + 1), db_type, moment(:,:,m), counts, displs, db_type, 0, MPI_COMM_WORLD, ierr)
END DO
!$OMP END MASTER
!$OMP BARRIER
END SUBROUTINE moments_gather
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Solves Poisson equation using FEM. Distributes the result on all MPI workers.
!---------------------------------------------------------------------------
SUBROUTINE poisson
USE basic, ONLY: rhs, nrank, nlend, nthet, rgrid, nr, nz
USE bsplines, ONLY: spline2d, gridval
USE mumps_bsplines, ONLY: bsolve, vmx
USE futils
USE geometry
INTEGER:: ierr, i, j, m, iend
real(kind=db), allocatable::reducedrhsreal(:,:),reducedrhsimag(:,:)
real(kind=db), allocatable::reducedsolreal(:,:), reducedsolimag(:,:), tempcol(:)
! FFT specific variables
TYPE(C_PTR) :: plan_theta2m ! transform from (r,theta,z) to (r,m,z)
TYPE(C_PTR) :: plan_m2theta ! transform from (r,m,z) to (r,theta,z)
REAL(kind=db), DIMENSION(:,:), ALLOCATABLE :: rhstemp, phitemp
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE :: dftrhstemp, dftphitemp
! double complex, allocatable:: rhstemp(:), soltemp(:)
! TYPE(zmumps_mat) :: mattemp
! ALLOCATE(rhstemp(nbreducedspline),soltemp(nbreducedspline))
! DO i=1,nbreducedspline
! rhstemp(i)=cmplx(reducedrhsreal(i),reducedrhsimag(i))
! END DO
!mattemp=0
! rhstemp=0
! soltemp=0
! CALL bsolve(mattemp, rhstemp, soltemp)
- ! Compute DFT of RHS along azimuthal direction
ALLOCATE(rhstemp(nthet,nrank(1)*nrank(2)))
ALLOCATE(dftrhstemp(nthet/2+1,nrank(1)*nrank(2)))
-
- rhstemp=transpose(rhs)
- CALL dfftw_plan_many_dft_r2c(plan_theta2m,1,nthet,nrank(1)*nrank(2),rhstemp,nthet,1,nthet,dftrhstemp,nthet/2+1,1,nthet/2+1,FFTW_ESTIMATE)
- CALL dfftw_execute_dft_r2c(plan_theta2m,rhstemp,dftrhstemp)
- CALL dfftw_destroy_plan(plan_theta2m)
- dftrhsreal=real(transpose(dftrhstemp))
- dftrhsimag=aimag(transpose(dftrhstemp))
-
- DEALLOCATE(rhstemp)
- DEALLOCATE(dftrhstemp)
-
- ! Solves Poisson equation
ALLOCATE(reducedrhsreal(nrank(1)*nrank(2),nthet/2+1),reducedrhsimag(nrank(1)*nrank(2),nthet/2+1))
ALLOCATE(reducedsolreal(nbreducedspline,nthet/2+1),reducedsolimag(nbreducedspline,nthet/2+1))
ALLOCATE(tempcol(nrank(1)*nrank(2)))
+ ALLOCATE(dftphitemp(nthet/2+1,nrank(1)*nrank(2)))
+ ALLOCATE(phitemp(nthet,nrank(1)*nrank(2)))
- DO m=1,nthet/2+1
- IF(mpirank .eq. 0) THEN
+ IF(mpirank .eq. 0) THEN
+ ! Compute DFT of RHS along azimuthal direction
+ rhstemp=transpose(rhs)
+ CALL dfftw_plan_many_dft_r2c(plan_theta2m,1,nthet,nrank(1)*nrank(2),rhstemp,nthet,1,nthet,dftrhstemp,nthet/2+1,1,nthet/2+1,FFTW_ESTIMATE)
+ CALL dfftw_execute_dft_r2c(plan_theta2m,rhstemp,dftrhstemp)
+ CALL dfftw_destroy_plan(plan_theta2m)
+ dftrhsreal=real(transpose(dftrhstemp))
+ dftrhsimag=aimag(transpose(dftrhstemp))
+
+ ! Solves Poisson equation
+ DO m=1,nthet/2+1
IF(nlweb) THEN ! We use the web-spline reduction for stability
reducedrhsreal(:,m) = vmx(etilde, dftrhsreal(:,m))
reducedrhsimag(:,m) = vmx(etilde, dftrhsimag(:,m))
CALL bsolve(reduccedmat(m), reducedrhsreal(1:nbreducedspline,m), reducedsolreal(:,m))
CALL bsolve(reduccedmat(m), reducedrhsimag(1:nbreducedspline,m), reducedsolimag(:,m))
+ tempcol = 0
+ tempcol(1:nbreducedspline) = reducedsolreal(:,m)
+ dftphireal(:,m) = vmx(etildet, tempcol)
+ tempcol(1:nbreducedspline) = reducedsolimag(:,m)
+ dftphiimag(:,m) = vmx(etildet, tempcol)
ELSE
CALL bsolve(femat(m), dftrhsreal(:,m), dftphireal(:,m))
CALL bsolve(femat(m), dftrhsimag(:,m), dftphiimag(:,m))
END IF
- END IF
-
- !$OMP MASTER
- ! Distributes the result on all MPI workers
- IF(nlweb) THEN
- CALL MPI_Bcast(reducedsolreal(:,m), nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr)
- CALL MPI_Bcast(reducedsolimag(:,m), nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr)
- ELSE
- CALL MPI_Bcast(dftphireal(:,m), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
- CALL MPI_Bcast(dftphiimag(:,m), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
- END IF
- IF(nlweb) THEN
- tempcol = 0
- tempcol(1:nbreducedspline) = reducedsolreal(:,m)
- dftphireal(:,m) = vmx(etildet, tempcol)
- tempcol(1:nbreducedspline) = reducedsolimag(:,m)
- dftphiimag(:,m) = vmx(etildet, tempcol)
- END IF
- ! DEALLOCATE(reducedsolreal,reducedsolimag)
- !$OMP END MASTER
- !$OMP BARRIER
- END DO
-
- ! Compute the inverse DFT of the potential
- ALLOCATE(dftphitemp(nthet/2+1,nrank(1)*nrank(2)))
- ALLOCATE(phitemp(nthet,nrank(1)*nrank(2)))
-
- DO m=1,nthet/2+1
- DO i=1,nrank(1)*nrank(2)
- dftphitemp(m,i)=cmplx(dftphireal(i,m),dftphiimag(i,m))
END DO
- END DO
- CALL dfftw_plan_many_dft_c2r(plan_m2theta,1,nthet,nrank(1)*nrank(2),dftphitemp,nthet/2+1,1,nthet/2+1,phitemp,nthet,1,nthet,FFTW_ESTIMATE)
- CALL dfftw_execute_dft_c2r(plan_m2theta,dftphitemp,phitemp)
- CALL dfftw_destroy_plan(plan_m2theta)
- phi=transpose(phitemp)
- ! Compute the inverse DFT of azimuthal derivative of potential
- DO m=1,nthet/2+1
- DO i=1,nrank(1)*nrank(2)
- dftphitemp(m,i)=(m-1)*cmplx(0.0,1.0)*dftphitemp(m,i)
+ ! Compute the inverse DFT of the potential
+ DO m=1,nthet/2+1
+ DO i=1,nrank(1)*nrank(2)
+ dftphitemp(m,i)=cmplx(dftphireal(i,m),dftphiimag(i,m))
+ END DO
END DO
- END DO
- CALL dfftw_plan_many_dft_c2r(plan_m2theta,1,nthet,nrank(1)*nrank(2),dftphitemp,nthet/2+1,1,nthet/2+1,phitemp,nthet,1,nthet,FFTW_ESTIMATE)
- CALL dfftw_execute_dft_c2r(plan_m2theta,dftphitemp,phitemp)
- CALL dfftw_destroy_plan(plan_m2theta)
- dphidthet=transpose(phitemp)
+ CALL dfftw_plan_many_dft_c2r(plan_m2theta,1,nthet,nrank(1)*nrank(2),dftphitemp,nthet/2+1,1,nthet/2+1,phitemp,nthet,1,nthet,FFTW_ESTIMATE)
+ CALL dfftw_execute_dft_c2r(plan_m2theta,dftphitemp,phitemp)
+ CALL dfftw_destroy_plan(plan_m2theta)
+ phi=transpose(phitemp)
+
+ ! Compute the inverse DFT of azimuthal derivative of potential
+ DO m=1,nthet/2+1
+ DO i=1,nrank(1)*nrank(2)
+ dftphitemp(m,i)=(m-1)*cmplx(0.0,1.0)*dftphitemp(m,i)
+ END DO
+ END DO
+ CALL dfftw_plan_many_dft_c2r(plan_m2theta,1,nthet,nrank(1)*nrank(2),dftphitemp,nthet/2+1,1,nthet/2+1,phitemp,nthet,1,nthet,FFTW_ESTIMATE)
+ CALL dfftw_execute_dft_c2r(plan_m2theta,dftphitemp,phitemp)
+ CALL dfftw_destroy_plan(plan_m2theta)
+ dphidthet=transpose(phitemp)
+ END IF
- DEALLOCATE(phitemp)
- DEALLOCATE(dftphitemp)
+ ! Distributes the result on all MPI workers
+ DO m=1,nthet
+ CALL MPI_Bcast(phi(:,m), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
+ CALL MPI_Bcast(dphidthet(:,m), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
+ END DO
END SUBROUTINE poisson
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Solves Poisson equation using FEM. Distributes the result on all MPI workers.
!---------------------------------------------------------------------------
SUBROUTINE poisson2D
USE basic, ONLY: rhs, nrank, nlend
USE bsplines, ONLY: spline2d, gridval
USE mumps_bsplines, ONLY: bsolve, vmx
USE futils
USE geometry
INTEGER:: ierr, i, j, iend
real(kind=db), allocatable::reducedrhs(:), reducedsol(:), tempcol(:)
ALLOCATE(reducedrhs(nrank(1)*nrank(2)))
ALLOCATE(reducedsol(nbreducedspline))
ALLOCATE(tempcol(nrank(1)*nrank(2)))
! Solves Poisson equation
IF(nlweb .and. mpirank .eq. 0) THEN ! We use the web-spline reduction for stability
reducedrhs = vmx(etilde, rhs(:,1))
CALL bsolve(reduccedmat(1), reducedrhs(1:nbreducedspline), reducedsol)
+ tempcol = 0
+ tempcol(1:nbreducedspline) = reducedsol
+ phiext = vmx(etildet, tempcol)
ELSE IF(mpirank.eq.0) THEN
- CALL bsolve(femat(1), rhs(:,1), phi(:,1))
+ CALL bsolve(femat(1), rhs(:,1), phiext)
END IF
+ DEALLOCATE(reducedsol,reducedrhs,tempcol)
! Distributes the result on all MPI workers
- !$OMP MASTER
- IF(nlweb) THEN
- CALL MPI_Bcast(reducedsol, nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr)
- ELSE
- CALL MPI_Bcast(phi(:,1), nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
- END IF
- IF(nlweb) THEN
- tempcol = 0
- tempcol(1:nbreducedspline) = reducedsol
- phi(:,1) = vmx(etildet, tempcol)
- END IF
- DEALLOCATE(reducedsol)
- !$OMP END MASTER
- !$OMP BARRIER
+ CALL MPI_Bcast(phiext, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
END SUBROUTINE poisson2D
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Computes the electric fields and potential at the particles position
!> and on the grid for diagnostics purposes
!---------------------------------------------------------------------------
SUBROUTINE Efield_comp(splinevar,p)
USE basic, ONLY: thetgrid, invdthet, nthet, nlend
USE geometry
USE splinebound
USE bsplines, ONLY: spline2d, gridval
USE beam, ONLY: particles
TYPE(spline2d) :: splinevar(:)
TYPE(particles), INTENT(INOUT) :: p
INTEGER :: i,j,k,m,iend,nbunch
INTEGER :: thetindex
REAL(kind=db) :: wthet, pot1, pot2, Er1, Er2, Ez1, Ez2, dphidthet1, dphidthet2, potext, Erext, Ezext
REAL(kind=db), ALLOCATABLE :: gtildeloc(:,:) !,potext(:), Erext(:), Ezext(:)
INTEGER:: Evalpot(2), EvalEr(2), EvalEz(2)
Evalpot=(/0,0/)
EvalEr=(/1,0/)
EvalEz=(/0,1/)
nbunch = 64 ! Particle bunch size
ALLOCATE(gtildeloc(0:2,0:nbunch - 1))
! ALLOCATE(potext(nbunch), Erext(nbunch), Ezext(nbunch))
! Update the splines with phi
DO m=1,nthet
!$OMP DO SIMD collapse(2)
DO j=1,nrank(2)
DO i=1,nrank(1)
- matcoef(i,j) = phi((j-1)*nrank(1)+i,m)
+ matcoef(i,j) = phi((j-1)*nrank(1)+i,m) + phiext((j-1)*nrank(1)+i)
END DO
END DO
!$OMP END DO SIMD
! Update the ppform coefficients
CALL updt_ppform2d(splinevar(m), matcoef)
!$OMP BARRIER
END DO
! Evaluate pot, Er, Ez at the particle position
!$OMP DO SIMD
DO i=1,p%Nploc,nbunch
! Avoid segmentation fault by accessing non relevant data
iend=min(i+nbunch-1,p%Nploc)
- CALL speval(splrz_ext, p%pos(1,i:iend), p%pos(3,i:iend),p%cellindex(1,i:iend),p%cellindex(3,i:iend), p%pot(i:iend),p%E(1,i:iend),p%E(3,i:iend)) ! external field
+ ! CALL speval(splrz_ext, p%pos(1,i:iend), p%pos(3,i:iend),p%cellindex(1,i:iend),p%cellindex(3,i:iend), p%pot(i:iend),p%E(1,i:iend),p%E(3,i:iend)) ! external field
DO k=1,(iend-i+1)
thetindex=p%cellindex(2,i+k-1)
wthet = (p%pos(2,i+k-1)-thetgrid(thetindex))*invdthet
IF(thetindex .eq. nthet-1) THEN
CALL speval1(splrz(thetindex+1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), pot1, Er1, Ez1)
CALL speval1(splrz(1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), pot2, Er2, Ez2)
ELSE
CALL speval1(splrz(thetindex+1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), pot1, Er1, Ez1)
CALL speval1(splrz(thetindex+2), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), pot2, Er2, Ez2)
END IF
!CALL speval1(splrz_ext, p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), potext, Erext, Ezext)
- p%pot(i+k-1)=p%pot(i+k-1)+(1-wthet)*pot1+wthet*pot2 !+potext
- p%E(1,i+k-1)=p%E(1,i+k-1)+(1-wthet)*Er1+wthet*Er2 !+Erext
- p%E(3,i+k-1)=p%E(3,i+k-1)+(1-wthet)*Ez1+wthet*Ez2 !+Ezext
+ p%pot(i+k-1)=(1-wthet)*pot1+wthet*pot2 !+potext+p%pot(i+k-1)
+ p%E(1,i+k-1)=(1-wthet)*Er1+wthet*Er2 !+Erext+p%E(1,i+k-1)
+ p%E(3,i+k-1)=(1-wthet)*Ez1+wthet*Ez2 !+Ezext+p%E(3,i+k-1)
END DO
CALL total_gtilde(p%pos(1,i:iend), p%pos(3,i:iend), gtildeloc(:,0:iend - i), p%geomweight(:,i:iend))
p%E(3,i:iend) = -p%E(3,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(1,i:iend) - gtildeloc(1,0:iend - i)
p%E(1,i:iend) = -p%E(1,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(2,i:iend) - gtildeloc(2,0:iend - i)
p%pot(i:iend) = p%geomweight(0,i:iend)*p%pot(i:iend) + gtildeloc(0,0:iend - i)
END DO
!$OMP END DO SIMD NOWAIT
! On the root process, evaluate pot, Er, Ez on the grid for diagnostic purposes
IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN
DO m=1,nthet
!$OMP DO
DO i=1,size(pot,1),16
iend=min(size(pot,1),i+15)
CALL gridval(splinevar(m), vec1(i:iend), vec2(i:iend), pot(i:iend,m), Evalpot)
CALL gridval(splinevar(m), vec1(i:iend), vec2(i:iend), Er(i:iend,m), EvalEr)
CALL gridval(splinevar(m), vec1(i:iend), vec2(i:iend), Ez(i:iend,m), EvalEz)
- Ez(i:iend,m) = - pot(i:iend,m)*gridwdir(1,i:iend) - Ez(i:iend,m)*gridwdir(0,i:iend) !- gtilde(1,i:iend)
- Er(i:iend,m) = - pot(i:iend,m)*gridwdir(2,i:iend) - Er(i:iend,m)*gridwdir(0,i:iend) !- gtilde(2,i:iend)
- pot(i:iend,m) = pot(i:iend,m)*gridwdir(0,i:iend) !+ gtilde(0,i:iend)
+ Ez(i:iend,m) = - pot(i:iend,m)*gridwdir(1,i:iend) - Ez(i:iend,m)*gridwdir(0,i:iend) - gtilde(1,i:iend)
+ Er(i:iend,m) = - pot(i:iend,m)*gridwdir(2,i:iend) - Er(i:iend,m)*gridwdir(0,i:iend) - gtilde(2,i:iend)
+ pot(i:iend,m) = pot(i:iend,m)*gridwdir(0,i:iend) + gtilde(0,i:iend)
! Add the external field
- Ez(i:iend,m) = Ez(i:iend,m) + Ezxt(i:iend)
- Er(i:iend,m) = Er(i:iend,m) + Erxt(i:iend)
- pot(i:iend,m) = pot(i:iend,m) + potxt(i:iend)
+ ! Ez(i:iend,m) = Ez(i:iend,m) + Ezxt(i:iend)
+ ! Er(i:iend,m) = Er(i:iend,m) + Erxt(i:iend)
+ ! pot(i:iend,m) = pot(i:iend,m) + potxt(i:iend)
END DO
!$OMP END DO NOWAIT
END DO
END IF
! Update the splines with dphi/dthet
DO m=1,nthet
!$OMP DO SIMD collapse(2)
DO j=1,nrank(2)
DO i=1,nrank(1)
matcoef(i,j) = dphidthet((j-1)*nrank(1)+i,m)
END DO
END DO
!$OMP END DO SIMD
! Update the ppform coefficients
CALL updt_ppform2d(splinevar(m), matcoef)
!$OMP BARRIER
END DO
! Evaluate Ethet at the particle position
!$OMP DO SIMD
DO i=1,p%Nploc,nbunch
! Avoid segmentation fault by accessing non relevant data
iend=min(i+nbunch-1,p%Nploc)
DO k=1,(iend-i+1)
thetindex=p%cellindex(2,i+k-1)
wthet = (p%pos(2,i+k-1)-thetgrid(thetindex))*invdthet
IF(thetindex .eq. nthet-1) THEN
CALL speval1(splrz(thetindex+1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), dphidthet1)
CALL speval1(splrz(1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), dphidthet2)
ELSE
CALL speval1(splrz(thetindex+1), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), dphidthet1)
CALL speval1(splrz(thetindex+2), p%pos(1,i+k-1), p%pos(3,i+k-1),p%cellindex(1,i+k-1),p%cellindex(3,i+k-1), dphidthet2)
END IF
- p%E(2,i+k-1)=p%geomweight(0,i+k-1)*((1-wthet)*dphidthet1+wthet*dphidthet2)/p%pos(1,i+k-1)
+ p%E(2,i+k-1)=-p%geomweight(0,i+k-1)*((1-wthet)*dphidthet1+wthet*dphidthet2)/p%pos(1,i+k-1)
END DO
END DO
!$OMP END DO SIMD NOWAIT
! On the root process, evaluate Ethet on the grid for diagnostic purposes
IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN
DO m=1,nthet
!$OMP DO
DO i=1,size(Ethet,1),16
iend=min(size(pot,1),i+15)
CALL gridval(splinevar(m), vec1(i:iend), vec2(i:iend), Ethet(i:iend,m), Evalpot)
DO k=i,iend
IF(rgrid(modulo(k-1,nr+1)) .eq. 0) THEN
Ethet(k,m) = 0
ELSE
- Ethet(k,m) = 1/rgrid(modulo(k-1,nr+1))*Ethet(k,m)*gridwdir(0,k)
+ Ethet(k,m) = -1/rgrid(modulo(k-1,nr+1))*Ethet(k,m)*gridwdir(0,k)
END IF
END DO
END DO
!$OMP END DO NOWAIT
END DO
END IF
END SUBROUTINE Efield_comp
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Constucts the FEM matrix using bsplines initialized in fields_init
!---------------------------------------------------------------------------
SUBROUTINE fematrix(mat,m)
USE bsplines
USE geometry
USE omp_lib
USE sparse
type(mumps_mat) :: mat
INTEGER :: m
REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :)
INTEGER, ALLOCATABLE :: f(:, :), aux(:)
REAL(kind=db), ALLOCATABLE :: coefs(:)
REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :)
REAL(kind=db) :: contrib
INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :)
integer,allocatable:: iid(:),jid(:)
INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize
kterms=8
If (allocated(fun)) deallocate (fun)
If (allocated(fun2)) deallocate (fun2)
ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2)))
If (allocated(wgeom)) deallocate (wgeom)
ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2))) !Arrays keeping values of b-splines at gauss node
ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines
ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2))
ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !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
CALL coefeq(splrz(1)%sp1%knots(0:1), idert, iderw, iderg, coefs, kterms)
! Assemble FEM matrix
!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, fun, fun2,iid,jid), collapse(2)
DO j = 1, nz ! Loop on z position
DO i = 1, nr ! Loop on r position
! Computation of gauss weight and position in r and z direction for gaussian integration
CALL calc_gauss(splrz(1), ngauss, i, j, xgauss, wgauss, gausssize)
iid=i
jid=j
IF(gausssize .gt. 1) THEN
CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize))
CALL basfun(xgauss(1:gausssize, 1), splrz(1)%sp1, fun(:,:,1:gausssize), iid(1:gausssize))
CALL basfun(xgauss(1:gausssize, 2), splrz(1)%sp2, fun2(:,:,1:gausssize), jid(1:gausssize))
END IF
DO jt = 1, (1 + femorder(1))*(femorder(2) + 1)
irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1
mu = irow + (jcol - 1)*nrank(1)
DO iw = 1, (1 + femorder(1))*(femorder(2) + 1)
irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1
mu2 = irow2 + (jcol2 - 1)*nrank(1)
contrib=0.0_db
DO igauss = 1, gausssize ! Loop on gaussian weights and positions
DO iterm = 1, kterms ! Loop on the two integration dimensions
contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* &
& fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* &
& fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* &
& wgauss(igauss)*xgauss(igauss, 1)
END DO
contrib = contrib + (m-1)**2*wgeom(0,igauss)*wgeom(0,igauss)*fun(f(jt, 1),0,igauss)*fun(f(iw, 1),0,igauss)*fun2(f(jt, 2),0,igauss)*fun2(f(iw, 2),0,igauss)*wgauss(igauss)/xgauss(igauss,1)
END DO
CALL omp_set_lock(mu_lock(mu))
CALL updt_sploc(mat%mat%row(mu), mu2, contrib)
CALL omp_unset_lock(mu_lock(mu))
END DO
END DO
END DO
END DO
!$OMP END PARALLEL DO
DEALLOCATE (f, aux)
DEALLOCATE (idert, iderw, coefs, fun, fun2)
END SUBROUTINE fematrix
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> @brief
!> Computes the volume of the splines cells needed to display the density in post-processing
!---------------------------------------------------------------------------
SUBROUTINE comp_volume
USE bsplines
USE geometry
USE basic, ONLY: Volume,nrank,rnorm,nthet
REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :)
INTEGER, ALLOCATABLE :: f(:, :), aux(:)
REAL(kind=db), ALLOCATABLE :: coefs(:)
REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), ftestpt(:, :)
Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw
INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms
Real(kind=db)::newcontrib
CALL timera(0, "comp_volume")
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(f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) ! Auxiliary arrays ordering bsplines
nterms = 4
ALLOCATE(idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms))
! 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
volume = 0
volume = 0
if (walltype .lt. 0) fverif = 0
volume = 0
if (walltype .lt. 0) fverif = 0
! Assemble Volume matrix
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib), collapse(2)
DO j = 1, nz ! Loop on z position
DO i = 1, nr ! Loop on r position
! Computation of gauss weight and position in r and z direction for gaussian integration
Call calc_gauss(splrz(1), ngauss, i, j, xgauss, wgauss, gausssize)
IF(allocated(wgeom)) DEALLOCATE(wgeom)
IF(gausssize .gt. 0) THEN
ALLOCATE(wgeom(0:2,size(xgauss, 1)))
CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom)
END IF
DO igauss = 1, gausssize ! Loop on gaussian weights and positions
CALL basfun(xgauss(igauss, 1), splrz(1)%sp1, fun, i)
CALL basfun(xgauss(igauss, 2), splrz(1)%sp2, fun2, j)
DO jt = 1, (1 + femorder(1))*(femorder(2) + 1)
irow = i + f(jt, 1) - 1;
jcol = j + f(jt, 2) - 1
mu = irow + (jcol - 1)*nrank(1)
newcontrib = 2*pi*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)*wgauss(igauss)*xgauss(igauss, 1)!*wgeom(igauss,0)
!$OMP ATOMIC UPDATE
volume(mu) = volume(mu) + newcontrib*rnorm**3
!$OMP END ATOMIC
END DO
END DO
END DO
END DO
!$OMP END PARALLEL DO
DEALLOCATE (f, aux)
DEALLOCATE (fun, fun2)
volume=volume/nthet
CALL timera(1, "comp_volume")
END SUBROUTINE comp_volume
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the gradient of the gtilde function for the web-spline method needed to correctly apply the dirichlet boundary conditions
!---------------------------------------------------------------------------
SUBROUTINE comp_gradgtilde
USE bsplines
USE geometry
REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :)
INTEGER, ALLOCATABLE :: f(:, :), aux(:)
REAL(kind=db), ALLOCATABLE :: coefs(:)
REAL(kind=db), ALLOCATABLE :: fun(:, :,:), fun2(:, :,:), gtildeintegr(:, :)
Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw
integer,allocatable:: iid(:),jid(:)
INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms
Real(kind=db)::newcontrib
!call timera(0, "comp_gradgtilde")
ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node
ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node
ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines
nterms = 4
Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms))
ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2)))
If (allocated(gtildeintegr)) deallocate (gtildeintegr)
ALLOCATE (gtildeintegr(0:2,3*ngauss(1)*ngauss(2)))
! 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
CALL coefeqext(splrz(1)%sp1%knots(0:1), idt, idw, idg, idp, coefs)
!$OMP DO SIMD
do j=1,size(gradgtilde)
gradgtilde(j) = 0
END DO
!$OMP END DO SIMD
!$OMP BARRIER
! Assemble gradgtilde matrix
!$OMP DO collapse(2)
DO j = 1, nz ! Loop on z position
DO i = 1, nr ! Loop on r position
! Computation of gauss weight and position in r and z direction for gaussian integration
Call calc_gauss(splrz(1), ngauss, i, j, xgauss, wgauss, gausssize)
iid=i
jid=j
if (gausssize .gt. 1) then
!If (allocated(wgeom)) deallocate (wgeom)
!ALLOCATE (wgeom(0:2,gausssize))
CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize))
CALL basfun(xgauss(1:gausssize, 1), splrz(1)%sp1, fun(:,:,1:gausssize), iid(1:gausssize))
CALL basfun(xgauss(1:gausssize, 2), splrz(1)%sp2, fun2(:,:,1:gausssize), jid(1:gausssize))
Call total_gtilde(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), gtildeintegr(:,1:gausssize),wgeom(:,1:gausssize))
End if
DO jt = 1, (1 + femorder(1))*(femorder(2) + 1)
irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1
mu = irow + (jcol - 1)*nrank(1)
newcontrib = 0.0_db
DO igauss = 1, gausssize ! Loop on gaussian weights and positions
Do iterm = 1, nterms
newcontrib = newcontrib + wgeom( idg(iterm),igauss)*gtildeintegr( idp(iterm),igauss)* &
& fun(f(jt, 1), idt(iterm),igauss)*fun2(f(jt, 2), idw(iterm),igauss)* &
& wgauss(igauss)*xgauss(igauss, 1)
End do
end do
!$OMP ATOMIC UPDATE
gradgtilde(mu) = gradgtilde(mu) + newcontrib
!$OMP END ATOMIC
END DO
END DO
END DO
!$OMP END DO
!DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg)
DEALLOCATE (f, aux)
DEALLOCATE (fun, fun2)
!call timera(1, "comp_gradgtilde")
END SUBROUTINE comp_gradgtilde
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Imposes the dirichlet boundary conditions on the FEM matrix for the case where we use regular splines ( not web-splines).
!---------------------------------------------------------------------------
SUBROUTINE fe_dirichlet
REAL(kind=db), ALLOCATABLE :: arr(:)
INTEGER :: i
ALLOCATE (arr(nrank(1)*nrank(2)))
DO i = 1, nrank(1)
IF (rgrid(0) .ne. 0.0_db) THEN
arr = 0; arr(i) = 1;
CALL putrow(femat(1), i, arr)
END IF
arr = 0; arr(nrank(1)*nrank(2) + 1 - i) = 1;
CALL putrow(femat(1), nrank(1)*nrank(2) + 1 - i, arr)
END DO
DEALLOCATE (arr)
END SUBROUTINE fe_dirichlet
!________________________________________________________________________________
SUBROUTINE coefeq(x, idt, idw, idg, c, kterms)
REAL(kind=db), INTENT(in) :: x(:)
INTEGER, INTENT(out) :: idt(:, :), idw(:, :), idg(:, :),kterms
REAL(kind=db), INTENT(out) :: c(:)
kterms=8
c = x(2)
idt(1, 1) = 0
idt(1, 2) = 0
idw(1, 1) = 0
idw(1, 2) = 0
idg(1, 1) = 1
idg(1, 2) = 1
idt(2, 1) = 0
idt(2, 2) = 0
idw(2, 1) = 0
idw(2, 2) = 1
idg(2, 1) = 1
idg(2, 2) = 0
idt(3, 1) = 0
idt(3, 2) = 0
idw(3, 1) = 1
idw(3, 2) = 0
idg(3, 1) = 0
idg(3, 2) = 1
idt(4, 1) = 0
idt(4, 2) = 0
idw(4, 1) = 1
idw(4, 2) = 1
idg(4, 1) = 0
idg(4, 2) = 0
idt(5, 1) = 0
idt(5, 2) = 0
idw(5, 1) = 0
idw(5, 2) = 0
idg(5, 1) = 2
idg(5, 2) = 2
idt(6, 1) = 0
idt(6, 2) = 1
idw(6, 1) = 0
idw(6, 2) = 0
idg(6, 1) = 2
idg(6, 2) = 0
idt(7, 1) = 1
idt(7, 2) = 0
idw(7, 1) = 0
idw(7, 2) = 0
idg(7, 1) = 0
idg(7, 2) = 2
idt(8, 1) = 1
idt(8, 2) = 1
idw(8, 1) = 0
idw(8, 2) = 0
idg(8, 1) = 0
idg(8, 2) = 0
END SUBROUTINE coefeq
SUBROUTINE coefeqext(x, idt, idw, idg, idp, c)
REAL(kind=db), INTENT(in) :: x(:)
INTEGER, INTENT(out) :: idp(:), idt(:), idw(:), idg(:)
REAL(kind=db), INTENT(out) :: c(:)
c(1) = x(2)
idp(1) = 1
idg(1) = 1
idt(1) = 0
idw(1) = 0
c(2) = x(2)
idp(2) = 1
idg(2) = 0
idt(2) = 0
idw(2) = 1
c(3) = x(2)
idp(3) = 2
idg(3) = 2
idt(3) = 0
idw(3) = 0
c(4) = x(2)
idp(4) = 2
idg(4) = 0
idt(4) = 1
idw(4) = 0
END SUBROUTINE coefeqext
!---------------------------------------------------------------------------
!> @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
USE basic, ONLY: rhs
INTEGER:: i,m
DO i = 1, nrank(1)*nrank(2)
CALL omp_destroy_lock(mu_lock(i))
END DO
DEALLOCATE(mu_lock)
DEALLOCATE(matcoef)
DEALLOCATE(pot)
DEALLOCATE(rhs)
DEALLOCATE(dftrhsreal)
DEALLOCATE(dftrhsimag)
DEALLOCATE(loc_rhs)
DEALLOCATE(loc_moments)
DEALLOCATE(phi)
DEALLOCATE(dphidthet)
DEALLOCATE(dftphireal)
DEALLOCATE(dftphiimag)
DEALLOCATE(Br,Bz)
DEALLOCATE(Er,Ethet,Ez)
DEALLOCATE(vec1,vec2)
DO m=1,nthet
CALL DESTROY_SP(splrz(m))
END DO
DEALLOCATE(splrz)
CALL DESTROY_SP(splrz_ext)
END SUBROUTINE clean_fields
SUBROUTINE updt_sploc(arow, j, val)
!
! Update element j of row arow or insert it in an increasing "index"
!
USE sparse
TYPE(sprow), TARGET :: arow
INTEGER, INTENT(in) :: j
DOUBLE PRECISION, INTENT(in) :: val
!
TYPE(elt), TARGET :: pre_root
TYPE(elt), POINTER :: t, p
!
if(val.eq.0) return
pre_root%next => arow%row0 ! pre_root is linked to the head of the list.
t => pre_root
DO WHILE (ASSOCIATED(t%next))
p => t%next
IF (p%index .EQ. j) THEN
p%val = p%val + val
RETURN
END IF
IF (p%index .GT. j) EXIT
t => t%next
END DO
ALLOCATE (p)
p = elt(j, val, t%next)
t%next => p
!
arow%nnz = arow%nnz + 1
arow%row0 => pre_root%next ! In case the head is altered
END SUBROUTINE updt_sploc
SUBROUTINE updt_ppform2d(sp,c)
USE bsplines
TYPE(spline2d), INTENT(inout) :: sp
DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: c
!DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:)
INTEGER:: i,m,mm
INTEGER :: d1, d2, k1, k2, n1, n2
d1 = sp%sp1%dim
d2 = sp%sp2%dim
k1 = sp%sp1%order
k2 = sp%sp2%order
n1 = sp%sp1%nints
n2 = sp%sp2%nints
!$OMP SINGLE
IF( ASSOCIATED(sp%bcoefs) ) DEALLOCATE(sp%bcoefs)
ALLOCATE(sp%bcoefs(SIZE(c,1),SIZE(c,2)))
!$OMP END SINGLE
!ALLOCATE(work(d2,k1,n1))
!$OMP DO
DO m=1,SIZE(c,2)
CALL topp0(sp%sp1, c(:,m), ppformwork(m,:,:))
sp%bcoefs(:,m)=c(:,m)
END DO
!$OMP END DO NOWAIT
!$OMP SINGLE
IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform)
ALLOCATE(sp%ppform(k1,n1,k2,n2))
!$OMP END SINGLE
!$OMP DO
DO mm=1,SIZE(ppformwork,3)
DO m=1,SIZE(ppformwork,2)
CALL topp0(sp%sp2, ppformwork(:,m,mm), sp%ppform(m,mm,:,:))
END DO
END DO
!$OMP END DO
!DEALLOCATE(work)
END SUBROUTINE updt_ppform2d
!===========================================================================
SUBROUTINE topp0(sp, c, ppform)
!
! Compute PPFORM of a fuction defined by the spline SP
! and spline coefficients C(1:d)
!
use bsplines
TYPE(spline1d), INTENT(in) :: sp
DOUBLE PRECISION, INTENT(in) :: c(:)
DOUBLE PRECISION, INTENT(out) :: ppform(0:,:)
INTEGER :: p, nints, i, j, k
!
p = sp%order - 1
nints = sp%nints
!
ppform = 0.0d0
DO i=1,nints ! on each knot interval
DO j=1,p+1 ! all spline in interval i
DO k=0,p ! k_th derivatives
ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1)
END DO
END DO
END DO
!
END SUBROUTINE topp0
!+
END MODULE fields
diff --git a/src/mpihelper_mod.f90_old b/src/mpihelper_mod.f90_old
deleted file mode 100644
index 764c19d..0000000
--- a/src/mpihelper_mod.f90_old
+++ /dev/null
@@ -1,401 +0,0 @@
-!------------------------------------------------------------------------------
-! EPFL/Swiss Plasma Center
-!------------------------------------------------------------------------------
-!
-! MODULE: mpihelper
-!
-!> @author
-!> Guillaume Le Bars EPFL/SPC
-!
-! DESCRIPTION:
-!> Module responsible for setting up the MPI variables used in the communications.
-!------------------------------------------------------------------------------
-
-MODULE mpihelper
- USE constants
- use mpi
- USE particletypes
-
- IMPLICIT NONE
-
- INTEGER, SAVE :: basicdata_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for communicating basicdata
- INTEGER, SAVE :: particle_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for particles communication between nodes
- !INTEGER, SAVE :: particles_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for particles gathering on node 0 and broadcast from node 0
- INTEGER, SAVE :: rhsoverlap_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for the communication of a rhs column
- INTEGER, SAVE :: densityoverlap_type=MPI_DATATYPE_NULL
- INTEGER, SAVE :: db_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for the communication of a REAL(kind=db)
- INTEGER, SAVE :: momentsoverlap_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for the communication of a column of a grid variable
- INTEGER, SAVE :: rcvrhsoverlap_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for the receive communication of a rhs column
- INTEGER, SAVE :: rcvdensityoverlap_type=MPI_DATATYPE_NULL
- INTEGER, SAVE :: rcvmomentsoverlap_type=MPI_DATATYPE_NULL !< Stores the MPI data type used for the receive communication of a column of a grid variable
- INTEGER, SAVE:: db_sum_op !< Store the MPI sum operation for db_type
- REAL(kind=db), ALLOCATABLE, SAVE:: rhsoverlap_buffer(:) !< buffer used for storing the rhs ghost cells
- !< received from the left or right MPI process
- REAL(kind=db), ALLOCATABLE, SAVE:: densityoverlap_buffer(:)
- REAL(kind=db), ALLOCATABLE, SAVE:: momentsoverlap_buffer(:) !< buffer used for storing the moments ghost cells
- !< received from the left or right MPI process
- !INTEGER, SAVE:: momentsoverlap_requests(2) = MPI_REQUEST_NULL
- !INTEGER, SAVE:: rhsoverlap_requests(2) = MPI_REQUEST_NULL
- INTEGER:: rhsoverlap_tag= 200
- INTEGER:: densityoverlap_tag=400
- INTEGER:: momentsoverlap_tag= 300
- INTEGER:: partsgather_tag= 500
- INTEGER:: partsexchange_tag=600
- INTEGER:: nbpartsexchange_tag=700
-
-
- CONTAINS
-
- !---------------------------------------------------------------------------
- !> @author
- !> Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !>
- !> @brief
- !> Initialize the MPI types used for inter process communications
- !
- !---------------------------------------------------------------------------
- SUBROUTINE mpitypes_init
- IMPLICIT NONE
- INTEGER:: ierr
- ! Initialize db_type to use real(kind=db) in MPI and the sum operator for reduce
- CALL MPI_TYPE_CREATE_F90_REAL(dprequestedprec,MPI_UNDEFINED,db_type,ierr)
- CALL MPI_Type_commit(db_type,ierr)
- CALL MPI_Op_Create(DB_sum, .true., db_sum_op, ierr)
-
-
- CALL init_particlempi
-
- END SUBROUTINE mpitypes_init
-
- !---------------------------------------------------------------------------
- !> @author
- !> Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !>
- !> @brief
- !> Computes the sum in MPI_Reduce operations involving Real(kinc=db)
- !
- !---------------------------------------------------------------------------
-
- SUBROUTINE DB_sum(INVEC, INOUTVEC, LEN, TYPE)bind(c)
- !REAL(kind=db):: INVEC(0:LEN-1), INOUTVEC(0:LEN-1)
- use, intrinsic:: iso_c_binding, ONLY: c_ptr,c_f_pointer
- use mpi
- implicit none
- TYPE(C_PTR), VALUE:: invec, inoutvec
- real(kind=db),pointer:: ivec(:), iovec(:)
- INTEGER:: LEN
- INTEGER:: TYPE
- INTEGER:: i
-
- call c_f_pointer(INVEC,ivec, (/len/))
- call c_f_pointer(inoutvec,iovec, (/len/))
- Do i=1,LEN
- IOVEC(i)=IVEC(i)+IOVEC(i)
- END DO
-
-
-
- END SUBROUTINE DB_sum
-
- !--------------------------------------------------------------------------
- !> @author
- !> Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !>
- !> @brief
- !> Initialize the MPI communicators used for allreduce between neighbors
- !
- !> @param[in] nrank ranks of the FEM array in (1) z direction and (2) r direction
- !> @param[in] femorder finite element method order in z and r direction
- !> @param[in] zlimleft z index delimiting the mpi local left boundary
- !> @param[in] zlimright z index delimiting the mpi local right boundary
- !> @param[in] nbmoments number of moments calculated and stored.
- !
- !---------------------------------------------------------------------------
- SUBROUTINE init_overlaps(nrank, femorder, zlimleft, zlimright, nbmoments)
- INTEGER, INTENT(IN):: nrank(:), femorder(:), zlimright, zlimleft, nbmoments
-
- IF(ALLOCATED(rhsoverlap_buffer)) DEALLOCATE(rhsoverlap_buffer)
- IF(ALLOCATED(densityoverlap_buffer)) DEALLOCATE(densityoverlap_buffer)
- IF(ALLOCATED(momentsoverlap_buffer)) DEALLOCATE(momentsoverlap_buffer)
- ALLOCATE(rhsoverlap_buffer(nrank(1)*femorder(2)))
- ALLOCATE(densityoverlap_buffer(nrank(1)*femorder(2)))
- ALLOCATE(momentsoverlap_buffer(nbmoments*nrank(1)*femorder(2)))
-
- ! Initialize the MPI column overlap type for rhs
- CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(2), 1, 1, db_type, rhsoverlap_type)
- CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(2), 1, 1, db_type, densityoverlap_type)
- ! Initialize the MPI grid col type
- CALL init_coltypempi(nrank(2), zlimright-zlimleft+femorder(2), nbmoments, 1, db_type, momentsoverlap_type)
- ! Initialize the MPI receive column overlap type for rhs
- CALL init_coltypempi(nrank(1), nrank(2), 1, 1, db_type, rcvrhsoverlap_type)
- CALL init_coltypempi(nrank(1), nrank(2), 1, 1, db_type, rcvdensityoverlap_type)
- ! Initialize the MPI receive grid col type
- CALL init_coltypempi(nrank(1), nrank(2), nbmoments, 1, db_type, rcvmomentsoverlap_type)
-
- END SUBROUTINE init_overlaps
-
-
- SUBROUTINE start_persistentcomm(requests, mpirank, leftproc, rightproc)
- INTEGER,INTENT(INOUT):: requests(:)
- INTEGER, INTENT(IN):: mpirank, leftproc, rightproc
- INTEGER:: ierr
- INTEGER:: stats(MPI_STATUS_SIZE,2)
- LOGICAL:: completed=.false.
-
- IF(leftproc .lt. mpirank) THEN
- ! Start to receive
- CALL MPI_START(requests(2),ierr)
- IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in recv_init"
- END IF
-
- IF(rightproc .gt. mpirank) THEN
- ! Start to send
- CALL MPI_START(requests(1),ierr)
- IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in send_init"
- END IF
-
- IF(leftproc .lt. mpirank) THEN
- ! Start to receive
- completed=.FALSE.
- DO WHILE(.not. completed)
- CALL MPI_TEST(requests(2), completed,stats(:,2),ierr)
- END DO
- WRITE(*,*)"status 2", completed, stats(:,2)
- !CALL MPI_WAIT(requests(2),stats(:,2),ierr)
- !WRITE(*,*)"status 2", stats(:,2)
- IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in recv_init"
- END IF
-
- IF(rightproc .gt. mpirank) THEN
- ! Start to send
- completed=.FALSE.
- DO WHILE(.not. completed)
- CALL MPI_TEST(requests(1), completed,stats(:,1),ierr)
- END DO
- !CALL MPI_WAIT(requests(1),stats(:,1),ierr)
- IF(IERR .ne. MPI_SUCCESS) WRITE(*,*) "error in send_init"
- END IF
- END SUBROUTINE start_persistentcomm
-
-
- SUBROUTINE rhsoverlapcomm(mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright)
- INTEGER, INTENT(IN):: mpirank, leftproc, rightproc
- REAL(kind=db), DIMENSION(:), INTENT(INOUT):: moments
- INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright
- INTEGER, SAVE:: rhsoverlap_requests(2) = MPI_REQUEST_NULL
- INTEGER:: ierr
- INTEGER:: stats(MPI_STATUS_SIZE,2)
- rhsoverlap_requests=MPI_REQUEST_NULL
- rhsoverlap_buffer=0
- IF(rightproc .gt. mpirank .and. rightproc .ge. 0) THEN
- CALL MPI_ISEND(moments(zlimright*nrank(1)+1), femorder(2)*nrank(1), db_type, rightproc, rhsoverlap_tag, &
- & MPI_COMM_WORLD, rhsoverlap_requests(1), ierr )
- END IF
-
- ! If the processor on the left has actually lower z positions
- IF(leftproc .lt. mpirank .and. leftproc .ge. 0) THEN
- CALL MPI_IRECV(rhsoverlap_buffer, nrank(1)*(femorder(2)), db_type, leftproc, rhsoverlap_tag, &
- & MPI_COMM_WORLD, rhsoverlap_requests(2), ierr )
- END IF
-
- CALL MPI_WAITALL(2,rhsoverlap_requests,stats, ierr)
- END SUBROUTINE rhsoverlapcomm
-
- SUBROUTINE densityoverlapcomm(mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright)
- INTEGER, INTENT(IN):: mpirank, leftproc, rightproc
- REAL(kind=db), DIMENSION(:), INTENT(INOUT):: moments
- INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright
- INTEGER, SAVE:: densityoverlap_requests(2) = MPI_REQUEST_NULL
- INTEGER:: ierr
- INTEGER:: stats(MPI_STATUS_SIZE,2)
- densityoverlap_requests=MPI_REQUEST_NULL
- densityoverlap_buffer=0
- IF(rightproc .gt. mpirank .and. rightproc .ge. 0) THEN
- CALL MPI_ISEND(moments(zlimright*nrank(1)+1), femorder(2)*nrank(1), db_type, rightproc, densityoverlap_tag, &
- & MPI_COMM_WORLD, densityoverlap_requests(1), ierr )
- END IF
-
- ! If the processor on the left has actually lower z positions
- IF(leftproc .lt. mpirank .and. leftproc .ge. 0) THEN
- CALL MPI_IRECV(densityoverlap_buffer, nrank(1)*(femorder(2)), db_type, leftproc, densityoverlap_tag, &
- & MPI_COMM_WORLD, densityoverlap_requests(2), ierr )
- END IF
-
- CALL MPI_WAITALL(2,densityoverlap_requests,stats, ierr)
- END SUBROUTINE densityoverlapcomm
-
- SUBROUTINE momentsoverlapcomm(mpirank, leftproc, rightproc, moments, nrank, femorder, zlimright)
- INTEGER, INTENT(IN):: mpirank, leftproc, rightproc
- REAL(kind=db), DIMENSION(:,:), INTENT(INOUT):: moments
- INTEGER, INTENT(IN):: nrank(2), femorder(2), zlimright
- INTEGER, SAVE:: momentsoverlap_requests(2) = MPI_REQUEST_NULL
- INTEGER:: ierr
- INTEGER:: stats(MPI_STATUS_SIZE,2)
- momentsoverlap_requests=MPI_REQUEST_NULL
- momentsoverlap_buffer=0
-
- IF(rightproc .gt. mpirank .and. rightproc .ge. 0) THEN
- CALL MPI_ISEND(moments(1,zlimright*nrank(1)+1), femorder(2)*nrank(1)*10, db_type, rightproc, momentsoverlap_tag, &
- & MPI_COMM_WORLD, momentsoverlap_requests(1), ierr )
- END IF
-
- ! If the processor on the left has actually lower z positions
- IF(leftproc .lt. mpirank .and. leftproc .ge. 0) THEN
- CALL MPI_IRECV(momentsoverlap_buffer, 10*nrank(1)*(femorder(2)), db_type, leftproc, momentsoverlap_tag, &
- & MPI_COMM_WORLD, momentsoverlap_requests(2), ierr )
- END IF
-
- CALL MPI_WAITALL(2,momentsoverlap_requests,stats, ierr)
- END SUBROUTINE momentsoverlapcomm
-
- !---------------------------------------------------------------------------
- !> @author
- !> Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !>
- !> @brief
- !> Initialize the particle MPI type used for inter process communications and publish it to
- !> the process in the communicator
- !
- !---------------------------------------------------------------------------
- SUBROUTINE init_particlempi()
- INTEGER :: nblock = 7
- INTEGER:: blocklength(7)
- INTEGER(kind=MPI_ADDRESS_KIND):: displs(7), displ0
- INTEGER:: types(7)
- TYPE(particle) :: part
- INTEGER:: ierr
-
- CALL mpi_get_address(part%partindex, displs(1), ierr)
- CALL mpi_get_address(part%cellindex, displs(2), ierr)
- types(1:2)=MPI_INTEGER
- CALL mpi_get_address(part%pos, displs(3), ierr)
- CALL mpi_get_address(part%U, displs(4), ierr)
- CALL mpi_get_address(part%geomweight, displs(5), ierr)
- CALL mpi_get_address(part%GAMMA, displs(6), ierr)
- CALL mpi_get_address(part%pot, displs(7), ierr)
- types(3:7)=db_type
- blocklength(1:7) = 1
- blocklength(2:5) =3
- CALL mpi_get_address(part, displ0, ierr)
- displs=displs-displ0
-
-
- CALL MPI_Type_create_struct(nblock, blocklength, displs, types, particle_type, ierr)
- CALL MPI_Type_commit(particle_type,ierr)
-
- END SUBROUTINE init_particlempi
-
-
- !---------------------------------------------------------------------------
- !> @author
- !> Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !>
- !> @brief
- !> Initialize the particles MPI type used for gathering particles to the root and broadcast them and publish it to
- !> the process in the communicator
- !
- !---------------------------------------------------------------------------
- SUBROUTINE init_particles_gather_mpi(p,idstart,nsend,mpi_particles_type)
- INTEGER:: mpi_particles_type
- INTEGER:: nsend
- INTEGER:: idstart
- INTEGER :: nblock = 8
- INTEGER:: blocklength(8)
- INTEGER(kind=MPI_ADDRESS_KIND):: displs(8), displ0
- INTEGER:: types(8)
- TYPE(particles), INTENT(INOUT):: p
- INTEGER:: ierr
- INTEGER:: temptype
-
- IF(nsend .lt. 1) RETURN
-
- temptype=MPI_DATATYPE_NULL
- IF( mpi_particles_type .ne. MPI_DATATYPE_NULL) CALL MPI_TYPE_FREE(mpi_particles_type,ierr)
-
- !CALL mpi_get_address(p%Rindex(idstart), displs(1), ierr)
- !CALL mpi_get_address(p%Zindex(idstart), displs(2), ierr)
- CALL mpi_get_address(p%cellindex(1,idstart), displs(1), ierr)
- CALL mpi_get_address(p%partindex(idstart), displs(2), ierr)
- types(1:2)=MPI_INTEGER
- CALL mpi_get_address(p%pos(1,idstart), displs(3), ierr)
- !CALL mpi_get_address(p%Z(idstart), displs(5), ierr)
- !CALL mpi_get_address(p%THET(idstart), displs(6), ierr)
- CALL mpi_get_address(p%pot(idstart), displs(4), ierr)
- CALL mpi_get_address(p%U(1,idstart), displs(5), ierr)
- CALL mpi_get_address(p%Uold(1,idstart), displs(6), ierr)
- !CALL mpi_get_address(p%UTHET(idstart), displs(10), ierr)
- !CALL mpi_get_address(p%UTHETold(idstart), displs(11), ierr)
- !CALL mpi_get_address(p%UZ(idstart), displs(12), ierr)
- !CALL mpi_get_address(p%UZold(idstart), displs(13), ierr)
-
- CALL mpi_get_address(p%GAMMA(idstart), displs(7), ierr)
- CALL mpi_get_address(p%GAMMAold(idstart), displs(8), ierr)
- types(3:8)=db_type
- blocklength = nsend
- blocklength(1)=3*nsend
- blocklength(3)=3*nsend
- blocklength(5:6)=3*nsend
- CALL mpi_get_address(p, displ0, ierr)
- displs=displs-displ0
-
- CALL MPI_Type_create_struct(nblock, blocklength, displs, types, mpi_particles_type, ierr)
- CALL MPI_TYPE_COMMIT(mpi_particles_type, ierr)
-
- END SUBROUTINE init_particles_gather_mpi
-
-
- !---------------------------------------------------------------------------
- !> @author Guillaume Le Bars EPFL/SPC
- !
- ! DESCRIPTION:
- !> @brief Initialize the column MPI type used for inter process communications and publish it to
- !> the processes in the communicator (can be rhs or grid quantities)
- !
- !> @param[in] nr number of elements in the r direction
- !> @param[in] nz number of elements in the z direction
- !> @param[in] init_type MPI type of the initial data
- !> @param[inout] mpi_coltype final type usable in communications
- !---------------------------------------------------------------------------
- SUBROUTINE init_coltypempi(nr, nz, block_size, stride, init_type, mpi_coltype)
- INTEGER, INTENT(IN) :: nr
- INTEGER, INTENT(IN) :: nz
- INTEGER, INTENT(IN) :: block_size
- INTEGER, INTENT(IN) :: stride
- INTEGER, INTENT(IN) :: init_type
- INTEGER, INTENT(OUT) :: mpi_coltype
- INTEGER :: temp_mpi_coltype
- INTEGER:: ierr
- INTEGER(KIND=MPI_ADDRESS_KIND):: init_type_lb, init_type_extent
- !(nrank(2), nrank(1), 1, 10, db_type, rhsoverlap_type)
- ! if mpi_coltype was used, we free it first
- IF( mpi_coltype .ne. MPI_DATATYPE_NULL) CALL MPI_TYPE_FREE(mpi_coltype,ierr)
- ! Create vector type of length nx
- CALL MPI_TYPE_VECTOR(nr, block_size, stride*block_size*nz, init_type, temp_mpi_coltype, ierr)
- CALL MPI_TYPE_COMMIT(temp_mpi_coltype, ierr)
-
- ! Get the size in bytes of the initial type
- CALL MPI_TYPE_GET_EXTENT(init_type, init_type_lb, init_type_extent, ierr)
-
- if(mpi_coltype .ne. MPI_DATATYPE_NULL) CALL MPI_TYPE_FREE(mpi_coltype,ierr)
- ! Resize temp_mpi_coltype such that the next item to read is at j+1
- CALL MPI_TYPE_CREATE_RESIZED(temp_mpi_coltype, init_type_lb, stride*block_size*init_type_extent ,&
- & mpi_coltype, ierr)
- CALL MPI_TYPE_COMMIT(mpi_coltype, ierr)
-
- CALL MPI_TYPE_FREE(temp_mpi_coltype,ierr)
-
-
- END SUBROUTINE init_coltypempi
-
- END MODULE mpihelper
-
\ No newline at end of file