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

Reducing overhead when no collisions are present

Browse source 
Particles are added to lists only if there are MCC collisions. Hopefully
this will reduce overhead when OpenMP is used and no collisions are
active.
This commit is contained in:
Jorge Gonzalez 2023-01-07 12:12:37 +01:00
commit 7ce1b7a4dd
13 changed files with 142 additions and 114 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/fpakc.f90
View file

@ -74,7 +74,10 @@ PROGRAM fpakc
tColl = omp_get_wtime() tColl = omp_get_wtime()
!$OMP END SINGLE !$OMP END SINGLE
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions(t) IF (doMCC) THEN
CALL meshForMCC%doCollisions(t)
END IF
!$OMP SINGLE !$OMP SINGLE
tColl = omp_get_wtime() - tColl tColl = omp_get_wtime() - tColl
@ -83,7 +86,10 @@ PROGRAM fpakc
tCoul = omp_get_wTime() tCoul = omp_get_wTime()
!$OMP END SINGLE !$OMP END SINGLE
IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb() IF (ASSOCIATED(mesh%doCoulomb)) THEN
CALL mesh%doCoulomb()
END IF
!$OMP SINGLE !$OMP SINGLE
tCoul = omp_get_wTime() - tCoul tCoul = omp_get_wTime() - tCoul

37
src/modules/init/moduleInput.f90
View file

