Homec4science
Diffusion FENNECS lite (master) Browse

fennecslite/src/fields_mod.f90
6e6d693148f1master

Last Change
Raw File

fields_mod.f90
View Options

!------------------------------------------------------------------------------
! 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.f'
REAL(kind=db), allocatable, SAVE :: matcoef(:,:), 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(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 Update_phi2D(splrz_ext)
!$OMP BARRIER
!$OMP DO
! On the root process, compute the electric field for diagnostic purposes
DO i=1,size(pot,1),16
iend=min(size(pot,1),i+15)
potxt(i:iend) = pot(i:iend,1)
Erxt(i:iend) = Er(i:iend,1)
Ezxt(i:iend) = Ez(i:iend,1)
END DO
!$OMP END DO NOWAIT
!$OMP END PARALLEL
END SUBROUTINE vacuum_field
!---------------------------------------------------------------------------
!> @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
! FFT specific variables
TYPE(C_PTR) :: plan_theta2m ! transform from (r,theta,z) to (r,m,z)
REAL(kind=db), DIMENSION(:,:), ALLOCATABLE :: rhstemp
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE :: dftrhstemp
ALLOCATE(rhstemp(nthet,nrank(1)*nrank(2)))
ALLOCATE(dftrhstemp(nthet/2+1,nrank(1)*nrank(2)))
IF (nlphis) THEN ! We calculate the self-consistent field
!$OMP DO SIMD
DO i=1,size(loc_rhs,1)
DO j=1,size(loc_rhs,2)
loc_rhs(i,j)=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(rhs,1)
! loc_rhs(i,:)=loc_rhs(i,:)-gradgtilde(i)
! 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)
! rhs(i,:)=loc_rhs(i,:)
! END DO
! !$OMP END DO
! END IF
! In serial copy loc_rhs to rhs
!$OMP DO
DO i=1,size(loc_rhs,1)
DO j=1,size(loc_rhs,2)
rhs(i,j)=loc_rhs(i,j)
END DO
END DO
!$OMP END DO
ELSE ! We only consider the externally imposed field
!$OMP DO
DO i=1,size(rhs,1)
DO j=1,size(rhs,2)
rhs(i,j)=0
END DO
END DO
!$OMP END DO
END IF
! 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_FORWARD,FFTW_PATIENT)
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)
END SUBROUTINE rhscon
!---------------------------------------------------------------------------
!> @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
!> 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
! INTEGER:: ierr, i, j
! !$OMP MASTER
! !WRITE(*,*) mpirank, "wE communicate overlap rhs"
! CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs, nrank, femorder, loc_zspan - femorder(2))
! !$OMP END MASTER
! !$OMP BARRIER
! IF (mpirank .gt. 0) THEN
! !$OMP DO SIMD
! DO i = 1, femorder(2)*nrank(1)
! loc_rhs(i) = loc_rhs(i)&
! & + rhsoverlap_buffer(i)
! END DO
! !$OMP END DO SIMD
! END IF
! !!$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
! REAL(kind=db), DIMENSION(:), INTENT(INOUT):: rhs
! INTEGER:: ierr, i, j
! 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
! CALL MPI_GATHERV(loc_rhs, counts(mpirank + 1), db_type, &
! & rhs, counts, displs, db_type, 0, MPI_COMM_WORLD, ierr)
! !$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
! REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: moment
! INTEGER:: ierr, i, j
! 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
! CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments, 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) = loc_moments(1:nbmoments, i)&
! & + momentsoverlap_buffer((i-1)*nbmoments+1:i*nbmoments)
! END DO
! !$OMP END DO SIMD
! END IF
! !$OMP MASTER
! ! Set communication vector type
! IF (mpirank .eq. 0) THEN
! moment = 0
! END IF
! CALL MPI_GATHERV(loc_moments, counts(mpirank + 1), db_type, &
! & moment, counts, displs, db_type, 0, MPI_COMM_WORLD, ierr)
! !$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
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(:)
! 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)
ALLOCATE(reducedrhsreal(nrank(1)*nrank(2)),reducedrhsimag(nrank(1)*nrank(2)))
ALLOCATE(reducedsolreal(nbreducedspline),reducedsolimag(nbreducedspline))
ALLOCATE(tempcol(nrank(1)*nrank(2)))
DO m=1,nthet/2+1
! Solves Poisson equation
IF(nlweb .and. mpirank .eq. 0) THEN ! We use the web-spline reduction for stability
reducedrhsreal = vmx(etilde, dftrhsreal(:,m))
reducedrhsimag = vmx(etilde, dftrhsimag(:,m))
CALL bsolve(reduccedmat(m), reducedrhsreal(1:nbreducedspline), reducedsolreal)
CALL bsolve(reduccedmat(m), reducedrhsimag(1:nbreducedspline), reducedsolimag)
ELSE IF(mpirank.eq.0) THEN
CALL bsolve(femat(m), dftrhsreal(:,m), dftphireal(:,m))
CALL bsolve(femat(m), dftrhsimag(:,m), dftphiimag(:,m))
END IF
! Distributes the result on all MPI workers
IF(nlweb) THEN
CALL MPI_Bcast(reducedsolreal, nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr)
CALL MPI_Bcast(reducedsolimag, 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
dftphireal(:,m) = vmx(etildet, tempcol)
tempcol(1:nbreducedspline) = reducedsolimag
dftphiimag(:,m) = vmx(etildet, tempcol)
END IF
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)
ELSE IF(mpirank.eq.0) THEN
CALL bsolve(femat(1), rhs(:,1), phi(:,1))
END IF
! Distributes the result on all MPI workers
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)
END SUBROUTINE poisson2D
!---------------------------------------------------------------------------
!> @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 Update_phi
USE basic, ONLY: nrank, pot, nlend, gradthet, rgrid, nr, nz
USE futils
USE geometry
INTEGER:: i,m
! FFT specific variables
TYPE(C_PTR) :: plan_m2theta ! transform from (r,m,z) to (r,theta,z)
REAL(kind=db), DIMENSION(:,:), ALLOCATABLE :: phitemp
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE :: dftphitemp
ALLOCATE(dftphitemp(nthet/2+1,nrank(1)*nrank(2)))
ALLOCATE(phitemp(nthet,nrank(1)*nrank(2)))
! 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
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_BACKWARD,FFTW_PATIENT)
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_BACKWARD,FFTW_PATIENT)
CALL dfftw_execute_dft_c2r(plan_m2theta,dftphitemp,phitemp)
CALL dfftw_destroy_plan(plan_m2theta)
dphidthet=transpose(phitemp)
DEALLOCATE(phitemp)
DEALLOCATE(dftphitemp)
END SUBROUTINE Update_phi
!---------------------------------------------------------------------------
!> @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 Update_phi2D(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)
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 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN
!$OMP DO
DO i=1,size(pot,1),16
iend=min(size(pot,1),i+15)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), pot(i:iend,1), Evalpot)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Er(i:iend,1), EvalEr)
CALL gridval(splinevar, vec1(i:iend), vec2(i:iend), Ez(i:iend,1), EvalEz)
Ez(i:iend,1) = -pot(i:iend,1)*gridwdir(1,i:iend) - Ez(i:iend,1)*gridwdir(0,i:iend) - gtilde(1,i:iend)
Er(i:iend,1) = -pot(i:iend,1)*gridwdir(2,i:iend) - Er(i:iend,1)*gridwdir(0,i:iend) - gtilde(2,i:iend)
pot(i:iend,1) = pot(i:iend,1)*gridwdir(0,i:iend) + gtilde(0,i:iend)
END DO
!$OMP END DO NOWAIT
END IF
END SUBROUTINE Update_phi2D
!---------------------------------------------------------------------------
!> @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 Efieldscomp(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, Ethet1, Ethet2
REAL(kind=db), ALLOCATABLE :: gtildeloc(:,:)
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))
! Update the splines
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)
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)
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
p%pot(i+k-1)=(1-wthet)*pot1+wthet*pot2
p%E(1,i+k-1)=(1-wthet)*Er1+wthet*Er2
p%E(3,i+k-1)=(1-wthet)*Ez1+wthet*Ez2
END DO
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%potxt(i:iend))
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)
p%potxt(i:iend) = p%geomweight(0,i:iend)*p%potxt(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)
END DO
!$OMP END DO NOWAIT
END DO
END IF
! Update the splines
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), Ethet1)
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), Ethet2)
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), Ethet1)
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), Ethet2)
END IF
p%E(2,i+k-1)=(1-wthet)*Ethet1+wthet*Ethet2
p%E(2,i+k-1)=0
END DO
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(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
Ethet(k,m) = 1/rgrid(modulo(k-1,nr+1))*Ethet(k,m)*gridwdir(0,k)
END DO
END DO
!$OMP END DO NOWAIT
END DO
END IF
END SUBROUTINE Efieldscomp
! !---------------------------------------------------------------------------
! !> @author
! !> Guillaume Le Bars EPFL/SPC
! !
! ! DESCRIPTION:
! !> @brief
! !> Computes the electric fields and potential at the particles position for particles
! !> between positions nstart and nend in the list
! !
! !> @param[in] p the particles type storing the desired specie parameters
! !> @param[in] nstart starting index for the particle list
! !> @param[in] nend ending index for the particle list
! !---------------------------------------------------------------------------
! SUBROUTINE Efieldscompatparts(p, nstart, nend)
! USE basic, ONLY: thetgrid, invdthet
! USE beam, ONLY: particles
! USE geometry
! USE splinebound
! TYPE(particles), INTENT(INOUT) :: p
! INTEGER, OPTIONAL ::nstart, nend
! INTEGER :: i, iend, k, nst, nnd, nbunch
! INTEGER :: thetindex
! REAL(kind=db) :: wthet, pot1, pot2, Er1, Er2, Ez1, Ez2
! REAL(kind=db), ALLOCATABLE :: gtildeloc(:,:)
! IF(.not. present(nstart)) nst = 1
! IF(.not. present(nend)) nnd = p%Nploc
! nbunch = 64 ! Particle bunch size
! ALLOCATE(gtildeloc(0:2,0:nbunch - 1))
! ! Evaluate the electric potential and field at the particles position
! !$OMP DO SIMD
! DO i=nst,nnd,nbunch
! ! Avoid segmentation fault by accessing non relevant data
! iend=min(i+nbunch-1,nnd)
! 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
! p%pot(i+k-1)=(1-wthet)*pot1+wthet*pot2
! p%E(1,i+k-1)=(1-wthet)*Er1+wthet*Er2
! p%E(2,i+k-1)=0
! p%E(3,i+k-1)=(1-wthet)*Ez1+wthet*Ez2
! END DO
! 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%potxt(i:iend))
! 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)
! p%potxt(i:iend) = p%geomweight(0,i:iend)*p%potxt(i:iend) + gtildeloc(0,0:iend - i)
! END DO
! !$OMP END DO SIMD NOWAIT
! END SUBROUTINE Efieldscompatparts
!---------------------------------------------------------------------------
!> @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
c4science · Help