!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) ! Surface of edge REAL(8):: surface = 0.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(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS !GENERIC PROCEDURES PROCEDURE, PASS:: constructGlobalK PROCEDURE, PASS:: doCoulomb 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:: 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 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 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%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 (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, 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 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 MODULE moduleMesh