diff --git a/runs/Argon_Expansion/base_case.json b/runs/Argon_Expansion/base_case.json index 64b0e21..78e16ed 100644 --- a/runs/Argon_Expansion/base_case.json +++ b/runs/Argon_Expansion/base_case.json @@ -41,7 +41,7 @@ "radius": 1.88e-10 }, "case": { - "tau": [1.0e-4, 1.0e-6], + "tau": [1.0e-6, 1.0e-6], "time": 4.0e-3, "pusher": ["2DCylNeutral", "2DCylNeutral"], "WeightingScheme": "Volume" diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index f239fbc..bde3aaf 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -92,11 +92,11 @@ MODULE moduleInject SELECT CASE(units) CASE ("sccm") !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tau(self%sp)*ti_ref/species(sp)%obj%weight) + self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight) CASE ("A") !Input current in Ampers - self%nParticles = INT(flow*tau(self%sp)*ti_ref/(qe*species(sp)%obj%weight)) + self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight)) CASE DEFAULT CALL criticalError("No support for units: " // units, 'initInject') @@ -106,7 +106,6 @@ MODULE moduleInject IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') !Gets the edge elements from which particles are injected - !TODO: Improve this A LOT DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface @@ -114,7 +113,6 @@ MODULE moduleInject self%nEdges = COUNT(phSurface == physicalSurface) ALLOCATE(inject(i)%edges(1:self%nEdges)) - ! ALLOCATE(inject(i)%weight(1:self%nEdges)) et = 0 DO e=1, mesh%numEdges IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN @@ -124,7 +122,8 @@ MODULE moduleInject END IF END DO - ! self%sumWeight = SUM(self%weight) + + nPartInj = nPartInj + self%nParticles END SUBROUTINE initInject @@ -137,17 +136,12 @@ MODULE moduleInject INTEGER:: i !$OMP SINGLE - nPartInj = 0 - DO i = 1, nInject - IF (solver%pusher(inject(i)%sp)%pushSpecies) nPartInj = nPartInj + inject(i)%nParticles - - END DO IF (ALLOCATED(partInj)) DEALLOCATE(partInj) ALLOCATE(partInj(1:nPartInj)) !$OMP END SINGLE DO i=1, nInject - IF (solver%pusher(inject(i)%sp)%pushSpecies) CALL inject(i)%addParticles() + CALL inject(i)%addParticles() END DO END SUBROUTINE doInjects @@ -206,19 +200,13 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(in):: self INTEGER:: randomX - INTEGER:: i!, j INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array INTEGER:: n CLASS(meshEdge), POINTER:: randomEdge !Insert particles !$OMP SINGLE - nMin = 0 - DO i = 1, self%id - 1 - IF (solver%pusher(inject(i)%sp)%pushSpecies) nMin = nMin + inject(i)%nParticles - - END DO - nMin = nMin + 1 + nMin = SUM(inject(1:(self%id-1))%nParticles) + 1 nMax = nMin + self%nParticles - 1 !Assign particle type partInj(nMin:nMax)%sp = self%sp @@ -254,8 +242,8 @@ MODULE moduleInject self%v(2)%obj%randomVel(), & self%v(3)%obj%randomVel() /) - !Push new particle - CALL solver%pusher(self%sp)%pushParticle(partInj(n)) + !Push new particle with the minimum time step + CALL solver%pusher(self%sp)%pushParticle(partInj(n), tauMin) !Assign cell to new particle CALL solver%updateParticleCell(partInj(n)) diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index ba146d4..dfd5d4a 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -27,10 +27,11 @@ MODULE moduleSolver INTERFACE !Push a particle - PURE SUBROUTINE push_interafece(part) + PURE SUBROUTINE push_interafece(part, tauIn) USE moduleSpecies TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn END SUBROUTINE push_interafece @@ -134,7 +135,7 @@ MODULE moduleSolver !Checks if the species sp is update this iteration IF (solver%pusher(sp)%pushSpecies) THEN !Push particle - CALL solver%pusher(sp)%pushParticle(partOld(n)) + CALL solver%pusher(sp)%pushParticle(partOld(n), tau(sp)) !Find cell in wich particle reside CALL solver%updateParticleCell(partOld(n)) @@ -146,27 +147,25 @@ MODULE moduleSolver END SUBROUTINE doPushes !Push one particle. Boris pusher for 2D Cyl Neutral particle - PURE SUBROUTINE pushCylNeutral(part) + PURE SUBROUTINE pushCylNeutral(part, tauIn) USE moduleSpecies IMPLICIT NONE TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn TYPE(particle):: part_temp - REAL(8):: tauSp REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha REAL(8):: v_p_oh_star(2:3) part_temp = part - !Time step for the species - tauSp = tau(part_temp%sp) !z part_temp%v(1) = part%v(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn !r,theta v_p_oh_star(2) = part%v(2) - x_new = part%r(2) + v_p_oh_star(2)*tauSp + x_new = part%r(2) + v_p_oh_star(2)*tauIn v_p_oh_star(3) = part%v(3) - y_new = v_p_oh_star(3)*tauSp + y_new = v_p_oh_star(3)*tauIn r = DSQRT(x_new**2+y_new**2) part_temp%r(2) = r IF (r > 0.D0) THEN @@ -185,31 +184,29 @@ MODULE moduleSolver END SUBROUTINE pushCylNeutral !Push one particle. Boris pusher for 2D Cyl Charged particle - PURE SUBROUTINE pushCylCharged(part) + PURE SUBROUTINE pushCylCharged(part, tauIn) USE moduleSpecies USE moduleEM IMPLICIT NONE TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn REAL(8):: v_p_oh_star(2:3) TYPE(particle):: part_temp REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha - REAL(8):: tauSp - REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass + REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass part_temp = part - !Time step for the species - tauSp = tau(part_temp%sp) !Get electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp + qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn !z part_temp%v(1) = part%v(1) + qmEFt(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn !r,theta v_p_oh_star(2) = part%v(2) + qmEFt(2) - x_new = part%r(2) + v_p_oh_star(2)*tauSp + x_new = part%r(2) + v_p_oh_star(2)*tauIn v_p_oh_star(3) = part%v(3) + qmEFt(3) - y_new = v_p_oh_star(3)*tauSp + y_new = v_p_oh_star(3)*tauIn r = DSQRT(x_new**2+y_new**2) part_temp%r(2) = r IF (r > 0.D0) THEN @@ -228,25 +225,23 @@ MODULE moduleSolver END SUBROUTINE pushCylCharged !Push charged particles in 1D cartesian coordinates - PURE SUBROUTINE push1DCartCharged(part) + PURE SUBROUTINE push1DCartCharged(part, tauIn) USE moduleSPecies USE moduleEM IMPLICIT NONE TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn TYPE(particle):: part_temp - REAL(8):: tauSp REAL(8):: qmEFt(1:3) part_temp = part - !Time step for particle species - tauSp = tau(part_temp%sp) !Get the electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp + qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn !x part_temp%v(1) = part%v(1) + qmEFt(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn part_temp%n_in = .FALSE. @@ -255,28 +250,27 @@ MODULE moduleSolver END SUBROUTINE push1DCartCharged !Push one particle. Boris pusher for 1D Radial Charged particle - PURE SUBROUTINE push1DRadCharged(part) + PURE SUBROUTINE push1DRadCharged(part, tauIn) USE moduleSpecies USE moduleEM IMPLICIT NONE TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn REAL(8):: v_p_oh_star(1:2) TYPE(particle):: part_temp REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha - REAL(8):: tauSp - REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass + REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass part_temp = part !Time step for the species - tauSp = tau(part_temp%sp) !Get electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp + qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin !r,theta v_p_oh_star(1) = part%v(1) + qmEFt(1) - x_new = part%r(1) + v_p_oh_star(1)*tauSp + x_new = part%r(1) + v_p_oh_star(1)*tauIn v_p_oh_star(2) = part%v(2) + qmEFt(2) - y_new = v_p_oh_star(2)*tauSp + y_new = v_p_oh_star(2)*tauIn r = DSQRT(x_new**2+y_new**2) part_temp%r(1) = r IF (r > 0.D0) THEN