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

First attempt at Coulomb collisions #46

Merged
JorgeGonz merged 21 commits from feature/CoulombLinear into development 2023-07-16 14:47:59 +02:00
Conversation 2 Commits 21 Files changed 12 +551 -50
12 changed files with 551 additions and 50 deletions
Download patch file Download diff file Expand all files Collapse all files

23
doc/user-manual/bibliography.bib
View file

@ -62,4 +62,27 @@
publisher = {Taylor \& Francis},
}
@Article{sherlock2008monte,
author = {Sherlock, Mark},
journal = {Journal of Computational Physics},
title = {A Monte-Carlo method for Coulomb collisions in hybrid plasma models},
year = {2008},
number = {4},
pages = {2286--2292},
volume = {227},
groups = {Particle-in-cell},
publisher = {Elsevier},
}
@article{lemons2009small,
title={Small-angle Coulomb collision model for particle-in-cell simulations},
author={Lemons, Don S and Winske, Dan and Daughton, William and Albright, Brian},
journal={Journal of Computational Physics},
volume={228},
number={5},
pages={1391--1403},
year={2009},
publisher={Elsevier}
}
@Comment{jabref-meta: databaseType:bibtex;}

BIN
doc/user-manual/fpakc_UserManual.pdf
View file

Binary file not shown.

23
doc/user-manual/fpakc_UserManual.tex
View file

@ -223,8 +223,15 @@
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{\acrlong{cs}}
% Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline.
\subsection{\acrlong{cs}}
A simple linearization of the Coulomb operator based on Ref.~\cite{sherlock2008monte} is implemented.
This method assumes that the species involved in the scattering process have a Maxwellian distribution.
The method is made to conserve momentum and kinetic energy based on the approach in Ref.~\cite{lemons2009small} for self (same species) and intra (different species) collisions.
The user must specify the charged species that will interact together.
The Coulomb logarithm involved in these processes is currently set to a fix value of $10$.
This method is not valid for situations in which the distribution functions are far from Maxwellian.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Reset of particle array}
@ -619,6 +626,8 @@ make
Type of distribution function used to obtain injected particle velocity:
\begin{itemize}
\item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
\item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction.
Only takes into account the positive part of the half-Maxwellian.
\item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction.
\end{itemize}
\item \textbf{T}: Real.
@ -757,6 +766,16 @@ make
Only valid for \textbf{ionization} and \textbf{recombination} processes.
\end{itemize}
\end{itemize}
\item \textbf{Coulomb}: Array of objects.
Contains the information about which species must use the Coulomb linear scattering.
This method assumes a Maxwellian distribution for all species involved.
Each object in the array is defined by:
\begin{itemize}
\item \textbf{species\_i}, \textbf{species\_j}: Character.
Define the two species involved in the collision processes.
Order is indiferent.
\end{itemize}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{parallel}

4
src/fpakc.f90
View file

@ -74,7 +74,7 @@ PROGRAM fpakc
tColl = omp_get_wtime()
!$OMP END SINGLE
IF (doMCC) THEN
IF (doMCCollisions) THEN
CALL meshForMCC%doCollisions(t)
END IF
@ -86,7 +86,7 @@ PROGRAM fpakc
tCoul = omp_get_wTime()
!$OMP END SINGLE
IF (ASSOCIATED(mesh%doCoulomb)) THEN
IF (doCoulombScattering) THEN
CALL mesh%doCoulomb()
END IF

2
src/makefile
View file

@ -4,7 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \
$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o $(OBJDIR)/moduleCoulomb.o \
$(OBJDIR)/moduleMeshInoutCommon.o \
$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \

51
src/modules/init/moduleInput.f90
View file

