diff --git a/src/fpakc.f90 b/src/fpakc.f90 index cb68754..feee89a 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -1,39 +1,15 @@ ! FPAKC main program PROGRAM fpakc - use moduleCompTime, only: tStep, & - tEMField, & - tCoul, & - tColl, & - tPush, & - tReset, & - tWeight - use omp_lib, only: omp_get_wtime - use moduleErrors, only: criticalError, & - verboseError - use moduleInput, only: readCOnfig, & - initOutput - use moduleCaseParam, only: timeStep, & - tInitial, & - tFinal - use moduleProbe, only: resetProbes - use moduleSolver, only: doScatter, & - doEMField, & - doOutput, & - solver, & - doPushes, & - doReset, & - doAverage, & - doInjects - use moduleMesh, only: boundariesEM_update, & - boundariesEM_update, & - boundariesParticle_update, & - mesh, & - meshForMCC, & - doMCCollisions, & - doCollisions, & - doCoulombScattering - use moduleInject, only: updateInjects - implicit none + USE moduleCompTime + USE moduleCaseParam + USE moduleInput + USE moduleInject + USE moduleSolver + USE moduleMesh + USE moduleProbe + USE moduleErrors + USE OMP_LIB + IMPLICIT NONE ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index a1bc3fc..f4dcc46 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -67,20 +67,19 @@ MODULE moduleRandom end function randomMaxwellian !Returns a random number in a Maxwellian distribution of mean 0 and width 1 - function randomHalfMaxwellian() result(rnd) + FUNCTION randomHalfMaxwellian() RESULT(rnd) IMPLICIT NONE - real(8):: rnd - real(8):: v1 + REAL(8):: rnd + REAL(8):: x rnd = 0.D0 - v1 = 0.D0 - do while (v1 == 0.D0) - v1 = random() - - end do + x = 0.D0 + DO WHILE (x == 0.D0) + CALL RANDOM_NUMBER(x) + END DO - rnd = sqrt(-2.d0*log(v1)) + rnd = DSQRT(-DLOG(x)) END FUNCTION randomHalfMaxwellian diff --git a/src/modules/common/velocityDistribution.f90 b/src/modules/common/velocityDistribution.f90 index 359eebd..5b6c342 100644 --- a/src/modules/common/velocityDistribution.f90 +++ b/src/modules/common/velocityDistribution.f90 @@ -53,11 +53,12 @@ module velocityDistribution function randomVelMaxwellian(self) result (v) use moduleRandom implicit none + class(velDistMaxwellian), intent(in):: self real(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian() + v = self%vTh*randomMaxwellian()/sqrt(2.d0) end function randomVelMaxwellian diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index eb8bde5..bcf4dcc 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1431,6 +1431,7 @@ MODULE moduleInput CALL config%info('inject', found, n_children = nInject) ALLOCATE(injects(1:nInject)) + nPartInj = 0 DO i = 1, nInject WRITE(iString, '(i2)') i object = 'inject(' // trim(iString) // ')' @@ -1441,6 +1442,12 @@ MODULE moduleInput END DO + !Allocate array for injected particles + IF (nPartInj > 0) THEN + ALLOCATE(partInj(1:nPartInj)) + + END IF + END SUBROUTINE readInject SUBROUTINE readAverage(config) diff --git a/src/modules/makefile b/src/modules/makefile index 5789276..9f62c85 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -14,7 +14,7 @@ output.o: moduleSpecies.o common.o mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o $(MAKE) -C mesh all -solver.o: moduleInject.o moduleSpecies.o moduleProbe.o common.o output.o mesh.o +solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o $(MAKE) -C solver all init.o: common.o solver.o moduleInject.o diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 942189c..39de1ad 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -12,7 +12,7 @@ MODULE moduleInject REAL(8):: temperature(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n - INTEGER:: nParticles !Number of particles to inject each time step + INTEGER:: nParticles !Number of particles to introduce each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) @@ -72,6 +72,7 @@ MODULE moduleInject USE moduleRefParam USE moduleConstParam USE moduleSpecies + USE moduleSolver USE moduleErrors use json_module IMPLICIT NONE @@ -210,6 +211,8 @@ MODULE moduleInject END DO + self%nParticles = SUM(self%particlesPerEdge) + ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight @@ -217,14 +220,13 @@ MODULE moduleInject self%particlesPerEdge(e) = max(1,FLOOR(fluxPerStep*self%edges(e)%obj%surface / self%species%weight)) END DO + self%nParticles = SUM(self%particlesPerEdge) + !Rescale weight to match flux self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) END IF - ! Total number of particles to inject - self%nParticles = SUM(self%particlesPerEdge) - !Scale particles for different species steps IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') @@ -312,85 +314,126 @@ MODULE moduleInject end subroutine updateQuasiNeutral + !Injection of particles + SUBROUTINE doInjects() + USE moduleSpecies + USE moduleSolver + IMPLICIT NONE + + INTEGER:: i + + !$OMP SINGLE + nPartInj = 0 + DO i = 1, nInject + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + nPartInj = nPartInj + injects(i)%obj%nParticles + + END IF + + END DO + + IF (ALLOCATED(partInj)) DEALLOCATE(partInj) + ALLOCATE(partInj(1:nPartInj)) + !$OMP END SINGLE + + DO i=1, nInject + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + CALL injects(i)%obj%addParticles() + + END IF + END DO + + END SUBROUTINE doInjects + !Add particles for the injection SUBROUTINE addParticles(self) USE moduleSpecies + USE moduleSolver USE moduleMesh USE moduleRandom USE moduleErrors - use moduleList, only: partInj IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - type(particle), pointer:: part + INTEGER, SAVE:: nMin INTEGER:: i, e - INTEGER:: sp - CLASS(meshEdge), POINTER:: edge - + INTEGER:: n, sp + CLASS(meshEdge), POINTER:: randomEdge REAL(8):: direction(1:3) !Insert particles + !$OMP SINGLE + nMin = 0 + DO i = 1, self%id -1 + IF (solver%pusher(injects(i)%obj%species%n)%pushSpecies) THEN + nMin = nMin + injects(i)%obj%nParticles + + END IF + + END DO + nMin = nMin + 1 + !$OMP END SINGLE + !$OMP DO DO e = 1, self%nEdges ! Select edge for injection - edge => self%edges(e)%obj + randomEdge => self%edges(e)%obj ! Inject particles in edge DO i = 1, self%particlesPerEdge(e) - allocate(part) - ! Index in the partInj array - !Particle is considered to be inside the domain - part%n_in = .true. + ! Index in the global partInj array + n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i + !Particle is considered to be outside the domain + partInj(n)%n_in = .FALSE. !Random position in edge - part%r = edge%randPos() + partInj(n)%r = randomEdge%randPos() !Assign weight to particle. - part%weight = self%weightPerEdge(e) + partInj(n)%weight = self%weightPerEdge(e) !Volume associated to the edge: - IF (ASSOCIATED(edge%e1)) THEN - part%cell = edge%e1%n + IF (ASSOCIATED(randomEdge%e1)) THEN + partInj(n)%cell = randomEdge%e1%n - ELSEIF (ASSOCIATED(edge%e2)) THEN - part%cell = edge%e2%n + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN + partInj(n)%cell = randomEdge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - part%cellColl = edge%eColl%n + partInj(n)%cellColl = randomEdge%eColl%n sp = self%species%n !Assign particle type - part%species => self%species + partInj(n)%species => self%species if (all(self%n == 0.D0)) then - direction = edge%normal + direction = randomEdge%normal else direction = self%n end if - part%v = 0.D0 + partInj(n)%v = 0.D0 - ! do while(dot_product(partInj(n)%v, direction) <= 0.d0) - part%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + do while(dot_product(partInj(n)%v, direction) <= 0.d0) + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) !If injecting a no-drift distribution and velocity is negative, reflect if ((self%vMod == 0.D0) .and. & - (dot_product(part%v, direction) <= 0.D0)) then - part%v = - part%v + (dot_product(partInj(n)%v, direction) <= 0.D0)) then + partInj(n)%v = - partInj(n)%v end if - ! end do + end do !Obtain natural coordinates of particle in cell - part%Xi = mesh%cells(part%cell)%obj%phy2log(part%r) - - ! Add particle to global list - call partInj%setLock() - call partInj%add(part) - call partInj%unsetLock() + 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 + CALL solver%updateParticleCell(partInj(n)) END DO @@ -457,7 +500,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -467,7 +510,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index 44524c6..b1dafdc 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -23,7 +23,6 @@ MODULE moduleList END TYPE listNode - TYPE(listNode):: partInj !Particles comming from injections TYPE(listNode):: partWScheme !Particles comming from the nonAnalogue scheme TYPE(listNode):: partCollisions !Particles created in collisional process TYPE(listNode):: partSurfaces !Particles created in surface interactions diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 index 8c575da..ab08f08 100644 --- a/src/modules/moduleSpecies.f90 +++ b/src/modules/moduleSpecies.f90 @@ -48,8 +48,11 @@ MODULE moduleSpecies !Number of old particles INTEGER:: nPartOld + !Number of injected particles + INTEGER:: nPartInj !Arrays that contain the particles TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration + TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles CONTAINS FUNCTION speciesName2Index(speciesName) RESULT(sp) diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index e0ca06c..f951839 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -235,21 +235,23 @@ MODULE moduleSolver INTEGER:: nn, n, e INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nOldIn, nInj, nWScheme, nCollisions, nSurfaces + INTEGER, SAVE:: nInjIn, nOldIn, nWScheme, nCollisions, nSurfaces TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) INTEGER:: s !$OMP SECTIONS !$OMP SECTION + nInjIn = 0 + IF (ALLOCATED(partInj)) THEN + nInjIn = COUNT(partInj%n_in) + + END IF + !$OMP SECTION nOldIn = 0 IF (ALLOCATED(partOld)) THEN nOldIn = COUNT(partOld%n_in) END IF - - !$OMP SECTION - nInj = partInj%amount - !$OMP SECTION nWScheme = partWScheme%amount @@ -265,14 +267,30 @@ MODULE moduleSolver !$OMP SINGLE CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInj + nOldIn + nWScheme + nCollisions + nSurfaces + nPartNew = nInjIn + nOldIn + nWScheme + nCollisions + nSurfaces ALLOCATE(partOld(1:nPartNew)) !$OMP END SINGLE !$OMP SECTIONS !$OMP SECTION - !Reset particles from previous iteration + !Reset particles from injection nn = 0 + DO n = 1, nPartInj + IF (partInj(n)%n_in) THEN + nn = nn + 1 + partOld(nn) = partInj(n) + IF (nProbes > 0) THEN + CALL doProbes(partOld(nn)) + + END IF + + END IF + + END DO + + !$OMP SECTION + !Reset particles from previous iteration + nn = nInjIn DO n = 1, nPartOld IF (partTemp(n)%n_in) THEN nn = nn + 1 @@ -286,24 +304,19 @@ MODULE moduleSolver END DO - !$OMP SECTION - !Reset particles from injection - nn = nOldIn - call resetList(partInj, partOld, nn) - !$OMP SECTION !Reset particles from weighting scheme - nn = nOldIn + nInj + nn = nInjIn + nOldIn CALL resetList(partWScheme, partOld, nn) !$OMP SECTION !Reset particles from collisional process - nn = nOldIn + nInj + nWScheme + nn = nInjIn + nOldIn + nWScheme CALL resetList(partCollisions, partOld, nn) !$OMP SECTION !Reset particles from surface process - nn = nOldIn + nInj + nWScheme + nCollisions + nn = nInjIn + nOldIn + nWScheme + nCollisions CALL resetList(partSurfaces, partOld, nn) !$OMP SECTION @@ -460,28 +473,6 @@ MODULE moduleSolver END SUBROUTINE splitParticle - !Injection of particles - SUBROUTINE doInjects() - USE moduleSpecies - use moduleInject - IMPLICIT NONE - - INTEGER:: i, n, sp - - DO i=1, nInject - associate(inject => injects(i)%obj) - sp = inject%species%n - IF (solver%pusher(sp)%pushSpecies) THEN - - CALL inject%addParticles() - - END IF - end associate - - END DO - - END SUBROUTINE doInjects - SUBROUTINE updateParticleCell(self, part) USE moduleSpecies USE moduleMesh