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

I'm dying with hay fever but I have to commit

Browse source 
I'm feeling awful but I have work in my desktop that I need to commit so
I can work with my laptop while I'm at the IEPC 2022 in Boston.
This commit is contained in:
Jorge Gonzalez 2022-06-10 16:07:14 +02:00
commit 23e2fe9bae
17 changed files with 313 additions and 140 deletions
Download patch file Download diff file Expand all files Collapse all files

3
data/collisions/EL_e-Ar.dat Normal file
View file

@ -0,0 +1,3 @@
# Relative energy (eV) cross section (m^2)
0.0 1e-20
1000.0 1e-20

148
data/collisions/IO_e-Li.dat Normal file
View file

@ -0,0 +1,148 @@
# Rusudan I., et al., Atoms 9(4):90 2021
# Relative energy (eV) cross section (m^2)
5.393 1.061E-23
6.000 5.885E-21
7.000 1.352E-20
8.000 1.927E-20
9.000 2.364E-20
10.000 2.698E-20
11.000 2.956E-20
12.000 3.156E-20
13.000 3.309E-20
14.000 3.427E-20
15.000 3.517E-20
16.000 3.585E-20
17.000 3.635E-20
18.000 3.670E-20
19.000 3.694E-20
20.000 3.708E-20
21.000 3.714E-20
22.000 3.714E-20
23.000 3.709E-20
24.000 3.700E-20
25.000 3.686E-20
26.000 3.671E-20
27.000 3.652E-20
28.000 3.632E-20
29.000 3.610E-20
30.000 3.588E-20
31.000 3.564E-20
32.000 3.539E-20
33.000 3.514E-20
34.000 3.488E-20
35.000 3.462E-20
36.000 3.435E-20
37.000 3.409E-20
38.000 3.383E-20
39.000 3.356E-20
40.000 3.330E-20
41.000 3.304E-20
42.000 3.277E-20
43.000 3.252E-20
44.000 3.226E-20
45.000 3.201E-20
46.000 3.175E-20
47.000 3.151E-20
48.000 3.126E-20
49.000 3.102E-20
50.000 3.078E-20
51.000 3.054E-20
52.000 3.031E-20
53.000 3.008E-20
54.000 2.985E-20
55.000 2.963E-20
56.000 2.941E-20
57.000 2.919E-20
58.000 2.898E-20
59.000 2.877E-20
60.000 2.856E-20
61.000 2.835E-20
62.000 2.815E-20
63.000 2.795E-20
64.000 2.776E-20
65.000 2.756E-20
66.000 2.737E-20
67.000 2.719E-20
68.000 2.700E-20
69.000 2.682E-20
70.000 2.664E-20
71.000 2.646E-20
72.000 2.629E-20
73.000 2.612E-20
74.000 2.595E-20
75.000 2.578E-20
76.000 2.562E-20
77.000 2.546E-20
78.000 2.530E-20
79.000 2.514E-20
80.000 2.499E-20
81.000 2.483E-20
82.000 2.468E-20
83.000 2.453E-20
84.000 2.439E-20
85.000 2.424E-20
86.000 2.410E-20
87.000 2.396E-20
88.000 2.382E-20
89.000 2.369E-20
90.000 2.355E-20
91.000 2.342E-20
92.000 2.329E-20
93.000 2.316E-20
94.000 2.303E-20
95.000 2.290E-20
96.000 2.278E-20
97.000 2.266E-20
98.000 2.253E-20
99.000 2.241E-20
100.000 2.230E-20
101.000 2.218E-20
102.000 2.206E-20
103.000 2.195E-20
104.000 2.184E-20
105.000 2.173E-20
106.000 2.162E-20
107.000 2.151E-20
108.000 2.140E-20
109.000 2.129E-20
110.000 2.119E-20
111.000 2.109E-20
112.000 2.098E-20
113.000 2.088E-20
114.000 2.078E-20
115.000 2.068E-20
116.000 2.059E-20
117.000 2.049E-20
118.000 2.039E-20
119.000 2.030E-20
120.000 2.021E-20
121.000 2.011E-20
122.000 2.002E-20
123.000 1.993E-20
124.000 1.984E-20
125.000 1.976E-20
126.000 1.967E-20
127.000 1.958E-20
128.000 1.950E-20
129.000 1.941E-20
130.000 1.933E-20
131.000 1.924E-20
132.000 1.916E-20
133.000 1.908E-20
134.000 1.900E-20
135.000 1.892E-20
136.000 1.884E-20
137.000 1.877E-20
138.000 1.869E-20
139.000 1.861E-20
140.000 1.854E-20
141.000 1.846E-20
142.000 1.839E-20
143.000 1.831E-20
144.000 1.824E-20
145.000 1.817E-20
146.000 1.810E-20
147.000 1.803E-20
148.000 1.796E-20
149.000 1.789E-20
150.000 1.782E-20

2
src/modules/mesh/0D/moduleMesh0D.f90
View file

@ -77,8 +77,6 @@ MODULE moduleMesh0D
self%volume = 1.D0 self%volume = 1.D0
self%n1%v = 1.D0 self%n1%v = 1.D0
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

2
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -216,8 +216,6 @@ MODULE moduleMesh1DCart
self%n1%v = self%n1%v + self%arNodes(1) self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2) self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

2
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -218,8 +218,6 @@ MODULE moduleMesh1DRad
self%n1%v = self%n1%v + self%arNodes(1) self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2) self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

4
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -307,8 +307,6 @@ MODULE moduleMesh2DCart
self%n3%v = self%n3%v + self%arNodes(3) self%n3%v = self%n3%v + self%arNodes(3)
self%n4%v = self%n4%v + self%arNodes(4) self%n4%v = self%n4%v + self%arNodes(4)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
@ -637,8 +635,6 @@ MODULE moduleMesh2DCart
self%n2%v = self%n2%v + self%arNodes(2) self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3) self%n3%v = self%n3%v + self%arNodes(3)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

4
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -295,8 +295,6 @@ MODULE moduleMesh2DCyl
self%n3%v = self%n3%v + self%arNodes(3) self%n3%v = self%n3%v + self%arNodes(3)
self%n4%v = self%n4%v + self%arNodes(4) self%n4%v = self%n4%v + self%arNodes(4)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))
@ -658,8 +656,6 @@ MODULE moduleMesh2DCyl
self%n2%v = self%n2%v + self%arNodes(2) self%n2%v = self%n2%v + self%arNodes(2)
self%n3%v = self%n3%v + self%arNodes(3) self%n3%v = self%n3%v + self%arNodes(3)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

2
src/modules/mesh/3DCart/moduleMesh3DCart.f90
View file

@ -281,8 +281,6 @@ MODULE moduleMesh3DCart
self%n3%v = self%n3%v + volNodes(3) self%n3%v = self%n3%v + volNodes(3)
self%n4%v = self%n4%v + volNodes(4) self%n4%v = self%n4%v + volNodes(4)
self%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%listPart_in(1:nSpecies))

3
src/modules/mesh/inout/0D/moduleMeshOutput0D.f90
View file

@ -45,6 +45,7 @@ MODULE moduleMeshOutput0D
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
INTEGER:: k
fileName='OUTPUT_Collisions.dat' fileName='OUTPUT_Collisions.dat'
IF (t == 0) THEN IF (t == 0) THEN
@ -56,7 +57,7 @@ MODULE moduleMeshOutput0D
END IF END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
WRITE(20, "(ES20.6E3, I20)") REAL(t)*tauMin*ti_ref, self%vols(1)%obj%nColl WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%vols(1)%obj%tallyColl(k)%tally, k=1,nCollPairs)
CLOSE(20) CLOSE(20)
END SUBROUTINE printColl0D END SUBROUTINE printColl0D

9
src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90
View file