@ -338,7 +338,7 @@ MODULE moduleInput
!Mean velocity and temperature at particle position !Mean velocity and temperature at particle position
REAL(8):: velocityXi(1:3), temperatureXi REAL(8):: velocityXi(1:3), temperatureXi
INTEGER:: nNewPart = 0.D0 INTEGER:: nNewPart = 0.D0
CLASS(meshCell), POINTER:: vol CLASS(meshCell), POINTER:: cell
TYPE(particle), POINTER:: partNew TYPE(particle), POINTER:: partNew
REAL(8):: vTh REAL(8):: vTh
TYPE(lNode), POINTER:: partCurr, partNext TYPE(lNode), POINTER:: partCurr, partNext
@ -394,12 +394,12 @@ MODULE moduleInput
partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian()
partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian()
partNew%vol = e partNew%cell = e
IF (doubleMesh) THEN IF (doubleMesh) THEN
partNew%volColl = findCellBrute(meshColl, partNew%r) partNew%cellColl = findCellBrute(meshColl, partNew%r)
ELSE ELSE
partNew%volColl = partNew%vol partNew%cellColl = partNew%cell
END IF END IF
@ -411,11 +411,14 @@ MODULE moduleInput
CALL partInitial%add(partNew) CALL partInitial%add(partNew)
!Assign particle to list in volume !Assign particle to list in volume
vol => meshforMCC%cells(partNew%volColl)%obj IF (doMCC) THEN
CALL OMP_SET_LOCK(vol%lock) cell => meshforMCC%cells(partNew%cellColl)%obj
CALL vol%listPart_in(sp)%add(partNew) CALL OMP_SET_LOCK(cell%lock)
vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight CALL cell%listPart_in(sp)%add(partNew)
CALL OMP_UNSET_LOCK(vol%lock) cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight
CALL OMP_UNSET_LOCK(cell%lock)
END IF
END DO END DO
@ -628,7 +631,7 @@ MODULE moduleInput
REAL(8):: energyThreshold, energyBinding REAL(8):: energyThreshold, energyBinding
CHARACTER(:), ALLOCATABLE:: electron CHARACTER(:), ALLOCATABLE:: electron
INTEGER:: e INTEGER:: e
CLASS(meshCell), POINTER:: vol CLASS(meshCell), POINTER:: cell
!Firstly, check if the object 'interactions' exists !Firstly, check if the object 'interactions' exists
CALL config%info('interactions', found) CALL config%info('interactions', found)
@ -725,17 +728,17 @@ MODULE moduleInput
!Init the required arrays in each volume to account for MCC. !Init the required arrays in each volume to account for MCC.
DO e = 1, meshForMCC%numCells DO e = 1, meshForMCC%numCells
vol => meshForMCC%cells(e)%obj cell => meshForMCC%cells(e)%obj
!Allocate Maximum cross section per collision pair and assign the initial collision rate !Allocate Maximum cross section per collision pair and assign the initial collision rate
ALLOCATE(vol%sigmaVrelMax(1:nCollPairs)) ALLOCATE(cell%sigmaVrelMax(1:nCollPairs))
vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref)
IF (collOutput) THEN IF (collOutput) THEN
ALLOCATE(vol%tallyColl(1:nCollPairs)) ALLOCATE(cell%tallyColl(1:nCollPairs))
DO k = 1, nCollPairs DO k = 1, nCollPairs
ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount))
vol%tallyColl(k)%tally = 0 cell%tallyColl(k)%tally = 0
END DO END DO
@ -892,6 +895,8 @@ MODULE moduleInput
END IF END IF
doMCC = ASSOCIATED(meshForMCC)
!Get the dimension of the geometry !Get the dimension of the geometry
CALL config%get(object // '.dimension', mesh%dimen, found) CALL config%get(object // '.dimension', mesh%dimen, found)
IF (.NOT. found) THEN IF (.NOT. found) THEN

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

@ -59,7 +59,7 @@ MODULE moduleMesh1DCart
PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: neighbourElement => neighbourElementSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm
!PARTICLUAR PROCEDURES !PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeSegm PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm
END TYPE meshCell1DCartSegm END TYPE meshCell1DCartSegm
@ -193,7 +193,7 @@ MODULE moduleMesh1DCart
self%x = (/ r1(1), r2(1) /) self%x = (/ r1(1), r2(1) /)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
@ -419,7 +419,7 @@ MODULE moduleMesh1DCart
END SUBROUTINE neighbourElementSegm END SUBROUTINE neighbourElementSegm
!Compute element vol !Compute element volume
PURE SUBROUTINE volumeSegm(self) PURE SUBROUTINE volumeSegm(self)
IMPLICIT NONE IMPLICIT NONE
@ -478,10 +478,10 @@ MODULE moduleMesh1DCart
INTEGER:: e, et INTEGER:: e, et
DO e = 1, self%numCells DO e = 1, self%numCells
!Connect Vol-Vol !Connect Cell-Cell
DO et = 1, self%numCells DO et = 1, self%numCells
IF (e /= et) THEN IF (e /= et) THEN
CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj) CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)
END IF END IF
@ -489,9 +489,9 @@ MODULE moduleMesh1DCart
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
!Connect Vol-Edge !Connect Cell-Edge
DO et = 1, self%numEdges DO et = 1, self%numEdges
CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj) CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)
END DO END DO
@ -501,7 +501,7 @@ MODULE moduleMesh1DCart
END SUBROUTINE connectMesh1DCart END SUBROUTINE connectMesh1DCart
SUBROUTINE connectVolVol(elemA, elemB) SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
@ -517,7 +517,7 @@ MODULE moduleMesh1DCart
END SELECT END SELECT
END SUBROUTINE connectVolVol END SUBROUTINE connectCellCell
SUBROUTINE connectSegmSegm(elemA, elemB) SUBROUTINE connectSegmSegm(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
@ -542,7 +542,7 @@ MODULE moduleMesh1DCart
END SUBROUTINE connectSegmSegm END SUBROUTINE connectSegmSegm
SUBROUTINE connectVolEdge(elemA, elemB) SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
@ -558,7 +558,7 @@ MODULE moduleMesh1DCart
END SELECT END SELECT
END SUBROUTINE connectVolEdge END SUBROUTINE connectCellEdge
SUBROUTINE connectSegmEdge(elemA, elemB) SUBROUTINE connectSegmEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE

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

@ -59,7 +59,7 @@ MODULE moduleMesh1DRad
PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: neighbourElement => neighbourElementSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm
!PARTICLUAR PROCEDURES !PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeSegm PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm
END TYPE meshCell1DRadSegm END TYPE meshCell1DRadSegm
@ -193,7 +193,7 @@ MODULE moduleMesh1DRad
self%r = (/ r1(1), r2(1) /) self%r = (/ r1(1), r2(1) /)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
@ -427,7 +427,7 @@ MODULE moduleMesh1DRad
END SUBROUTINE neighbourElementSegm END SUBROUTINE neighbourElementSegm
!Compute element vol !Compute element volume
PURE SUBROUTINE volumeSegm(self) PURE SUBROUTINE volumeSegm(self)
USE moduleConstParam, ONLY: PI4 USE moduleConstParam, ONLY: PI4
IMPLICIT NONE IMPLICIT NONE
@ -493,10 +493,10 @@ MODULE moduleMesh1DRad
INTEGER:: e, et INTEGER:: e, et
DO e = 1, self%numCells DO e = 1, self%numCells
!Connect Vol-Vol !Connect Cell-Cell
DO et = 1, self%numCells DO et = 1, self%numCells
IF (e /= et) THEN IF (e /= et) THEN
CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj) CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)
END IF END IF
@ -504,9 +504,9 @@ MODULE moduleMesh1DRad
SELECT TYPE(self) SELECT TYPE(self)
TYPE IS(meshParticles) TYPE IS(meshParticles)
!Connect Vol-Edge !Connect Cell-Edge
DO et = 1, self%numEdges DO et = 1, self%numEdges
CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj) CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)
END DO END DO
@ -516,7 +516,7 @@ MODULE moduleMesh1DRad
END SUBROUTINE connectMesh1DRad END SUBROUTINE connectMesh1DRad
SUBROUTINE connectVolVol(elemA, elemB) SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
@ -532,7 +532,7 @@ MODULE moduleMesh1DRad
END SELECT END SELECT
END SUBROUTINE connectVolVol END SUBROUTINE connectCellCell
SUBROUTINE connectSegmSegm(elemA, elemB) SUBROUTINE connectSegmSegm(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
@ -557,7 +557,7 @@ MODULE moduleMesh1DRad
END SUBROUTINE connectSegmSegm END SUBROUTINE connectSegmSegm
SUBROUTINE connectVolEdge(elemA, elemB) SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA CLASS(meshCell), INTENT(inout):: elemA
@ -573,7 +573,7 @@ MODULE moduleMesh1DRad
END SELECT END SELECT
END SUBROUTINE connectVolEdge END SUBROUTINE connectCellEdge
SUBROUTINE connectSegmEdge(elemA, elemB) SUBROUTINE connectSegmEdge(elemA, elemB)
IMPLICIT NONE IMPLICIT NONE

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

@ -66,7 +66,7 @@ MODULE moduleMesh2DCart
PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: neighbourElement => neighbourElementQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad
!PARTICLUAR PROCEDURES !PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeQuad PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad
END TYPE meshCell2DCartQuad END TYPE meshCell2DCartQuad
@ -97,7 +97,7 @@ MODULE moduleMesh2DCart
PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: neighbourElement => neighbourElementTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria
!PARTICULAR PROCEDURES !PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeTria PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria
END TYPE meshCell2DCartTria END TYPE meshCell2DCartTria
@ -267,7 +267,7 @@ MODULE moduleMesh2DCart
self%y = (/r1(2), r2(2), r3(2), r4(2)/) self%y = (/r1(2), r2(2), r3(2), r4(2)/)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
@ -609,7 +609,7 @@ MODULE moduleMesh2DCart
self%x = (/r1(1), r2(1), r3(1)/) self%x = (/r1(1), r2(1), r3(1)/)
self%y = (/r1(2), r2(2), r3(2)/) self%y = (/r1(2), r2(2), r3(2)/)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)

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

