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

Small changes

Browse source 
Just some small changes to the code to improve its quality.

Nothing regarding conservation was done yet.
This commit is contained in:
Jorge Gonzalez 2023-07-11 11:22:19 +02:00
commit 63fc2842be
1 changed files with 18 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -971,9 +971,10 @@ MODULE moduleMesh
TYPE(outputFormat):: output
REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node
REAL(8):: density, velocity(1:3), temperature!values at particle position
REAL(8):: C(1:3), W(1:3) !relative velocity and velocity in the relative frame of reference
REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference
REAL(8):: l, lW, l2
REAL(8):: normC, normW
REAL(8):: GlW, HlW
REAL(8):: normC
REAL(8):: theta, phi !angles between w and c
REAL(8):: cosThe, sinThe
REAL(8):: cosPhi, sinPhi
@ -1038,11 +1039,17 @@ MODULE moduleMesh
END IF
theta = ACOS(C(3) / normC)
cosThe = COS(theta)
sinThe = SIN(theta)
phi = SIGN(1.D0, C(1)) * ACOS(C(1) / SQRT(C(1)**2 + C(2)**2))
C_per = SQRT(C(1)**2 + C(2)**2)
IF (C_per > eps) THEN
phi = SIGN(1.D0, C(1)) * ACOS(C(1) / C_per)
ELSE
phi = 0.D0
END IF
cosPhi = COS(phi)
sinPhi = SIN(phi)
@ -1052,14 +1059,18 @@ MODULE moduleMesh
!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
deltaW_par = - A * coulombMatrix(k)%one_plus_massRatio_ij * l2 * G(lW) * tauMin
deltaW_par_square = SQRT(AW * G(lW) * tauMin)*randomMaxwellian()
deltaW_per_square = SQRT(AW * H(lW) * tauMin)*randomMaxwellian()
!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()
theta_per = PI2*random()
!Change W
W(1) = deltaW_per_square * COS(theta_per)
W(2) = deltaW_per_square * SIN(theta_per)