JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

First attempt at Coulomb collisions #46

Merged
JorgeGonz merged 21 commits from feature/CoulombLinear into development 2023-07-16 14:47:59 +02:00
Conversation 2 Commits 21 Files changed 12 +59 -47
2 changed files with 59 additions and 47 deletions
Download patch file Download diff file Expand all files Collapse all files
Showing only changes of commit 6113ac3305 - Show all commits

Correction with conservation

Prev Next 
Now the method is much better in conserving total energy.
However, still there is an issue with collisions between species of
dispaprate mass.
Jorge Gonzalez 2023-03-06 16:16:17 +01:00
commit 6113ac3305

91
src/modules/mesh/moduleMesh.f90
View file

@ -966,7 +966,6 @@ MODULE moduleMesh
INTEGER:: i, j
INTEGER:: n
TYPE(lNode), POINTER:: partTemp
REAL(8):: W(3), dW(2) !Relative velocity between particle and species and its increment
INTEGER(8), ALLOCATABLE:: cellNodes(:)
CLASS(meshNode), POINTER:: node
TYPE(outputFormat):: output
@ -974,7 +973,8 @@ MODULE moduleMesh
REAL(8):: density, velocity(1:3), temperature!values at particle position
REAL(8), DIMENSION(1:3):: e1, e2, e3
REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square
REAL(8):: l, lW, AW
REAL(8):: W(3), dW(2), normW !Relative velocity between particle and species and its increment
REAL(8):: l2, l, lW, AW
REAL(8):: rnd
@ -1007,12 +1007,25 @@ MODULE moduleMesh
density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
l = coulombMatrix(k)%l_j/SQRT(temperature)
l2 = coulombMatrix(k)%l2_j/temperature
l = SQRT(l2)
W = partTemp%part%v - velocity
lW = l * NORM2(W)
AW = coulombMatrix(k)%A_ij/NORM2(W)
!Axis of the relative velocity
normW = NORM2(W)
lW = l * normW
AW = coulombMatrix(k)%A_i/normW
delta_par = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW)
delta_par_square = AW*density*G(lW)
delta_per_square = AW*density*H(lW)
dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
!System of reference for the velocity change
!First one is parallel to the relative velocity
e1 = normalize(W)
!Second one is perpendicular to it
@ -1024,18 +1037,12 @@ MODULE moduleMesh
e3 = crossProduct(e2, e1)
e3 = normalize(e3)
delta_par = -coulombMatrix(k)%A_ij*coulombMatrix(k)%one_plus_massRatio_ij*density*l**2*G(lW)
!Random number for direction
rnd = PI2*random()
delta_par_square = AW*density*G(lW)
delta_per_square = AW*density*H(lW)
dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
rnd = random()
partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + &
SIN(PI2*rnd)*e3)
!Change particle velocity
partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
SIN(rnd)*e3)
partTemp => partTemp%next
@ -1049,20 +1056,33 @@ MODULE moduleMesh
densityNodes(n) = output%density/n_ref
velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
temperatureNodes(n) = output%temperature/T_ref
END DO
partTemp => cell%listPart_in(j)%head
DO WHILE(ASSOCIATED(partTemp))
density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
l = coulombMatrix(k)%l_i/SQRT(temperature)
l2 = coulombMatrix(k)%l2_i/temperature
l = SQRT(l2)
W = partTemp%part%v - velocity
lW = l * NORM2(W)
AW = coulombMatrix(k)%A_ji/NORM2(W)
!Axis of the relative velocity
normW = NORM2(W)
lW = l * normW
AW = coulombMatrix(k)%A_j/normW
delta_par = -coulombMatrix(k)%A_j*coulombMatrix(k)%one_plus_massRatio_ji*density*l2*G(lW)
delta_par_square = AW*density*G(lW)
delta_per_square = AW*density*H(lW)
dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
!System of reference for the velocity change
!First one is parallel to the relative velocity
e1 = normalize(W)
!Second one is perpendicular to it
@ -1073,24 +1093,15 @@ MODULE moduleMesh
!Third one is perpendicular to the other two
e3 = crossProduct(e2, e1)
e3 = normalize(e3)
delta_par = -coulombMatrix(k)%A_ji*coulombMatrix(k)%one_plus_massRatio_ji*density*l**2*G(lW)
delta_par_square = AW*density*G(lW)
delta_per_square = AW*density*H(lW)
dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
rnd = random()
partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + &
SIN(PI2*rnd)*e3)
rnd = PI2*random()
partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + &
SIN(rnd)*e3)
partTemp => partTemp%next
END DO
END IF
END DO

15
src/modules/moduleCoulomb.f90
View file

@ -10,8 +10,8 @@ MODULE moduleCoulomb
CLASS(speciesGeneric), POINTER:: sp_j
REAL(8):: one_plus_massRatio_ij, one_plus_massRatio_ji
REAL(8):: lnCoulomb !This can be done a function in the future
REAL(8):: A_ij, A_ji
REAL(8):: l_j, l_i
REAL(8):: A_i, A_j
REAL(8):: l2_j, l2_i
CONTAINS
PROCEDURE, PASS:: init => initInteractionCoulomb
@ -79,12 +79,13 @@ MODULE moduleCoulomb
END SELECT
self%lnCoulomb = 12.0
self%A_ij = 8.D0*PI*Z_i**2*Z_j**2*e_ref**4*self%lnCoulomb / self%sp_i%m
self%A_ji = 8.D0*PI*Z_j**2*Z_i**2*e_ref**4*self%lnCoulomb / self%sp_j%m
self%lnCoulomb = 12.0 !Make this function dependent
self%l_j = SQRT(self%sp_j%m / 2.D0) !Missing temperature because it's cell dependent
self%l_i = SQRT(self%sp_i%m / 2.D0) !Missing temperature because it's cell dependent
self%A_i = 8.D0*PI*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2
self%A_j = 8.D0*PI*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2
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