@ -155,7 +155,7 @@ MODULE moduleInput
sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity
ELSE
sigmaVrel_ref = 0.D0 !Assume no collisions
sigmaVrel_ref = L_ref**2 * v_ref
END IF
@ -412,7 +412,7 @@ MODULE moduleInput
CALL partInitial%add(partNew)
!Assign particle to list in volume
IF (doMCC) THEN
IF (listInCells) THEN
cell => meshforMCC%cells(partNew%cellColl)%obj
CALL OMP_SET_LOCK(cell%lock)
CALL cell%listPart_in(sp)%add(partNew)
@ -613,6 +613,7 @@ MODULE moduleInput
USE moduleSpecies
USE moduleList
USE moduleCollisions
USE moduleCoulomb
USE moduleErrors
USE moduleMesh
USE moduleCaseParam
@ -640,10 +641,11 @@ MODULE moduleInput
!Firstly, check if the object 'interactions' exists
CALL config%info('interactions', found)
IF (found) THEN
!Checks if MC collisions have been defined
!Check if MC collisions have been defined
CALL config%info('interactions.collisions', found)
IF (found) THEN
!Reads collision time step
doMCCollisions = .TRUE.
!Read collision time step
CALL config%info('interactions.timeStep', found)
IF (found) THEN
CALL config%get('interactions.timeStep', tauColl, found)
@ -752,8 +754,35 @@ MODULE moduleInput
END IF
!Check if Coulomb scattering is defined
CALL config%info('interactions.Coulomb', found)
IF (found) THEN
CALL config%info('interactions.Coulomb', found, n_children = nPairs)
IF (nPairs > 0) THEN
nCoulombPairs = nPairs
doCoulombScattering = .TRUE.
ALLOCATE(coulombMatrix(1:nPairs))
DO i = 1, nPairs
WRITE(iString, '(I2)') i
object = 'interactions.Coulomb(' // TRIM(iString) // ')'
CALL config%get(object // '.species_i', species_i, found)
pt_i = speciesName2Index(species_i)
CALL config%get(object // '.species_j', species_j, found)
pt_j = speciesName2Index(species_j)
CALL coulombMatrix(i)%init(pt_i, pt_j)
END DO
END IF
END IF
END IF
listInCells = doMCCollisions .OR. doCoulombScattering
END SUBROUTINE readInteractions
!Reads boundary conditions for the mesh
@ -906,8 +935,6 @@ MODULE moduleInput
END IF
doMCC = ASSOCIATED(meshForMCC)
!Get the dimension of the geometry
CALL config%get(object // '.dimension', mesh%dimen, found)
IF (.NOT. found) THEN
@ -1195,7 +1222,6 @@ MODULE moduleInput
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v
REAL(8), ALLOCATABLE:: T(:), normal(:)
LOGICAL:: fixDirection
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
INTEGER:: physicalSurface
@ -1216,10 +1242,7 @@ MODULE moduleInput
CALL config%get(object // '.v', v, found)
CALL config%get(object // '.T', T, found)
CALL config%get(object // '.n', normal, found)
IF (found) THEN
fixDirection = .TRUE.
ELSE
fixDirection = .FALSE.
IF (.NOT. found) THEN
ALLOCATE(normal(1:3))
normal = 0.D0
END IF
@ -1227,7 +1250,7 @@ MODULE moduleInput
CALL config%get(object // '.units', units, found)
CALL config%get(object // '.physicalSurface', physicalSurface, found)
CALL inject(i)%init(i, v, normal, fixDirection, T, flow, units, sp, physicalSurface)
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
CALL readVelDistr(config, inject(i), object)
@ -1306,6 +1329,10 @@ MODULE moduleInput
T = inj%T(i)
CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
CASE ("Half-Maxwellian")
T = inj%T(i)
CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m)
CASE ("Delta")
v = inj%vMod*inj%n(i)
CALL initVelDistDelta(inj%v(i)%obj)

7
src/modules/makefile
View file

@ -1,6 +1,6 @@
OBJS = common.o output.o mesh.o solver.o init.o \
moduleBoundary.o moduleCollisions.o moduleInject.o \
moduleList.o moduleProbe.o \
moduleList.o moduleProbe.o moduleCoulomb.o \
moduleSpecies.o
all: $(OBJS)
@ -11,7 +11,7 @@ common.o:
output.o: moduleSpecies.o common.o
$(MAKE) -C output all
mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o
mesh.o: moduleCollisions.o moduleCoulomb.o moduleBoundary.o output.o common.o
$(MAKE) -C mesh all
solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
@ -26,6 +26,9 @@ moduleBoundary.o: common.o moduleBoundary.f90
moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleCoulomb.o: moduleSpecies.o common.o moduleCoulomb.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleList.o: common.o moduleSpecies.o moduleList.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

2
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl
r2 = self%n2%getCoordinates()
self%z = (/r1(1), r2(1)/)
self%r = (/r1(2), r2(2)/)
self%weight = SUM(self%r)*5.D-1
self%weight = r2(2)**2 - r1(2)**2
!Normal vector
self%normal = (/ -(self%r(2)-self%r(1)), &
self%z(2)-self%z(1) , &

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

@ -300,7 +300,7 @@ MODULE moduleMesh
END FUNCTION inside_interface
PURE FUNCTION phy2log_interface(self,r) RESULT(Xi)
PURE FUNCTION phy2log_interface(self,r) RESULT(Xi)
IMPORT:: meshCell
CLASS(meshCell), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
@ -387,17 +387,17 @@ MODULE moduleMesh
!Array of boundary elements
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
!Global stiffness matrix
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Permutation matrix for P L U factorization
INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
!PROCEDURES SPECIFIC OF FILE TYPE
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL()
PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
CONTAINS
!GENERIC PROCEDURES
PROCEDURE, PASS:: constructGlobalK
PROCEDURE, PASS:: doCoulomb
END TYPE meshParticles
@ -424,7 +424,7 @@ MODULE moduleMesh
CLASS(meshParticles), INTENT(inout):: self
END SUBROUTINE doCoulomb_interface
!Prints average values
SUBROUTINE printAverage_interface(self)
IMPORT meshParticles
@ -481,7 +481,11 @@ MODULE moduleMesh
!Logical to indicate if an specific mesh for MC Collisions is used
LOGICAL:: doubleMesh
!Logical to indicate if MCC collisions are performed
LOGICAL:: doMCC
LOGICAL:: doMCCollisions = .FALSE.
!Logical to indicate if Coulomb scattering is performed
LOGICAL:: doCoulombScattering = .FALSE.
!Logica to indicate if particles have to be listed in list inside the cells
LOGICAL:: listInCells = .FALSE.
!Complete path for the two meshes
CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle
@ -511,7 +515,7 @@ MODULE moduleMesh
END DO
END DO
DEALLOCATE(n, localK)
END DO
@ -616,7 +620,7 @@ MODULE moduleMesh
tensorS = outerProduct(part%v, part%v)
sp = part%species%n
DO i = 1, nNodes
node => mesh%nodes(cellNodes(i))%obj
CALL OMP_SET_LOCK(node%lock)
@ -650,7 +654,7 @@ MODULE moduleMesh
part%Xi = Xi
part%n_in = .TRUE.
!Assign particle to listPart_in
IF (doMCC) THEN
IF (listInCells) THEN
CALL OMP_SET_LOCK(self%lock)
sp = part%species%n
CALL self%listPart_in(sp)%add(part)
@ -673,7 +677,7 @@ MODULE moduleMesh
CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part)
!If particle is still inside the domain, call findCell
IF (part%n_in) THEN
IF (part%n_in) THEN
IF(PRESENT(oldCell)) THEN
CALL self%findCell(part, oldCell)
@ -719,13 +723,13 @@ MODULE moduleMesh
INTEGER:: sp
found = .FALSE.
cell => meshColl%cells(part%cellColl)%obj
DO WHILE(.NOT. found)
Xi = cell%phy2log(part%r)
IF (cell%inside(Xi)) THEN
part%cellColl = cell%n
IF (doMCC) THEN
IF (listInCells) THEN
CALL OMP_SET_LOCK(cell%lock)
sp = part%species%n
CALL cell%listPart_in(sp)%add(part)
@ -923,7 +927,7 @@ MODULE moduleMesh
END DO
END IF
!Deallocate arrays for next collision
DEALLOCATE(sigmaVrel, probabilityColl)
@ -946,9 +950,308 @@ MODULE moduleMesh
END SUBROUTINE doCollisions
SUBROUTINE doCoulomb(self)
USE moduleCoulomb
USE moduleRandom
USE moduleOutput
USE moduleList
USE moduleMath
USE moduleRefParam
USE moduleConstParam
IMPLICIT NONE
CLASS(meshParticles), INTENT(inout):: self
CLASS(meshParticles), INTENT(in), TARGET:: self
CLASS(meshCell), POINTER:: cell
TYPE(interactionsCoulomb):: pair
INTEGER:: e
INTEGER:: k
INTEGER:: i, j
INTEGER:: n
INTEGER:: p
TYPE(lNode), POINTER:: partTemp
INTEGER(8), ALLOCATABLE:: cellNodes(:)
CLASS(meshNode), POINTER:: node
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), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference
REAL(8):: l, lW, l2
REAL(8):: GlW, HlW
REAL(8):: normC
REAL(8):: cosThe, sinThe
REAL(8):: cosPhi, sinPhi
REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame
REAL(8):: A, AW
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), 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)
DO e = 1, self%numCells
cell => self%cells(e)%obj
cellNodes = cell%getNodes(cell%nNodes)
ALLOCATE(densityNodes(1:cell%nNodes), &
velocityNodes(1:cell%nNodes, 1:3), &
temperatureNodes(1:cell%nNodes))
DO k=1, nCoulombPairs
pair = coulombMatrix(k)
i = pair%sp_i%n
j = pair%sp_j%n
!Do scattering of particles from species_i due to species j
!Compute background properties of species_j
DO n = 1, cell%nNodes
node => self%nodes(cellNodes(n))%obj
CALL calculateOutput(node%output(j), output, node%v, pair%sp_j)
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_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)
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_j/temperature
l = SQRT(l2)
ELSE
partTemp => partTemp%next
CYCLE
END IF
A = pair%A_i*density
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)
!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_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
!Move to the next particle in the list
partTemp => partTemp%next
p = p + 1
END DO
END IF
!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
DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)
END DO
!$OMP END DO
END SUBROUTINE doCoulomb

98
src/modules/moduleCoulomb.f90 Normal file
View file

@ -0,0 +1,98 @@
MODULE moduleCoulomb
USE moduleSpecies
IMPLICIT NONE
INTEGER:: nCoulombPairs = 0
!Type for Coulomb iteraction matrix
TYPE:: interactionsCoulomb
CLASS(speciesGeneric), POINTER:: sp_i
CLASS(speciesGeneric), POINTER:: sp_j
REAL(8):: one_plus_massRatio_ij
REAL(8):: lnCoulomb !This can be done a function in the future
REAL(8):: A_i
REAL(8):: A_j
REAL(8):: l2_j
REAL(8):: l2_i
CONTAINS
PROCEDURE, PASS:: init => initInteractionCoulomb
END TYPE interactionsCoulomb
!Coulomb collision 'matrix'
TYPE(interactionsCoulomb), ALLOCATABLE:: coulombMatrix(:)
CONTAINS
PURE REAL(8) FUNCTION G(x)
USE moduleConstParam
IMPLICIT NONE
REAL(8), INTENT(in):: x
G = 0.D0
IF (x /= 0.D0) THEN
G = (ERF(x) - x*2.D0/SQRT(PI)*EXP(-x**2))/(2.D0*x**2)
END IF
END FUNCTION G
PURE REAL(8) FUNCTION H(x)
IMPLICIT NONE
REAL(8), INTENT(in):: x
H = ERF(x) - G(x)
END FUNCTION H
SUBROUTINE initInteractionCoulomb(self, i, j)
USE moduleSpecies
USE moduleErrors
USE moduleConstParam
USE moduleRefParam
IMPLICIT NONE
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
self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m
Z_i = 0.D0
Z_j = 0.D0
SELECT TYPE(sp => self%sp_i)
TYPE IS (speciesCharged)
Z_i = sp%q
CLASS DEFAULT
CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
END SELECT
SELECT TYPE(sp => self%sp_j)
TYPE IS (speciesCharged)
Z_j = sp%q
CLASS DEFAULT
CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb')
END SELECT
self%lnCoulomb = 10.0 !Make this function dependent
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 !Missing density because it's cell dependent
self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%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_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent
END SUBROUTINE initInteractionCoulomb
END MODULE moduleCoulomb

57
src/modules/moduleInject.f90
View file

@ -35,6 +35,13 @@ MODULE moduleInject
END TYPE velDistMaxwellian
TYPE, EXTENDS(velDistGeneric):: velDistHalfMaxwellian
REAL(8):: vTh !Thermal Velocity
CONTAINS
PROCEDURE, PASS:: randomVel => randomVelHalfMaxwellian
END TYPE velDistHalfMaxwellian
!Dirac's delta distribution function
TYPE, EXTENDS(velDistGeneric):: velDistDelta
CONTAINS
@ -68,7 +75,7 @@ MODULE moduleInject
CONTAINS
!Initialize an injection of particles
SUBROUTINE initInject(self, i, v, n, fixDirection, T, flow, units, sp, physicalSurface)
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
@ -80,7 +87,6 @@ MODULE moduleInject
CLASS(injectGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
LOGICAL, INTENT(in):: fixDirection
INTEGER, INTENT(in):: sp, physicalSurface
REAL(8):: tauInject
REAL(8), INTENT(in):: flow
@ -91,12 +97,7 @@ MODULE moduleInject
self%id = i
self%vMod = v / v_ref
IF (.NOT. fixDirection) THEN
self%n = n / NORM2(n)
ELSE
self%n = 0.D0
END IF
self%fixDirection = fixDirection
self%n = n / NORM2(n)
self%T = T / T_ref
self%species => species(sp)%obj
tauInject = tau(self%species%n)
@ -212,6 +213,16 @@ MODULE moduleInject
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: T, m
velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m))
END SUBROUTINE initVelDistHalfMaxwellian
SUBROUTINE initVelDistDelta(velDist)
IMPLICIT NONE
@ -234,6 +245,22 @@ MODULE moduleInject
END FUNCTION randomVelMaxwellian
!Random velocity from Half Maxwellian distribution
FUNCTION randomVelHalfMaxwellian(self) RESULT (v)
USE moduleRandom
IMPLICIT NONE
CLASS(velDistHalfMaxwellian), INTENT(in):: self
REAL(8):: v
v = 0.D0
DO WHILE (v <= 0.D0)
v = self%vTh*randomMaxwellian()
END DO
END FUNCTION randomVelHalfMaxwellian
!Random velocity from Dirac's delta distribution
PURE FUNCTION randomVelDelta(self) RESULT(v)
IMPLICIT NONE
@ -305,18 +332,20 @@ MODULE moduleInject
!Assign particle type
partInj(n)%species => self%species
IF (self%fixDirection) THEN
direction = self%n
direction = self%n
ELSE
direction = randomEdge%normal
END IF
partInj(n)%v = 0.D0
partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
!If velocity is not in the right direction, invert it
IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN
partInj(n)%v = - partInj(n)%v
END IF
!Obtain natural coordinates of particle in cell
partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)
!Push new particle with the minimum time step

5
src/modules/solver/moduleSolver.f90
View file

@ -70,7 +70,6 @@ MODULE moduleSolver
CHARACTER(:), ALLOCATABLE:: pusherType
REAL(8):: tau, tauSp
!TODO: Reorganize if Cart pushers are combined
SELECT CASE(mesh%dimen)
CASE(0)
self%pushParticle => push0D
@ -322,7 +321,7 @@ MODULE moduleSolver
!$OMP SECTION
!Erase the list of particles inside the cell if particles have been pushed
IF (doMCC) THEN
IF (listInCells) THEN
DO s = 1, nSpecies
DO e = 1, mesh%numCells
IF (solver%pusher(s)%pushSpecies) THEN
@ -456,7 +455,7 @@ MODULE moduleSolver
CALL partWScheme%unsetLock()
!Add particle to cell list
sp = part%species%n
IF (doMCC) THEN
IF (listInCells) THEN
CALL OMP_SET_lock(cell%lock)
CALL cell%listPart_in(sp)%add(newPart)
CALL OMP_UNSET_lock(cell%lock)