Jorge Gonzalez
de1d4567f3
Due to a high convergence value (1.0e-2) in phy2logQuad (variable conv),
particles were being stuck in some elements, reaching a segmentation
fault. The new limit (1.0e-4) should avoid this.
955 lines
29 KiB
Fortran
955 lines
29 KiB
Fortran
!moduleMesh: General module for Finite Element mesh
|
|
MODULE moduleMesh
|
|
USE moduleList
|
|
USE moduleOutput
|
|
USE moduleBoundary
|
|
USE moduleCollisions
|
|
IMPLICIT NONE
|
|
|
|
!Generic mesh element
|
|
TYPE, PUBLIC, ABSTRACT:: meshElement
|
|
!Index
|
|
INTEGER:: n = 0
|
|
CONTAINS
|
|
|
|
END TYPE meshElement
|
|
|
|
!Parent of Node element
|
|
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshNode
|
|
!Node volume
|
|
REAL(8):: v = 0.D0
|
|
!Output values
|
|
TYPE(outputNode), ALLOCATABLE:: output(:)
|
|
TYPE(emNode):: emData
|
|
!Lock indicator for scattering
|
|
INTEGER(KIND=OMP_LOCK_KIND):: lock
|
|
CONTAINS
|
|
!DEFERED PROCEDURES
|
|
PROCEDURE(initNode_interface), DEFERRED, PASS:: init
|
|
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
|
|
!GENERIC PROCEDURES
|
|
PROCEDURE, PASS:: resetOutput
|
|
|
|
END TYPE meshNode
|
|
|
|
ABSTRACT INTERFACE
|
|
!Interface of init a node (3D generic coordinates)
|
|
SUBROUTINE initNode_interface(self, n, r)
|
|
IMPORT:: meshNode
|
|
CLASS(meshNode), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
|
|
END SUBROUTINE initNode_interface
|
|
|
|
!Interface to get coordinates from node
|
|
PURE FUNCTION getCoord_interface(self) RESULT(r)
|
|
IMPORT:: meshNode
|
|
CLASS(meshNode), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
END FUNCTION getCoord_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Containers for nodes in the mesh
|
|
TYPE:: meshNodeCont
|
|
CLASS(meshNode), ALLOCATABLE:: obj
|
|
CONTAINS
|
|
|
|
END TYPE meshNodeCont
|
|
|
|
!Type for array of boundary functions (one per species)
|
|
TYPE, PUBLIC:: fBoundaryGeneric
|
|
PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
|
|
CONTAINS
|
|
|
|
END TYPE
|
|
|
|
!Parent of Edge element
|
|
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge
|
|
!Nomber of nodes in the edge
|
|
INTEGER:: nNodes
|
|
!Connectivity to cells
|
|
CLASS(meshCell), POINTER:: e1 => NULL(), e2 => NULL()
|
|
!Connectivity to cells in meshColl
|
|
CLASS(meshCell), POINTER:: eColl => NULL()
|
|
!Normal vector
|
|
REAL(8):: normal(1:3)
|
|
!Weight for random injection of particles
|
|
REAL(8):: weight = 1.D0
|
|
!Pointer to boundary type
|
|
TYPE(boundaryCont), POINTER:: boundary
|
|
!Array of functions for boundary conditions
|
|
TYPE(fBoundaryGeneric), ALLOCATABLE:: fBoundary(:)
|
|
!Physical surface for the edge
|
|
INTEGER:: physicalSurface
|
|
CONTAINS
|
|
!DEFERED PROCEDURES
|
|
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
|
|
PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
|
|
PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
|
|
PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos
|
|
|
|
END TYPE meshEdge
|
|
|
|
ABSTRACT INTERFACE
|
|
!Inits the edge parameters
|
|
SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
|
|
IMPORT:: meshEdge
|
|
|
|
CLASS(meshEdge), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
INTEGER, INTENT(in):: p(:)
|
|
INTEGER, INTENT(in):: bt
|
|
INTEGER, INTENT(in):: physicalSurface
|
|
|
|
END SUBROUTINE initEdge_interface
|
|
|
|
!Get nodes index from node
|
|
PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n)
|
|
IMPORT:: meshEdge
|
|
CLASS(meshEdge), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
INTEGER:: n(1:nNodes)
|
|
|
|
END FUNCTION getNodesEdge_interface
|
|
|
|
!Returns the intersecction between an edge and a line defined by point r0
|
|
PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r)
|
|
IMPORT:: meshEdge
|
|
CLASS(meshEdge), INTENT(in):: self
|
|
REAL(8), INTENT(in), DIMENSION(1:3):: r0
|
|
REAL(8):: r(1:3)
|
|
|
|
END FUNCTION intersectionEdge_interface
|
|
|
|
!Returns a random position in the edge
|
|
FUNCTION randPosEdge_interface(self) RESULT(r)
|
|
IMPORT:: meshEdge
|
|
CLASS(meshEdge), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
END FUNCTION randPosEdge_interface
|
|
|
|
END INTERFACE
|
|
|
|
INTERFACE
|
|
SUBROUTINE boundary_interface(edge, part)
|
|
USE moduleSpecies
|
|
|
|
IMPORT:: meshEdge
|
|
CLASS (meshEdge), INTENT(inout):: edge
|
|
CLASS (particle), INTENT(inout):: part
|
|
|
|
END SUBROUTINE
|
|
|
|
END INTERFACE
|
|
|
|
!Containers for edges in the mesh
|
|
TYPE:: meshEdgeCont
|
|
CLASS(meshEdge), ALLOCATABLE:: obj
|
|
|
|
END TYPE meshEdgeCont
|
|
|
|
!Parent of cell element
|
|
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshCell
|
|
!Number of nodes in the cell
|
|
INTEGER:: nNodes
|
|
!Maximum collision rate
|
|
REAL(8), ALLOCATABLE:: sigmaVrelMax(:)
|
|
!Arrays for counting number of collisions
|
|
TYPE(tallyCollisions), ALLOCATABLE:: tallyColl(:)
|
|
!Volume
|
|
REAL(8):: volume = 0.D0
|
|
!List of particles inside the volume
|
|
TYPE(listNode), ALLOCATABLE:: listPart_in(:)
|
|
!Lock indicator for listPart_in
|
|
INTEGER(KIND=OMP_LOCK_KIND):: lock
|
|
!Total weight of particles inside cell
|
|
REAL(8), ALLOCATABLE:: totalWeight(:)
|
|
CONTAINS
|
|
!DEFERRED PROCEDURES
|
|
!Init the cell
|
|
PROCEDURE(initCell_interface), DEFERRED, PASS:: init
|
|
!Get the index of the nodes
|
|
PROCEDURE(getNodesCell_interface), DEFERRED, PASS:: getNodes
|
|
!Calculate random position on the cell
|
|
PROCEDURE(randPosCell_interface), DEFERRED, PASS:: randPos
|
|
!Obtain functions and values of cell natural functions
|
|
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
|
|
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
|
|
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
|
|
PROCEDURE(detJac_interface), DEFERRED, NOPASS:: detJac
|
|
PROCEDURE(invJac_interface), DEFERRED, NOPASS:: invJac
|
|
!Procedures to get specific values in the node
|
|
PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherElectricField
|
|
PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherMagneticField
|
|
!Compute K and F to solve PDE on the mesh
|
|
PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK
|
|
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
|
|
!Check if particle is inside the cell
|
|
PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside
|
|
!Convert physical coordinates (r) into logical coordinates (Xi)
|
|
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
|
|
!Returns the neighbour element based on particle position outside the cell
|
|
PROCEDURE(neighbourElement_interface), DEFERRED, PASS:: neighbourElement
|
|
!Scatter properties of particles on cell nodes
|
|
PROCEDURE, PASS:: scatter
|
|
!Subroutine to find in which cell a particle is located
|
|
PROCEDURE, PASS:: findCell
|
|
!Gather value and spatial derivative on the nodes at position Xi
|
|
PROCEDURE, PASS, PRIVATE:: gatherF_scalar
|
|
PROCEDURE, PASS, PRIVATE:: gatherF_array
|
|
PROCEDURE, PASS, PRIVATE:: gatherDF_scalar
|
|
GENERIC:: gatherF => gatherF_scalar, gatherF_array
|
|
GENERIC:: gatherDF => gatherDF_scalar
|
|
|
|
END TYPE meshCell
|
|
|
|
ABSTRACT INTERFACE
|
|
SUBROUTINE initCell_interface(self, n, p, nodes)
|
|
IMPORT:: meshCell
|
|
IMPORT meshNodeCont
|
|
CLASS(meshCell), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
INTEGER, INTENT(in):: p(:)
|
|
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
|
|
|
|
END SUBROUTINE initCell_interface
|
|
|
|
PURE FUNCTION getNodesCell_interface(self, nNodes) RESULT(n)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
INTEGER:: n(1:nNodes)
|
|
|
|
END FUNCTION getNodesCell_interface
|
|
|
|
FUNCTION randPosCell_interface(self) RESULT(r)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
END FUNCTION randPosCell_interface
|
|
|
|
PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi)
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8):: fPsi(1:nNodes)
|
|
|
|
END FUNCTION fPsi_interface
|
|
|
|
PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi)
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8):: dPsi(1:3, 1:nNodes)
|
|
|
|
END FUNCTION dPsi_interface
|
|
|
|
PURE FUNCTION partialDer_interface(self, nNodes, dPsi) RESULT(pDer)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes)
|
|
REAL(8):: pDer(1:3, 1:3)
|
|
|
|
END FUNCTION partialDer_interface
|
|
|
|
PURE FUNCTION detJac_interface(pDer) RESULT(dJ)
|
|
REAL(8), INTENT(in):: pDer(1:3,1:3)
|
|
REAL(8):: dJ
|
|
|
|
END FUNCTION detJac_interface
|
|
|
|
PURE FUNCTION invJac_interface(pDer) RESULT(invJ)
|
|
REAL(8), INTENT(in):: pDer(1:3,1:3)
|
|
REAL(8):: invJ(1:3,1:3)
|
|
|
|
END FUNCTION invJac_interface
|
|
|
|
PURE FUNCTION gatherArray_interface(self, Xi) RESULT(array)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
REAL(8):: array(1:3)
|
|
|
|
END FUNCTION gatherArray_interface
|
|
|
|
PURE FUNCTION elemK_interface(self, nNodes) RESULT(localK)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8):: localK(1:nNodes,1:nNodes)
|
|
|
|
END FUNCTION elemK_interface
|
|
|
|
PURE FUNCTION elemF_interface(self, nNodes, source) RESULT(localF)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: source(1:nNodes)
|
|
REAL(8):: localF(1:nNodes)
|
|
|
|
END FUNCTION elemF_interface
|
|
|
|
PURE FUNCTION inside_interface(Xi) RESULT(ins)
|
|
IMPORT:: meshCell
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
LOGICAL:: ins
|
|
|
|
END FUNCTION inside_interface
|
|
|
|
PURE FUNCTION phy2log_interface(self,r) RESULT(Xi)
|
|
IMPORT:: meshCell
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
REAL(8):: Xi(1:3)
|
|
|
|
END FUNCTION phy2log_interface
|
|
|
|
SUBROUTINE neighbourElement_interface(self, Xi, neighbourElement)
|
|
IMPORT:: meshCell, meshElement
|
|
CLASS(meshCell), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
|
|
|
|
END SUBROUTINE neighbourElement_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Containers for cells in the mesh
|
|
TYPE:: meshCellCont
|
|
CLASS(meshCell), ALLOCATABLE:: obj
|
|
|
|
END TYPE meshCellCont
|
|
|
|
!Generic mesh type
|
|
TYPE, ABSTRACT:: meshGeneric
|
|
!Dimension of the mesh
|
|
INTEGER:: dimen
|
|
!Geometry of the mesh
|
|
CHARACTER(:), ALLOCATABLE:: geometry
|
|
!Number of elements
|
|
INTEGER:: numNodes, numCells
|
|
!Array of nodes
|
|
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
|
|
!Array of cell elements
|
|
TYPE(meshCellCont), ALLOCATABLE:: cells(:)
|
|
!PROCEDURES SPECIFIC OF FILE TYPE
|
|
PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL()
|
|
PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL()
|
|
PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL()
|
|
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
|
|
CONTAINS
|
|
!GENERIC PROCEDURES
|
|
PROCEDURE, PASS:: doCollisions
|
|
|
|
END TYPE meshGeneric
|
|
|
|
ABSTRACT INTERFACE
|
|
!Reads the mesh from a file
|
|
SUBROUTINE readMesh_interface(self, filename)
|
|
IMPORT meshGeneric
|
|
CLASS(meshGeneric), INTENT(inout):: self
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
|
|
|
|
END SUBROUTINE readMesh_interface
|
|
|
|
SUBROUTINE readInitial_interface(filename, density, velocity, temperature)
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
|
|
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
|
|
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
|
|
REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
|
|
|
|
END SUBROUTINE readInitial_interface
|
|
|
|
!Connects cell and edges to the mesh
|
|
SUBROUTINE connectMesh_interface(self)
|
|
IMPORT meshGeneric
|
|
CLASS(meshGeneric), INTENT(inout):: self
|
|
|
|
END SUBROUTINE connectMesh_interface
|
|
|
|
!Prints number of collisions in each cell
|
|
SUBROUTINE printColl_interface(self, t)
|
|
IMPORT meshGeneric
|
|
CLASS(meshGeneric), INTENT(in):: self
|
|
INTEGER, INTENT(in):: t
|
|
|
|
END SUBROUTINE printColl_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Particle mesh
|
|
TYPE, EXTENDS(meshGeneric), PUBLIC:: meshParticles
|
|
INTEGER:: numEdges
|
|
!Array of boundary elements
|
|
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
|
|
!Global stiffness matrix
|
|
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
|
|
|
|
END TYPE meshParticles
|
|
|
|
ABSTRACT INTERFACE
|
|
!Prints Species data
|
|
SUBROUTINE printOutput_interface(self, t)
|
|
IMPORT meshParticles
|
|
CLASS(meshParticles), INTENT(in):: self
|
|
INTEGER, INTENT(in):: t
|
|
|
|
END SUBROUTINE printOutput_interface
|
|
|
|
!Prints EM info
|
|
SUBROUTINE printEM_interface(self, t)
|
|
IMPORT meshParticles
|
|
CLASS(meshParticles), INTENT(in):: self
|
|
INTEGER, INTENT(in):: t
|
|
|
|
END SUBROUTINE printEM_interface
|
|
|
|
!Perform Coulomb Scattering
|
|
SUBROUTINE doCoulomb_interface(self)
|
|
IMPORT meshParticles
|
|
CLASS(meshParticles), INTENT(inout):: self
|
|
|
|
END SUBROUTINE doCoulomb_interface
|
|
|
|
!Prints average values
|
|
SUBROUTINE printAverage_interface(self)
|
|
IMPORT meshParticles
|
|
CLASS(meshParticles), INTENT(in):: self
|
|
|
|
END SUBROUTINE printAverage_interface
|
|
|
|
END INTERFACE
|
|
|
|
TYPE(meshParticles), TARGET:: mesh
|
|
|
|
!Collision (MCC) mesh
|
|
TYPE, EXTENDS(meshGeneric):: meshCollisions
|
|
CONTAINS
|
|
!GENERIC PROCEDURES
|
|
|
|
END TYPE meshCollisions
|
|
|
|
TYPE(meshCollisions), TARGET:: meshColl
|
|
|
|
ABSTRACT INTERFACE
|
|
SUBROUTINE readMeshColl_interface(self, filename)
|
|
IMPORT meshCollisions
|
|
CLASS(meshCollisions), INTENT(inout):: self
|
|
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
|
|
|
|
END SUBROUTINE readMeshColl_interface
|
|
|
|
SUBROUTINE connectMeshColl_interface(self)
|
|
IMPORT meshParticles
|
|
|
|
CLASS(meshParticles), INTENT(inout):: self
|
|
|
|
END SUBROUTINE connectMeshColl_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Pointer to mesh used for MC collisions
|
|
CLASS(meshGeneric), POINTER:: meshForMCC => NULL()
|
|
|
|
!Procedure to find a cell for a particle in meshColl
|
|
PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL()
|
|
|
|
ABSTRACT INTERFACE
|
|
SUBROUTINE findCellColl_interface(part)
|
|
USE moduleSpecies
|
|
|
|
TYPE(particle), INTENT(inout):: part
|
|
|
|
END SUBROUTINE findCellColl_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Logical to indicate if an specific mesh for MC Collisions is used
|
|
LOGICAL:: doubleMesh
|
|
!Logical to indicate if MCC collisions are performed
|
|
LOGICAL:: doMCC
|
|
!Complete path for the two meshes
|
|
CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle
|
|
|
|
CONTAINS
|
|
!Constructs the global K matrix
|
|
PURE SUBROUTINE constructGlobalK(self)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshParticles), INTENT(inout):: self
|
|
INTEGER:: e
|
|
INTEGER:: nNodes
|
|
INTEGER, ALLOCATABLE:: n(:)
|
|
REAL(8), ALLOCATABLE:: localK(:,:)
|
|
INTEGER:: i, j
|
|
|
|
DO e = 1, self%numCells
|
|
nNodes = self%cells(e)%obj%nNodes
|
|
ALLOCATE(n(1:nNodes))
|
|
ALLOCATE(localK(1:nNodes, 1:nNodes))
|
|
n = self%cells(e)%obj%getNodes(nNodes)
|
|
localK = self%cells(e)%obj%elemK(nNodes)
|
|
|
|
DO i = 1, nNodes
|
|
DO j = 1, nNodes
|
|
self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
|
|
|
|
END DO
|
|
|
|
END DO
|
|
|
|
DEALLOCATE(n, localK)
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE constructGlobalK
|
|
|
|
!Reset the output of node
|
|
PURE 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
|
|
|
|
!Gather the value of valNodes (scalar) at position Xi
|
|
PURE FUNCTION gatherF_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_scalar
|
|
|
|
!Gather the value of valNodes (array) at position Xi
|
|
PURE FUNCTION gatherF_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_array
|
|
|
|
!Gather the spatial derivative of valNodes (scalar) at position Xi
|
|
PURE FUNCTION gatherDF_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_scalar
|
|
|
|
!Scatters particle properties into cell nodes
|
|
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
|
|
|
|
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
|
|
CALL OMP_SET_LOCK(node%lock)
|
|
node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i)
|
|
node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:)
|
|
node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*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) = 0.D0
|
|
CLASS(meshElement), POINTER:: neighbourElement => NULL()
|
|
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 (doMCC) 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%fBoundary(part%species%n)%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 (doMCC) 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, t)
|
|
USE moduleCollisions
|
|
USE moduleSpecies
|
|
USE moduleList
|
|
use moduleRefParam
|
|
USE moduleRandom
|
|
USE moduleOutput
|
|
USE moduleMath
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshGeneric), INTENT(inout), TARGET:: self
|
|
INTEGER, INTENT(in):: t
|
|
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(t, 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 no lose these collisions. Maybe check new 'k' and use that for the collision, maybe?
|
|
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
|
|
CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
|
|
|
|
!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)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshParticles), INTENT(inout):: self
|
|
|
|
END SUBROUTINE doCoulomb
|
|
|
|
END MODULE moduleMesh
|