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)