From 23bba3100555218705995662374045ed8be6cad4 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 4 Feb 2026 13:52:24 +0100 Subject: [PATCH] Fixed an issue with the random Maxwellian subroutines --- src/modules/common/moduleRandom.f90 | 23 +++++++++++------------ src/modules/moduleInject.f90 | 14 ++++---------- 2 files changed, 15 insertions(+), 22 deletions(-) diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 156f5e4..1b8ec67 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -48,24 +48,23 @@ MODULE moduleRandom END FUNCTION randomIntAB !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method - FUNCTION randomMaxwellian() RESULT(rnd) - USE moduleConstParam, ONLY: PI - IMPLICIT NONE + function randomMaxwellian() result(rnd) + USE moduleConstParam, only: pi + implicit none - REAL(8):: rnd - REAL(8):: v1, v2, Rsquare + real(8):: rnd + real(8):: v1, v2, Rsquare - Rsquare = 1.D0 - do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) - v1 = 2.D0 * random() - 1.D0 - v2 = 2.D0 * random() - 1.D0 - Rsquare = v1**2 + v2**2 + v1 = 0.d0 + do while (v1 <= 0.d0) + v1 = random() 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 FUNCTION randomHalfMaxwellian() RESULT(rnd) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index cb1c7fb..ff8a694 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -257,7 +257,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -267,7 +267,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian @@ -289,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = sqrt(2.0)*self%vTh*randomMaxwellian() + v = self%vTh*randomMaxwellian()/sqrt(2.d0) END FUNCTION randomVelMaxwellian @@ -302,7 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() + v = self%vTh*randomHalfMaxwellian()/sqrt(2.d0) END FUNCTION randomVelHalfMaxwellian @@ -393,12 +393,6 @@ MODULE moduleInject self%v(3)%obj%randomVel() /) 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 partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)