diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index e30bf97..bd59e32 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -1236,6 +1236,7 @@ MODULE moduleInput REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface + INTEGER:: particlesPerEdge INTEGER:: sp CALL config%info('inject', found, n_children = nInject) @@ -1260,8 +1261,10 @@ MODULE moduleInput CALL config%get(object // '.flow', flow, found) CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) + particlesPerEdge = 0 + CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) + CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge) CALL readVelDistr(config, inject(i), object) diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 529b687..59748eb 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -61,6 +61,8 @@ MODULE moduleInject CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges + INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge + REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge REAL(8):: surface ! Total surface of injection TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS @@ -74,7 +76,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) + SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -86,52 +88,28 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i REAL(8), INTENT(in):: v, n(1:3), T(1:3) - INTEGER, INTENT(in):: sp, physicalSurface + INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) INTEGER:: nVolColl + REAL(8):: fluxPerStep = 0.D0 - self%id = i - self%vMod = v / v_ref - self%n = n / NORM2(n) - self%T = T / T_ref - self%species => species(sp)%obj - tauInject = tau(self%species%n) - SELECT CASE(units) - CASE ("sccm") - !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight) - - CASE ("A") - !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight)) - - CASE ("Am2") - !Input current in Ampers per square meter - self%nParticles = INT(flow*tauInject*ti_ref*L_ref**2/(qe*species(sp)%obj%weight)) - - CASE ("part/s") - !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight) - - CASE DEFAULT - CALL criticalError("No support for units: " // units, 'initInject') - - END SELECT - !Scale particles for different species steps - IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') - + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%T = T / T_ref !Gets the edge elements from which particles are injected DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface END DO - self%nEdges = COUNT(phSurface == physicalSurface) - ALLOCATE(inject(i)%edges(1:self%nEdges)) + ALLOCATE(self%edges(1:self%nEdges)) + ALLOCATE(self%particlesPerEdge(1:self%nEdges)) + ALLOCATE(self%weightPerEdge(1:self%nEdges)) et = 0 DO e=1, mesh%numEdges IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN @@ -170,6 +148,60 @@ MODULE moduleInject END DO + ! Information about species and flux + self%species => species(sp)%obj + tauInject = tau(self%species%n) + SELECT CASE(units) + CASE ("sccm") + !Standard cubic centimeter per minute + fluxPerStep = flow*sccm2atomPerS*tauInject*ti_ref + + CASE ("A") + !Current in Ampers + fluxPerStep = flow*tauInject*ti_ref/qe + + CASE ("Am2") + !Input current in Ampers per square meter + fluxPerStep = flow*tauInject*ti_ref*self%surface*L_ref**2/qe + + CASE ("part/s") + !Input current in Ampers + fluxPerStep = flow*tauInject*ti_ref + + CASE DEFAULT + CALL criticalError("No support for units: " // units, 'initInject') + + END SELECT + fluxPerStep = fluxPerStep / self%surface + + !Assign particles per edge + IF (particlesPerEdge > 0) THEN + ! Particles per edge defined by the user + self%particlesPerEdge = particlesPerEdge + DO et = 1, self%nEdges + self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge) + + END DO + + ELSE + ! No particles assigned per edge, use the species weight + self%weightPerEdge = self%species%weight + DO et = 1, self%nEdges + self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight) + + END DO + + END IF + + print *, self%particlesPerEdge + print *, self%weightPerEdge + + self%nParticles = SUM(self%particlesPerEdge) + print *, self%nParticles + + !Scale particles for different species steps + IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') + END SUBROUTINE initInject !Injection of particles @@ -285,9 +317,10 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(in):: self INTEGER:: randomX INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array - INTEGER:: i + INTEGER:: i, e INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge + INTEGER:: particlesPerEdge REAL(8):: direction(1:3) !Insert particles @@ -300,58 +333,62 @@ MODULE moduleInject END IF END DO - nMin = nMin + 1 - nMax = nMin + self%nParticles - 1 - !Particle is considered to be outside the domain - partInj(nMin:nMax)%n_in = .FALSE. !$OMP END SINGLE !$OMP DO - DO n = nMin, nMax - randomX = random(1, self%nEdges) - randomEdge => mesh%edges(self%edges(randomX))%obj - !Random position in edge - partInj(n)%r = randomEdge%randPos() - !Assign weight to particle. - partInj(n)%weight = self%species%weight * randomEdge%surface / self%surface - !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + DO e = 1, self%nEdges + ! Select edge for injection + randomEdge => mesh%edges(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. + !Random position in edge + partInj(n)%r = randomEdge%randPos() + !Assign weight to particle. + partInj(n)%weight = self%weightPerEdge(e) + !Volume associated to the edge: + IF (ASSOCIATED(randomEdge%e1)) THEN + partInj(n)%cell = randomEdge%e1%n - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN + partInj(n)%cell = randomEdge%e2%n - ELSE - CALL criticalError("No Volume associated to edge", 'addParticles') + ELSE + CALL criticalError("No Volume associated to edge", 'addParticles') - END IF - partInj(n)%cellColl = randomEdge%eColl%n - sp = self%species%n + END IF + partInj(n)%cellColl = randomEdge%eColl%n + sp = self%species%n - !Assign particle type - partInj(n)%species => self%species + !Assign particle type + partInj(n)%species => self%species - direction = self%n + direction = self%n - partInj(n)%v = 0.D0 + partInj(n)%v = 0.D0 - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) - !If velocity is not in the right direction, invert it - IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN - partInj(n)%v = - partInj(n)%v + !If velocity is not in the right direction, invert it + IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN + partInj(n)%v = - partInj(n)%v - END IF + END IF - !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)) + !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)) + + END DO END DO !$OMP END DO