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