JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Still unsure, but things fixed

Browse source 
There was an issue with the calculation of theta and phi for the
rotation from W to C. This was causing some velocities not being
correct.

Now the angles are properly computed. Still unsure about the e-i
collisions as they seem to be quite small. Probably a numerical issue
with the mass ratios still exists.
This commit is contained in:
Jorge Gonzalez 2023-07-12 11:38:12 +02:00
commit a891360b7a
2 changed files with 15 additions and 20 deletions
Download patch file Download diff file Expand all files Collapse all files

31
src/modules/mesh/moduleMesh.f90
View file

@ -982,7 +982,7 @@ MODULE moduleMesh
REAL(8):: A, AW REAL(8):: A, AW
REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W
REAL(8):: theta_per !Random angle for perpendicular direction REAL(8):: theta_per !Random angle for perpendicular direction
REAL(8):: eps = 1.D-10 REAL(8):: eps = 1.D-12
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
@ -1039,23 +1039,18 @@ MODULE moduleMesh
END IF END IF
theta = ACOS(C(3) / normC) !C_3 = z; C_1, C2 = x, y (per)
cosThe = COS(theta) C_per = NORM2(C(1:2))
sinThe = SIN(theta) cosPhi = C(1) / C_per
C_per = SQRT(C(1)**2 + C(2)**2) sinPhi = C(2) / C_per
IF (C_per > eps) THEN cosThe = C(3) / normC
phi = SIGN(1.D0, C(1)) * ACOS(C(1) / C_per) sinThe = C_per / normC
ELSE !Rotation matrix to go from W to C
phi = 0.D0 rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column
-sinPhi, cosPhi, 0.D0, & !Second column
END IF sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column
cosPhi = COS(phi) (/ 3, 3 /))
sinPhi = SIN(phi)
rotation(1, 1:3) = (/ cosThe*cosPhi, -sinPhi, sinThe*cosPhi /)
rotation(2, 1:3) = (/ cosThe*sinPhi, cosPhi, sinThe*sinPhi /)
rotation(3, 1:3) = (/-sinThe, 0.D0, cosThe /)
!W at start is = (0, 0, normC), so normW = normC !W at start is = (0, 0, normC), so normW = normC
lW = l * normC lW = l * normC
@ -1078,7 +1073,7 @@ MODULE moduleMesh
W(3) = normC + deltaW_par + deltaW_par_square W(3) = normC + deltaW_par + deltaW_par_square
!Update particle velocity and return to laboratory frame !Update particle velocity and return to laboratory frame
partTemp%part%v = velocity + MATMUL(rotation, W) partTemp%part%v = MATMUL(rotation, W) + velocity
!Move to the next particle in the list !Move to the next particle in the list
partTemp => partTemp%next partTemp => partTemp%next

4
src/modules/moduleCoulomb.f90
View file

@ -83,9 +83,9 @@ MODULE moduleCoulomb
self%lnCoulomb = 10.0 !Make this function dependent self%lnCoulomb = 10.0 !Make this function dependent
scaleFactor = (n_ref * qe**4) / (eps_0**2 * m_ref**2 * v_ref**3) * ti_ref scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3)
self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent
self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent