From e05c0d4635c3b2ae362334381de37569317cbfee Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 16 Jul 2023 14:28:07 +0200 Subject: [PATCH] Coulomb Scattering fully conservative Coulomb scattering is now fully conservative thanks to the method in lemos2009small. The trick was to conserve the momentum and energy of ALL particles involved in the scattering in each cell. The substeps in Coulomb collisions have been removed as they are no longer necessary. Still some issues with e-i, but I don't know right now. --- src/modules/init/moduleInput.f90 | 8 +- src/modules/mesh/moduleMesh.f90 | 274 +++++++++++++++++++------------ src/modules/moduleCoulomb.f90 | 8 +- 3 files changed, 175 insertions(+), 115 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index e3927e3..9891806 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -637,7 +637,6 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: electron INTEGER:: e CLASS(meshCell), POINTER:: cell - INTEGER:: subSteps !Firstly, check if the object 'interactions' exists CALL config%info('interactions', found) @@ -771,13 +770,8 @@ MODULE moduleInput pt_i = speciesName2Index(species_i) CALL config%get(object // '.species_j', species_j, found) pt_j = speciesName2Index(species_j) - CALL config%get(object // '.subSteps', subSteps, found) - IF (.NOT. found) THEN - subSteps = 1 - END IF - - CALL coulombMatrix(i)%init(pt_i, pt_j, subSteps) + CALL coulombMatrix(i)%init(pt_i, pt_j) END DO diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 3a65809..7e30326 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -966,7 +966,6 @@ MODULE moduleMesh INTEGER:: k INTEGER:: i, j INTEGER:: n - INTEGER:: t INTEGER:: p TYPE(lNode), POINTER:: partTemp INTEGER(8), ALLOCATABLE:: cellNodes(:) @@ -978,7 +977,6 @@ MODULE moduleMesh REAL(8):: l, lW, l2 REAL(8):: GlW, HlW REAL(8):: normC - REAL(8):: theta, phi !angles between w and c REAL(8):: cosThe, sinThe REAL(8):: cosPhi, sinPhi REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame @@ -986,8 +984,13 @@ 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(:,:) + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij + REAL(8):: massSum_ij + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji + REAL(8):: massSum_ji + REAL(8):: alpha_num, alpha_den, alpha, beta(1:3) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) @@ -1015,9 +1018,12 @@ MODULE moduleMesh END DO - totalP_ij = 0.D0 + ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) !Loop over particles of species_i partTemp => cell%listPart_in(i)%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) @@ -1037,8 +1043,89 @@ MODULE moduleMesh A = pair%A_i*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 * tauMin + 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 + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + !Compute changes in velocity for each particle + deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ij(p) = pair%sp_i%m*partTemp%part%weight + p_ij(p,1:3) = mass_ij(p)*partTemp%part%v + + !Move to the next particle in the list + partTemp => partTemp%next + p = p + 1 + + 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 + + ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) + !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 + C = partTemp%part%v - velocity normC = NORM2(C) @@ -1062,9 +1149,9 @@ MODULE moduleMesh 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() + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin + 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() @@ -1074,105 +1161,90 @@ MODULE moduleMesh W(2) = deltaW_per_square * SIN(theta_per) 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) + !Compute changes in velocity for each particle + deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ji(p) = pair%sp_j%m*partTemp%part%weight + p_ji(p,1:3) = mass_ji(p)*partTemp%part%v - END DO - - !Move to the next particle in the list - partTemp => partTemp%next - - 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 + p = p + 1 END DO END IF - print *, k, NORM2(totalP_ij), NORM2(totalP_ji) + !Calculate correction + !Total mass + massSum_ij = SUM(mass_ij) + massSum_ji = 0.D0 + + !Beta + beta = 0.D0 + DO p = 1, cell%listPart_in(i)%amount + beta = beta + mass_ij(p) * deltaV_ij(p,1:3) + + END DO + + IF (i /= j) THEN + massSum_ji = SUM(mass_ji) + DO p = 1, cell%listPart_in(j)%amount + beta = beta + mass_ji(p) * deltaV_ji(p,1:3) + + END DO + + END IF + + beta = beta / (massSum_ij + massSum_ji) + + !Alpha + alpha_num = 0.D0 + alpha_den = 0.D0 + DO p =1, cell%listPart_in(i)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2 + + END DO + + IF (i /= j) THEN + DO p = 1, cell%listPart_in(j)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2 + + END DO + + END IF + + alpha = -2.D0*alpha_num / alpha_den + + !Apply correction to particles velocity + partTemp => cell%listPart_in(i)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + IF (i /= j) THEN + partTemp => cell%listPart_in(j)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + END IF + + DEALLOCATE(deltaV_ij, p_ij, mass_ij) + + IF (i /= j) THEN + DEALLOCATE(deltaV_ji, p_ji, mass_ji) + + END IF END DO diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 5faa277..166b6b2 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -14,8 +14,6 @@ MODULE moduleCoulomb REAL(8):: A_j REAL(8):: l2_j REAL(8):: l2_i - REAL(8):: tauSub - INTEGER:: nSubSteps CONTAINS PROCEDURE, PASS:: init => initInteractionCoulomb @@ -48,7 +46,7 @@ MODULE moduleCoulomb END FUNCTION H - SUBROUTINE initInteractionCoulomb(self, i, j, subSteps) + SUBROUTINE initInteractionCoulomb(self, i, j) USE moduleSpecies USE moduleErrors USE moduleConstParam @@ -56,14 +54,10 @@ MODULE moduleCoulomb IMPLICIT NONE CLASS(interactionsCoulomb), INTENT(out):: self - INTEGER, INTENT(in):: subSteps INTEGER, INTENT(in):: i, j REAL(8):: Z_i, Z_j REAL(8):: scaleFactor - self%nSubSteps = subSteps - self%tauSub = tauMin / REAL(self%nSubSteps) - self%sp_i => species(i)%obj self%sp_j => species(j)%obj