Subroutine that samples a thermal distribution function of electrons for generating the initial condition of a set of simulated particles.
This subroutine uses the Inverse Transform Sampling Method along with the Metropolis-Hastings algorithm to generate an initial condition of the velocity distribution that follows a 3-dimensional (in velocity space) thermal distribution.
Check that the gyro-distribution is initialized right in this function.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(KORC_PARAMS), | intent(in) | :: | params | Core KORC simulation parameters. |
||
| type(SPECIES), | intent(inout) | :: | spp | An instance of the derived type SPECIES containing all the parameters and simulation variables of the different species in the simulation. |
subroutine thermal_distribution(params,spp)
!! @note Subroutine that samples a thermal distribution function
!! of electrons for generating the initial condition of a set of
!! simulated particles. @endnote
!! This subroutine uses the Inverse Transform Sampling Method along
!! with the Metropolis-Hastings algorithm to generate an
!! initial condition of the velocity distribution that follows a
!! 3-dimensional (in velocity space) thermal distribution.
!! @todo Check that the gyro-distribution is initialized right in
!! this function.
TYPE(KORC_PARAMS), INTENT(IN) :: params
!! Core KORC simulation parameters.
TYPE(SPECIES), INTENT(INOUT) :: spp
!! An instance of the derived type SPECIES containing all the
!! parameters and simulation variables of the different species
!! in the simulation.
REAL(rp) :: Vmax
!! Velocity cutoff where we stop sampling the tail of the thermal
!! distribution.
REAL(rp) :: vsq
!!
REAL(rp) :: Vth
!! Thermal velocity of the sampled distribution \(v_{T_e} =
!! \sqrt{T_e/m_e}\).
REAL(rp) :: sv
!! Step to sample the velocity space using the Metropolis-Hastings
!! algorithm.
REAL(rp) :: ratio
!! Ratio used to accept or reject a sampling in the Metropolis-Hastings
!! algorithm.
REAL(rp) :: rand_unif
!! Uniform random deviate in the interval \([0,1]\).
REAL(rp), DIMENSION(3) :: V
!! Sampled velocity.
REAL(rp), DIMENSION(3) :: U
!! Sampled velocity.
REAL(rp), DIMENSION(3) :: b = (/1.0_rp,0.0_rp,0.0_rp/)
!! Temporary variable representing a unit vector along the \(x\)-axis.
INTEGER :: ii
!! Iterator.
INTEGER :: ppp
!! Number of particles per species.
Vmax = 0.9_rp
Vth = SQRT(spp%Eo*ABS(spp%q)/spp%m)
ppp = spp%ppp
V = (/0.0_rp,0.0_rp,0.0_rp/)
sv = Vth/10.0_rp
ii=2_idef
do while (ii .LE. 1000_idef)
U(1) = V(1) + random_norm(0.0_rp,sv)
do while (ABS(U(1)) .GT. Vmax)
U(1) = V(1) + random_norm(0.0_rp,sv)
end do
U(2) = V(2) + random_norm(0.0_rp,sv)
do while (ABS(U(2)) .GT. Vmax)
U(2) = V(2) + random_norm(0.0_rp,sv)
end do
U(3) = V(3) + random_norm(0.0_rp,sv)
do while (ABS(U(3)) .GT. Vmax)
U(3) = V(3) + random_norm(0.0_rp,sv)
end do
ratio = fth_3V(Vth,U)/fth_3V(Vth,V)
if (ratio .GE. 1.0_rp) then
V = U
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (ratio .GT. rand_unif) then
V = U
ii = ii + 1_idef
end if
end if
end do
spp%vars%V(1,1) = V(1)
spp%vars%V(1,2) = V(2)
spp%vars%V(1,3) = V(3)
ii=2_idef
do while (ii .LE. ppp)
U(1) = spp%vars%V(ii-1,1) + random_norm(0.0_rp,sv)
do while (ABS(U(1)) .GT. Vmax)
U(1) = spp%vars%V(ii-1,1) + random_norm(0.0_rp,sv)
end do
U(2) = spp%vars%V(ii-1,2) + random_norm(0.0_rp,sv)
do while (ABS(U(2)) .GT. Vmax)
U(2) = spp%vars%V(ii-1,2) + random_norm(0.0_rp,sv)
end do
U(3) = spp%vars%V(ii-1,3) + random_norm(0.0_rp,sv)
do while (ABS(U(3)) .GT. Vmax)
U(3) = spp%vars%V(ii-1,3) + random_norm(0.0_rp,sv)
end do
ratio = fth_3V(Vth,U)/fth_3V(Vth,spp%vars%V(ii-1,:))
if (ratio .GE. 1.0_rp) then
spp%vars%V(ii,1) = U(1)
spp%vars%V(ii,2) = U(2)
spp%vars%V(ii,3) = U(3)
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (ratio .GT. rand_unif) then
spp%vars%V(ii,1) = U(1)
spp%vars%V(ii,2) = U(2)
spp%vars%V(ii,3) = U(3)
ii = ii + 1_idef
end if
end if
end do
do ii=1_idef,ppp
vsq = spp%vars%V(ii,1)*spp%vars%V(ii,1) &
+ spp%vars%V(ii,2)*spp%vars%V(ii,2) &
+ spp%vars%V(ii,3)*spp%vars%V(ii,3)
spp%vars%g(ii) = 1.0_rp/SQRT(1.0_rp - vsq)
spp%vars%eta(ii) = ACOS(spp%vars%V(ii,1)/SQRT(vsq))
end do
spp%go = spp%Eo/(spp%m*C_C**2)
spp%etao = 90.0_rp
end subroutine thermal_distribution