@ -66,7 +66,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: phy2log => phy2logQuad
PROCEDURE, PASS:: neighbourElement => neighbourElementQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad
!PARTICLUAR PROCEDURES !PARTICLUAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeQuad PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad
END TYPE meshCell2DCylQuad END TYPE meshCell2DCylQuad
@ -97,7 +97,7 @@ MODULE moduleMesh2DCyl
PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: phy2log => phy2logTria
PROCEDURE, PASS:: neighbourElement => neighbourElementTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria
!PARTICULAR PROCEDURES !PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeTria PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria
END TYPE meshCell2DCylTria END TYPE meshCell2DCylTria
@ -275,7 +275,7 @@ MODULE moduleMesh2DCyl
self%r = (/r1(2), r2(2), r3(2), r4(2)/) self%r = (/r1(2), r2(2), r3(2), r4(2)/)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)
@ -636,7 +636,7 @@ MODULE moduleMesh2DCyl
self%z = (/r1(1), r2(1), r3(1)/) self%z = (/r1(1), r2(1), r3(1)/)
self%r = (/r1(2), r2(2), r3(2)/) self%r = (/r1(2), r2(2), r3(2)/)
!Assign node volume !Assign node volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)

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

@ -60,7 +60,7 @@ MODULE moduleMesh3DCart
PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: phy2log => phy2logTetra
PROCEDURE, PASS:: neighbourElement => neighbourElementTetra PROCEDURE, PASS:: neighbourElement => neighbourElementTetra
!PARTICULAR PROCEDURES !PARTICULAR PROCEDURES
PROCEDURE, PASS, PRIVATE:: vol => volumeTetra PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTetra
END TYPE meshCell3DCartTetra END TYPE meshCell3DCartTetra
@ -255,7 +255,7 @@ MODULE moduleMesh3DCart
self%z = (/r1(3), r2(3), r3(3), r4(3)/) self%z = (/r1(3), r2(3), r3(3), r4(3)/)
!Computes the element volume !Computes the element volume
CALL self%vol() CALL self%calculateVolume()
CALL OMP_INIT_LOCK(self%lock) CALL OMP_INIT_LOCK(self%lock)

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

