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

Finally, the right values for velocity distributions

Browse source
This commit is contained in:
Jorge Gonzalez 2026-04-04 20:03:20 +02:00
commit cf4488976a
3 changed files with 44 additions and 90 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -67,19 +67,20 @@ MODULE moduleRandom
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)
IMPLICIT NONE IMPLICIT NONE
REAL(8):: rnd real(8):: rnd
REAL(8):: x real(8):: v1
rnd = 0.D0 rnd = 0.D0
x = 0.D0 v1 = 0.D0
DO WHILE (x == 0.D0) do while (v1 == 0.D0)
CALL RANDOM_NUMBER(x) v1 = random()
END DO
end do
rnd = DSQRT(-DLOG(x)) rnd = sqrt(-2.d0*log(v1))
END FUNCTION randomHalfMaxwellian END FUNCTION randomHalfMaxwellian

3
src/modules/common/velocityDistribution.f90
View file

@ -53,12 +53,11 @@ module velocityDistribution
function randomVelMaxwellian(self) result (v) function randomVelMaxwellian(self) result (v)
use moduleRandom use moduleRandom
implicit none implicit none
class(velDistMaxwellian), intent(in):: self class(velDistMaxwellian), intent(in):: self
real(8):: v real(8):: v
v = 0.D0 v = 0.D0
v = self%vTh*randomMaxwellian()/sqrt(2.d0) v = self%vTh*randomMaxwellian()
end function randomVelMaxwellian end function randomVelMaxwellian

114
src/modules/moduleInject.f90
View file

