From fe94615a275abd17ff23cc62c2eb4ecd4cfbaa32 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Wed, 8 Mar 2023 16:37:45 +0100 Subject: [PATCH] First conservative implementation of Coulomb I am doing a trick in which I ensure that energy is conserved for Coulomb collisions. This was not happening and what an issue for different mass ratios. Still, this can cause an issue on getting the right relaxation rates, still necessary to check it. --- src/modules/common/moduleRefParam.f90 | 2 +- src/modules/init/moduleInput.f90 | 1 - src/modules/mesh/moduleMesh.f90 | 29 ++++++++++++++++++--------- src/modules/moduleCoulomb.f90 | 7 +++++-- 4 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/modules/common/moduleRefParam.f90 b/src/modules/common/moduleRefParam.f90 index e8819ce..051190b 100644 --- a/src/modules/common/moduleRefParam.f90 +++ b/src/modules/common/moduleRefParam.f90 @@ -3,7 +3,7 @@ MODULE moduleRefParam !Parameters that define the problem (inputs) REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref !Reference parameters for non-dimensional problem - REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref, e_ref + REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref END MODULE moduleRefParam diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 7af67db..c8beb0f 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -144,7 +144,6 @@ MODULE moduleInput EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field Volt_ref = EF_ref*L_ref !reference voltage B_ref = m_ref / (ti_ref * qe) !reference magnetic field - e_ref = qe/(m_ref*v_ref**2) !reference charge !If a reference cross section is given, it is used CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index c80eec3..a62cfb2 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -975,10 +975,11 @@ 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):: rnd - !$OMP SINGLE + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) DO e = 1, self%numCells cell => self%cells(e)%obj cellNodes = cell%getNodes(cell%nNodes) @@ -1003,6 +1004,7 @@ 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) @@ -1036,16 +1038,18 @@ MODULE moduleMesh e2(3) = 0.D0 e2 = normalize(e2) !Third one is perpendicular to the other two - e3 = crossProduct(e2, e1) + e3 = crossProduct(e1, e2) e3 = normalize(e3) !Random number for direction rnd = PI2*random() - !Change particle velocity - partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + & - SIN(rnd)*e3) + 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 + partTemp => partTemp%next END DO @@ -1062,6 +1066,8 @@ MODULE moduleMesh END DO + !Divide total momentum exchanged among all the particles of species j + normDeltaV = totalDeltaV_ij / REAL(cell%listPart_in(j)%amount) * (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m) !Loop over particles of species_j partTemp => cell%listPart_in(j)%head DO WHILE(ASSOCIATED(partTemp)) @@ -1086,6 +1092,10 @@ MODULE moduleMesh 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, it should not be neccessary but I don't lnow why it is. + dW = normDeltaV*dW/NORM2(dW) + !System of reference for the velocity change !First one is parallel to the relative velocity e1 = normalize(W) @@ -1095,15 +1105,16 @@ MODULE moduleMesh e2(3) = 0.D0 e2 = normalize(e2) !Third one is perpendicular to the other two - e3 = crossProduct(e2, e1) + 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 + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + & - SIN(rnd)*e3) + partTemp%part%v = partTemp%part%v + deltaV partTemp => partTemp%next @@ -1114,7 +1125,7 @@ MODULE moduleMesh END DO END DO - !$OMP END SINGLE + !$OMP END DO END SUBROUTINE doCoulomb diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 85f4f49..ab931e0 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -54,6 +54,7 @@ MODULE moduleCoulomb CLASS(interactionsCoulomb), INTENT(out):: self INTEGER, INTENT(in):: i, j REAL(8):: Z_i, Z_j + REAL(8):: scaleFactor self%sp_i => species(i)%obj self%sp_j => species(j)%obj @@ -81,8 +82,10 @@ MODULE moduleCoulomb self%lnCoulomb = 10.0 !Make this function dependent - self%A_i = 8.D0*PI*Z_i**2*Z_j**2*self%lnCoulomb / self%sp_i%m**2 - self%A_j = 8.D0*PI*Z_j**2*Z_i**2*self%lnCoulomb / self%sp_j%m**2 + 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