@ -480,6 +480,8 @@ MODULE moduleMesh
!Logical to indicate if an specific mesh for MC Collisions is used !Logical to indicate if an specific mesh for MC Collisions is used
LOGICAL:: doubleMesh LOGICAL:: doubleMesh
!Logical to indicate if MCC collisions are performed
LOGICAL:: doMCC
!Complete path for the two meshes !Complete path for the two meshes
CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle
@ -644,15 +646,18 @@ MODULE moduleMesh
Xi = self%phy2log(part%r) Xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell !Checks if particle is inside 'self' cell
IF (self%inside(Xi)) THEN IF (self%inside(Xi)) THEN
part%vol = self%n part%cell = self%n
part%Xi = Xi part%Xi = Xi
part%n_in = .TRUE. part%n_in = .TRUE.
!Assign particle to listPart_in !Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock) IF (doMCC) THEN
sp = part%species%n CALL OMP_SET_LOCK(self%lock)
CALL self%listPart_in(sp)%add(part) sp = part%species%n
self%totalWeight(sp) = self%totalWeight(sp) + part%weight CALL self%listPart_in(sp)%add(part)
CALL OMP_UNSET_LOCK(self%lock) self%totalWeight(sp) = self%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(self%lock)
END IF
ELSE ELSE
!If not, searches for a neighbour and repeats the process. !If not, searches for a neighbour and repeats the process.
@ -688,14 +693,14 @@ MODULE moduleMesh
END SUBROUTINE findCell END SUBROUTINE findCell
!If Coll and Particle are the same, simply copy the part%vol into part%volColl !If Coll and Particle are the same, simply copy the part%cell into part%cellColl
SUBROUTINE findCellSameMesh(part) SUBROUTINE findCellSameMesh(part)
USE moduleSpecies USE moduleSpecies
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
part%volColl = part%vol part%cellColl = part%cell
END SUBROUTINE findCellSameMesh END SUBROUTINE findCellSameMesh
@ -715,16 +720,19 @@ MODULE moduleMesh
found = .FALSE. found = .FALSE.
cell => meshColl%cells(part%volColl)%obj cell => meshColl%cells(part%cellColl)%obj
DO WHILE(.NOT. found) DO WHILE(.NOT. found)
Xi = cell%phy2log(part%r) Xi = cell%phy2log(part%r)
IF (cell%inside(Xi)) THEN IF (cell%inside(Xi)) THEN
part%volColl = cell%n part%cellColl = cell%n
CALL OMP_SET_LOCK(cell%lock) IF (doMCC) THEN
sp = part%species%n CALL OMP_SET_LOCK(cell%lock)
CALL cell%listPart_in(sp)%add(part) sp = part%species%n
cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight CALL cell%listPart_in(sp)%add(part)
CALL OMP_UNSET_LOCK(cell%lock) cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight
CALL OMP_UNSET_LOCK(cell%lock)
END IF
found = .TRUE. found = .TRUE.
ELSE ELSE

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

@ -156,10 +156,10 @@ MODULE moduleMeshBoundary
newElectron%r = edge%randPos() newElectron%r = edge%randPos()
newIon%r = newElectron%r newIon%r = newElectron%r
newElectron%vol = part%vol newElectron%cell = part%cell
newIon%vol = part%vol newIon%cell = part%cell
newElectron%Xi = mesh%cells(part%vol)%obj%phy2log(newElectron%r) newElectron%Xi = mesh%cells(part%cell)%obj%phy2log(newElectron%r)
newIon%Xi = newElectron%Xi newIon%Xi = newElectron%Xi
newElectron%weight = part%weight newElectron%weight = part%weight

8
src/modules/moduleInject.f90
View file

@ -285,16 +285,16 @@ MODULE moduleInject
partInj(n)%r = randomEdge%randPos() partInj(n)%r = randomEdge%randPos()
!Volume associated to the edge: !Volume associated to the edge:
IF (ASSOCIATED(randomEdge%e1)) THEN IF (ASSOCIATED(randomEdge%e1)) THEN
partInj(n)%vol = randomEdge%e1%n partInj(n)%cell = randomEdge%e1%n
ELSEIF (ASSOCIATED(randomEdge%e2)) THEN ELSEIF (ASSOCIATED(randomEdge%e2)) THEN
partInj(n)%vol = randomEdge%e2%n partInj(n)%cell = randomEdge%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)%volColl = randomEdge%eColl%n partInj(n)%cellColl = randomEdge%eColl%n
sp = self%species%n sp = self%species%n
!Assign particle type !Assign particle type
@ -305,7 +305,7 @@ MODULE moduleInject
self%v(3)%obj%randomVel() /) self%v(3)%obj%randomVel() /)
!Obtain natural coordinates of particle in cell !Obtain natural coordinates of particle in cell
partInj(n)%Xi = mesh%cells(partInj(n)%vol)%obj%phy2log(partInj(n)%r) partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
!Push new particle with the minimum time step !Push new particle with the minimum time step
CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp))
!Assign cell to new particle !Assign cell to new particle

4
src/modules/moduleSpecies.f90
View file

@ -38,8 +38,8 @@ MODULE moduleSpecies
REAL(8):: r(1:3) !Position REAL(8):: r(1:3) !Position
REAL(8):: v(1:3) !Velocity REAL(8):: v(1:3) !Velocity
CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle
INTEGER:: vol !Index of element in which the particle is located INTEGER:: cell !Index of element in which the particle is located
INTEGER:: volColl !Index of element in which the particle is located in the Collision Mesh INTEGER:: cellColl !Index of element in which the particle is located in the Collision Mesh
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

79
src/modules/solver/moduleSolver.f90
View file

