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

Trying to implement Lemos Coulomb Scatering

Browse source 
I was having tones of issues with the previous implementation. I think
the problem was the velocity vector and how it was returning to the
normal reference frame.

I hope this new implementation works better.
This commit is contained in:
Jorge Gonzalez 2023-07-11 07:51:49 +02:00
commit d75af4bda7
1 changed files with 20 additions and 29 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -973,9 +973,12 @@ 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):: W(3), dW(2), normW !Relative velocity between particle and species and its increment
REAL(8):: l2, l, lW, AW
REAL(8):: deltaV(1:3)
REAL(8):: normW
REAL(8):: l2, l, lW, AW, AW2, AW3
REAL(8):: deltaW, deltaThe, deltaPhi
REAL(8):: cosDeltaPhi, sinDeltaPhi
REAL(8):: cosDeltaThe, sinDeltaThe
REAL(8):: deltaV(1:3) !Relative particle velocity increment
REAL(8):: rnd
REAL(8):: eps = 1.D-10
@ -1023,8 +1026,7 @@ MODULE moduleMesh
END IF
W = partTemp%part%v - velocity
normW = NORM2(W)
normW = NORM2(partTemp%part%v - velocity)
!If relative velocity is too low, skip collision to avoid division by zero and move to next particle
IF (normW < eps) THEN
@ -1035,33 +1037,22 @@ MODULE moduleMesh
END IF
lW = l * normW
AW = coulombMatrix(k)%A_i/normW
AW = coulombMatrix(k)%A_i/normW
AW2 = coulombMatrix(k)%A_i/normW**2 / 2.D0
AW3 = coulombMatrix(k)%A_i/normW**3 / 2.D0
delta_par = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW)
deltaW = (-coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*l2*G(lW) - AW2*G(lW) + AW2*ERF(lw)) * tauMin + &
SQRT(AW * G(lW) * tauMin) * ABS(randomMaxwellian())
deltaThe = SQRT(2.D0 * AW3 * H(lW) * tauMin)*ABS(randomMaxwellian())
deltaPhi = PI2 * random()
delta_par_square = AW*density*G(lW)
cosDeltaThe = COS(deltaThe)
sinDeltaThe = SIN(deltaThe)
cosDeltaPhi = COS(deltaPhi)
sinDeltaPhi = SIN(deltaPhi)
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
e2(1) = -e1(2)
e2(2) = e1(1)
e2(3) = 0.D0
e2 = normalize(e2)
!Third one is perpendicular to the other two
e3 = crossProduct(e1, e2)
e3 = normalize(e3)
!Random number for direction
rnd = PI2*random()
deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
!Rotate velocity frame assuming theta = phi = 0
deltaV = (normW + deltaW) * (/ sinDeltaThe * cosDeltaPhi, sinDeltaThe * sinDeltaPhi, cosDeltaPhi /)
!Change particle velocity
partTemp%part%v = partTemp%part%v + deltaV