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

Change in calculation of reduced mass and relative energy

Browse source 
Now reduced mass and relative energy are calculated on the fly per
collision.
This commit is contained in:
Jorge Gonzalez 2022-12-09 09:54:44 +01:00
commit ae8aa9075e
2 changed files with 6 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

6
src/modules/mesh/moduleMesh.f90
View file

@ -660,6 +660,7 @@ MODULE moduleMesh
use moduleRefParam use moduleRefParam
USE moduleRandom USE moduleRandom
USE moduleOutput USE moduleOutput
USE moduleMath
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self CLASS(meshGeneric), INTENT(inout), TARGET:: self
@ -673,7 +674,7 @@ MODULE moduleMesh
TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:) TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
TYPE(particle), POINTER:: part_i, part_j TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n, c INTEGER:: n, c
REAL(8):: vRel, eRel REAL(8):: vRel, rMass, eRel
REAL(8):: sigmaVrelTotal REAL(8):: sigmaVrelTotal
REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
REAL(8):: rnd !Random number for collision REAL(8):: rnd !Random number for collision
@ -749,7 +750,8 @@ MODULE moduleMesh
!Obtain the cross sections for the different processes !Obtain the cross sections for the different processes
!TODO: From here it might be a procedure in interactionMatrix !TODO: From here it might be a procedure in interactionMatrix
vRel = NORM2(part_i%v-part_j%v) vRel = NORM2(part_i%v-part_j%v)
eRel = interactionMatrix(k)%rMass*vRel**2 rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
eRel = rMass*vRel**2
CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
!Update maximum sigma*v_rel !Update maximum sigma*v_rel

8
src/modules/moduleCollisions.f90
View file

@ -70,7 +70,6 @@ MODULE moduleCollisions
CLASS(speciesGeneric), POINTER:: sp_i CLASS(speciesGeneric), POINTER:: sp_i
CLASS(speciesGeneric), POINTER:: sp_j CLASS(speciesGeneric), POINTER:: sp_j
INTEGER:: amount INTEGER:: amount
REAL(8):: rMass !Reduced mass
TYPE(collisionCont), ALLOCATABLE:: collisions(:) TYPE(collisionCont), ALLOCATABLE:: collisions(:)
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initInteractionBinary PROCEDURE, PASS:: init => initInteractionBinary
@ -133,7 +132,6 @@ MODULE moduleCollisions
ALLOCATE(interactionMatrix(1:nCollPairs)) ALLOCATE(interactionMatrix(1:nCollPairs))
interactionMatrix(:)%amount = 0 interactionMatrix(:)%amount = 0
interactionMatrix(:)%rMass = 0.D0
END SUBROUTINE initInteractionMatrix END SUBROUTINE initInteractionMatrix
@ -169,8 +167,6 @@ MODULE moduleCollisions
mass_i = species(i)%obj%m mass_i = species(i)%obj%m
mass_j = species(j)%obj%m mass_j = species(j)%obj%m
self%rMass = reducedMass(mass_i, mass_j)
ALLOCATE(self%collisions(1:self%amount)) ALLOCATE(self%collisions(1:self%amount))
END SUBROUTINE initInteractionBinary END SUBROUTINE initInteractionBinary
@ -234,7 +230,7 @@ MODULE moduleCollisions
m_i = part_i%species%m m_i = part_i%species%m
m_j = part_j%species%m m_j = part_j%species%m
!Applies the collision !Applies the collision
vCM = velocityCM(m_i, part_i%v, m_j, part_j%v) vCM = velocityCM(part_i%weight*m_i, part_i%v, part_j%weight*m_j, part_j%v)
vp = vRel*randomDirectionVHS() vp = vRel*randomDirectionVHS()
!Assign velocities to particles !Assign velocities to particles
@ -311,7 +307,7 @@ MODULE moduleCollisions
REAL(8), DIMENSION(1:3):: vChange REAL(8), DIMENSION(1:3):: vChange
TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newElectron
rMass = reducedMass(part_i%species%m, part_j%species%m) rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
eRel = rMass*vRel**2 eRel = rMass*vRel**2
!Relative energy must be higher than threshold !Relative energy must be higher than threshold
IF (eRel > self%eThreshold) THEN IF (eRel > self%eThreshold) THEN