diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 9891806..e3927e3 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -637,6 +637,7 @@ 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) @@ -770,8 +771,13 @@ 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 - CALL coulombMatrix(i)%init(pt_i, pt_j) + END IF + + CALL coulombMatrix(i)%init(pt_i, pt_j, subSteps) END DO diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 234ae3c..81f53c3 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -961,10 +961,12 @@ MODULE moduleMesh CLASS(meshParticles), INTENT(in), TARGET:: self CLASS(meshCell), POINTER:: cell + TYPE(interactionsCoulomb):: pair INTEGER:: e INTEGER:: k INTEGER:: i, j INTEGER:: n + INTEGER:: t TYPE(lNode), POINTER:: partTemp INTEGER(8), ALLOCATABLE:: cellNodes(:) CLASS(meshNode), POINTER:: node @@ -995,14 +997,15 @@ MODULE moduleMesh temperatureNodes(1:cell%nNodes)) DO k=1, nCoulombPairs - i = coulombMatrix(k)%sp_i%n - j = coulombMatrix(k)%sp_j%n + pair = coulombMatrix(k) + i = pair%sp_i%n + j = pair%sp_j%n !Do scattering of particles from species_i due to species j !Compute background properties of species_j DO n = 1, cell%nNodes node => self%nodes(cellNodes(n))%obj - CALL calculateOutput(node%output(j), output, node%v, coulombMatrix(k)%sp_j) + CALL calculateOutput(node%output(j), output, node%v, pair%sp_j) densityNodes(n) = output%density/n_ref velocityNodes(n,1:3) = output%velocity(1:3)/v_ref temperatureNodes(n) = output%temperature/T_ref @@ -1018,7 +1021,7 @@ MODULE moduleMesh !If cell temperature is too low, skip particle to avoid division by zero IF (temperature>eps) THEN - l2 = coulombMatrix(k)%l2_j/temperature + l2 = pair%l2_j/temperature l = SQRT(l2) ELSE @@ -1028,52 +1031,49 @@ MODULE moduleMesh END IF - C = partTemp%part%v - velocity - normC = NORM2(C) + A = pair%A_i*density - !If relative velocity is too low, skip collision to avoid division by zero and move to next particle - IF (normC < eps) THEN - partTemp => partTemp%next + !Do the required substeps + DO t = 1, pair%nSubSteps + C = partTemp%part%v - velocity + normC = NORM2(C) - CYCLE + !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 - END IF + !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 /)) - !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 + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / 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 /)) + !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() - !W at start is = (0, 0, normC), so normW = normC - lW = l * normC - GlW = G(lW) - HlW = H(lW) - A = coulombMatrix(k)%A_i*density - AW = A / normC + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() - !Calculate changes in W due to collision process - deltaW_par = - A * coulombMatrix(k)%one_plus_massRatio_ij * l2 * GlW * tauMin - deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() - deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + !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 - !Random angle to distribute perpendicular change in velocity - theta_per = PI2*random() + !Update particle velocity and return to laboratory frame + partTemp%part%v = MATMUL(rotation, W) + velocity - !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 - - !Update particle velocity and return to laboratory frame - partTemp%part%v = MATMUL(rotation, W) + velocity + END DO !Move to the next particle in the list partTemp => partTemp%next diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 653eaa4..40e85c2 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -12,6 +12,8 @@ MODULE moduleCoulomb REAL(8):: lnCoulomb !This can be done a function in the future REAL(8):: A_i REAL(8):: l2_j + REAL(8):: tauSub + INTEGER:: nSubSteps CONTAINS PROCEDURE, PASS:: init => initInteractionCoulomb @@ -44,7 +46,7 @@ MODULE moduleCoulomb END FUNCTION H - SUBROUTINE initInteractionCoulomb(self, i, j) + SUBROUTINE initInteractionCoulomb(self, i, j, subSteps) USE moduleSpecies USE moduleErrors USE moduleConstParam @@ -52,10 +54,14 @@ 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