@ -12,7 +12,7 @@ MODULE moduleInject
REAL(8):: temperature(1:3) !Temperature REAL(8):: temperature(1:3) !Temperature
REAL(8):: n(1:3) !Direction of injection REAL(8):: n(1:3) !Direction of injection
LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
INTEGER:: nParticles !Number of particles to introduce each time step INTEGER:: nParticles !Number of particles to inject each time step
CLASS(speciesGeneric), POINTER:: species !Species of injection CLASS(speciesGeneric), POINTER:: species !Species of injection
INTEGER:: nEdges INTEGER:: nEdges
type(meshEdgePointer), allocatable:: edges(:) type(meshEdgePointer), allocatable:: edges(:)
@ -72,7 +72,6 @@ MODULE moduleInject
USE moduleRefParam USE moduleRefParam
USE moduleConstParam USE moduleConstParam
USE moduleSpecies USE moduleSpecies
USE moduleSolver
USE moduleErrors USE moduleErrors
use json_module use json_module
IMPLICIT NONE IMPLICIT NONE
@ -211,8 +210,6 @@ MODULE moduleInject
END DO END DO
self%nParticles = SUM(self%particlesPerEdge)
ELSE ELSE
! No particles assigned per edge, use the species weight ! No particles assigned per edge, use the species weight
self%weightPerEdge = self%species%weight self%weightPerEdge = self%species%weight
@ -220,13 +217,14 @@ MODULE moduleInject
self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight)) self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight))
END DO END DO
self%nParticles = SUM(self%particlesPerEdge)
!Rescale weight to match flux !Rescale weight to match flux
self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles))
END IF END IF
! Total number of particles to inject
self%nParticles = SUM(self%particlesPerEdge)
!Scale particles for different species steps !Scale particles for different species steps
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
@ -314,126 +312,82 @@ MODULE moduleInject
end subroutine updateQuasiNeutral end subroutine updateQuasiNeutral
!Injection of particles
SUBROUTINE doInjects()
USE moduleSpecies
USE moduleSolver
IMPLICIT NONE
INTEGER:: i
!$OMP SINGLE
nPartInj = 0
DO i = 1, nInject
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
nPartInj = nPartInj + injects(i)%obj%nParticles
END IF
END DO
IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
ALLOCATE(partInj(1:nPartInj))
!$OMP END SINGLE
DO i=1, nInject
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
CALL injects(i)%obj%addParticles()
END IF
END DO
END SUBROUTINE doInjects
!Add particles for the injection !Add particles for the injection
SUBROUTINE addParticles(self) SUBROUTINE addParticles(self, partArray)
USE moduleSpecies USE moduleSpecies
USE moduleSolver
USE moduleMesh USE moduleMesh
USE moduleRandom USE moduleRandom
USE moduleErrors USE moduleErrors
IMPLICIT NONE IMPLICIT NONE
CLASS(injectGeneric), INTENT(in):: self CLASS(injectGeneric), INTENT(in):: self
INTEGER, SAVE:: nMin type(particle), allocatable, intent(inout):: partArray(:)
INTEGER:: i, e INTEGER:: i, e
INTEGER:: n, sp INTEGER:: n, sp
CLASS(meshEdge), POINTER:: randomEdge CLASS(meshEdge), POINTER:: edge
REAL(8):: direction(1:3) REAL(8):: direction(1:3)
!$omp single
n = 0
!$omp end single
!Insert particles !Insert particles
!$OMP SINGLE
nMin = 0
DO i = 1, self%id -1
IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN
nMin = nMin + injects(i)%obj%nParticles
END IF
END DO
nMin = nMin + 1
!$OMP END SINGLE
!$OMP DO !$OMP DO
DO e = 1, self%nEdges DO e = 1, self%nEdges
! Select edge for injection ! Select edge for injection
randomEdge => self%edges(e)%obj edge => self%edges(e)%obj
! Inject particles in edge ! Inject particles in edge
DO i = 1, self%particlesPerEdge(e) DO i = 1, self%particlesPerEdge(e)
! Index in the global partInj array ! Index in the partInj array
n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i n = n + 1
!Particle is considered to be outside the domain !Particle is considered to be outside the domain for now
partInj(n)%n_in = .FALSE. partArray(n)%n_in = .FALSE.
!Random position in edge !Random position in edge
partInj(n)%r = randomEdge%randPos() partArray(n)%r = edge%randPos()
!Assign weight to particle. !Assign weight to particle.
partInj(n)%weight = self%weightPerEdge(e) partArray(n)%weight = self%weightPerEdge(e)
!Volume associated to the edge: !Volume associated to the edge:
IF (ASSOCIATED(randomEdge%e1)) THEN IF (ASSOCIATED(edge%e1)) THEN
partInj(n)%cell = randomEdge%e1%n partArray(n)%cell = edge%e1%n
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN ELSEIF (ASSOCIATED(edge%e2)) THEN
partInj(n)%cell = randomEdge%e2%n partArray(n)%cell = edge%e2%n
ELSE ELSE
CALL criticalError("No Volume associated to edge", 'addParticles') CALL criticalError("No Volume associated to edge", 'addParticles')
END IF END IF
partInj(n)%cellColl = randomEdge%eColl%n partArray(n)%cellColl = edge%eColl%n
sp = self%species%n sp = self%species%n
!Assign particle type !Assign particle type
partInj(n)%species => self%species partArray(n)%species => self%species
if (all(self%n == 0.D0)) then if (all(self%n == 0.D0)) then
direction = randomEdge%normal direction = edge%normal
else else
direction = self%n direction = self%n
end if end if
partInj(n)%v = 0.D0 partArray(n)%v = 0.D0
do while(dot_product(partInj(n)%v, direction) <= 0.d0) ! do while(dot_product(partInj(n)%v, direction) <= 0.d0)
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & partArray(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), & self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /) self%v(3)%obj%randomVel() /)
!If injecting a no-drift distribution and velocity is negative, reflect !If injecting a no-drift distribution and velocity is negative, reflect
if ((self%vMod == 0.D0) .and. & if ((self%vMod == 0.D0) .and. &
(dot_product(partInj(n)%v, direction) <= 0.D0)) then (dot_product(partArray(n)%v, direction) <= 0.D0)) then
partInj(n)%v = - partInj(n)%v partArray(n)%v = - partArray(n)%v
end if end if
end do ! end do
!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) partArray(n)%Xi = mesh%cells(partArray(n)%cell)%obj%phy2log(partArray(n)%r)
!Push new particle with the minimum time step
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle
CALL solver%updateParticleCell(partInj(n))
END DO END DO
@ -500,7 +454,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(2.d0*temperature/m)) velDist = velDistMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistMaxwellian END SUBROUTINE initVelDistMaxwellian
@ -510,7 +464,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(2.d0*temperature/m)) velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m))
END SUBROUTINE initVelDistHalfMaxwellian END SUBROUTINE initVelDistHalfMaxwellian