941 lines
29 KiB
Fortran
941 lines
29 KiB
Fortran
submodule(moduleMesh) elements
|
|
CONTAINS
|
|
!Reset the output of node
|
|
PURE module SUBROUTINE resetOutput(self)
|
|
USE moduleSpecies
|
|
USE moduleOutput
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshNode), INTENT(inout):: self
|
|
INTEGER:: k
|
|
|
|
DO k = 1, nSpecies
|
|
self%output(k)%den = 0.D0
|
|
self%output(k)%mom = 0.D0
|
|
self%output(k)%tensorS = 0.D0
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE resetOutput
|
|
|
|
module subroutine meshNodePointer_add(self, node)
|
|
implicit none
|
|
|
|
class(meshNodePointer), allocatable, intent(inout):: self(:)
|
|
integer, intent(in):: node
|
|
integer:: n
|
|
logical:: inArray
|
|
type(meshNodePointer):: nodeToAdd
|
|
|
|
|
|
nodeToAdd%obj => mesh%nodes(node)%obj
|
|
inArray = .false.
|
|
! I cannot use the procedure 'any' for this
|
|
do n = 1 , size(self)
|
|
if (self(n) == nodeToAdd) then
|
|
! Node already in array
|
|
inArray = .true.
|
|
exit
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
if (.not. inArray) then
|
|
! If not, add element to array of nodes
|
|
self = [self, nodeToAdd]
|
|
|
|
end if
|
|
|
|
end subroutine meshNodePointer_add
|
|
|
|
module function meshNodePointer_equal_type_type(self, other) result(isEqual)
|
|
implicit none
|
|
|
|
class(meshNodePointer), intent(in):: self, other
|
|
logical:: isEqual
|
|
|
|
isEqual = self%obj%n == other%obj%n
|
|
|
|
end function meshNodePointer_equal_type_type
|
|
|
|
module function meshNodePointer_equal_type_int(self, other) result(isEqual)
|
|
implicit none
|
|
|
|
class(meshNodePointer), intent(in):: self
|
|
integer, intent(in):: other
|
|
logical:: isEqual
|
|
|
|
isEqual = self%obj%n == other
|
|
|
|
end function meshNodePointer_equal_type_int
|
|
|
|
module function meshEdgePointer_equal_type_type(self, other) result(isEqual)
|
|
implicit none
|
|
|
|
class(meshEdgePointer), intent(in):: self, other
|
|
logical:: isEqual
|
|
|
|
isEqual = self%obj%n == other%obj%n
|
|
|
|
end function meshEdgePointer_equal_type_type
|
|
|
|
module subroutine meshEdgePointer_add(self, edge)
|
|
implicit none
|
|
|
|
class(meshEdgePointer), allocatable, intent(inout):: self(:)
|
|
integer, intent(in):: edge
|
|
integer:: n
|
|
logical:: inArray
|
|
type(meshEdgePointer):: edgeToAdd
|
|
|
|
|
|
edgeToAdd%obj => mesh%edges(edge)%obj
|
|
inArray = .false.
|
|
! I cannot use the procedure 'any' for this
|
|
do n = 1 , size(self)
|
|
if (self(n) == edgeToAdd) then
|
|
! Node already in array
|
|
inArray = .true.
|
|
exit
|
|
|
|
end if
|
|
|
|
end do
|
|
|
|
if (.not. inArray) then
|
|
! If not, add element to array of nodes
|
|
self = [self, edgeToAdd]
|
|
|
|
end if
|
|
|
|
end subroutine meshEdgePointer_add
|
|
|
|
module function meshEdgePointer_equal_type_int(self, other) result(isEqual)
|
|
implicit none
|
|
|
|
class(meshEdgePointer), intent(in):: self
|
|
integer, intent(in):: other
|
|
logical:: isEqual
|
|
|
|
isEqual = self%obj%n == other
|
|
|
|
end function meshEdgePointer_equal_type_int
|
|
|
|
!Constructs the global K matrix
|
|
module SUBROUTINE constructGlobalK(self)
|
|
use moduleErrors, only: criticalError
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshParticles), INTENT(inout):: self
|
|
INTEGER:: c
|
|
INTEGER, ALLOCATABLE:: nodes(:)
|
|
REAL(8), ALLOCATABLE:: localK(:,:)
|
|
INTEGER:: i, j
|
|
integer:: n, b, ni
|
|
INTEGER:: info
|
|
EXTERNAL:: dgetrf
|
|
|
|
DO c = 1, self%numCells
|
|
associate(nNodes => self%cells(c)%obj%nNodes)
|
|
ALLOCATE(nodes(1:nNodes))
|
|
ALLOCATE(localK(1:nNodes, 1:nNodes))
|
|
nodes = self%cells(c)%obj%getNodes(nNodes)
|
|
localK = self%cells(c)%obj%elemK(nNodes)
|
|
|
|
DO i = 1, nNodes
|
|
DO j = 1, nNodes
|
|
self%K(nodes(i), nodes(j)) = self%K(nodes(i), nodes(j)) + localK(i, j)
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
DEALLOCATE(nodes, localK)
|
|
|
|
end associate
|
|
|
|
END DO
|
|
|
|
! Modify K matrix due to EM boundary conditions
|
|
DO b = 1, nBoundariesEM
|
|
SELECT TYPE(boundary => boundariesEM(b)%obj)
|
|
TYPE IS(boundaryEMDirichlet)
|
|
DO n = 1, boundary%nNodes
|
|
ni = boundary%nodes(n)%obj%n
|
|
self%K(ni, :) = 0.D0
|
|
self%K(ni, ni) = 1.D0
|
|
|
|
END DO
|
|
|
|
TYPE IS(boundaryEMDirichletTime)
|
|
DO n = 1, boundary%nNodes
|
|
ni = boundary%nodes(n)%obj%n
|
|
self%K(ni, :) = 0.D0
|
|
self%K(ni, ni) = 1.D0
|
|
|
|
END DO
|
|
|
|
TYPE IS(boundaryEMFloating)
|
|
DO n = 1, boundary%nNodes
|
|
ni = boundary%nodes(n)%obj%n
|
|
self%K(ni, :) = 0.D0
|
|
self%K(ni, ni) = 1.D0
|
|
|
|
END DO
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
!Compute the PLU factorization of K once boundary conditions have been read
|
|
CALL dgetrf(self%numNodes, self%numNodes, self%K, self%numNodes, self%IPIV, info)
|
|
IF (info /= 0) THEN
|
|
CALL criticalError('Factorization of K matrix failed', 'constructGlobalK')
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE constructGlobalK
|
|
|
|
! Gather the value of valNodes at position Xi of an edge
|
|
pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f)
|
|
implicit none
|
|
|
|
class(meshEdge), intent(in):: self
|
|
real(8), intent(in):: Xi(1:3)
|
|
integer, intent(in):: nNodes
|
|
real(8), intent(in):: valNodes(1:nNodes)
|
|
real(8):: f
|
|
real(8):: fPsi(1:nNodes)
|
|
|
|
fPsi = self%fPsi(Xi, nNodes)
|
|
f = dot_product(fPsi, valNodes)
|
|
|
|
end function gatherF_edge_scalar
|
|
|
|
!Gather the value of valNodes (scalar) at position Xi
|
|
PURE module FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: valNodes(1:nNodes)
|
|
REAL(8):: f
|
|
REAL(8):: fPsi(1:nNodes)
|
|
|
|
fPsi = self%fPsi(Xi, nNodes)
|
|
f = DOT_PRODUCT(fPsi, valNodes)
|
|
|
|
END FUNCTION gatherF_cell_scalar
|
|
|
|
!Gather the value of valNodes (array) at position Xi
|
|
PURE module FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3)
|
|
REAL(8):: f(1:3)
|
|
REAL(8):: fPsi(1:nNodes)
|
|
|
|
fPsi = self%fPsi(Xi, nNodes)
|
|
f = MATMUL(fPsi, valNodes)
|
|
|
|
END FUNCTION gatherF_cell_array
|
|
|
|
!Gather the spatial derivative of valNodes (scalar) at position Xi
|
|
PURE module FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: valNodes(1:nNodes)
|
|
REAL(8):: df(1:3)
|
|
REAL(8):: dPsi(1:3, 1:nNodes)
|
|
REAL(8):: pDer(1:3,1:3)
|
|
REAL(8):: dPsiR(1:3, 1:nNodes)
|
|
REAL(8):: invJ(1:3, 1:3), detJ
|
|
|
|
dPsi = self%dPsi(Xi, nNodes)
|
|
pDer = self%partialDer(nNodes, dPsi)
|
|
detJ = self%detJac(pDer)
|
|
invJ = self%invJac(pDer)
|
|
dPsiR = MATMUL(invJ, dPsi)/detJ
|
|
df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), &
|
|
DOT_PRODUCT(dPsiR(2,:), valNodes), &
|
|
DOT_PRODUCT(dPsiR(3,:), valNodes) /)
|
|
|
|
END FUNCTION gatherDF_cell_scalar
|
|
|
|
!Scatters particle properties into cell nodes
|
|
module SUBROUTINE scatter(self, nNodes, part)
|
|
USE moduleMath
|
|
USE moduleSpecies
|
|
USE OMP_LIB
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(inout):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
CLASS(particle), INTENT(in):: part
|
|
REAL(8):: fPsi(1:nNodes)
|
|
INTEGER:: cellNodes(1:nNodes)
|
|
REAL(8):: tensorS(1:3, 1:3)
|
|
INTEGER:: sp
|
|
INTEGER:: i
|
|
CLASS(meshNode), POINTER:: node
|
|
REAL(8):: pFraction !Particle fraction
|
|
|
|
cellNodes = self%getNodes(nNodes)
|
|
fPsi = self%fPsi(part%Xi, nNodes)
|
|
|
|
tensorS = outerProduct(part%v, part%v)
|
|
|
|
sp = part%species%n
|
|
|
|
DO i = 1, nNodes
|
|
node => mesh%nodes(cellNodes(i))%obj
|
|
pFraction = fPsi(i)*part%weight
|
|
CALL OMP_SET_LOCK(node%lock)
|
|
node%output(sp)%den = node%output(sp)%den + pFraction
|
|
node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:)
|
|
node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
|
|
CALL OMP_UNSET_LOCK(node%lock)
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE scatter
|
|
|
|
!Find next cell for particle
|
|
RECURSIVE SUBROUTINE findCell(self, part, oldCell)
|
|
USE moduleSpecies
|
|
USE moduleErrors
|
|
USE OMP_LIB
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(inout):: self
|
|
CLASS(particle), INTENT(inout), TARGET:: part
|
|
CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell
|
|
REAL(8):: Xi(1:3)
|
|
CLASS(meshElement), POINTER:: neighbourElement
|
|
INTEGER:: sp
|
|
|
|
Xi = self%phy2log(part%r)
|
|
!Checks if particle is inside 'self' cell
|
|
IF (self%inside(Xi)) THEN
|
|
part%cell = self%n
|
|
part%Xi = Xi
|
|
part%n_in = .TRUE.
|
|
!Assign particle to listPart_in
|
|
IF (listInCells) THEN
|
|
CALL OMP_SET_LOCK(self%lock)
|
|
sp = part%species%n
|
|
CALL self%listPart_in(sp)%add(part)
|
|
self%totalWeight(sp) = self%totalWeight(sp) + part%weight
|
|
CALL OMP_UNSET_LOCK(self%lock)
|
|
|
|
END IF
|
|
|
|
ELSE
|
|
!If not, searches for a neighbour and repeats the process.
|
|
CALL self%neighbourElement(Xi, neighbourElement)
|
|
!Defines the next step
|
|
SELECT TYPE(neighbourElement)
|
|
CLASS IS(meshCell)
|
|
!Particle moved to new cell, repeat find procedure
|
|
CALL neighbourElement%findCell(part, self)
|
|
|
|
CLASS IS (meshEdge)
|
|
!Particle encountered a surface, apply boundary
|
|
CALL neighbourElement%boundariesParticle(part%species%n)%obj%apply(neighbourElement,part)
|
|
|
|
!If particle is still inside the domain, call findCell
|
|
IF (part%n_in) THEN
|
|
IF(PRESENT(oldCell)) THEN
|
|
CALL self%findCell(part, oldCell)
|
|
|
|
ELSE
|
|
CALL self%findCell(part)
|
|
|
|
END IF
|
|
END IF
|
|
|
|
CLASS DEFAULT
|
|
WRITE (*, "(A, I6)") "Element = ", self%n
|
|
CALL criticalError("No connectivity found for element", "findCell")
|
|
|
|
END SELECT
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE findCell
|
|
|
|
!If Coll and Particle are the same, simply copy the part%cell into part%cellColl
|
|
SUBROUTINE findCellSameMesh(part)
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
TYPE(particle), INTENT(inout):: part
|
|
|
|
part%cellColl = part%cell
|
|
|
|
END SUBROUTINE findCellSameMesh
|
|
|
|
!TODO: try to combine this with the findCell for a regular mesh
|
|
!Find the volume in which particle reside in the mesh for collisions
|
|
!No boundary interaction taken into account
|
|
SUBROUTINE findCellCollMesh(part)
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
TYPE(particle), INTENT(inout):: part
|
|
LOGICAL:: found
|
|
CLASS(meshCell), POINTER:: cell
|
|
REAL(8), DIMENSION(1:3):: Xi
|
|
CLASS(meshElement), POINTER:: neighbourElement
|
|
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 (listInCells) THEN
|
|
CALL OMP_SET_LOCK(cell%lock)
|
|
sp = part%species%n
|
|
CALL cell%listPart_in(sp)%add(part)
|
|
cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight
|
|
CALL OMP_UNSET_LOCK(cell%lock)
|
|
|
|
END IF
|
|
found = .TRUE.
|
|
|
|
ELSE
|
|
CALL cell%neighbourElement(Xi, neighbourElement)
|
|
SELECT TYPE(neighbourElement)
|
|
CLASS IS(meshCell)
|
|
!Try next element
|
|
cell => neighbourElement
|
|
|
|
CLASS DEFAULT
|
|
!Should never happend, but just in case, stops loops
|
|
found = .TRUE.
|
|
|
|
END SELECT
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE findCellCollMesh
|
|
|
|
!Returns index of volume associated to a position (if any)
|
|
!If no voulme is found, returns 0
|
|
!WARNING: This function is slow and should only be used in initialization phase
|
|
FUNCTION findCellBrute(self, r) RESULT(nVol)
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshGeneric), INTENT(in):: self
|
|
REAL(8), DIMENSION(1:3), INTENT(in):: r
|
|
INTEGER:: nVol
|
|
INTEGER:: e
|
|
REAL(8), DIMENSION(1:3):: Xi
|
|
|
|
!Inits RESULT
|
|
nVol = 0
|
|
|
|
DO e = 1, self%numCells
|
|
Xi = self%cells(e)%obj%phy2log(r)
|
|
IF(self%cells(e)%obj%inside(Xi)) THEN
|
|
nVol = self%cells(e)%obj%n
|
|
EXIT
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END FUNCTION findCellBrute
|
|
|
|
!Computes collisions in element
|
|
SUBROUTINE doCollisions(self)
|
|
USE moduleCollisions
|
|
USE moduleSpecies
|
|
USE moduleList
|
|
use moduleRefParam
|
|
USE moduleRandom
|
|
USE moduleOutput
|
|
USE moduleMath
|
|
USE moduleCaseParam, ONLY: timeStep
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshGeneric), INTENT(inout), TARGET:: self
|
|
INTEGER:: e
|
|
CLASS(meshCell), POINTER:: cell
|
|
INTEGER:: k, i, j
|
|
INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
|
|
REAL(8):: pMax !Maximum probability of collision
|
|
INTEGER:: nColl
|
|
TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
|
|
TYPE(particle), POINTER:: part_i, part_j
|
|
INTEGER:: n, c
|
|
REAL(8):: vRel, rMass, eRel
|
|
REAL(8):: sigmaVrelTotal
|
|
REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
|
|
REAL(8):: rnd_real !Random number for collision
|
|
INTEGER:: rnd_int !Random number for collision
|
|
|
|
IF (MOD(timeStep, everyColl) == 0) THEN
|
|
!Collisions need to be performed in this iteration
|
|
!$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
|
|
DO e=1, self%numCells
|
|
|
|
cell => self%cells(e)%obj
|
|
|
|
!TODO: Simplify this, to many sublevels
|
|
!Iterate over the number of pairs
|
|
DO k = 1, nCollPairs
|
|
!Reset tally of collisions
|
|
IF (collOutput) THEN
|
|
cell%tallyColl(k)%tally = 0
|
|
|
|
END IF
|
|
|
|
IF (interactionMatrix(k)%amount > 0) THEN
|
|
!Select the species for the collision pair
|
|
i = interactionMatrix(k)%sp_i%n
|
|
j = interactionMatrix(k)%sp_j%n
|
|
|
|
!Number of particles per species in the collision pair
|
|
nPart_i = cell%listPart_in(i)%amount
|
|
nPart_j = cell%listPart_in(j)%amount
|
|
|
|
IF (nPart_i > 0 .AND. nPart_j > 0) THEN
|
|
!Total number of particles for the collision pair
|
|
nPart = nPart_i + nPart_j
|
|
|
|
!Resets the number of collisions in the cell
|
|
nColl = 0
|
|
|
|
!Probability of collision for pair i-j
|
|
pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume
|
|
|
|
!Number of collisions in the cell
|
|
nColl = NINT(REAL(nPart)*pMax*0.5D0)
|
|
|
|
!Converts the list of particles to an array for easy access
|
|
IF (nColl > 0) THEN
|
|
partTemp_i = cell%listPart_in(i)%convert2Array()
|
|
partTemp_j = cell%listPart_in(j)%convert2Array()
|
|
|
|
END IF
|
|
|
|
DO n = 1, nColl
|
|
!Select random particles
|
|
part_i => NULL()
|
|
part_j => NULL()
|
|
rnd_int = random(1, nPart_i)
|
|
part_i => partTemp_i(rnd_int)%part
|
|
rnd_int = random(1, nPart_j)
|
|
part_j => partTemp_j(rnd_int)%part
|
|
!If they are the same particle, skip
|
|
!TODO: Maybe try to improve this
|
|
IF (ASSOCIATED(part_i, part_j)) THEN
|
|
CYCLE
|
|
|
|
END IF
|
|
|
|
!If particles do not belong to the species, skip collision
|
|
!This can happen, for example, if particle has been previously ionized or removed
|
|
!TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision?
|
|
IF (part_i%species%n /= i .OR. &
|
|
part_j%species%n /= j) THEN
|
|
CYCLE
|
|
|
|
END IF
|
|
!Obtain the cross sections for the different processes
|
|
!TODO: From here it might be a procedure in interactionMatrix
|
|
vRel = NORM2(part_i%v-part_j%v)
|
|
rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
|
|
eRel = rMass*vRel**2
|
|
CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
|
|
|
|
!Update maximum sigma*v_rel
|
|
IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN
|
|
cell%sigmaVrelMax(k) = sigmaVrelTotal
|
|
|
|
END IF
|
|
|
|
ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
|
|
probabilityColl = 0.0
|
|
DO c = 1, interactionMatrix(k)%amount
|
|
probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1))
|
|
|
|
END DO
|
|
|
|
!Selects random number between 0 and 1
|
|
rnd_real = random()
|
|
|
|
!If the random number is below the total probability of collision, collide particles
|
|
IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN
|
|
|
|
!Loop over collisions
|
|
DO c = 1, interactionMatrix(k)%amount
|
|
IF (rnd_real <= probabilityColl(c)) THEN
|
|
!$OMP CRITICAL
|
|
CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
|
|
!$OMP END CRITICAL
|
|
|
|
!If collisions are gonna be output, count the collision
|
|
IF (collOutput) THEN
|
|
cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1
|
|
|
|
END IF
|
|
|
|
!A collision has ocurred, exit the loop
|
|
EXIT
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
END IF
|
|
|
|
!Deallocate arrays for next collision
|
|
DEALLOCATE(sigmaVrel, probabilityColl)
|
|
|
|
!End loop collisions in cell
|
|
END DO
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!End loop collision pairs
|
|
END DO
|
|
|
|
!End loop volumes
|
|
END DO
|
|
!$OMP END DO
|
|
|
|
END IF
|
|
|
|
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(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))
|
|
deltaV_ij = 0.D0
|
|
p_ij = 0.D0
|
|
mass_ij = 0.D0
|
|
!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))
|
|
deltaV_ji = 0.D0
|
|
p_ji = 0.D0
|
|
mass_ji = 0.D0
|
|
!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
|
|
|
|
end submodule elements
|