diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index 6a6c93b..b451cf4 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,1318 +1,1318 @@ !------------------------------------------------------------------------------ ! 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 mag_init 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 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, ONLY: pot, nlperiod, nrank, rhs, 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 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, volume, rgrid implicit none ! set up the geometry module for setting up non-conforming boundary conditions call timera(0, "geom_init") call geom_init(splrz, vec1, vec2) call timera(1, "geom_init") ! 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 Call comp_gradgtilde 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, 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 (plist(i)%is_field) 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) 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(.not. allocated(p%moments))THEN if(mpirank.eq.0)THEN Allocate(p%moments(nbmoments,nrank(1)*nrank(2))) else Allocate(p%moments(0,0)) end if end if ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL moments_gather(p%moments) ELSE p%moments = 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,i) 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(1:iend-i+1), erext(1:iend-i+1)) !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(:, :) 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 :: 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(j,i,xgauss,wgauss,gausssize,wgeom, igauss, iterm,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, iderw, idert, iderg, coefs, fun, fun2) 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(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: 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 if (walltype .lt. 0) fverif = 0 ! Assemble Volume matrix !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, gtildeintegr, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib) 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 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 !> 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(:, :), 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_gradgtilde") + !call timera(0, "comp_gradgtilde") 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 gradgtilde = 0 ! Assemble Volume matrix !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, gtildeintegr, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib) 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 = 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_gradgtilde") + !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, 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(:) 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(:) 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(mat, xi, yj, xcoeff, ycoeff) REAL(kind=db):: interp2 Real(kind=db), INTENT(IN):: mat(:, :), xcoeff, ycoeff INTEGER, INTENT(IN):: xi, yj interp2 = mat(xi, yj)*xcoeff*ycoeff + mat(xi + 1, yj)*(1 - xcoeff)*ycoeff & & + mat(xi, yj + 1)*xcoeff*(1 - ycoeff) + mat(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 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 (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 diff --git a/src/psupply_mod.f90 b/src/psupply_mod.f90 index 8727d76..7989750 100644 --- a/src/psupply_mod.f90 +++ b/src/psupply_mod.f90 @@ -1,309 +1,309 @@ module psupply use constants implicit none type power_supply logical :: active = .false. ! is the power supply active real(kind=db):: geomcapacitor = 1 ! capacitance of the metalic vessel normalised to the neutral density in vessel used in the neutcol module real(kind=db):: PSresistor = 1 ! internal resistance of the power supply normalised to the actual neutral density in vessel real(kind=db):: targetbias = 0 ! Set voltage on the power supply real(kind=db):: bias = 0 ! current voltage on the power supply integer :: nbbounds = 2 ! number of boundaries defined in the geometry integer :: lststp = 0 ! previous step on which the bias was updated real(kind=db):: current(3) = 0 ! current collected on the boundaries normalised to the simulated collision neutral density set in neutcol module ! 1 is at time i-2nbhdt, 2 is at time i-nbhdt and 3 is at time i integer, allocatable:: bdpos(:) ! sign of each boundary for collected charge to determine direction of current real(kind=db),allocatable:: charge(:) ! Charge collected on each boundary and per nbdt real(kind=db),allocatable:: biases(:) ! Actual potentials at each boundary integer :: nbhdt = 10 ! half of the number of time steps between each calls to RK4 real(kind=db):: expdens ! [m-3] experimental neutral density real(kind=db):: neutcoldens ! [m-3] neutral density used in neutcol module end type type(power_supply):: the_ps contains ! read the input parameters from the input file and setup the necesary variables for the module to work subroutine psupply_init(fileid,cstep,nbbounds,neutcoldens,rstbias) use splinebound use basic, only: phinorm, tnorm, mpirank, qnorm, potinn, potout use constants use mpi use geometry use weighttypes use fields integer:: fileid, cstep, nbbounds, istat, nbhdt, ierr,i real(kind=db),OPTIONAL, INTENT(IN):: rstbias real(kind=db):: neutcoldens ! [m-3] real(kind=db):: expneutdens = 1 ! [m-3] real(kind=db):: PsResistor = 1 ! [Ohm] real(kind=db):: geomcapacitor = 1 ! [F] real(kind=db):: targetbias = 0 ! [V] integer, allocatable:: bdpos(:) character(len=1000) :: line logical :: active = .false. NAMELIST /psupplyparams/ expneutdens, PsResistor, geomcapacitor, targetbias, nbhdt, active, bdpos the_ps%lststp=cstep the_ps%nbbounds=nbbounds allocate(the_ps%bdpos(nbbounds),bdpos(nbbounds)) allocate(the_ps%charge(nbbounds),the_ps%biases(nbbounds)) the_ps%bdpos=0 bdpos=0 the_ps%charge=0 the_ps%biases=0 the_ps%current=0 ! read the input parameters from file Rewind(fileid) READ(fileid, psupplyparams, iostat=istat) if (istat.gt.0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in pssupplyparams: '//trim(line) call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if ! save the parameters on output IF(mpirank .eq. 0) THEN WRITE(*, psupplyparams) END IF ! rescale the targetbias set on the power supply the_ps%targetbias=abs(targetbias)/phinorm the_ps%bdpos=bdpos ! save the experimental neutral density the_ps%expdens=expneutdens ! save the neutral collision density the_ps%neutcoldens=neutcoldens if(present(rstbias))then ! Initialize the current bias from the restart value the_ps%bias=rstbias else ! initialize with the file input parameters if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds if(the_ps%bdpos(i) .lt. 0)then the_ps%bias=-the_domain%boundaries(i)%Dirichlet_val exit end if end do else the_ps%bias=(potout-potinn) end if end if ! set the initial bias where(the_ps%bdpos .lt. 0) the_ps%biases=the_ps%bias*the_ps%bdpos end where ! Normalise resistor and capacitor to adapt to experimental pressure the_ps%PSresistor = PSresistor*the_ps%expdens/the_ps%neutcoldens*qnorm/(tnorm*phinorm) the_ps%geomcapacitor = geomcapacitor*phinorm/qnorm the_ps%nbhdt = nbhdt the_ps%active = active if( .not. the_ps%active) return ! Initialize the biases if (the_domain%nbsplines.gt.0) then do i=1,the_ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=the_ps%biases(i) end do else potinn=the_ps%biases(1) Potout=0 end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwgeom) call comp_gradgtilde end subroutine ! save to the result file the parameters of this module read from the input Subroutine psupply_diag(File_handle, str) use mpi Use futils use basic, only: tnorm, phinorm, qnorm implicit none Integer:: File_handle Character(len=*):: str CHARACTER(len=256):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0 .and. the_ps%active) THEN Write(grpname,'(a,a)') trim(str),"/psupply" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "expdens", the_ps%expdens) Call attach(File_handle, trim(grpname), "targetbias", the_ps%targetbias*phinorm) Call attach(File_handle, trim(grpname), "PSresistor", the_ps%PSresistor/the_ps%expdens*the_ps%neutcoldens/qnorm*(tnorm*phinorm)) Call attach(File_handle, trim(grpname), "geomcapacitor", the_ps%geomcapacitor/phinorm*qnorm) Call attach(File_handle, trim(grpname), "nbhdt", the_ps%nbhdt) Call putarr(File_handle,trim(grpname)//"/bdpos", the_ps%bdpos) END IF End subroutine psupply_diag ! gneral routine called from stepon to update the psupply bias subroutine psupply_step(ps,p,cstep) use particletypes use geometry use weighttypes use fields use basic, only: Potinn, potout type(power_supply):: ps type(particles):: p(:) integer:: cstep, i if (.not. ps%active ) return ! calculate the charge collected on each boundary due to the contribution of each specie call add_charge(ps,p,cstep) ! calculate the current flowing between the electrodes due to the cloud call calc_current(ps,cstep) ! calculate the bias at the new time step call updt_bias(ps,cstep) if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! update the bias on the geometry for the Dirichlet b.c. if (the_domain%nbsplines.gt.0) then do i=1,ps%nbbounds the_domain%boundaries(i)%Dirichlet_val=ps%biases(i) end do else potinn=ps%biases(1) Potout=0 end if ! recalculate gtilde to adapt for the new biases CALL total_gtilde(vec1, vec2, gtilde, gridwgeom) call comp_gradgtilde end subroutine ! calculates the current flowing between the electrodes due to the cloud subroutine calc_current(ps,cstep) use geometry use basic, only: phinorm, dt use fields type(power_supply):: ps integer:: cstep if(mod(cstep-ps%lststp,ps%nbhdt).eq.0) then ! communicate the charge accumulation in this timestep call reduce_charge(ps) if(mod(cstep-ps%lststp,ps%nbhdt*2).eq.0)then ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(3)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if ps%lststp=cstep else ! calculates the current by adding the contribution of each boundary if (mpirank .eq. 0)then ps%current(2)=sum(-ps%charge*ps%bdpos)/(ps%nbhdt*dt) end if end if ps%charge=0 end if end subroutine ! calculate the charge deposited by each specie on the electrodes (used to calculate the resulting current) subroutine add_charge(ps,p, cstep) use particletypes use basic, only: qnorm type(power_supply):: ps type(particles):: p(:) integer::cstep, i do i=1,size(p,1) if(.not. p(i)%is_field) cycle !Add the normalised contribution of each specie ps%charge=ps%charge+p(i)%nblost(5:)*p(i)%weight*p(i)%q/qnorm end do end subroutine ! Time integrate the ODE of the actual bias between the accelerating electrodes ! and broadcast it to all the workers subroutine updt_bias(ps,cstep) use basic, only: dt, phinorm implicit none type(power_supply):: ps integer:: cstep real(kind=db):: bias,k1,k2,k3,k4, hdeltat if(mod(cstep-ps%lststp,2*ps%nbhdt) .ne. 0) return ! half delta t hdeltat=dt*ps%nbhdt/(ps%PSresistor*ps%geomcapacitor) bias=ps%bias ! we update the bias using RK4 k1=-(bias+ps%current(1)*ps%PSresistor-ps%targetbias) k2=-(bias+hdeltat*k1+ps%current(2)*ps%PSresistor-ps%targetbias) k3=-(bias+hdeltat*k2+ps%current(2)*ps%PSresistor-ps%targetbias) k4=-(bias+2*hdeltat*k3+ps%current(3)*ps%PSresistor-ps%targetbias) ps%bias=bias+(k1+2*k2+2*k3+k4)*2*hdeltat/6 - Write(*,*) " new bias ", ps%bias*phinorm + !Write(*,*) " new bias ", ps%bias*phinorm where (ps%bdpos .lt. 0) ps%biases=-ps%bias end where ! broadcast the bias to all the mpi processes call bcast_bias(ps) ps%current(1)=ps%current(3) end subroutine updt_bias ! gather on node 0 the collected charge on each metallic boundary subroutine reduce_charge(ps) use mpi use mpihelper use basic, ONLY: mpirank type(power_supply):: ps integer:: ierr if(mpirank .eq. 0) then call MPI_REDUCE(MPI_IN_PLACE,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) - Write(*,*) "curr charge ", ps%charge + !Write(*,*) "curr charge ", ps%charge else call MPI_REDUCE(ps%charge,ps%charge,ps%nbbounds,db_type,db_sum_op,0,MPI_COMM_WORLD,ierr) end if end subroutine ! broadcast to all the nodes the new bias imposed by the power supply on the electrodes subroutine bcast_bias(ps) use mpi use mpihelper type(power_supply):: ps integer:: ierr call MPI_BCAST(ps%biases,ps%nbbounds,db_type,0,MPI_COMM_WORLD,ierr) end subroutine end module psupply