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

fennecs/src/fields_mod.f90
e4669700ed0dgeom_testmumps

Last Change
Raw File

fields_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: beam
!
!> @author
!> Patryk Kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for initializing the magnetic field, solving the Poisson equation and computing the moments of the particles distribution function
!------------------------------------------------------------------------------
MODULE fields
USE constants
USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, 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
IMPLICIT NONE
REAL(kind=db), allocatable, SAVE :: matcoef(:, :), phi_spline(:), vec1(:), vec2(:)
REAL(kind=db), allocatable, SAVE :: loc_moments(:, :), loc_rhs(:), gradgtilde(:), fverif(:)
INTEGER, SAVE:: loc_zspan
TYPE(mumps_mat), SAVE :: femat !< Finite Element Method matrix
TYPE(mumps_mat), SAVE :: reduccedmat !< Finite Element Method matrix
REAL(kind=db), SAVE:: potextA, potextB !< Variables used for fast calculation of the external potential
INTEGER :: nbmoments = 10 !< number of moments to be calculated and stored
INTEGER(kind=omp_lock_kind), Allocatable:: mu_lock(:) !< Stores the lock for fields parallelism
CONTAINS
SUBROUTINE init_mag
USE basic, ONLY: magnetfile, nr, nz
USE bsplines
USE mumps_bsplines
USE mpihelper
USE geometry
ALLOCATE (Br((nr + 1)*(nz + 1)), Bz((nr + 1)*(nz + 1)))
ALLOCATE (Athet((nr + 1)*(nz + 1)))
! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed)
! or load it from magnetfile if present
CALL magnet(magnetfile)
end subroutine init_mag
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid
!
!---------------------------------------------------------------------------
SUBROUTINE fields_init
USE basic, ONLY: pot, nlperiod, nrank, rhs, moments, volume, rgrid
USE bsplines
USE geometry
USE mumps_bsplines
USE mpihelper
INTEGER :: nrz(2), i
! Set up 2d spline splrz used in the FEM
CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod)
! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential
CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz_ext, nlppform=nlppform, period=nlperiod)
! Calculate dimension of splines
nrz(1) = nz
nrz(2) = nr
CALL get_dim(splrz, nrank, nrz, femorder)
! Allocate necessary variables
ALLOCATE (matcoef(nrank(1), nrank(2)))
ALLOCATE (pot((nr + 1)*(nz + 1)))
ALLOCATE (potxt((nr + 1)*(nz + 1)))
ALLOCATE (Erxt((nr + 1)*(nz + 1)))
ALLOCATE (Ezxt((nr + 1)*(nz + 1)))
ALLOCATE (rhs(nrank(1)*nrank(2)))
ALLOCATE (gradgtilde(nrank(1)*nrank(2)))
gradgtilde = 0
IF (mpirank .eq. 0) THEN
ALLOCATE (moments(nbmoments, nrank(1)*nrank(2)))
ELSE
ALLOCATE (moments(0, 0))
END IF
ALLOCATE (phi_spline(nrank(1)*nrank(2)))
ALLOCATE (volume(nrank(1)*nrank(2)))
volume = 0
ALLOCATE (Er((nr + 1)*(nz + 1)), Ez((nr + 1)*(nz + 1)))
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
SUBROUTINE fields_start
USE geometry
USE basic, ONLY: pot, nlperiod, nrank, rhs, moments, volume, rgrid
implicit none
! set up the geometry module for setting up non-conforming boundary conditions
call timera(0, "init_geom")
call init_geom(splrz, vec1, vec2)
call timera(1, "init_geom")
! Calculate and factorise FEM matrix (depends only on mesh)
CALL fematrix
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
if (nlweb) then
! Calculate reduced matrix for use of web splines
call timera(0, "reduce femat")
call Reducematrix(femat, reduccedmat)
call factor(reduccedmat)
call timera(1, "reduce femat")
else
CALL factor(femat)
end if
! Computes the externally imposed electric field
rhs = -gradgtilde
if (walltype .lt. 0) rhs = rhs + fverif
call poisson(splrz_ext)
rhs = 0
potxt = pot
erxt = Er
Ezxt = Ez
END SUBROUTINE fields_start
!---------------------------------------------------------------------------
!> @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
USE mpihelper
INTEGER:: Zbounds(0:)
loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(1)
if (allocated(loc_moments)) deallocate (loc_moments)
ALLOCATE (loc_moments(nbmoments, loc_zspan*nrank(2)))
if (allocated(loc_rhs)) deallocate (loc_rhs)
ALLOCATE (loc_rhs(loc_zspan*nrank(2)))
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: rnorm, rhs, moments, nlend, step, it2d
USE beam, ONLY: particles
USE mpihelper
Use geometry
type(particles), INTENT(INOUT):: plist(:)
INTEGER:: i
IF (nlphis) then ! We calculate the self-consistent field
loc_rhs = 0 ! Reset the moments matrix
! Assemble rhs for each specie
Do i = 1, size(plist, 1)
if (.not. plist(i)%is_test) CALL deposit_charge(plist(i), loc_rhs)
END Do
! If we are using MPI parallelism, reduce the rhs on the root process
IF (mpisize .gt. 1) THEN
CALL rhs_gather(rhs)
ELSE
rhs = loc_rhs
END IF
ELSE ! We only consider the externally imposed field
rhs = 0
END IF
IF (mpirank .eq. 0) THEN
rhs = rhs - gradgtilde
if (walltype .lt. 0) rhs = rhs + fverif
END IF
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, moment)
USE bsplines
USE mpi
USE basic, ONLY: rnorm, nlend, step, it2d
USE beam, ONLY: particles
USE mpihelper
Use geometry
type(particles), INTENT(INOUT):: p
REAL(kind=db), INTENT(INOUT):: moment(:, :)
loc_moments = 0 ! Reset the moments matrix
! Assemble rhs
IF (p%Nploc .ne. 0) THEN
CALL deposit_moments(p, loc_moments)
END IF
! If we are using MPI parallelism, reduce the rhs on the root process
IF (mpisize .gt. 1) THEN
CALL moments_gather(moment)
ELSE
moment = loc_moments
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
!> @param[in] p_loc_moments local tensor used to store the moments of the given specie
!---------------------------------------------------------------------------
SUBROUTINE deposit_moments(p, p_loc_moments)
USE bsplines
USE mpi
USE basic, ONLY: Zbounds
USE beam, ONLY: particles
USE mpihelper
USE geometry
USE omp_lib
TYPE(particles), INTENT(IN):: p
REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: p_loc_moments
INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch
INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft
REAL(kind=db) :: vr, vthet, vz, coeff
REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :), wgeom(:, :)
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
ALLOCATE (zleft(nbunch), rleft(nbunch))
ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node
ALLOCATE (wgeom(nbunch, 0:0))
! Assemble rhs
IF (p%Nploc .ne. 0) THEN
!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, iend, irow, jcol, mu, k, vr, vz, vthet, coeff, fun, fun2)
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
!CALL locintv(splrz%sp2, p%R(i:iend), rleft(1:k))
!CALL locintv(splrz%sp1, p%Z(i:iend), zleft(1:k))
rleft(1:k) = p%rindex(i:iend)
zleft(1:k) = p%zindex(i:iend)
! Compute the value of the splines at the particles positions
CALL basfun(p%Z(i:iend), splrz%sp1, fun(:, :, 1:k), zleft(1:k) + 1)
CALL basfun(p%R(i:iend), splrz%sp2, fun2(:, :, 1:k), rleft(1:k) + 1)
!CALL geom_weight(p%Z(i:iend),p%R(i:iend),wgeom)
DO k = 1, (iend - i + 1)
DO jw = 1, (femorder(2) + 1)
DO it = 1, (femorder(1) + 1)
irow = zleft(k) + it - Zbounds(mpirank)
jcol = rleft(k) + jw
mu = irow + (jcol - 1)*(loc_zspan)
coeff = p%weight*fun(it, 0, k)*fun2(jw, 0, k)
! Add contribution of particle nbunch to rhs grid point mu
vr = 0.5*(p%UR(i + k - 1)/p%Gamma(i + k - 1) + p%URold(i + k - 1)/p%Gammaold(i + k - 1))
vz = 0.5*(p%UZ(i + k - 1)/p%Gamma(i + k - 1) + p%UZold(i + k - 1)/p%Gammaold(i + k - 1))
vthet = 0.5*(p%UTHET(i + k - 1)/p%Gamma(i + k - 1) + p%UTHETold(i + k - 1)/p%Gammaold(i + k - 1))
call omp_set_lock(mu_lock(mu))
!!$OMP ATOMIC UPDATE
p_loc_moments(1, mu) = p_loc_moments(1, mu) + coeff
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(2, mu) = p_loc_moments(2, mu) + coeff*vr
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(3, mu) = p_loc_moments(3, mu) + coeff*vthet
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(4, mu) = p_loc_moments(4, mu) + coeff*vz
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(5, mu) = p_loc_moments(5, mu) + coeff*vr*vr
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(6, mu) = p_loc_moments(6, mu) + coeff*vr*vthet
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(7, mu) = p_loc_moments(7, mu) + coeff*vr*vz
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(8, mu) = p_loc_moments(8, mu) + coeff*vthet*vthet
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(9, mu) = p_loc_moments(9, mu) + coeff*vthet*vz
!!$OMP END ATOMIC
!!$OMP ATOMIC UPDATE
p_loc_moments(10, mu) = p_loc_moments(10, mu) + coeff*vz*vz
!!$OMP END ATOMIC
call omp_unset_lock(mu_lock(mu))
END DO
END DO
END DO
END DO
!$OMP END PARALLEL DO
END IF
DEALLOCATE (fun, fun2, zleft, rleft)
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
!> @param[in] p_loc_moments local tensor used to store the moments of the given specie
!---------------------------------------------------------------------------
SUBROUTINE deposit_charge(p, p_loc_moments)
USE bsplines
USE mpi
USE constants
USE basic, ONLY: Zbounds, rnorm, phinorm
USE beam, ONLY: particles
USE mpihelper
USE geometry
USE omp_lib
TYPE(particles), INTENT(IN):: p
REAL(kind=db), DIMENSION(:), INTENT(INOUT):: p_loc_moments
INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch
INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft
REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :)
INTEGER:: num_threads
real(kind=db):: contrib, chargecoeff
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
chargecoeff = p%weight*p%q/(2*pi*eps_0*phinorm*rnorm) ! Normalized charge density simulated by each macro particle
! Assemble rhs
IF (p%Nploc .ne. 0) THEN
!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(zleft, rleft, jw, it, iend, irow, jcol, mu, k, fun, fun2, contrib)
ALLOCATE (zleft(nbunch), rleft(nbunch))
ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node
!$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
rleft(1:k) = p%rindex(i:iend)
zleft(1:k) = p%zindex(i:iend)
! Compute the value of the splines at the particles positions
CALL basfun(p%Z(i:iend), splrz%sp1, fun, zleft(1:k) + 1)
CALL basfun(p%R(i:iend), splrz%sp2, fun2, rleft(1:k) + 1)
!CALL geom_weight(p%Z(i:iend),p%R(i:iend),wgeom)
DO k = 1, (iend - i + 1)
DO jw = 1, (femorder(2) + 1)
DO it = 1, (femorder(1) + 1)
irow = zleft(k) + it - Zbounds(mpirank)
jcol = rleft(k) + jw
mu = irow + (jcol - 1)*(loc_zspan)
! Add contribution of particle k to rhs grid point mu
contrib = fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(i + k - 1, 0)*chargecoeff
!$OMP ATOMIC UPDATE
p_loc_moments(mu) = p_loc_moments(mu) + contrib
!$OMP END ATOMIC
END DO
END DO
END DO
END DO
!$OMP END DO
DEALLOCATE (fun, fun2, zleft, rleft)
!$OMP END PARALLEL
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
!> and reduce the result on the host
!
!---------------------------------------------------------------------------
SUBROUTINE rhs_gather(rhs)
USE mpihelper
USE Basic, ONLY: Zbounds, mpirank, nlend, it2d, step, 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)
counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1)
counts(mpisize) = counts(mpisize) + femorder(1)
CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs, nrank, femorder, loc_zspan - femorder(1))
IF (mpirank .gt. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i)
DO j = 1, femorder(1)
DO i = 1, nrank(2)
loc_rhs((i - 1)*loc_zspan + j) = loc_rhs((i - 1)*loc_zspan + j)&
& + rhsoverlap_buffer(nrank(2)*(j - 1) + i)
END DO
END DO
!$OMP END PARALLEL DO SIMD
END IF
! Set communication vector type
overlap_type = rhsoverlap_type
rcvoverlap_type = rcvrhsoverlap_type
IF (mpirank .eq. 0) THEN
rhs = 0
END IF
CALL MPI_GATHERV(loc_rhs, counts(mpirank + 1), overlap_type, &
& rhs, counts, displs, rcvoverlap_type, 0, MPI_COMM_WORLD, ierr)
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, nlend, it2d, step, leftproc, rightproc
REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: moment
INTEGER:: ierr, i, j
INTEGER:: displs(mpisize), counts(mpisize)
displs = Zbounds(0:mpisize - 1)
counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1)
counts(mpisize) = counts(mpisize) + femorder(1)
CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan - femorder(1))
IF (mpirank .gt. 0) THEN
!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i)
DO j = 1, femorder(1)
DO i = 1, nrank(2)
loc_moments(1:nbmoments, (i - 1)*loc_zspan + j) = loc_moments(1:nbmoments, (i - 1)*loc_zspan + j)&
& + momentsoverlap_buffer(nbmoments*(nrank(2)*(j - 1) + i - 1) + 1:nbmoments*(nrank(2)*(j - 1) + i))
END DO
END DO
!$OMP END PARALLEL DO SIMD
END IF
! Set communication vector type
IF (mpirank .eq. 0) THEN
moment = 0
END IF
CALL MPI_GATHERV(loc_moments, counts(mpirank + 1), momentsoverlap_type, &
& moment, counts, displs, rcvmomentsoverlap_type, 0, MPI_COMM_WORLD, ierr)
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 and interpolate the electric forces
!> for each particle.
!
!---------------------------------------------------------------------------
SUBROUTINE poisson(splinevar)
USE basic, ONLY: rhs, nrank, pot, nlend
USE bsplines, ONLY: spline2d, gridval
USE mumps_bsplines, ONLY: bsolve, vmx
USE futils
Use geometry
type(spline2d):: splinevar
INTEGER:: ierr, jder(2)
real(kind=db), allocatable::reducedrhs(:)
real(kind=db), allocatable:: reducedsol(:), tempcol(:)
jder(1) = 0
jder(2) = 0
! Only the root process solves Poisson
IF (mpirank .eq. 0) then
if (nlweb) then ! we use the web-spline reduction for stability
allocate (reducedrhs(nrank(1)*nrank(2)))
allocate (reducedsol(nbreducedspline), tempcol(nrank(1)*nrank(2)))
reducedrhs = vmx(etilde, rhs)
Call bsolve(reduccedmat, reducedrhs(1:nbreducedspline), reducedsol)
tempcol = 0
tempcol(1:nbreducedspline) = reducedsol
!phi_spline = 0
phi_spline = vmx(etildet, tempcol)
else
CALL bsolve(femat, rhs, phi_spline)
end if
end if
! Broadcast the solution to all MPI workers
CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr)
matcoef = reshape(phi_spline, (/nrank(1), nrank(2)/))
! Evaluate the electric potential on the grid points
CALL gridval(splinevar, vec1(1:2), vec2(1:2), pot(1:4), jder, matcoef)
IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN
! On the root process, compute the electric field for diagnostic purposes
CALL gridval(splinevar, vec1, vec2, pot, (/0, 0/))
CALL gridval(splinevar, vec1, vec2, Ez, (/1, 0/))
CALL gridval(splinevar, vec1, vec2, Er, (/0, 1/))
Ez = -pot*gridwgeom(:, 1) - Ez*gridwgeom(:, 0) - gtilde(:, 1)
Er = -pot*gridwgeom(:, 2) - Er*gridwgeom(:, 0) - gtilde(:, 2)
pot = pot*gridwgeom(:, 0) + gtilde(:, 0)
END IF
END SUBROUTINE poisson
!---------------------------------------------------------------------------
!> @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 beam, ONLY: particles
Use geometry
TYPE(particles), INTENT(INOUT):: p
INTEGER, OPTIONAL::nstart, nend
INTEGER:: i, iend, nst, nnd
INTEGER:: nbunch
INTEGER:: num_threads
Real(kind=db), ALLOCATABLE:: erext(:), ezext(:), gtildeloc(:, :)
if (.not. present(nstart)) nst = 1
if (.not. present(nend)) nnd = p%Nploc
num_threads = omp_get_max_threads()
nbunch = (nnd - nst + 1)/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
Allocate (erext(nbunch), ezext(nbunch), gtildeloc(0:nbunch - 1, 0:2))
! Evaluate the electric potential and field at the particles position
!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(iend) firstprivate(erext,ezext,gtildeloc)
DO i = nst, nnd, nbunch
! Avoid segmentation fault by accessing non relevant data
iend = min(i + nbunch - 1, nnd)
CALL getgrad(splrz, p%Z(i:iend), p%R(i:iend), p%pot(i:iend), p%EZ(i:iend), p%ER(i:iend))
CALL getgrad(splrz_ext, p%Z(i:iend), p%R(i:iend), p%potxt(i:iend), EZext(:), erext(:))
!Call geom_weight(p%Z(i:iend), p%R(i:iend), wgeom(0:iend-i,:))
Call total_gtilde(p%Z(i:iend), p%R(i:iend), gtildeloc(0:iend - i, :),p%geomweight(i:iend,:))
p%EZ(i:iend) = -p%Ez(i:iend)*p%geomweight(i:iend, 0) - p%pot(i:iend)*p%geomweight(i:iend, 1) - gtildeloc(0:iend - i, 1)
p%ER(i:iend) = -p%ER(i:iend)*p%geomweight(i:iend, 0) - p%pot(i:iend)*p%geomweight(i:iend, 2) - gtildeloc(0:iend - i, 2)
p%pot(i:iend) = p%geomweight(i:iend, 0)*p%pot(i:iend) + gtildeloc(0:iend - i, 0)
p%potxt(i:iend) = p%geomweight(i:iend, 0)*p%potxt(i:iend) + gtildeloc(0:iend - i, 0)
END DO
!$OMP END PARALLEL DO
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
USE bsplines
USE geometry
USE omp_lib
USE sparse
REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :)
REAL(kind=db), ALLOCATABLE :: f(:, :), aux(:)
REAL(kind=db), ALLOCATABLE :: coefs(:)
REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :)
REAL(kind=db) :: contrib
INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :)
INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize
kterms=8
ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node
!ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays
ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines
ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2))
!Pointers on the order of derivatives
call timera(0, "fematrix")
! Constuction of auxiliary array ordering bsplines in given interval
DO i = 1, (femorder(1) + 1)
aux(i) = i
END DO
DO i = 1, (femorder(2) + 1)
f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux
f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i
END DO
! Initialisation of FEM matrix
CALL init(nrank(1)*nrank(2), 2, femat)
! Assemble FEM matrix
!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(i,xgauss,wgauss,gausssize,wgeom, igauss, iterm,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, iderw, idert, iderg, coefs, fun, fun2), firstprivate(ngauss)
DO j = 1, nr ! Loop on r position
DO i = 1, nz ! Loop on z position
!! Computation of gauss weight and position in r and z direction for gaussian integration
Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize)
if (gausssize .gt. 0) then
If (allocated(wgeom)) deallocate (wgeom)
ALLOCATE (wgeom(gausssize, 0:2))
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%sp1, fun, i)
CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j)
CALL coefeq(xgauss(igauss, :), idert, iderw, iderg, coefs, kterms)
DO iterm = 1, kterms ! Loop on the two integration dimensions
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 = wgeom(igauss, iderg(iterm, 1))*wgeom(igauss, iderg(iterm, 2))* &
& fun(f(jt, 1), idert(iterm, 1))*fun(f(iw, 1), idert(iterm, 2))* &
& fun2(f(jt, 2), iderw(iterm, 1))*fun2(f(iw, 2), iderw(iterm, 2))* &
& wgauss(igauss)*coefs(iterm)
call omp_set_lock(mu_lock(mu))
CALL updt_sploc(femat%mat%row(mu), mu2, contrib)
call omp_unset_lock(mu_lock(mu))
END DO
END DO
END DO
END DO
END DO
END DO
!$OMP End parallel do
DEALLOCATE (f, aux)
DEALLOCATE (idert, iderw, coefs, fun, fun2)
call timera(1, "fematrix")
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
REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :)
REAL(kind=db), ALLOCATABLE :: f(:, :), aux(:), coefs(:)
REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), gtildeintegr(:, :), 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(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays
ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines
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
gradgtilde = 0
if (walltype .lt. 0) fverif = 0
! Assemble Volume matrix
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,xgauss,wgauss,gausssize,wgeom, igauss, gtildeintegr, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib), firstprivate(ngauss)
DO j = 1, nr ! Loop on r position
DO i = 1, nz ! Loop on z position
! Computation of gauss weight and position in r and z direction for gaussian integration
Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize)
If (allocated(wgeom)) deallocate (wgeom)
if (gausssize .gt. 0) then
ALLOCATE (wgeom(size(xgauss, 1), 0:2))
CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom)
End if
If (allocated(gtildeintegr)) deallocate (gtildeintegr)
ALLOCATE (gtildeintegr(size(xgauss, 1), 0:2))
Call total_gtilde(xgauss(:, 1), xgauss(:, 2), gtildeintegr,wgeom)
if (walltype .lt. 0) then
If (allocated(ftestpt)) deallocate (ftestpt)
ALLOCATE (ftestpt(size(xgauss, 1), 0:0))
CALL ftest(xgauss(:, 1), xgauss(:, 2), ftestpt)
end if
DO igauss = 1, gausssize ! Loop on gaussian weights and positions
CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i)
CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j)
CALL coefeqext(xgauss(igauss, :), idt, idw, idg, idp, coefs)
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, 2)!*wgeom(igauss,0)
!$OMP ATOMIC UPDATE
volume(mu) = volume(mu) + newcontrib
!$OMP END ATOMIC
if (walltype .lt. 0) THEN
newcontrib = ftestpt(igauss, 0)*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)&
&*wgeom(igauss, 0)*wgauss(igauss)*xgauss(igauss, 2)
!$OMP ATOMIC UPDATE
fverif(mu) = fverif(mu) + newcontrib
!$OMP END ATOMIC
end if
newcontrib = 0
Do iterm = 1, nterms
newcontrib = newcontrib + wgeom(igauss, idg(iterm))*gtildeintegr(igauss, idp(iterm))* &
& fun(f(jt, 1), idt(iterm))*fun2(f(jt, 2), idw(iterm))* &
& wgauss(igauss)*coefs(iterm)
End do
!$OMP ATOMIC UPDATE
gradgtilde(mu) = gradgtilde(mu) + newcontrib
!$OMP END ATOMIC
END DO
END DO
END DO
END DO
!$OMP END PARALLEL DO
!DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg)
DEALLOCATE (f, aux)
DEALLOCATE (fun, fun2)
call timera(1, "comp_volume")
END SUBROUTINE comp_volume
!---------------------------------------------------------------------------
!> @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, i, arr)
END IF
arr = 0; arr(nrank(1)*nrank(2) + 1 - i) = 1;
CALL putrow(femat, nrank(1)*nrank(2) + 1 - i, arr)
END DO
DEALLOCATE (arr)
END SUBROUTINE fe_dirichlet
!________________________________________________________________________________
SUBROUTINE coefeq(x, idt, idw, idg, c, kterms)
REAL(kind=db), INTENT(in) :: x(2)
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) = 1
idw(2, 1) = 0
idw(2, 2) = 0
idg(2, 1) = 1
idg(2, 2) = 0
idt(3, 1) = 1
idt(3, 2) = 0
idw(3, 1) = 0
idw(3, 2) = 0
idg(3, 1) = 0
idg(3, 2) = 1
idt(4, 1) = 1
idt(4, 2) = 1
idw(4, 1) = 0
idw(4, 2) = 0
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) = 0
idw(6, 1) = 0
idw(6, 2) = 1
idg(6, 1) = 2
idg(6, 2) = 0
idt(7, 1) = 0
idt(7, 2) = 0
idw(7, 1) = 1
idw(7, 2) = 0
idg(7, 1) = 0
idg(7, 2) = 2
idt(8, 1) = 0
idt(8, 2) = 0
idw(8, 1) = 1
idw(8, 2) = 1
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(2)
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) = 1
idw(2) = 0
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) = 0
idw(4) = 1
END SUBROUTINE coefeqext
!---------------------------------------------------------------------------
!> @author
!> Patryk kaminski EPFL/SPC
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the magnetic field on the grid according to a magnetic mirror,
!> or according to the linear interpolation of the values on the
!> grid saved in h5 file stored at magfile.
!> @param[in] magfile filname of .h5 file containing the definitions of A and B
!---------------------------------------------------------------------------
SUBROUTINE magnet(magfile)
USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm
USE constants, ONLY: Pi
CHARACTER(LEN=*), INTENT(IN), OPTIONAL:: magfile
REAL(kind=db) :: rg, zg, halfLz, MirrorRatio
INTEGER :: i, rindex
IF (len_trim(magfile) .lt. 1) THEN
halfLz = (zgrid(nz) + zgrid(0))/2
MirrorRatio = (Rcurv - 1)/(Rcurv + 1)
DO i = 1, (nr + 1)*(nz + 1)
rindex = (i - 1)/(nz + 1)
rg = rgrid(rindex)
zg = zgrid(i - rindex*(nz + 1) - 1) - halfLz
Br(i) = -B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm
Bz(i) = B0*(1 - MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm
Athet(i) = 0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm))
END DO
ELSE
CALL load_mag_from_h5(magfile)
END IF
END SUBROUTINE magnet
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Loads the magnetic field defined in the .h5 file at location magfile
!> @param[in] magfile filname of .h5 file containing the definitions of A and B
!---------------------------------------------------------------------------
SUBROUTINE load_mag_from_h5(magfile)
USE basic, ONLY: B0, rgrid, zgrid, rnorm, nr, nz, bnorm, bscaling
USE constants, ONLY: Pi
USE futils
CHARACTER(LEN=*), INTENT(IN):: magfile
REAL(kind=db), ALLOCATABLE :: magr(:), magz(:)
REAL(kind=db), ALLOCATABLE :: tempBr(:, :), tempBz(:, :), tempAthet(:, :)
REAL(kind=db) :: zi, rj, rcoeff, zcoeff, maxB
INTEGER :: i, j, l, m, magfid
LOGICAL:: B_is_saved
INTEGER :: magn(2), magrank
CALL openf(trim(magfile), magfid, 'r', real_prec='d')
CALL getdims(magfid, '/mag/Athet', magrank, magn)
ALLOCATE (magr(magn(2)), magz(magn(1)))
ALLOCATE (tempAthet(magn(1), magn(2)), tempBr(magn(1), magn(2)), tempBz(magn(1), magn(2)))
! Read r and z coordinates for the definition of A_\thet, and B
CALL getarr(magfid, '/mag/r', magr)
CALL getarr(magfid, '/mag/z', magz)
CALL getarr(magfid, '/mag/Athet', tempAthet)
IF (isdataset(magfid, '/mag/Br') .and. isdataset(magfid, '/mag/Bz')) THEN
CALL getarr(magfid, '/mag/Br', tempBr)
CALL getarr(magfid, '/mag/Bz', tempBz)
IF(bscaling .gt. 0) then
maxB=sqrt(maxval(tempBr**2+tempBz**2))
tempBr=tempBr/maxB*B0
tempBz=tempBz/maxB*B0
end if
B_is_saved = .true.
ELSE
B_is_saved = .false.
END IF
! Interpolate the magnetic field and potential at the grid points
m = 1
Do j = 0, nr
rj = rgrid(j)*rnorm
Do while (magr(m) .lt. rj)
m = m + 1
END DO
rcoeff = (rj - magr(m - 1))/(magr(m) - magr(m - 1))
l = 1
Do i = 0, nz
zi = zgrid(i)*rnorm
Do while (magz(l) .lt. zi)
l = l + 1
END DO
zcoeff = (zi - magz(l - 1))/(magz(l) - magz(l - 1))
Athet(i + j*(nz + 1) + 1) = interp2(tempAthet, l, m, zcoeff, rcoeff)
IF (B_is_saved) THEN
Br(i + j*(nz + 1) + 1) = interp2(tempBr, l, m, zcoeff, rcoeff)
Bz(i + j*(nz + 1) + 1) = interp2(tempBz, l, m, zcoeff, rcoeff)
ELSE
Br(i + j*(nz + 1) + 1) = -(tempAthet(l + 1, m + 1) + tempAthet(l + 1, m) - tempAthet(l, m + 1) - tempAthet(l, m))/2 &
& /(magz(l) - magz(l - 1))
Bz(i + j*(nz + 1) + 1) = ((tempAthet(l + 1, m + 1) + tempAthet(l, m + 1))*magr(m + 1) &
& - (tempAthet(l + 1, m) - tempAthet(l, m))*magr(m)) &
& /(magr(m) - magr(m - 1))/(magr(m) + magr(m + 1))
END IF
END DO
END DO
if( bscaling .lt. 0 ) then
maxB = maxval(sqrt(Bz**2 + Br**2))
Bz = Bz/maxB*B0
Br = Br/maxB*B0
end if
! We normalize
Br = Br/bnorm
Bz = Bz/bnorm
CALL closef(magfid)
END SUBROUTINE load_mag_from_h5
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the weighted average of the 4 elements of matrix defined by
!> xi, xi+1, yj, yj+1
!
!> @param[in] matrix The matrix containing the initial values to be interpolated
!> @param[in] xi left reference along the first dimension for the average
!> @param[in] yj left reference along the second dimension for the average
!> @param[in] xcoeff left weight along the first dimension for the average
!> @param[in] ycoeff left weight along the second dimension for the average
!> @param[in] interp2 weighted average value
!---------------------------------------------------------------------------
FUNCTION interp2(matrix, xi, yj, xcoeff, ycoeff)
REAL(kind=db):: interp2
Real(kind=db), INTENT(IN):: matrix(:, :), xcoeff, ycoeff
INTEGER, INTENT(IN):: xi, yj
interp2 = matrix(xi, yj)*xcoeff*ycoeff + matrix(xi + 1, yj)*(1 - xcoeff)*ycoeff &
& + matrix(xi, yj + 1)*xcoeff*(1 - ycoeff) + matrix(xi + 1, yj + 1)*(1 - xcoeff)*(1 - ycoeff)
END function
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the zero order
FUNCTION bessi0(x)
REAL(kind=db) :: bessi0, x
REAL(kind=db) :: ax
REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y
SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9
DATA p1, p2, p3, p4, p5, p6, p7/1.0d0, 3.5156229d0, 3.0899424d0, 1.2067492d0, 0.2659732d0, 0.360768d-1, 0.45813d-2/
DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2, 0.916281d-2, &
& -0.2057706d-1, 0.2635537d-1, -0.1647633d-1, 0.392377d-2/
if (abs(x) .lt. 3.75) then
y = (x/3.75)**2
bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))
else
ax = abs(x)
y = 3.75/ax
bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))))
end if
return
END FUNCTION bessi0
!________________________________________________________________________________
!Modified Bessel functions of the first kind of the first order
FUNCTION bessi1(x)
REAL(kind=db) :: bessi1, x
REAL(kind=db) :: ax
REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y
SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9
DATA p1, p2, p3, p4, p5, p6, p7/0.5d0, 0.87890594d0, 0.51498869d0, 0.15084934d0, 0.2658733d-1, 0.301532d-2, 0.32411d-3/
DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2, -0.1031555d-1, &
& 0.2282967d-1, -0.2895312d-1, 0.1787654d-1, -0.420059d-2/
if (abs(x) .lt. 3.75D0) then
y = (x/3.75D0)**2
bessi1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))))
else
ax = abs(x)
y = 3.75D0/ax
bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9))))))))
if (x .lt. 0.) bessi1 = -bessi1
end if
return
END FUNCTION bessi1
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Computes the value of the external electric field at position R,Z for the case of a coaxial configuration of constant radius
!
!> @param[in] R radial positon
!> @param[in] Z axial position
!> @param[out] external potential at position R,Z
!---------------------------------------------------------------------------
ELEMENTAL FUNCTION ext_potential(R)
REAL(kind=db), INTENT(IN):: R
REAL(kind=db):: ext_potential
ext_potential = potextA*log(R) + potextB
END FUNCTION ext_potential
!---------------------------------------------------------------------------
!> @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, moments
INTEGER:: i
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 (loc_rhs)
DEALLOCATE (loc_moments)
DEALLOCATE (phi_spline)
DEALLOCATE (moments)
DEALLOCATE (Br, Bz)
DEALLOCATE (Er, Ez)
DEALLOCATE (vec1, vec2)
Call DESTROY_SP(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
!
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
!+
END MODULE fields
c4science · Help