From c3a6f77ffc36f2ca5249b3a2c1290f10d30c409c Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Wed, 12 Jul 2023 15:17:26 +0200 Subject: [PATCH] Combining ij - ji collisions In an attempt to make the operator fully conservarive I have combined ij and ji collisions (when i/=j). Now the matter is to find a way that makes this conserve momentum and energy for intraspecies. --- src/modules/mesh/moduleMesh.f90 | 94 +++++++++++++++++++++++++++++++++ src/modules/moduleCoulomb.f90 | 4 ++ 2 files changed, 98 insertions(+) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 81f53c3..3a65809 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -967,6 +967,7 @@ MODULE moduleMesh INTEGER:: i, j INTEGER:: n INTEGER:: t + INTEGER:: p TYPE(lNode), POINTER:: partTemp INTEGER(8), ALLOCATABLE:: cellNodes(:) CLASS(meshNode), POINTER:: node @@ -985,6 +986,8 @@ MODULE moduleMesh REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W REAL(8):: theta_per !Random angle for perpendicular direction REAL(8):: eps = 1.D-12 + REAL(8):: preV(1:3), totalP_ij(1:3), totalP_ji(1:3) + REAL(8), ALLOCATABLE:: deltaV_ji(:,:) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) @@ -1012,6 +1015,7 @@ MODULE moduleMesh END DO + totalP_ij = 0.D0 !Loop over particles of species_i partTemp => cell%listPart_in(i)%head DO WHILE(ASSOCIATED(partTemp)) @@ -1071,7 +1075,9 @@ MODULE moduleMesh W(3) = normC + deltaW_par + deltaW_par_square !Update particle velocity and return to laboratory frame + preV = partTemp%part%v partTemp%part%v = MATMUL(rotation, W) + velocity + totalP_ij = totalP_ij + pair%sp_i%m*(partTemp%part%v - preV) END DO @@ -1080,6 +1086,94 @@ MODULE moduleMesh END DO + !Do corresponding collisions + 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 + CALL calculateOutput(node%output(i), output, node%v, pair%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 + + totalP_ji = 0.D0 + ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount,1:3)) + !Loop over particles of species_j + partTemp => cell%listPart_in(j)%head + p = 1 + 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 = pair%l2_i/temperature + l = SQRT(l2) + + ELSE + partTemp => partTemp%next + + CYCLE + + END IF + A = pair%A_j*density + + !Do the required substeps + DO t = 1, pair%nSubSteps + C = partTemp%part%v - velocity + normC = NORM2(C) + + !C_3 = z; C_1, C2 = x, y (per) + C_per = NORM2(C(1:2)) + cosPhi = C(1) / C_per + sinPhi = C(2) / C_per + cosThe = C(3) / normC + sinThe = C_per / normC + + !Rotation matrix to go from W to C + rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column + -sinPhi, cosPhi, 0.D0, & !Second column + sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column + (/ 3, 3 /)) + + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / normC + + !Calculate changes in W due to collision process + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * pair%tauSub + deltaW_par_square = SQRT(AW * GlW * pair%tauSub)*randomMaxwellian() + deltaW_per_square = SQRT(AW * HlW * pair%tauSub)*randomMaxwellian() + + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() + + !Change W + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + preV = partTemp%part%v + partTemp%part%v = MATMUL(rotation, W) + velocity + totalP_ji = totalP_ji + pair%sp_j%m*(partTemp%part%v - preV) + + END DO + + !Move to the next particle in the list + partTemp => partTemp%next + + END DO + + END IF + + print *, k, NORM2(totalP_ij), NORM2(totalP_ji) + END DO DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 40e85c2..5faa277 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -11,7 +11,9 @@ MODULE moduleCoulomb REAL(8):: one_plus_massRatio_ij REAL(8):: lnCoulomb !This can be done a function in the future REAL(8):: A_i + REAL(8):: A_j REAL(8):: l2_j + REAL(8):: l2_i REAL(8):: tauSub INTEGER:: nSubSteps CONTAINS @@ -92,8 +94,10 @@ MODULE moduleCoulomb scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3) self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent + self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent 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