From 96c563c1462c048e4ddf96ba73b142ab00dd70ff Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Thu, 11 Jul 2024 14:39:56 +0200 Subject: [PATCH] Finally, some progress I rewrote how particles are injected. Now the particles per edge and its weight are calculated in the initialization. There is the possibility for the user to select the particles per edge. TODO: Write documentation for new feature. TODO: Test in 2DCyl --- src/modules/init/moduleInput.f90 | 5 +- src/modules/moduleInject.f90 | 183 +++++++++++++++++++------------ 2 files changed, 114 insertions(+), 74 deletions(-) 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