@ -43,14 +43,14 @@ MODULE moduleSolver
END SUBROUTINE solveEM_interface END SUBROUTINE solveEM_interface
!Apply nonAnalogue scheme to a particle !Apply nonAnalogue scheme to a particle
SUBROUTINE weightingScheme_interface(part, volOld, volNew) SUBROUTINE weightingScheme_interface(part, cellOld, cellNew)
USE moduleSpecies USE moduleSpecies
USE moduleMesh USE moduleMesh
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshCell), POINTER, INTENT(in):: volOld CLASS(meshCell), POINTER, INTENT(in):: cellOld
CLASS(meshCell), POINTER, INTENT(inout):: volNew CLASS(meshCell), POINTER, INTENT(inout):: cellNew
END SUBROUTINE weightingScheme_interface END SUBROUTINE weightingScheme_interface
@ -322,31 +322,37 @@ MODULE moduleSolver
!$OMP SECTION !$OMP SECTION
!Erase the list of particles inside the cell if particles have been pushed !Erase the list of particles inside the cell if particles have been pushed
DO s = 1, nSpecies IF (doMCC) THEN
DO e = 1, mesh%numCells DO s = 1, nSpecies
IF (solver%pusher(s)%pushSpecies) THEN DO e = 1, mesh%numCells
CALL mesh%cells(e)%obj%listPart_in(s)%erase() IF (solver%pusher(s)%pushSpecies) THEN
mesh%cells(e)%obj%totalWeight(s) = 0.D0 CALL mesh%cells(e)%obj%listPart_in(s)%erase()
mesh%cells(e)%obj%totalWeight(s) = 0.D0
END IF END IF
END DO
END DO END DO
END DO END IF
!$OMP SECTION !$OMP SECTION
!Erase the list of particles inside the cell in coll mesh !Erase the list of particles inside the cell in coll mesh
DO s = 1, nSpecies IF (doubleMesh) THEN
DO e = 1, meshColl%numCells DO s = 1, nSpecies
IF (solver%pusher(s)%pushSpecies) THEN DO e = 1, meshColl%numCells
CALL meshColl%cells(e)%obj%listPart_in(s)%erase() IF (solver%pusher(s)%pushSpecies) THEN
meshColl%cells(e)%obj%totalWeight(s) = 0.D0 CALL meshColl%cells(e)%obj%listPart_in(s)%erase()
meshColl%cells(e)%obj%totalWeight(s) = 0.D0
END IF END IF
END DO
END DO END DO
END DO END IF
!$OMP END SECTIONS !$OMP END SECTIONS
@ -368,7 +374,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
cell => mesh%cells(partOld(n)%vol)%obj cell => mesh%cells(partOld(n)%cell)%obj
CALL cell%scatter(cell%nNodes, partOld(n)) CALL cell%scatter(cell%nNodes, partOld(n))
END DO END DO
@ -387,28 +393,28 @@ MODULE moduleSolver
END SUBROUTINE doEMField END SUBROUTINE doEMField
!Split particles as a function of cell volume and splits particle !Split particles as a function of cell volume and splits particle
SUBROUTINE volumeWScheme(part, volOld, volNew) SUBROUTINE volumeWScheme(part, cellOld, cellNew)
USE moduleSpecies USE moduleSpecies
USE moduleMesh USE moduleMesh
USE moduleRandom USE moduleRandom
IMPLICIT NONE IMPLICIT NONE
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshCell), POINTER, INTENT(in):: volOld CLASS(meshCell), POINTER, INTENT(in):: cellOld
CLASS(meshCell), POINTER, INTENT(inout):: volNew CLASS(meshCell), POINTER, INTENT(inout):: cellNew
REAL(8):: fractionVolume, pSplit REAL(8):: fractionVolume, pSplit
!If particle changes volume to smaller cell !If particle changes volume to smaller cell
IF (volOld%volume > volNew%volume .AND. & IF (cellOld%volume > cellNew%volume .AND. &
part%weight >= part%species%weight*1.0D-1) THEN part%weight >= part%species%weight*1.0D-1) THEN
fractionVolume = volOld%volume/volNew%volume fractionVolume = cellOld%volume/cellNew%volume
!Calculate probability of splitting particle !Calculate probability of splitting particle
pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1) pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1)
IF (random() < pSplit) THEN IF (random() < pSplit) THEN
!Split particle in two !Split particle in two
CALL splitParticle(part, 2, volNew) CALL splitParticle(part, 2, cellNew)
END IF END IF
@ -418,7 +424,7 @@ MODULE moduleSolver
!Subroutine to split the particle 'part' into a number 'nSplit' of particles. !Subroutine to split the particle 'part' into a number 'nSplit' of particles.
!'nSplit-1' particles are added to the partNAScheme list !'nSplit-1' particles are added to the partNAScheme list
SUBROUTINE splitParticle(part, nSplit, vol) SUBROUTINE splitParticle(part, nSplit, cell)
USE moduleSpecies USE moduleSpecies
USE moduleList USE moduleList
USE moduleMesh USE moduleMesh
@ -427,7 +433,7 @@ MODULE moduleSolver
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
INTEGER, INTENT(in):: nSplit INTEGER, INTENT(in):: nSplit
CLASS(meshCell), INTENT(inout):: vol CLASS(meshCell), INTENT(inout):: cell
REAL(8):: newWeight REAL(8):: newWeight
TYPE(particle), POINTER:: newPart TYPE(particle), POINTER:: newPart
INTEGER:: p INTEGER:: p
@ -449,10 +455,13 @@ MODULE moduleSolver
CALL partWScheme%add(newPart) CALL partWScheme%add(newPart)
CALL partWScheme%unsetLock() CALL partWScheme%unsetLock()
!Add particle to cell list !Add particle to cell list
CALL OMP_SET_lock(vol%lock)
sp = part%species%n sp = part%species%n
CALL vol%listPart_in(sp)%add(newPart) IF (doMCC) THEN
CALL OMP_UNSET_lock(vol%lock) CALL OMP_SET_lock(cell%lock)
CALL cell%listPart_in(sp)%add(newPart)
CALL OMP_UNSET_lock(cell%lock)
END IF
END DO END DO
@ -465,18 +474,18 @@ MODULE moduleSolver
CLASS(solverGeneric), INTENT(in):: self CLASS(solverGeneric), INTENT(in):: self
TYPE(particle), INTENT(inout):: part TYPE(particle), INTENT(inout):: part
CLASS(meshCell), POINTER:: volOld, volNew CLASS(meshCell), POINTER:: cellOld, cellNew
!Assume that particle is outside the domain !Assume that particle is outside the domain
part%n_in = .FALSE. part%n_in = .FALSE.
volOld => mesh%cells(part%vol)%obj cellOld => mesh%cells(part%cell)%obj
CALL volOld%findCell(part) CALL cellOld%findCell(part)
CALL findCellColl(part) CALL findCellColl(part)
!Call the NA shcme !Call the NA shcme
IF (ASSOCIATED(self%weightingScheme)) THEN IF (ASSOCIATED(self%weightingScheme)) THEN
volNew => mesh%cells(part%vol)%obj cellNew => mesh%cells(part%cell)%obj
CALL self%weightingScheme(part, volOld, volNew) CALL self%weightingScheme(part, cellOld, cellNew)
END IF END IF

