From 8d35123508d032ed63c9c53f737e0b77155fb175 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Tue, 4 Jul 2023 17:01:02 +0200 Subject: [PATCH] First Coulomb implementation that works 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. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 2 +- src/modules/mesh/moduleMesh.f90 | 73 +++++++++++++++------- src/modules/moduleInject.f90 | 2 +- 3 files changed, 51 insertions(+), 26 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index a2ffb7a..d4baedd 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = SUM(self%r)*5.D-1 + self%weight = r2(2)**2 - r1(2)**2 !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 03dd356..c5f301a 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -977,6 +977,7 @@ MODULE moduleMesh REAL(8):: l2, l, lW, AW REAL(8):: deltaV(1:3), totalDeltaV_ij, normDeltaV REAL(8):: rnd + REAL(8):: eps = 1.D-10 !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) @@ -1011,13 +1012,22 @@ MODULE moduleMesh velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) - l2 = coulombMatrix(k)%l2_j/temperature - l = SQRT(l2) + !If cell temperature is too low, skip particle to avoid division by zero + 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 normW = NORM2(W) - 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 to avoid division by zero and move to next particle + IF (normW < eps) THEN partTemp => partTemp%next CYCLE @@ -1027,7 +1037,6 @@ MODULE moduleMesh 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) @@ -1063,8 +1072,8 @@ MODULE moduleMesh END DO + !Do scattering of particles from species_j due to species i IF (i /= j) THEN - !Do scattering of particles from species_j due to species i !Compute background properties of species_i DO n = 1, cell%nNodes node => self%nodes(cellNodes(n))%obj @@ -1072,10 +1081,11 @@ 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 !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) * & (coulombMatrix(k)%sp_i%weight*coulombMatrix(k)%sp_i%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) temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) - l2 = coulombMatrix(k)%l2_i/temperature - l = SQRT(l2) - - W = partTemp%part%v - velocity - normW = NORM2(W) - IF (normW < 1.D-12) THEN - !If relative velocity is too low, skip collision and move to next particle + !If cell temperature is too low, skip particle to avoid division by zero + IF (temperature>eps) THEN + l2 = coulombMatrix(k)%l2_i/temperature + l = SQRT(l2) + + ELSE partTemp => partTemp%next CYCLE 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 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) @@ -1114,7 +1133,13 @@ MODULE moduleMesh !Normalize with average exchange per particle !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 !First one is parallel to the relative velocity @@ -1127,24 +1152,24 @@ MODULE moduleMesh !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) !Change particle velocity partTemp%part%v = partTemp%part%v + deltaV - + partTemp => partTemp%next - + END DO END IF END DO - DEALLOCATE(densityNodes, velocityNodes, temperatureNodes) + DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) END DO !$OMP END DO diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index e551f0f..4e57083 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -97,7 +97,7 @@ MODULE moduleInject self%id = i self%vMod = v / v_ref - self%n = n / NORM2(n) + self%n = n / NORM2(n) self%T = T / T_ref self%species => species(sp)%obj tauInject = tau(self%species%n)