diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index f4dcc46..a1bc3fc 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -67,19 +67,20 @@ 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):: x + real(8):: rnd + real(8):: v1 rnd = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO + v1 = 0.D0 + do while (v1 == 0.D0) + v1 = random() + + end do - rnd = DSQRT(-DLOG(x)) + rnd = sqrt(-2.d0*log(v1)) END FUNCTION randomHalfMaxwellian diff --git a/src/modules/common/velocityDistribution.f90 b/src/modules/common/velocityDistribution.f90 index 5b6c342..359eebd 100644 --- a/src/modules/common/velocityDistribution.f90 +++ b/src/modules/common/velocityDistribution.f90 @@ -53,12 +53,11 @@ 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()/sqrt(2.d0) + v = self%vTh*randomMaxwellian() end function randomVelMaxwellian diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 39de1ad..7a1c5d6 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 introduce each time step + INTEGER:: nParticles !Number of particles to inject each time step CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges type(meshEdgePointer), allocatable:: edges(:) @@ -72,7 +72,6 @@ MODULE moduleInject USE moduleRefParam USE moduleConstParam USE moduleSpecies - USE moduleSolver USE moduleErrors use json_module IMPLICIT NONE @@ -211,8 +210,6 @@ MODULE moduleInject END DO - self%nParticles = SUM(self%particlesPerEdge) - ELSE ! No particles assigned per edge, use the species weight self%weightPerEdge = self%species%weight @@ -220,13 +217,14 @@ 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') @@ -314,126 +312,82 @@ 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) + SUBROUTINE addParticles(self, partArray) USE moduleSpecies - USE moduleSolver USE moduleMesh USE moduleRandom USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - INTEGER, SAVE:: nMin + type(particle), allocatable, intent(inout):: partArray(:) INTEGER:: i, e INTEGER:: n, sp - CLASS(meshEdge), POINTER:: randomEdge + CLASS(meshEdge), POINTER:: edge REAL(8):: direction(1:3) + !$omp single + n = 0 + !$omp end single + !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 - randomEdge => self%edges(e)%obj + edge => self%edges(e)%obj ! Inject particles in edge DO i = 1, self%particlesPerEdge(e) - ! 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. + ! Index in the partInj array + n = n + 1 + !Particle is considered to be outside the domain for now + partArray(n)%n_in = .FALSE. !Random position in edge - partInj(n)%r = randomEdge%randPos() + partArray(n)%r = edge%randPos() !Assign weight to particle. - partInj(n)%weight = self%weightPerEdge(e) + partArray(n)%weight = self%weightPerEdge(e) !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + IF (ASSOCIATED(edge%e1)) THEN + partArray(n)%cell = edge%e1%n - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + ELSEIF (ASSOCIATED(edge%e2)) THEN + partArray(n)%cell = edge%e2%n ELSE CALL criticalError("No Volume associated to edge", 'addParticles') END IF - partInj(n)%cellColl = randomEdge%eColl%n + partArray(n)%cellColl = edge%eColl%n sp = self%species%n !Assign particle type - partInj(n)%species => self%species + partArray(n)%species => self%species if (all(self%n == 0.D0)) then - direction = randomEdge%normal + direction = edge%normal else direction = self%n end if - partInj(n)%v = 0.D0 + partArray(n)%v = 0.D0 - do while(dot_product(partInj(n)%v, direction) <= 0.d0) - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + ! do while(dot_product(partInj(n)%v, direction) <= 0.d0) + partArray(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(partInj(n)%v, direction) <= 0.D0)) then - partInj(n)%v = - partInj(n)%v + (dot_product(partArray(n)%v, direction) <= 0.D0)) then + partArray(n)%v = - partArray(n)%v end if - end do + ! end do !Obtain natural coordinates of particle in cell - 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)) + partArray(n)%Xi = mesh%cells(partArray(n)%cell)%obj%phy2log(partArray(n)%r) END DO @@ -500,7 +454,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistMaxwellian @@ -510,7 +464,7 @@ MODULE moduleInject CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(2.d0*temperature/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian