From 6113ac3305799e734b6032b671a696ffc9fe63ab Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Mon, 6 Mar 2023 16:16:17 +0100 Subject: [PATCH] Correction with conservation Now the method is much better in conserving total energy. However, still there is an issue with collisions between species of dispaprate mass. --- src/modules/mesh/moduleMesh.f90 | 91 ++++++++++++++++++--------------- src/modules/moduleCoulomb.f90 | 15 +++--- 2 files changed, 59 insertions(+), 47 deletions(-) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 6607fd0..478fd83 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -966,7 +966,6 @@ MODULE moduleMesh INTEGER:: i, j INTEGER:: n TYPE(lNode), POINTER:: partTemp - REAL(8):: W(3), dW(2) !Relative velocity between particle and species and its increment INTEGER(8), ALLOCATABLE:: cellNodes(:) CLASS(meshNode), POINTER:: node TYPE(outputFormat):: output @@ -974,7 +973,8 @@ MODULE moduleMesh REAL(8):: density, velocity(1:3), temperature!values at particle position REAL(8), DIMENSION(1:3):: e1, e2, e3 REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square - REAL(8):: l, lW, AW + REAL(8):: W(3), dW(2), normW !Relative velocity between particle and species and its increment + REAL(8):: l2, l, lW, AW REAL(8):: rnd @@ -1007,12 +1007,25 @@ MODULE moduleMesh 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) - l = coulombMatrix(k)%l_j/SQRT(temperature) + + l2 = coulombMatrix(k)%l2_j/temperature + l = SQRT(l2) W = partTemp%part%v - velocity - lW = l * NORM2(W) - AW = coulombMatrix(k)%A_ij/NORM2(W) - !Axis of the relative velocity + normW = NORM2(W) + lW = l * normW + AW = coulombMatrix(k)%A_i/normW + + delta_par = -coulombMatrix(k)%A_i*coulombMatrix(k)%one_plus_massRatio_ij*density*l2*G(lW) + + delta_par_square = AW*density*G(lW) + + delta_per_square = AW*density*H(lW) + + dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) + dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) + + !System of reference for the velocity change !First one is parallel to the relative velocity e1 = normalize(W) !Second one is perpendicular to it @@ -1024,18 +1037,12 @@ MODULE moduleMesh e3 = crossProduct(e2, e1) e3 = normalize(e3) - delta_par = -coulombMatrix(k)%A_ij*coulombMatrix(k)%one_plus_massRatio_ij*density*l**2*G(lW) + !Random number for direction + rnd = PI2*random() - delta_par_square = AW*density*G(lW) - - delta_per_square = AW*density*H(lW) - - dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) - dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) - - rnd = random() - partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + & - SIN(PI2*rnd)*e3) + !Change particle velocity + partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + & + SIN(rnd)*e3) partTemp => partTemp%next @@ -1049,20 +1056,33 @@ MODULE moduleMesh densityNodes(n) = output%density/n_ref velocityNodes(n,1:3) = output%velocity(1:3)/v_ref temperatureNodes(n) = output%temperature/T_ref - + END DO - + partTemp => cell%listPart_in(j)%head 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) - l = coulombMatrix(k)%l_i/SQRT(temperature) - + + l2 = coulombMatrix(k)%l2_i/temperature + l = SQRT(l2) + W = partTemp%part%v - velocity - lW = l * NORM2(W) - AW = coulombMatrix(k)%A_ji/NORM2(W) - !Axis of the relative velocity + normW = NORM2(W) + lW = l * normW + AW = coulombMatrix(k)%A_j/normW + + delta_par = -coulombMatrix(k)%A_j*coulombMatrix(k)%one_plus_massRatio_ji*density*l2*G(lW) + + delta_par_square = AW*density*G(lW) + + delta_per_square = AW*density*H(lW) + + dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) + dW(2) = ABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) + + !System of reference for the velocity change !First one is parallel to the relative velocity e1 = normalize(W) !Second one is perpendicular to it @@ -1073,24 +1093,15 @@ MODULE moduleMesh !Third one is perpendicular to the other two e3 = crossProduct(e2, e1) e3 = normalize(e3) - - delta_par = -coulombMatrix(k)%A_ji*coulombMatrix(k)%one_plus_massRatio_ji*density*l**2*G(lW) - - delta_par_square = AW*density*G(lW) - - delta_per_square = AW*density*H(lW) - - dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) - dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) - - rnd = random() - partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + & - SIN(PI2*rnd)*e3) - + + rnd = PI2*random() + partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(rnd)*e2 + & + SIN(rnd)*e3) + partTemp => partTemp%next - + END DO - + END IF END DO diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 index 734b898..775caad 100644 --- a/src/modules/moduleCoulomb.f90 +++ b/src/modules/moduleCoulomb.f90 @@ -10,8 +10,8 @@ MODULE moduleCoulomb CLASS(speciesGeneric), POINTER:: sp_j REAL(8):: one_plus_massRatio_ij, one_plus_massRatio_ji REAL(8):: lnCoulomb !This can be done a function in the future - REAL(8):: A_ij, A_ji - REAL(8):: l_j, l_i + REAL(8):: A_i, A_j + REAL(8):: l2_j, l2_i CONTAINS PROCEDURE, PASS:: init => initInteractionCoulomb @@ -79,12 +79,13 @@ MODULE moduleCoulomb END SELECT - self%lnCoulomb = 12.0 - self%A_ij = 8.D0*PI*Z_i**2*Z_j**2*e_ref**4*self%lnCoulomb / self%sp_i%m - self%A_ji = 8.D0*PI*Z_j**2*Z_i**2*e_ref**4*self%lnCoulomb / self%sp_j%m + self%lnCoulomb = 12.0 !Make this function dependent - self%l_j = SQRT(self%sp_j%m / 2.D0) !Missing temperature because it's cell dependent - self%l_i = SQRT(self%sp_i%m / 2.D0) !Missing temperature because it's cell 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 + + 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