10
src/modules/solver/pusher/modulePusher.f90
View file

@ -23,7 +23,7 @@ MODULE modulePusher
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 = mesh%cells(part%vol)%obj%gatherElectricField(part%Xi) qmEFT = mesh%cells(part%cell)%obj%gatherElectricField(part%Xi)
qmEFt = qmEFt*part%species%qm*tauMin qmEFt = qmEFt*part%species%qm*tauMin
!Update velocity !Update velocity
@ -50,14 +50,14 @@ MODULE modulePusher
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 = mesh%cells(part%vol)%obj%gatherElectricField(part%Xi) qmEFT = mesh%cells(part%cell)%obj%gatherElectricField(part%Xi)
qmEFt = qmEFt*part%species%qm*tauInHalf qmEFt = qmEFt*part%species%qm*tauInHalf
!Half step for electrostatic !Half step for electrostatic
v_minus = part%v + qmEFt v_minus = part%v + qmEFt
!Full step rotation !Full step rotation
B = mesh%cells(part%vol)%obj%gatherMagneticField(part%Xi) B = mesh%cells(part%cell)%obj%gatherMagneticField(part%Xi)
BNorm = NORM2(B) BNorm = NORM2(B)
IF (BNorm > 0.D0) THEN IF (BNorm > 0.D0) THEN
fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm
@ -126,7 +126,7 @@ MODULE modulePusher
part_temp = part part_temp = part
!Get electric field at particle position !Get electric field at particle position
qmEFT = mesh%cells(part_temp%vol)%obj%gatherElectricField(part_temp%Xi) qmEFT = mesh%cells(part_temp%cell)%obj%gatherElectricField(part_temp%Xi)
qmEFt = qmEFt*part_temp%species%qm*tauMin qmEFt = qmEFt*part_temp%species%qm*tauMin
!z !z
part_temp%v(1) = part%v(1) + qmEFt(1) part_temp%v(1) = part%v(1) + qmEFt(1)
@ -202,7 +202,7 @@ MODULE modulePusher
part_temp = part part_temp = part
!Get electric field at particle position !Get electric field at particle position
qmEFT = mesh%cells(part_temp%vol)%obj%gatherElectricField(part_temp%Xi) qmEFT = mesh%cells(part_temp%cell)%obj%gatherElectricField(part_temp%Xi)
qmEFt = qmEFt*part_temp%species%qm*tauMin qmEFt = qmEFt*part_temp%species%qm*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)