@ -98,6 +98,7 @@ MODULE moduleMeshOutputGmsh2
CLASS(meshGeneric), INTENT(inout):: self CLASS(meshGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
INTEGER:: numEdges INTEGER:: numEdges
INTEGER:: k, c
INTEGER:: n INTEGER:: n
REAL(8):: time REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName CHARACTER(:), ALLOCATABLE:: fileName
@ -125,9 +126,11 @@ MODULE moduleMeshOutputGmsh2
WRITE(60, "(A)") '$MeshFormat' WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8' WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat' WRITE(60, "(A)") '$EndMeshFormat'
DO k = 1, nCollPairs
DO c = 1, interactionMatrix(k)%amount
WRITE(60, "(A)") '$ElementData' WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1' WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"' WRITE(60, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c
WRITE(60, *) 1 WRITE(60, *) 1
WRITE(60, *) time WRITE(60, *) time
WRITE(60, *) 3 WRITE(60, *) 3
@ -135,9 +138,11 @@ MODULE moduleMeshOutputGmsh2
WRITE(60, *) 1 WRITE(60, *) 1
WRITE(60, *) self%numVols WRITE(60, *) self%numVols
DO n=1, self%numVols DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%nColl WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c)
END DO END DO
WRITE(60, "(A)") '$EndElementData' WRITE(60, "(A)") '$EndElementData'
END DO
END DO
CLOSE(60) CLOSE(60)

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

@ -3,6 +3,7 @@ MODULE moduleMesh
USE moduleList USE moduleList
USE moduleOutput USE moduleOutput
USE moduleBoundary USE moduleBoundary
USE moduleCollisions
IMPLICIT NONE IMPLICIT NONE
!Generic mesh element !Generic mesh element
@ -148,15 +149,15 @@ MODULE moduleMesh
!Parent of Volume element !Parent of Volume element
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshVol
!Maximum collision rate !Maximum collision rate
REAL(8):: sigmaVrelMax = 0.D0 REAL(8), ALLOCATABLE:: sigmaVrelMax(:)
!Arrays for counting number of collisions
TYPE(tallyCollisions), ALLOCATABLE:: tallyColl(:)
!Volume !Volume
REAL(8):: volume = 0.D0 REAL(8):: volume = 0.D0
!List of particles inside the volume !List of particles inside the volume
TYPE(listNode), ALLOCATABLE:: listPart_in(:) TYPE(listNode), ALLOCATABLE:: listPart_in(:)
!Lock indicator for listPart_in !Lock indicator for listPart_in
INTEGER(KIND=OMP_LOCK_KIND):: lock INTEGER(KIND=OMP_LOCK_KIND):: lock
!Number of collisions per volume
INTEGER:: nColl = 0
!Total weight of particles inside cell !Total weight of particles inside cell
REAL(8), ALLOCATABLE:: totalWeight(:) REAL(8), ALLOCATABLE:: totalWeight(:)
CONTAINS CONTAINS
@ -658,15 +659,17 @@ MODULE moduleMesh
USE moduleList USE moduleList
use moduleRefParam use moduleRefParam
USE moduleRandom USE moduleRandom
USE moduleOutput
IMPLICIT NONE IMPLICIT NONE
CLASS(meshGeneric), INTENT(inout), TARGET:: self CLASS(meshGeneric), INTENT(inout), TARGET:: self
INTEGER, INTENT(in):: t INTEGER, INTENT(in):: t
INTEGER:: e INTEGER:: e
CLASS(meshVol), POINTER:: vol CLASS(meshVol), POINTER:: vol
INTEGER:: k, nPairs, i, j INTEGER:: k, i, j
INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision REAL(8):: pMax !Maximum probability of collision
INTEGER:: nColl
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
@ -674,21 +677,23 @@ MODULE moduleMesh
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
INTEGER:: realCollisions
IF (MOD(t, everyColl) == 0) THEN IF (MOD(t, everyColl) == 0) THEN
!Collisions need to be performed in this iteration !Collisions need to be performed in this iteration
!$OMP DO SCHEDULE(DYNAMIC) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
DO e=1, self%numVols DO e=1, self%numVols
realCollisions = 0
vol => self%vols(e)%obj vol => self%vols(e)%obj
!TODO: Simplify this, to many sublevels !TODO: Simplify this, to many sublevels
!Iterate over the number of pairs !Iterate over the number of pairs
nPairs = SIZE(interactionMatrix) !TODO: This does not change, make a variable in a module DO k = 1, nCollPairs
DO k = 1, nPairs !Reset tally of collisions
IF (collOutput) THEN
vol%tallyColl(k)%tally = 0
END IF
IF (interactionMatrix(k)%amount > 0) THEN IF (interactionMatrix(k)%amount > 0) THEN
!Select the species for the collision pair !Select the species for the collision pair
i = interactionMatrix(k)%sp_i%n i = interactionMatrix(k)%sp_i%n
@ -703,22 +708,22 @@ MODULE moduleMesh
nPart = nPart_i + nPart_j nPart = nPart_i + nPart_j
!Resets the number of collisions in the cell !Resets the number of collisions in the cell
vol%nColl = 0 nColl = 0
!Probability of collision for pair i-j !Probability of collision for pair i-j
pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax*tauColl/vol%volume pMax = (vol%totalWeight(i) + vol%totalWeight(j))*vol%sigmaVrelMax(k)*tauColl/vol%volume
!Number of collisions in the cell !Number of collisions in the cell
vol%nColl = NINT(REAL(nPart)*pMax*0.5D0) nColl = NINT(REAL(nPart)*pMax*0.5D0)
!Converts the list of particles to an array for easy access !Converts the list of particles to an array for easy access
IF (vol%nColl > 0) THEN IF (nColl > 0) THEN
partTemp_i = vol%listPart_in(i)%convert2Array() partTemp_i = vol%listPart_in(i)%convert2Array()
partTemp_j = vol%listPart_in(j)%convert2Array() partTemp_j = vol%listPart_in(j)%convert2Array()
END IF END IF
DO n = 1, vol%nColl DO n = 1, nColl
!Select random particles !Select random particles
part_i => NULL() part_i => NULL()
part_j => NULL() part_j => NULL()
@ -735,13 +740,12 @@ MODULE moduleMesh
!If particles do not belong to the species, skip collision !If particles do not belong to the species, skip collision
!This can happen, for example, if particle has been previously ionized or removed !This can happen, for example, if particle has been previously ionized or removed
!TODO: Try to find a way to no lose these collisions. Check new 'k' and use that for the collision, maybe? !TODO: Try to find a way to no lose these collisions. Maybe check new 'k' and use that for the collision, maybe?
IF (part_i%species%n /= i .OR. & IF (part_i%species%n /= i .OR. &
part_j%species%n /= j) THEN part_j%species%n /= j) THEN
CYCLE CYCLE
END IF END IF
!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)
@ -749,26 +753,34 @@ MODULE moduleMesh
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
IF (sigmaVrelTotal > vol%sigmaVrelMax) THEN IF (sigmaVrelTotal > vol%sigmaVrelMax(k)) THEN
vol%sigmaVrelMax = sigmaVrelTotal vol%sigmaVrelMax(k) = sigmaVrelTotal
END IF END IF
ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
probabilityColl(0) = 0.0 probabilityColl = 0.0
probabilityColl(1:interactionMatrix(k)%amount) = sigmaVrel/vol%sigmaVrelMax DO c = 1, interactionMatrix(k)%amount
probabilityColl(c) = sigmaVrel(c)/vol%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1))
END DO
!Selects random number between 0 and 1 !Selects random number between 0 and 1
rnd = random() rnd = random()
!If the random number is below the total probability of collision, collide particles !If the random number is below the total probability of collision, collide particles
IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax) THEN IF (rnd < sigmaVrelTotal / vol%sigmaVrelMax(k)) THEN
!Loop over collisions !Loop over collisions
DO c = 1, interactionMatrix(k)%amount DO c = 1, interactionMatrix(k)%amount
IF (SUM(probabilityColl(0:c-1)) + rnd <= probabilityColl(c)) THEN IF (rnd <= probabilityColl(c)) THEN
CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
realCollisions = realCollisions + 1
!If collisions are gonna be output, count the collision
IF (collOutput) THEN
vol%tallyColl(k)%tally(c) = vol%tallyColl(k)%tally(c) + 1
END IF
!A collision has ocurred, exit the loop !A collision has ocurred, exit the loop
EXIT EXIT
@ -782,16 +794,17 @@ MODULE moduleMesh
!Deallocate arrays for next collision !Deallocate arrays for next collision
DEALLOCATE(sigmaVrel, probabilityColl) DEALLOCATE(sigmaVrel, probabilityColl)
!End loop collisions in cell
END DO END DO
END IF END IF
END IF END IF
!End loop collision pairs
END DO END DO
vol%nColl = realCollisions !End loop volumes
END DO END DO
!$OMP END DO !$OMP END DO

