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

First Coulomb implementation that works

Browse source 
After fixing all possible divisions by zero I was able to find in the
Coulomb collision I think that this is a first working implementation of
a Coulomb operator based on moments.

Still to test a few things, modify the manual but I would say that I'm
satisfiyed right now. This operator won't be used that often but maybe
improving efficiency is still needed.

In the future a binary operator is required to be able to study cases
out of Maxwellian equilibrium.
This commit is contained in:
Jorge Gonzalez 2023-07-04 17:01:02 +02:00
commit 8d35123508
3 changed files with 51 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl
r2 = self%n2%getCoordinates() r2 = self%n2%getCoordinates()
self%z = (/r1(1), r2(1)/) self%z = (/r1(1), r2(1)/)
self%r = (/r1(2), r2(2)/) self%r = (/r1(2), r2(2)/)
self%weight = SUM(self%r)*5.D-1 self%weight = r2(2)**2 - r1(2)**2
!Normal vector !Normal vector
self%normal = (/ -(self%r(2)-self%r(1)), & self%normal = (/ -(self%r(2)-self%r(1)), &
self%z(2)-self%z(1) , & self%z(2)-self%z(1) , &

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

@ -977,6 +977,7 @@ MODULE moduleMesh
REAL(8):: l2, l, lW, AW REAL(8):: l2, l, lW, AW
REAL(8):: deltaV(1:3), totalDeltaV_ij, normDeltaV REAL(8):: deltaV(1:3), totalDeltaV_ij, normDeltaV
REAL(8):: rnd REAL(8):: rnd
REAL(8):: eps = 1.D-10
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
@ -1011,13 +1012,22 @@ MODULE moduleMesh
velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
l2 = coulombMatrix(k)%l2_j/temperature !If cell temperature is too low, skip particle to avoid division by zero
l = SQRT(l2) IF (temperature>eps) THEN
l2 = coulombMatrix(k)%l2_j/temperature
l = SQRT(l2)
ELSE
partTemp => partTemp%next
CYCLE
END IF
W = partTemp%part%v - velocity W = partTemp%part%v - velocity
normW = NORM2(W) normW = NORM2(W)
IF (normW < 1.D-12) THEN !If relative velocity is too low, skip collision to avoid division by zero and move to next particle
!If relative velocity is too low, skip collision and move to next particle IF (normW < eps) THEN
partTemp => partTemp%next partTemp => partTemp%next
CYCLE CYCLE
@ -1027,7 +1037,6 @@ MODULE moduleMesh
lW = l * normW lW = l * normW
AW = coulombMatrix(k)%A_i/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 = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW)
delta_par_square = AW*density*G(lW) delta_par_square = AW*density*G(lW)
@ -1063,8 +1072,8 @@ MODULE moduleMesh
END DO END DO
!Do scattering of particles from species_j due to species i
IF (i /= j) THEN IF (i /= j) THEN
!Do scattering of particles from species_j due to species i
!Compute background properties of species_i !Compute background properties of species_i
DO n = 1, cell%nNodes DO n = 1, cell%nNodes
node => self%nodes(cellNodes(n))%obj node => self%nodes(cellNodes(n))%obj
@ -1072,10 +1081,11 @@ MODULE moduleMesh
densityNodes(n) = output%density/n_ref densityNodes(n) = output%density/n_ref
velocityNodes(n,1:3) = output%velocity(1:3)/v_ref velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
temperatureNodes(n) = output%temperature/T_ref temperatureNodes(n) = output%temperature/T_ref
END DO END DO
!Divide total momentum exchanged among all the particles of species j !Divide total momentum exchanged among all the particles of species j
!TODO: This is a dirty trick to ensure conservation between species
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)
@ -1087,26 +1097,35 @@ MODULE moduleMesh
velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
l2 = coulombMatrix(k)%l2_i/temperature !If cell temperature is too low, skip particle to avoid division by zero
l = SQRT(l2) IF (temperature>eps) THEN
l2 = coulombMatrix(k)%l2_i/temperature
W = partTemp%part%v - velocity l = SQRT(l2)
normW = NORM2(W)
IF (normW < 1.D-12) THEN ELSE
!If relative velocity is too low, skip collision and move to next particle
partTemp => partTemp%next partTemp => partTemp%next
CYCLE CYCLE
END IF END IF
W = partTemp%part%v - velocity
normW = NORM2(W)
!If relative velocity is too low, skip collision and move to next particle
IF (normW < eps) THEN
partTemp => partTemp%next
CYCLE
END IF
lW = l * normW lW = l * normW
AW = coulombMatrix(k)%A_j/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 = -coulombMatrix(k)%A_j*coulombMatrix(k)%one_plus_massRatio_ji*density*l2*G(lW)
delta_par_square = AW*density*G(lW) delta_par_square = AW*density*G(lW)
delta_per_square = AW*density*H(lW) delta_per_square = AW*density*H(lW)
dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin)
@ -1114,7 +1133,13 @@ MODULE moduleMesh
!Normalize with average exchange per particle !Normalize with average exchange per particle
!TODO: This is a dirty trick to ensure conservation between species !TODO: This is a dirty trick to ensure conservation between species
dW = normDeltaV*dW/NORM2(dW) IF (NORM2(dW) > eps) THEN
dW = normDeltaV*dW/NORM2(dW)
ELSE
dW = 0.D0
END IF
!System of reference for the velocity change !System of reference for the velocity change
!First one is parallel to the relative velocity !First one is parallel to the relative velocity
@ -1127,24 +1152,24 @@ MODULE moduleMesh
!Third one is perpendicular to the other two !Third one is perpendicular to the other two
e3 = crossProduct(e1, e2) e3 = crossProduct(e1, e2)
e3 = normalize(e3) e3 = normalize(e3)
!Random number for direction !Random number for direction
rnd = PI2*random() rnd = PI2*random()
deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3) deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3)
!Change particle velocity !Change particle velocity
partTemp%part%v = partTemp%part%v + deltaV partTemp%part%v = partTemp%part%v + deltaV
partTemp => partTemp%next partTemp => partTemp%next
END DO END DO
END IF END IF
END DO END DO
DEALLOCATE(densityNodes, velocityNodes, temperatureNodes) DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)
END DO END DO
!$OMP END DO !$OMP END DO

2
src/modules/moduleInject.f90
View file

@ -97,7 +97,7 @@ MODULE moduleInject
self%id = i self%id = i
self%vMod = v / v_ref self%vMod = v / v_ref
self%n = n / NORM2(n) self%n = n / NORM2(n)
self%T = T / T_ref self%T = T / T_ref
self%species => species(sp)%obj self%species => species(sp)%obj
tauInject = tau(self%species%n) tauInject = tau(self%species%n)