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 +6 -4
2 changed files with 6 additions and 4 deletions
Download patch file Download diff file Expand all files Collapse all files
Showing only changes of commit ed6c2c46e4 - Show all commits

I'm a fucking idiot

Prev Next 
The limit I set to avoid divisions by zero was wront and collisions were
being skipped. It is corrected now.
Jorge Gonzalez 2023-03-29 10:08:11 +02:00
commit ed6c2c46e4

2
src/modules/init/moduleInput.f90
View file

@ -155,7 +155,7 @@ MODULE moduleInput
sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
ELSE ELSE
sigmaVrel_ref = 0.D0 !Assume no collisions sigmaVrel_ref = L_ref**2 * v_ref
END IF END IF

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

@ -1016,7 +1016,7 @@ MODULE moduleMesh
W = partTemp%part%v - velocity W = partTemp%part%v - velocity
normW = NORM2(W) normW = NORM2(W)
IF (normW > 1.D-10) THEN IF (normW < 1.D-12) THEN
!If relative velocity is too low, skip collision and move to next particle !If relative velocity is too low, skip collision and move to next particle
partTemp => partTemp%next partTemp => partTemp%next
@ -1055,6 +1055,7 @@ MODULE moduleMesh
deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3) deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV) totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV)
!Change particle velocity !Change particle velocity
partTemp%part%v = partTemp%part%v + deltaV partTemp%part%v = partTemp%part%v + deltaV
@ -1078,6 +1079,7 @@ MODULE moduleMesh
normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) * & normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) * &
(coulombMatrix(k)%sp_i%weight*coulombMatrix(k)%sp_i%m) / & (coulombMatrix(k)%sp_i%weight*coulombMatrix(k)%sp_i%m) / &
(coulombMatrix(k)%sp_j%weight*coulombMatrix(k)%sp_j%m) (coulombMatrix(k)%sp_j%weight*coulombMatrix(k)%sp_j%m)
!Loop over particles of species_j !Loop over particles of species_j
partTemp => cell%listPart_in(j)%head partTemp => cell%listPart_in(j)%head
DO WHILE(ASSOCIATED(partTemp)) DO WHILE(ASSOCIATED(partTemp))
@ -1090,7 +1092,7 @@ MODULE moduleMesh
W = partTemp%part%v - velocity W = partTemp%part%v - velocity
normW = NORM2(W) normW = NORM2(W)
IF (normW > 1.D-10) THEN IF (normW < 1.D-12) THEN
!If relative velocity is too low, skip collision and move to next particle !If relative velocity is too low, skip collision and move to next particle
partTemp => partTemp%next partTemp => partTemp%next
@ -1111,7 +1113,7 @@ MODULE moduleMesh
dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin))
!Normalize with average exchange per particle !Normalize with average exchange per particle
!TODO: This is a dirty trick, it should not be neccessary but I don't lnow why it is. !TODO: This is a dirty trick to ensure conservation between species
dW = normDeltaV*dW/NORM2(dW) dW = normDeltaV*dW/NORM2(dW)
!System of reference for the velocity change !System of reference for the velocity change