7
src/modules/mesh/moduleMeshBoundary.f90
View file

@ -159,13 +159,6 @@ MODULE moduleMeshBoundary
newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r) newElectron%xi = mesh%vols(part%vol)%obj%phy2log(newElectron%r)
newIon%xi = newElectron%xi newIon%xi = newElectron%xi
newElectron%qm = part%qm
SELECT TYPE(spe => bound%species)
TYPE IS(speciesCharged)
newIon%qm = spe%qm
END SELECT
newElectron%weight = bound%species%weight newElectron%weight = bound%species%weight
newIon%weight = newElectron%weight newIon%weight = newElectron%weight

70
src/modules/moduleCollisions.f90
View file

@ -78,6 +78,14 @@ MODULE moduleCollisions
END TYPE interactionsBinary END TYPE interactionsBinary
!Type to count number of collisions
TYPE:: tallyCollisions
INTEGER, ALLOCATABLE:: tally(:)
END TYPE
!Number of collision pairs (nSpecies*(nSpecies+1)/2)
INTEGER:: nCollPairs = 0
!Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory !Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:) TYPE(interactionsBinary), ALLOCATABLE, TARGET:: interactionMatrix(:)
!Folder for collision cross section tables !Folder for collision cross section tables
@ -120,10 +128,9 @@ MODULE moduleCollisions
IMPLICIT NONE IMPLICIT NONE
TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:) TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:)
INTEGER:: nInteractions
nInteractions = (nSpecies*(nSpecies+1))/2 nCollPairs = (nSpecies*(nSpecies+1))/2
ALLOCATE(interactionMatrix(1:nInteractions)) ALLOCATE(interactionMatrix(1:nCollPairs))
interactionMatrix(:)%amount = 0 interactionMatrix(:)%amount = 0
interactionMatrix(:)%rMass = 0.D0 interactionMatrix(:)%rMass = 0.D0
@ -222,8 +229,7 @@ MODULE moduleCollisions
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
REAL(8), INTENT(in):: vRel REAL(8), INTENT(in):: vRel
REAL(8):: m_i, m_j REAL(8):: m_i, m_j
REAL(8), DIMENSION(1:3):: vCM REAL(8), DIMENSION(1:3):: vCM, vp
REAL(8):: vp(1:3)
m_i = part_i%species%m m_i = part_i%species%m
m_j = part_j%species%m m_j = part_j%species%m
@ -232,13 +238,8 @@ MODULE moduleCollisions
vp = vRel*randomDirectionVHS() vp = vRel*randomDirectionVHS()
!Assign velocities to particles !Assign velocities to particles
PRINT *, part_i%v
part_i%v = vCM + m_j*vp / (m_i + m_j) part_i%v = vCM + m_j*vp / (m_i + m_j)
PRINT *, part_i%v
PRINT *, part_j%v
part_j%v = vCM - m_i*vp / (m_i + m_j) part_j%v = vCM - m_i*vp / (m_i + m_j)
PRINT *, part_j%v
PRINT *
END SUBROUTINE collideBinaryElastic END SUBROUTINE collideBinaryElastic
@ -273,16 +274,21 @@ MODULE moduleCollisions
!Input energy is in eV. Convert to J with ev2J and then to !Input energy is in eV. Convert to J with ev2J and then to
!non-dimensional units. !non-dimensional units.
collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2) collision%eThreshold = energyThreshold*eV2J/(m_ref*v_ref**2)
!species for impacting electron
electronIndex = speciesName2Index(electron) electronIndex = speciesName2Index(electron)
SELECT TYPE(sp => species(electronIndex)%obj) SELECT TYPE(sp => species(electronIndex)%obj)
TYPE IS(speciesCharged) TYPE IS(speciesCharged)
collision%electron => sp collision%electron => sp
CLASS DEFAULT CLASS DEFAULT
CALL criticalError("Species " // sp%name // " chosen for ionization is not a charged species", 'initBinaryIonization') CALL criticalError("Species " // sp%name // " chosen for " // &
"secondary electron is not a charged species", 'initBinaryIonization')
END SELECT END SELECT
!momentum change per ionization process
collision%deltaV = sqrt(collision%eThreshold / collision%electron%m)
END SELECT END SELECT
END SUBROUTINE initBinaryIonization END SUBROUTINE initBinaryIonization
@ -302,8 +308,8 @@ MODULE moduleCollisions
REAL(8), INTENT(in):: vRel REAL(8), INTENT(in):: vRel
REAL(8):: rMass, eRel REAL(8):: rMass, eRel
TYPE(particle), POINTER:: electron => NULL(), neutral => NULL() TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
REAL(8), DIMENSION(1:3):: vChange
TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newElectron
REAL(8), DIMENSION(1:3):: vp_e, vp_n
rMass = reducedMass(part_i%species%m, part_j%species%m) rMass = reducedMass(part_i%species%m, part_j%species%m)
eRel = rMass*vRel**2 eRel = rMass*vRel**2
@ -322,35 +328,35 @@ MODULE moduleCollisions
END IF END IF
!Exchange energy between
vp_e = electron%v*(1.D0 - self%deltaV/NORM2(electron%v))
vp_n = neutral%v* (1.D0 + self%deltaV/NORM2(neutral%v) )
!Changes velocity of impacting electron
electron%v = vp_e
!Creates a new electron from ionization
ALLOCATE(newElectron)
newElectron%species => electron%species
newElectron%v = vp_n
newElectron%r = neutral%r
newElectron%xi = neutral%xi
newElectron%n_in = .TRUE.
newElectron%vol = neutral%vol
newElectron%volColl = neutral%volColl
newElectron%weight = neutral%weight
newElectron%qm = electron%qm
!Ionize neutral particle !Ionize neutral particle
SELECT TYPE(sp => neutral%species) SELECT TYPE(sp => neutral%species)
TYPE IS(speciesNeutral) TYPE IS(speciesNeutral)
CALL sp%ionize(neutral) CALL sp%ionize(neutral)
CLASS DEFAULT CLASS DEFAULT
CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization') ! CALL criticalError(sp%name // " is not a neutral", 'collideBinaryIonization')
RETURN
END SELECT END SELECT
!Exchange of velocity between particles
vChange = self%deltaV*randomDirectionVHS()
!Energy is loss by the primary electron
electron%v = electron%v - vChange
!Creates a new electron from ionization
ALLOCATE(newElectron)
newElectron%species => electron%species
!Secondary electorn gains energy from ionization
newElectron%v = vChange
newElectron%r = neutral%r
newElectron%xi = neutral%xi
newElectron%n_in = .TRUE.
newElectron%vol = neutral%vol
newElectron%volColl = neutral%volColl
newElectron%weight = neutral%weight
!Adds new electron to list of new particles from collisions !Adds new electron to list of new particles from collisions
CALL OMP_SET_LOCK(lockCollisions) CALL OMP_SET_LOCK(lockCollisions)
CALL partCollisions%add(newElectron) CALL partCollisions%add(newElectron)

