Jorge Gonzalez
e05c0d4635
Coulomb scattering is now fully conservative thanks to the method in lemos2009small. The trick was to conserve the momentum and energy of ALL particles involved in the scattering in each cell. The substeps in Coulomb collisions have been removed as they are no longer necessary. Still some issues with e-i, but I don't know right now.
98 lines
2.6 KiB
Fortran
98 lines
2.6 KiB
Fortran
MODULE moduleCoulomb
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
INTEGER:: nCoulombPairs = 0
|
|
|
|
!Type for Coulomb iteraction matrix
|
|
TYPE:: interactionsCoulomb
|
|
CLASS(speciesGeneric), POINTER:: sp_i
|
|
CLASS(speciesGeneric), POINTER:: sp_j
|
|
REAL(8):: one_plus_massRatio_ij
|
|
REAL(8):: lnCoulomb !This can be done a function in the future
|
|
REAL(8):: A_i
|
|
REAL(8):: A_j
|
|
REAL(8):: l2_j
|
|
REAL(8):: l2_i
|
|
CONTAINS
|
|
PROCEDURE, PASS:: init => initInteractionCoulomb
|
|
|
|
END TYPE interactionsCoulomb
|
|
|
|
!Coulomb collision 'matrix'
|
|
TYPE(interactionsCoulomb), ALLOCATABLE:: coulombMatrix(:)
|
|
|
|
CONTAINS
|
|
PURE REAL(8) FUNCTION G(x)
|
|
USE moduleConstParam
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: x
|
|
|
|
G = 0.D0
|
|
IF (x /= 0.D0) THEN
|
|
G = (ERF(x) - x*2.D0/SQRT(PI)*EXP(-x**2))/(2.D0*x**2)
|
|
|
|
END IF
|
|
|
|
END FUNCTION G
|
|
|
|
PURE REAL(8) FUNCTION H(x)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: x
|
|
|
|
H = ERF(x) - G(x)
|
|
|
|
END FUNCTION H
|
|
|
|
SUBROUTINE initInteractionCoulomb(self, i, j)
|
|
USE moduleSpecies
|
|
USE moduleErrors
|
|
USE moduleConstParam
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(interactionsCoulomb), INTENT(out):: self
|
|
INTEGER, INTENT(in):: i, j
|
|
REAL(8):: Z_i, Z_j
|
|
REAL(8):: scaleFactor
|
|
|
|
self%sp_i => species(i)%obj
|
|
self%sp_j => species(j)%obj
|
|
|
|
self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m
|
|
|
|
Z_i = 0.D0
|
|
Z_j = 0.D0
|
|
SELECT TYPE(sp => self%sp_i)
|
|
TYPE IS (speciesCharged)
|
|
Z_i = sp%q
|
|
|
|
CLASS DEFAULT
|
|
CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
|
|
|
|
END SELECT
|
|
|
|
SELECT TYPE(sp => self%sp_j)
|
|
TYPE IS (speciesCharged)
|
|
Z_j = sp%q
|
|
|
|
CLASS DEFAULT
|
|
CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
|
|
|
|
END SELECT
|
|
|
|
self%lnCoulomb = 10.0 !Make this function dependent
|
|
|
|
scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
|
|
|
|
self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
|
|
self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent
|
|
|
|
self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent
|
|
self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent
|
|
|
|
END SUBROUTINE initInteractionCoulomb
|
|
|
|
END MODULE moduleCoulomb
|