diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index c5f301a..b8d7501 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -975,7 +975,7 @@ MODULE moduleMesh 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), totalDeltaV_ij, normDeltaV + REAL(8):: deltaV(1:3) REAL(8):: rnd REAL(8):: eps = 1.D-10 @@ -1005,7 +1005,6 @@ MODULE moduleMesh END DO !Loop over particles of species_i - totalDeltaV_ij = 0.D0 partTemp => cell%listPart_in(i)%head DO WHILE(ASSOCIATED(partTemp)) density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) @@ -1026,6 +1025,7 @@ MODULE moduleMesh W = partTemp%part%v - velocity normW = NORM2(W) + !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 @@ -1063,8 +1063,6 @@ MODULE moduleMesh deltaV = dW(1)*e1 + dW(2)*(COS(rnd)*e2 + SIN(rnd)*e3) - totalDeltaV_ij = totalDeltaV_ij + NORM2(deltaV) - !Change particle velocity partTemp%part%v = partTemp%part%v + deltaV @@ -1072,101 +1070,6 @@ MODULE moduleMesh END DO - !Do scattering of particles from species_j due to species i - IF (i /= j) THEN - !Compute background properties of species_i - DO n = 1, cell%nNodes - node => self%nodes(cellNodes(n))%obj - CALL calculateOutput(node%output(i), output, node%v, coulombMatrix(k)%sp_i) - 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) - - !Loop over particles of species_j - partTemp => cell%listPart_in(j)%head - DO WHILE(ASSOCIATED(partTemp)) - density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) - velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) - temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) - - !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) - dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) - - !Normalize with average exchange per particle - !TODO: This is a dirty trick to ensure conservation between species - 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 - 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) - - !Change particle velocity - partTemp%part%v = partTemp%part%v + deltaV - - partTemp => partTemp%next - - END DO - - END IF - END DO DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index ab931e0..683e626 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -8,10 +8,10 @@ MODULE moduleCoulomb TYPE:: interactionsCoulomb CLASS(speciesGeneric), POINTER:: sp_i CLASS(speciesGeneric), POINTER:: sp_j - REAL(8):: one_plus_massRatio_ij, one_plus_massRatio_ji + REAL(8):: one_plus_massRatio_ij REAL(8):: lnCoulomb !This can be done a function in the future - REAL(8):: A_i, A_j - REAL(8):: l2_j, l2_i + REAL(8):: A_i + REAL(8):: l2_j CONTAINS PROCEDURE, PASS:: init => initInteractionCoulomb @@ -60,7 +60,6 @@ MODULE moduleCoulomb 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_ji = 1.D0 + (self%sp_j%weight*self%sp_j%m)/(self%sp_i%weight*self%sp_i%m) SELECT TYPE(sp => self%sp_i) TYPE IS (speciesCharged) @@ -85,10 +84,8 @@ MODULE moduleCoulomb scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref self%A_i = 2.D0*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2 * scaleFactor - self%A_j = 2.D0*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2 * scaleFactor self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent - self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent END SUBROUTINE initInteractionCoulomb