33
src/modules/moduleInject.f90
View file

@ -82,6 +82,7 @@ MODULE moduleInject
INTEGER, INTENT(in):: i INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3) REAL(8), INTENT(in):: v, n(1:3), T(1:3)
INTEGER, INTENT(in):: sp, physicalSurface INTEGER, INTENT(in):: sp, physicalSurface
REAL(8):: tauInject
REAL(8), INTENT(in):: flow REAL(8), INTENT(in):: flow
CHARACTER(:), ALLOCATABLE, INTENT(in):: units CHARACTER(:), ALLOCATABLE, INTENT(in):: units
INTEGER:: e, et INTEGER:: e, et
@ -93,18 +94,19 @@ MODULE moduleInject
self%n = n / NORM2(n) self%n = n / NORM2(n)
self%T = T / T_ref self%T = T / T_ref
self%species => species(sp)%obj self%species => species(sp)%obj
tauInject = tau(self%species%n)
SELECT CASE(units) SELECT CASE(units)
CASE ("sccm") CASE ("sccm")
!Standard cubic centimeter per minute !Standard cubic centimeter per minute
self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight) self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight)
CASE ("A") CASE ("A")
!Input current in Ampers !Input current in Ampers
self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight)) self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight))
CASE ("part/s") CASE ("part/s")
!Input current in Ampers !Input current in Ampers
self%nParticles = INT(flow*tauMin*ti_ref/species(sp)%obj%weight) self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight)
CASE DEFAULT CASE DEFAULT
CALL criticalError("No support for units: " // units, 'initInject') CALL criticalError("No support for units: " // units, 'initInject')
@ -161,8 +163,6 @@ MODULE moduleInject
END DO END DO
self%sumWeight = self%cumWeight(self%nEdges) self%sumWeight = self%cumWeight(self%nEdges)
nPartInj = nPartInj + self%nParticles
END SUBROUTINE initInject END SUBROUTINE initInject
!Injection of particles !Injection of particles
@ -174,12 +174,24 @@ MODULE moduleInject
INTEGER:: i INTEGER:: i
!$OMP SINGLE !$OMP SINGLE
nPartInj = 0
DO i = 1, nInject
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
nPartInj = nPartInj + inject(i)%nParticles
END IF
END DO
IF (ALLOCATED(partInj)) DEALLOCATE(partInj) IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
ALLOCATE(partInj(1:nPartInj)) ALLOCATE(partInj(1:nPartInj))
!$OMP END SINGLE !$OMP END SINGLE
DO i=1, nInject DO i=1, nInject
IF (solver%pusher(inject(i)%species%n)%pushSpecies) THEN
CALL inject(i)%addParticles() CALL inject(i)%addParticles()
END IF
END DO END DO
END SUBROUTINE doInjects END SUBROUTINE doInjects
@ -240,7 +252,7 @@ MODULE moduleInject
CLASS(injectGeneric), INTENT(in):: self CLASS(injectGeneric), INTENT(in):: self
INTEGER:: randomX INTEGER:: randomX
INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
INTEGER:: n INTEGER:: n, sp
CLASS(meshEdge), POINTER:: randomEdge CLASS(meshEdge), POINTER:: randomEdge
!Insert particles !Insert particles
@ -251,12 +263,6 @@ MODULE moduleInject
partInj(nMin:nMax)%weight = self%species%weight partInj(nMin:nMax)%weight = self%species%weight
!Particle is considered to be outside the domain !Particle is considered to be outside the domain
partInj(nMin:nMax)%n_in = .FALSE. partInj(nMin:nMax)%n_in = .FALSE.
!Assign charge/mass to charged particle.
SELECT TYPE(sp => self%species)
TYPE IS(speciesCharged)
partInj(nMin:nMax)%qm = sp%qm
END SELECT
!$OMP END SINGLE !$OMP END SINGLE
!$OMP DO !$OMP DO
@ -278,6 +284,7 @@ MODULE moduleInject
END IF END IF
partInj(n)%volColl = randomEdge%eColl%n partInj(n)%volColl = randomEdge%eColl%n
sp = self%species%n
!Assign particle type !Assign particle type
partInj(n)%species => self%species partInj(n)%species => self%species
@ -289,7 +296,7 @@ MODULE moduleInject
!Obtain natural coordinates of particle in cell !Obtain natural coordinates of particle in cell
partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r) partInj(n)%xi = mesh%vols(partInj(n)%vol)%obj%phy2log(partInj(n)%r)
!Push new particle with the minimum time step !Push new particle with the minimum time step
CALL solver%pusher(self%species%n)%pushParticle(partInj(n), tauMin) CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle !Assign cell to new particle
CALL solver%updateParticleCell(partInj(n)) CALL solver%updateParticleCell(partInj(n))

