JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Fixed an issue with the random Maxwellian subroutines

Browse source
This commit is contained in:
Jorge Gonzalez 2026-02-04 13:52:24 +01:00
commit 23bba31005
2 changed files with 15 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

23
src/modules/common/moduleRandom.f90
View file

@ -48,24 +48,23 @@ MODULE moduleRandom
END FUNCTION randomIntAB END FUNCTION randomIntAB
!Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method
FUNCTION randomMaxwellian() RESULT(rnd) function randomMaxwellian() result(rnd)
USE moduleConstParam, ONLY: PI USE moduleConstParam, only: pi
IMPLICIT NONE implicit none
REAL(8):: rnd real(8):: rnd
REAL(8):: v1, v2, Rsquare real(8):: v1, v2, Rsquare
Rsquare = 1.D0 v1 = 0.d0
do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) do while (v1 <= 0.d0)
v1 = 2.D0 * random() - 1.D0 v1 = random()
v2 = 2.D0 * random() - 1.D0
Rsquare = v1**2 + v2**2
end do end do
v2 = random()
rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) rnd = sqrt(-2.d0*log(v1))*cos(2*pi*v2)
END FUNCTION randomMaxwellian end function randomMaxwellian
!Returns a random number in a Maxwellian distribution of mean 0 and width 1 !Returns a random number in a Maxwellian distribution of mean 0 and width 1
FUNCTION randomHalfMaxwellian() RESULT(rnd) FUNCTION randomHalfMaxwellian() RESULT(rnd)

14
src/modules/moduleInject.f90
View file

@ -257,7 +257,7 @@ MODULE moduleInject
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m REAL(8), INTENT(in):: temperature, m
velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistMaxwellian END SUBROUTINE initVelDistMaxwellian
@ -267,7 +267,7 @@ MODULE moduleInject
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: temperature, m REAL(8), INTENT(in):: temperature, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian END SUBROUTINE initVelDistHalfMaxwellian
@ -289,7 +289,7 @@ MODULE moduleInject
REAL(8):: v REAL(8):: v
v = 0.D0 v = 0.D0
v = sqrt(2.0)*self%vTh*randomMaxwellian() v = self%vTh*randomMaxwellian()/sqrt(2.d0)
END FUNCTION randomVelMaxwellian END FUNCTION randomVelMaxwellian
@ -302,7 +302,7 @@ MODULE moduleInject
REAL(8):: v REAL(8):: v
v = 0.D0 v = 0.D0
v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0)
END FUNCTION randomVelHalfMaxwellian END FUNCTION randomVelHalfMaxwellian
@ -393,12 +393,6 @@ MODULE moduleInject
self%v(3)%obj%randomVel() /) self%v(3)%obj%randomVel() /)
end do end do
!If injecting a no-drift distribution and velocity is negative, reflect
if ((self%vMod == 0.D0) .and. &
(dot_product(direction, partInj(n)%v) < 0.D0)) then
partInj(n)%v = - partInj(n)%v
end if
!Obtain natural coordinates of particle in cell !Obtain natural coordinates of particle in cell
partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)