From 7ce1b7a4dd446919d6b5686d80493dab334a487d Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 7 Jan 2023 12:12:37 +0100 Subject: [PATCH] Reducing overhead when no collisions are present 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. --- src/fpakc.f90 | 10 ++- src/modules/init/moduleInput.f90 | 37 +++++---- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 22 +++--- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 22 +++--- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 8 +- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 8 +- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 4 +- src/modules/mesh/moduleMesh.f90 | 38 ++++++---- src/modules/mesh/moduleMeshBoundary.f90 | 6 +- src/modules/moduleInject.f90 | 8 +- src/modules/moduleSpecies.f90 | 4 +- src/modules/solver/moduleSolver.f90 | 79 +++++++++++--------- src/modules/solver/pusher/modulePusher.f90 | 10 +-- 13 files changed, 142 insertions(+), 114 deletions(-) diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 1224f53..3a8a033 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -74,7 +74,10 @@ PROGRAM fpakc tColl = omp_get_wtime() !$OMP END SINGLE - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%doCollisions(t) + IF (doMCC) THEN + CALL meshForMCC%doCollisions(t) + + END IF !$OMP SINGLE tColl = omp_get_wtime() - tColl @@ -83,7 +86,10 @@ PROGRAM fpakc tCoul = omp_get_wTime() !$OMP END SINGLE - IF (ASSOCIATED(mesh%doCoulomb)) CALL mesh%doCoulomb() + IF (ASSOCIATED(mesh%doCoulomb)) THEN + CALL mesh%doCoulomb() + + END IF !$OMP SINGLE tCoul = omp_get_wTime() - tCoul diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 60f33ec..81ababd 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -338,7 +338,7 @@ MODULE moduleInput !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi INTEGER:: nNewPart = 0.D0 - CLASS(meshCell), POINTER:: vol + CLASS(meshCell), POINTER:: cell TYPE(particle), POINTER:: partNew REAL(8):: vTh TYPE(lNode), POINTER:: partCurr, partNext @@ -394,12 +394,12 @@ MODULE moduleInput partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() - partNew%vol = e + partNew%cell = e IF (doubleMesh) THEN - partNew%volColl = findCellBrute(meshColl, partNew%r) + partNew%cellColl = findCellBrute(meshColl, partNew%r) ELSE - partNew%volColl = partNew%vol + partNew%cellColl = partNew%cell END IF @@ -411,11 +411,14 @@ MODULE moduleInput CALL partInitial%add(partNew) !Assign particle to list in volume - vol => meshforMCC%cells(partNew%volColl)%obj - CALL OMP_SET_LOCK(vol%lock) - CALL vol%listPart_in(sp)%add(partNew) - vol%totalWeight(sp) = vol%totalWeight(sp) + partNew%weight - CALL OMP_UNSET_LOCK(vol%lock) + IF (doMCC) THEN + cell => meshforMCC%cells(partNew%cellColl)%obj + CALL OMP_SET_LOCK(cell%lock) + CALL cell%listPart_in(sp)%add(partNew) + cell%totalWeight(sp) = cell%totalWeight(sp) + partNew%weight + CALL OMP_UNSET_LOCK(cell%lock) + + END IF END DO @@ -628,7 +631,7 @@ MODULE moduleInput REAL(8):: energyThreshold, energyBinding CHARACTER(:), ALLOCATABLE:: electron INTEGER:: e - CLASS(meshCell), POINTER:: vol + CLASS(meshCell), POINTER:: cell !Firstly, check if the object 'interactions' exists CALL config%info('interactions', found) @@ -725,17 +728,17 @@ MODULE moduleInput !Init the required arrays in each volume to account for MCC. 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(vol%sigmaVrelMax(1:nCollPairs)) - vol%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) + ALLOCATE(cell%sigmaVrelMax(1:nCollPairs)) + cell%sigmaVrelMax = sigmaVrel_ref/(L_ref**2 * v_ref) IF (collOutput) THEN - ALLOCATE(vol%tallyColl(1:nCollPairs)) + ALLOCATE(cell%tallyColl(1:nCollPairs)) DO k = 1, nCollPairs - ALLOCATE(vol%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) - vol%tallyColl(k)%tally = 0 + ALLOCATE(cell%tallyColl(k)%tally(1:interactionmatrix(k)%amount)) + cell%tallyColl(k)%tally = 0 END DO @@ -892,6 +895,8 @@ MODULE moduleInput END IF + doMCC = ASSOCIATED(meshForMCC) + !Get the dimension of the geometry CALL config%get(object // '.dimension', mesh%dimen, found) IF (.NOT. found) THEN diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 2f8b7b8..269f157 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -59,7 +59,7 @@ MODULE moduleMesh1DCart PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeSegm + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm END TYPE meshCell1DCartSegm @@ -193,7 +193,7 @@ MODULE moduleMesh1DCart self%x = (/ r1(1), r2(1) /) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) @@ -419,7 +419,7 @@ MODULE moduleMesh1DCart END SUBROUTINE neighbourElementSegm - !Compute element vol + !Compute element volume PURE SUBROUTINE volumeSegm(self) IMPLICIT NONE @@ -478,10 +478,10 @@ MODULE moduleMesh1DCart INTEGER:: e, et DO e = 1, self%numCells - !Connect Vol-Vol + !Connect Cell-Cell DO et = 1, self%numCells 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 @@ -489,9 +489,9 @@ MODULE moduleMesh1DCart SELECT TYPE(self) TYPE IS(meshParticles) - !Connect Vol-Edge + !Connect Cell-Edge 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 @@ -501,7 +501,7 @@ MODULE moduleMesh1DCart END SUBROUTINE connectMesh1DCart - SUBROUTINE connectVolVol(elemA, elemB) + SUBROUTINE connectCellCell(elemA, elemB) IMPLICIT NONE CLASS(meshCell), INTENT(inout):: elemA @@ -517,7 +517,7 @@ MODULE moduleMesh1DCart END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell SUBROUTINE connectSegmSegm(elemA, elemB) IMPLICIT NONE @@ -542,7 +542,7 @@ MODULE moduleMesh1DCart END SUBROUTINE connectSegmSegm - SUBROUTINE connectVolEdge(elemA, elemB) + SUBROUTINE connectCellEdge(elemA, elemB) IMPLICIT NONE CLASS(meshCell), INTENT(inout):: elemA @@ -558,7 +558,7 @@ MODULE moduleMesh1DCart END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge SUBROUTINE connectSegmEdge(elemA, elemB) IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 81a8864..d998267 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -59,7 +59,7 @@ MODULE moduleMesh1DRad PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: neighbourElement => neighbourElementSegm !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeSegm + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeSegm END TYPE meshCell1DRadSegm @@ -193,7 +193,7 @@ MODULE moduleMesh1DRad self%r = (/ r1(1), r2(1) /) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) @@ -427,7 +427,7 @@ MODULE moduleMesh1DRad END SUBROUTINE neighbourElementSegm - !Compute element vol + !Compute element volume PURE SUBROUTINE volumeSegm(self) USE moduleConstParam, ONLY: PI4 IMPLICIT NONE @@ -493,10 +493,10 @@ MODULE moduleMesh1DRad INTEGER:: e, et DO e = 1, self%numCells - !Connect Vol-Vol + !Connect Cell-Cell DO et = 1, self%numCells 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 @@ -504,9 +504,9 @@ MODULE moduleMesh1DRad SELECT TYPE(self) TYPE IS(meshParticles) - !Connect Vol-Edge + !Connect Cell-Edge 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 @@ -516,7 +516,7 @@ MODULE moduleMesh1DRad END SUBROUTINE connectMesh1DRad - SUBROUTINE connectVolVol(elemA, elemB) + SUBROUTINE connectCellCell(elemA, elemB) IMPLICIT NONE CLASS(meshCell), INTENT(inout):: elemA @@ -532,7 +532,7 @@ MODULE moduleMesh1DRad END SELECT - END SUBROUTINE connectVolVol + END SUBROUTINE connectCellCell SUBROUTINE connectSegmSegm(elemA, elemB) IMPLICIT NONE @@ -557,7 +557,7 @@ MODULE moduleMesh1DRad END SUBROUTINE connectSegmSegm - SUBROUTINE connectVolEdge(elemA, elemB) + SUBROUTINE connectCellEdge(elemA, elemB) IMPLICIT NONE CLASS(meshCell), INTENT(inout):: elemA @@ -573,7 +573,7 @@ MODULE moduleMesh1DRad END SELECT - END SUBROUTINE connectVolEdge + END SUBROUTINE connectCellEdge SUBROUTINE connectSegmEdge(elemA, elemB) IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index cba0cdf..c02078a 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -66,7 +66,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeQuad + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad END TYPE meshCell2DCartQuad @@ -97,7 +97,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeTria + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria END TYPE meshCell2DCartTria @@ -267,7 +267,7 @@ MODULE moduleMesh2DCart self%y = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) @@ -609,7 +609,7 @@ MODULE moduleMesh2DCart self%x = (/r1(1), r2(1), r3(1)/) self%y = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index c2b0674..307f71c 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -66,7 +66,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad !PARTICLUAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeQuad + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeQuad END TYPE meshCell2DCylQuad @@ -97,7 +97,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeTria + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTria END TYPE meshCell2DCylTria @@ -275,7 +275,7 @@ MODULE moduleMesh2DCyl self%r = (/r1(2), r2(2), r3(2), r4(2)/) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) @@ -636,7 +636,7 @@ MODULE moduleMesh2DCyl self%z = (/r1(1), r2(1), r3(1)/) self%r = (/r1(2), r2(2), r3(2)/) !Assign node volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index fcd4647..c451689 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -60,7 +60,7 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: neighbourElement => neighbourElementTetra !PARTICULAR PROCEDURES - PROCEDURE, PASS, PRIVATE:: vol => volumeTetra + PROCEDURE, PASS, PRIVATE:: calculateVolume => volumeTetra END TYPE meshCell3DCartTetra @@ -255,7 +255,7 @@ MODULE moduleMesh3DCart self%z = (/r1(3), r2(3), r3(3), r4(3)/) !Computes the element volume - CALL self%vol() + CALL self%calculateVolume() CALL OMP_INIT_LOCK(self%lock) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index f15f880..7482842 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -480,6 +480,8 @@ MODULE moduleMesh !Logical to indicate if an specific mesh for MC Collisions is used LOGICAL:: doubleMesh + !Logical to indicate if MCC collisions are performed + LOGICAL:: doMCC !Complete path for the two meshes CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle @@ -644,15 +646,18 @@ MODULE moduleMesh Xi = self%phy2log(part%r) !Checks if particle is inside 'self' cell IF (self%inside(Xi)) THEN - part%vol = self%n + part%cell = self%n part%Xi = Xi part%n_in = .TRUE. !Assign particle to listPart_in - CALL OMP_SET_LOCK(self%lock) - sp = part%species%n - CALL self%listPart_in(sp)%add(part) - self%totalWeight(sp) = self%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(self%lock) + IF (doMCC) THEN + CALL OMP_SET_LOCK(self%lock) + sp = part%species%n + CALL self%listPart_in(sp)%add(part) + self%totalWeight(sp) = self%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(self%lock) + + END IF ELSE !If not, searches for a neighbour and repeats the process. @@ -688,14 +693,14 @@ MODULE moduleMesh 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) USE moduleSpecies IMPLICIT NONE TYPE(particle), INTENT(inout):: part - part%volColl = part%vol + part%cellColl = part%cell END SUBROUTINE findCellSameMesh @@ -715,16 +720,19 @@ MODULE moduleMesh found = .FALSE. - cell => meshColl%cells(part%volColl)%obj + cell => meshColl%cells(part%cellColl)%obj DO WHILE(.NOT. found) Xi = cell%phy2log(part%r) IF (cell%inside(Xi)) THEN - part%volColl = cell%n - CALL OMP_SET_LOCK(cell%lock) - sp = part%species%n - CALL cell%listPart_in(sp)%add(part) - cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(cell%lock) + part%cellColl = cell%n + IF (doMCC) THEN + CALL OMP_SET_LOCK(cell%lock) + sp = part%species%n + CALL cell%listPart_in(sp)%add(part) + cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(cell%lock) + + END IF found = .TRUE. ELSE diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index 517835e..6120111 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -156,10 +156,10 @@ MODULE moduleMeshBoundary newElectron%r = edge%randPos() newIon%r = newElectron%r - newElectron%vol = part%vol - newIon%vol = part%vol + newElectron%cell = part%cell + 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 newElectron%weight = part%weight diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 050ad52..daa3846 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -285,16 +285,16 @@ MODULE moduleInject partInj(n)%r = randomEdge%randPos() !Volume associated to the edge: IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%vol = randomEdge%e1%n + partInj(n)%cell = randomEdge%e1%n ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%vol = randomEdge%e2%n + partInj(n)%cell = randomEdge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - partInj(n)%volColl = randomEdge%eColl%n + partInj(n)%cellColl = randomEdge%eColl%n sp = self%species%n !Assign particle type @@ -305,7 +305,7 @@ MODULE moduleInject self%v(3)%obj%randomVel() /) !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 CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) !Assign cell to new particle diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index ca7858c..ab08f08 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -38,8 +38,8 @@ MODULE moduleSpecies REAL(8):: r(1:3) !Position REAL(8):: v(1:3) !Velocity CLASS(speciesGeneric), POINTER:: species !Pointer to species associated with this particle - INTEGER:: vol !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:: cell !Index of element in which the particle is located + 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. LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 4ee6e7d..e557495 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -43,14 +43,14 @@ MODULE moduleSolver END SUBROUTINE solveEM_interface !Apply nonAnalogue scheme to a particle - SUBROUTINE weightingScheme_interface(part, volOld, volNew) + SUBROUTINE weightingScheme_interface(part, cellOld, cellNew) USE moduleSpecies USE moduleMesh IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshCell), POINTER, INTENT(in):: volOld - CLASS(meshCell), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: cellOld + CLASS(meshCell), POINTER, INTENT(inout):: cellNew END SUBROUTINE weightingScheme_interface @@ -322,31 +322,37 @@ MODULE moduleSolver !$OMP SECTION !Erase the list of particles inside the cell if particles have been pushed - DO s = 1, nSpecies - DO e = 1, mesh%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL mesh%cells(e)%obj%listPart_in(s)%erase() - mesh%cells(e)%obj%totalWeight(s) = 0.D0 + IF (doMCC) THEN + DO s = 1, nSpecies + DO e = 1, mesh%numCells + IF (solver%pusher(s)%pushSpecies) THEN + 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 IF !$OMP SECTION !Erase the list of particles inside the cell in coll mesh - DO s = 1, nSpecies - DO e = 1, meshColl%numCells - IF (solver%pusher(s)%pushSpecies) THEN - CALL meshColl%cells(e)%obj%listPart_in(s)%erase() - meshColl%cells(e)%obj%totalWeight(s) = 0.D0 + IF (doubleMesh) THEN + DO s = 1, nSpecies + DO e = 1, meshColl%numCells + IF (solver%pusher(s)%pushSpecies) THEN + 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 IF !$OMP END SECTIONS @@ -368,7 +374,7 @@ MODULE moduleSolver !Loops over the particles to scatter them !$OMP DO 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)) END DO @@ -387,28 +393,28 @@ MODULE moduleSolver END SUBROUTINE doEMField !Split particles as a function of cell volume and splits particle - SUBROUTINE volumeWScheme(part, volOld, volNew) + SUBROUTINE volumeWScheme(part, cellOld, cellNew) USE moduleSpecies USE moduleMesh USE moduleRandom IMPLICIT NONE TYPE(particle), INTENT(inout):: part - CLASS(meshCell), POINTER, INTENT(in):: volOld - CLASS(meshCell), POINTER, INTENT(inout):: volNew + CLASS(meshCell), POINTER, INTENT(in):: cellOld + CLASS(meshCell), POINTER, INTENT(inout):: cellNew REAL(8):: fractionVolume, pSplit !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 - fractionVolume = volOld%volume/volNew%volume + fractionVolume = cellOld%volume/cellNew%volume !Calculate probability of splitting particle pSplit = 1.D0 - DEXP(-fractionVolume*1.0D-1) IF (random() < pSplit) THEN !Split particle in two - CALL splitParticle(part, 2, volNew) + CALL splitParticle(part, 2, cellNew) END IF @@ -418,7 +424,7 @@ MODULE moduleSolver !Subroutine to split the particle 'part' into a number 'nSplit' of particles. !'nSplit-1' particles are added to the partNAScheme list - SUBROUTINE splitParticle(part, nSplit, vol) + SUBROUTINE splitParticle(part, nSplit, cell) USE moduleSpecies USE moduleList USE moduleMesh @@ -427,7 +433,7 @@ MODULE moduleSolver TYPE(particle), INTENT(inout):: part INTEGER, INTENT(in):: nSplit - CLASS(meshCell), INTENT(inout):: vol + CLASS(meshCell), INTENT(inout):: cell REAL(8):: newWeight TYPE(particle), POINTER:: newPart INTEGER:: p @@ -449,10 +455,13 @@ MODULE moduleSolver CALL partWScheme%add(newPart) CALL partWScheme%unsetLock() !Add particle to cell list - CALL OMP_SET_lock(vol%lock) sp = part%species%n - CALL vol%listPart_in(sp)%add(newPart) - CALL OMP_UNSET_lock(vol%lock) + IF (doMCC) THEN + CALL OMP_SET_lock(cell%lock) + CALL cell%listPart_in(sp)%add(newPart) + CALL OMP_UNSET_lock(cell%lock) + + END IF END DO @@ -465,18 +474,18 @@ MODULE moduleSolver CLASS(solverGeneric), INTENT(in):: self TYPE(particle), INTENT(inout):: part - CLASS(meshCell), POINTER:: volOld, volNew + CLASS(meshCell), POINTER:: cellOld, cellNew !Assume that particle is outside the domain part%n_in = .FALSE. - volOld => mesh%cells(part%vol)%obj - CALL volOld%findCell(part) + cellOld => mesh%cells(part%cell)%obj + CALL cellOld%findCell(part) CALL findCellColl(part) !Call the NA shcme IF (ASSOCIATED(self%weightingScheme)) THEN - volNew => mesh%cells(part%vol)%obj - CALL self%weightingScheme(part, volOld, volNew) + cellNew => mesh%cells(part%cell)%obj + CALL self%weightingScheme(part, cellOld, cellNew) END IF diff --git a/src/modules/solver/pusher/modulePusher.f90 b/src/modules/solver/pusher/modulePusher.f90 index c2aa46a..bc02912 100644 --- a/src/modules/solver/pusher/modulePusher.f90 +++ b/src/modules/solver/pusher/modulePusher.f90 @@ -23,7 +23,7 @@ MODULE modulePusher REAL(8):: qmEFt(1:3) !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 !Update velocity @@ -50,14 +50,14 @@ MODULE modulePusher tauInHalf = tauIn *0.5D0 !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 !Half step for electrostatic v_minus = part%v + qmEFt !Full step rotation - B = mesh%cells(part%vol)%obj%gatherMagneticField(part%Xi) + B = mesh%cells(part%cell)%obj%gatherMagneticField(part%Xi) BNorm = NORM2(B) IF (BNorm > 0.D0) THEN fn = DTAN(part%species%qm * tauInHalf*BNorm) / BNorm @@ -126,7 +126,7 @@ MODULE modulePusher part_temp = part !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 !z part_temp%v(1) = part%v(1) + qmEFt(1) @@ -202,7 +202,7 @@ MODULE modulePusher part_temp = part !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 !r,theta v_p_oh_star(1) = part%v(1) + qmEFt(1)