58
src/modules/moduleInput.f90
View file

@ -132,6 +132,15 @@ MODULE moduleInput
CALL config%get(object // '.temperature', T_ref, found) CALL config%get(object // '.temperature', T_ref, found)
IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference') IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
!Derived parameters
L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
ti_ref = L_ref/v_ref !reference time
Vol_ref = L_ref**3 !reference volume
EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field
Volt_ref = EF_ref*L_ref !reference voltage
B_ref = m_ref / (ti_ref * qe) !reference magnetic field
!If a reference cross section is given, it is used !If a reference cross section is given, it is used
CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found) CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found)
@ -148,15 +157,6 @@ MODULE moduleInput
END IF END IF
!Derived parameters
L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !reference length
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
ti_ref = L_ref/v_ref !reference time
Vol_ref = L_ref**3 !reference volume
EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field
Volt_ref = EF_ref*L_ref !reference voltage
B_ref = m_ref / (ti_ref * qe) !reference magnetic field
END SUBROUTINE readReference END SUBROUTINE readReference
!Reads the specific case parameters !Reads the specific case parameters
@ -417,15 +417,6 @@ MODULE moduleInput
END IF END IF
partNew%n_in = .TRUE. partNew%n_in = .TRUE.
partNew%weight = species(sp)%obj%weight partNew%weight = species(sp)%obj%weight
!If charged species, add qm to particle
SELECT TYPE(sp => species(sp)%obj)
TYPE IS (speciesCharged)
partNew%qm = sp%qm
CLASS DEFAULT
partNew%qm = 0.D0
END SELECT
!Assign particle to temporal list of particles !Assign particle to temporal list of particles
CALL partInitial%add(partNew) CALL partInitial%add(partNew)
@ -642,11 +633,13 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: crossSecFilePath CHARACTER(:), ALLOCATABLE:: crossSecFilePath
CHARACTER(:), ALLOCATABLE:: cType CHARACTER(:), ALLOCATABLE:: cType
LOGICAL:: found LOGICAL:: found
INTEGER:: nInteractions, nCollisions INTEGER:: nPairs, nCollisions
INTEGER:: i, k, ij INTEGER:: i, k, ij
INTEGER:: pt_i, pt_j INTEGER:: pt_i, pt_j
REAL(8):: energyThreshold, energyBinding REAL(8):: energyThreshold, energyBinding
CHARACTER(:), ALLOCATABLE:: electron CHARACTER(:), ALLOCATABLE:: electron
INTEGER:: e
CLASS(meshVol), POINTER:: vol
!Firstly, checks if the object 'interactions' exists !Firstly, checks if the object 'interactions' exists
CALL config%info('interactions', found) CALL config%info('interactions', found)
@ -681,8 +674,8 @@ MODULE moduleInput
!Inits lock for list of particles !Inits lock for list of particles
CALL OMP_INIT_LOCK(lockCollisions) CALL OMP_INIT_LOCK(lockCollisions)
CALL config%info('interactions.collisions', found, n_children = nInteractions) CALL config%info('interactions.collisions', found, n_children = nPairs)
DO i = 1, nInteractions DO i = 1, nPairs
WRITE(iString, '(I2)') i WRITE(iString, '(I2)') i
object = 'interactions.collisions(' // TRIM(iString) // ')' object = 'interactions.collisions(' // TRIM(iString) // ')'
CALL config%get(object // '.species_i', species_i, found) CALL config%get(object // '.species_i', species_i, found)
@ -691,6 +684,7 @@ MODULE moduleInput
pt_j = speciesName2Index(species_j) pt_j = speciesName2Index(species_j)
CALL config%info(object // '.cTypes', found, n_children = nCollisions) CALL config%info(object // '.cTypes', found, n_children = nCollisions)
ij = interactionIndex(pt_i,pt_j) ij = interactionIndex(pt_i,pt_j)
!Allocates the required number of collisions per each pair of species ij !Allocates the required number of collisions per each pair of species ij
CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j) CALL interactionMatrix(ij)%init(nCollisions, pt_i, pt_j)
@ -732,7 +726,7 @@ MODULE moduleInput
crossSecFilePath, energyBinding, electron) crossSecFilePath, energyBinding, electron)
CASE DEFAULT CASE DEFAULT
CALL criticalError('Collision type' // cType // 'not defined yet', 'readInteractions') CALL criticalError('Collision type' // cType // 'not defined', 'readInteractions')
END SELECT END SELECT
@ -740,6 +734,26 @@ MODULE moduleInput
END DO END DO
!Init the required arrays in each volume to account for MCC.
DO e = 1, meshForMCC%numVols
vol => meshForMCC%vols(e)%obj
!Allocate Maximum cross section per collision pair and assign the initial collision rate
ALLOCATE(vol%sigmaVrelMax(1:nCollPairs))
vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
IF (collOutput) THEN
ALLOCATE(vol%tallyColl(1:nCollPairs))
DO k = 1, nCollPairs
ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
vol%tallyColl(k)%tally = 0
END DO
END IF
END DO
END IF END IF
END IF END IF

14
src/modules/moduleSolver.f90
View file

@ -165,7 +165,7 @@ MODULE moduleSolver
INTEGER:: sp INTEGER:: sp
!$OMP DO !$OMP DO
DO n=1, nPartOld DO n = 1, nPartOld
!Select species type !Select species type
sp = partOld(n)%species%n sp = partOld(n)%species%n
!Checks if the species sp is update this iteration !Checks if the species sp is update this iteration
@ -204,7 +204,7 @@ MODULE moduleSolver
REAL(8):: qmEFt(1:3) REAL(8):: qmEFt(1:3)
!Get the electric field at particle position !Get the electric field at particle position
qmEFt = part%qm*gatherElecField(part)*tauIn qmEFt = part%species%qm*gatherElecField(part)*tauIn
!Update velocity !Update velocity
part%v = part%v + qmEFt part%v = part%v + qmEFt
@ -230,7 +230,7 @@ MODULE moduleSolver
tauInHalf = tauIn *0.5D0 tauInHalf = tauIn *0.5D0
!Half of the force o f the electric field !Half of the force o f the electric field
qmEFt = part%qm*gatherElecField(part)*tauInHalf qmEFt = part%species%qm*gatherElecField(part)*tauInHalf
!Half step for electrostatic !Half step for electrostatic
v_minus = part%v + qmEFt v_minus = part%v + qmEFt
@ -239,7 +239,7 @@ MODULE moduleSolver
B = gatherMagnField(part) B = gatherMagnField(part)
BNorm = NORM2(B) BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN IF (BNorm > 0.D0) THEN
fn = DTAN(part%qm * tauInHalf*BNorm) / Bnorm fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
v_prime = v_minus + fn * crossProduct(v_minus, B) v_prime = v_minus + fn * crossProduct(v_minus, B)
v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B) v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B)
@ -305,7 +305,7 @@ MODULE moduleSolver
part_temp = part part_temp = part
!Get electric field at particle position !Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauIn
!z !z
part_temp%v(1) = part%v(1) + qmEFt(1) part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn
@ -381,7 +381,7 @@ MODULE moduleSolver
part_temp = part part_temp = part
!Get electric field at particle position !Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin qmEFt = part_temp%species%qm*gatherElecField(part_temp)*tauMin
!r,theta !r,theta
v_p_oh_star(1) = part%v(1) + qmEFt(1) v_p_oh_star(1) = part%v(1) + qmEFt(1)
x_new = part%r(1) + v_p_oh_star(1)*tauIn x_new = part%r(1) + v_p_oh_star(1)*tauIn
@ -586,7 +586,7 @@ MODULE moduleSolver
!Loops over the particles to scatter them !Loops over the particles to scatter them
!$OMP DO !$OMP DO
DO n=1, nPartOld DO n = 1, nPartOld
CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n))
END DO END DO

