From f63e34e2664e97a7a64bd397546b2eca711a3375 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 11 Jul 2023 18:55:20 +0200 Subject: [PATCH] Not fully conservative but works The code is still not fully conservative in intra-species collisions (small error) but at least now is working. I have to test species with different weight. I have to implement a fully conservation for intra-species. --- src/modules/mesh/moduleMesh.f90 | 6 ++++-- src/modules/moduleCoulomb.f90 | 4 +++- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 0f1260e..3dea7c0 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -1055,7 +1055,7 @@ MODULE moduleMesh rotation(1, 1:3) = (/ cosThe*cosPhi, -sinPhi, sinThe*cosPhi /) rotation(2, 1:3) = (/ cosThe*sinPhi, cosPhi, sinThe*sinPhi /) - rotation(3, 1:3) = (/-sinThe, 0.D0, cosThe /) + rotation(3, 1:3) = (/-sinThe, 0.D0, cosThe /) !W at start is = (0, 0, normC), so normW = normC lW = l * normC @@ -1069,6 +1069,7 @@ MODULE moduleMesh deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + !Random angle to distribute perpendicular change in velocity theta_per = PI2*random() !Change W @@ -1076,9 +1077,10 @@ MODULE moduleMesh W(2) = deltaW_per_square * SIN(theta_per) W(3) = normC + deltaW_par + deltaW_par_square - !Change particle velocity + !Update particle velocity and return to laboratory frame partTemp%part%v = velocity + MATMUL(rotation, W) + !Move to the next particle in the list partTemp => partTemp%next END DO diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 7bd3f72..fe2a2f6 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -59,8 +59,10 @@ MODULE moduleCoulomb self%sp_i => species(i)%obj self%sp_j => species(j)%obj - self%one_plus_massRatio_ij = 1.D0 + (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m) + 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