5
src/modules/moduleSpecies.f90
View file

@ -6,7 +6,7 @@ MODULE moduleSpecies
TYPE, ABSTRACT:: speciesGeneric TYPE, ABSTRACT:: speciesGeneric
CHARACTER(:), ALLOCATABLE:: name CHARACTER(:), ALLOCATABLE:: name
REAL(8):: m=0.D0, weight=0.D0 REAL(8):: m=0.D0, weight=0.D0, qm=0.D0
INTEGER:: n=0 INTEGER:: n=0
END TYPE speciesGeneric END TYPE speciesGeneric
@ -18,7 +18,7 @@ MODULE moduleSpecies
END TYPE speciesNeutral END TYPE speciesNeutral
TYPE, EXTENDS(speciesGeneric):: speciesCharged TYPE, EXTENDS(speciesGeneric):: speciesCharged
REAL(8):: q=0.D0, qm=0.D0 REAL(8):: q=0.D0
CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL() CLASS(speciesGeneric), POINTER:: ion => NULL(), neutral => NULL()
CONTAINS CONTAINS
PROCEDURE, PASS:: ionize => ionizeCharged PROCEDURE, PASS:: ionize => ionizeCharged
@ -43,7 +43,6 @@ MODULE moduleSpecies
REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p. REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p.
LOGICAL:: n_in !Flag that indicates if a particle is in the domain LOGICAL:: n_in !Flag that indicates if a particle is in the domain
REAL(8):: weight=0.D0 !weight of particle REAL(8):: weight=0.D0 !weight of particle
REAL(8):: qm = 0.D0 !charge over mass
END TYPE particle END TYPE particle