fpakc/src/modules/mesh/moduleMesh.f90
Jorge Gonzalez ac2965621a Structure for 3D Cartesian Grid created.
Unification of boundary conditions into one file.

Some changes to input file for reference cases. This should have been
done in another branch but I wanto to commit to save progress and I
don't want to deal with tswitching branches right now, I'm very busy
watching Futurama.
2021-02-27 16:24:44 +01:00

649 lines
19 KiB
Fortran

!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
USE moduleList
USE moduleOutput
USE moduleBoundary
IMPLICIT NONE
!Parent of Node element
TYPE, PUBLIC, ABSTRACT:: meshNode
!Node index
INTEGER:: n = 0
!Node volume
REAL(8):: v = 0.D0
!Output values
TYPE(outputNode), ALLOCATABLE:: output(:)
TYPE(emNode):: emData
CONTAINS
PROCEDURE(initNode_interface), DEFERRED, PASS:: init
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
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
END TYPE meshNodeCont
!Type for array of boundary functions (one per species)
TYPE, PUBLIC:: fBoundaryGeneric
PROCEDURE(boundary_interface), POINTER, NOPASS:: apply => NULL()
END TYPE
!Parent of Edge element
TYPE, PUBLIC, ABSTRACT:: meshEdge
!Element index
INTEGER:: n = 0
!Connectivity to vols
CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL()
!Normal vector
REAL(8):: normal(1:3)
!Pointer to boundary element
TYPE(boundaryCont), POINTER:: boundary
!Array of functions for boundary conditions
TYPE(fBoundaryGeneric), ALLOCATABLE:: fBoundary(:)
!Physical surface for the edge
INTEGER:: physicalSurface
CONTAINS
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) RESULT(n)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION getNodesEdge_interface
!Returns the intersecction between an edge and a line defined by point r0 and vector v0
PURE FUNCTION intersectionEdge_interface(self, r0, v0) RESULT(r)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
REAL(8), INTENT(in), DIMENSION(1:3):: r0, v0
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 Volume element
TYPE, PUBLIC, ABSTRACT:: meshVol
!Volume index
INTEGER:: n = 0
!Maximum collision rate
REAL(8):: sigmaVrelMax = 0.D0
!Volume
REAL(8):: volume = 0.D0
!List of particles inside the volume
TYPE(listNode):: listPart_in
!Lock indicator for listPart_in
INTEGER(KIND=OMP_LOCK_KIND):: lock
!Number of collisions per volume
INTEGER:: nColl = 0
!Collisional fraction
REAL(8):: collFrac = 0.D0
!Total weight of particles inside cell
REAL(8):: totalWeight = 0.D0
CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos
PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
PROCEDURE, PASS:: findCell
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside
PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement
PROCEDURE, PASS:: collision
END TYPE meshVol
ABSTRACT INTERFACE
SUBROUTINE initVol_interface(self, n, p)
IMPORT:: meshVol
CLASS(meshVol), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
END SUBROUTINE initVol_interface
SUBROUTINE scatter_interface(self, part)
USE moduleSpecies
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
CLASS(particle), INTENT(in):: part
END SUBROUTINE scatter_interface
PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: EF(1:3)
END FUNCTION gatherEF_interface
PURE FUNCTION getNodesVol_interface(self) RESULT(n)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION getNodesVol_interface
PURE FUNCTION elemF_interface(self, source) RESULT(localF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
END FUNCTION elemF_interface
SUBROUTINE nextElement_interface(self, xi, nextElement)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
END SUBROUTINE nextElement_interface
PURE FUNCTION phy2log_interface(self,r) RESULT(xN)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
END FUNCTION phy2log_interface
PURE FUNCTION inside_interface(xi) RESULT(ins)
IMPORT:: meshVol
REAL(8), INTENT(in):: xi(1:3)
LOGICAL:: ins
END FUNCTION inside_interface
FUNCTION randPosVol_interface(self) RESULT(r)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8):: r(1:3)
END FUNCTION randPosVol_interface
END INTERFACE
!Containers for volumes in the mesh
TYPE:: meshVolCont
CLASS(meshVol), ALLOCATABLE:: obj
END TYPE meshVolCont
!Abstract type of mesh
TYPE, PUBLIC, ABSTRACT:: meshGeneric
INTEGER:: numEdges, numNodes, numVols
!Array of nodes
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of boundary elements
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
!Array of volume elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:)
!Global stiffness matrix
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Permutation matrix for P L U factorization
INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
CONTAINS
PROCEDURE(initMesh_interface), DEFERRED, PASS:: init
PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh
END TYPE meshGeneric
ABSTRACT INTERFACE
!Inits the mesh
SUBROUTINE initMesh_interface(self, meshFormat)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
END SUBROUTINE initMesh_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
!Prints Species data
SUBROUTINE printOutput_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints number of collisions
SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
!Prints EM info
SUBROUTINE printEM_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface
END INTERFACE
!Generic mesh
CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh
CONTAINS
!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
!Find next cell for particle
RECURSIVE SUBROUTINE findCell(self, part, oldCell)
USE moduleSpecies
USE OMP_LIB
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
CLASS(particle), INTENT(inout), TARGET:: part
REAL(8):: xi(1:3)
CLASS(*), POINTER:: nextElement
xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell
IF (self%inside(xi)) THEN
part%vol = self%n
part%xi = xi
part%n_in = .TRUE.
!Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock)
CALL self%listPart_in%add(part)
self%totalWeight = self%totalWeight + part%weight
CALL OMP_UNSET_LOCK(self%lock)
ELSE
!If not, searches for a neighbour and repeats the process.
CALL self%nextElement(xi, nextElement)
!Defines the next step
SELECT TYPE(nextElement)
CLASS IS(meshVol)
!Particle moved to new cell, repeat find procedure
CALL nextElement%findCell(part, self)
CLASS IS (meshEdge)
!Particle encountered a surface, apply boundary
CALL nextElement%fBoundary(part%sp)%apply(nextElement,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(*,*) "ERROR, CHECK findCell"
END SELECT
END IF
END SUBROUTINE findCell
!Computes collisions in element
SUBROUTINE collision(self)
USE moduleCollisions
USE moduleSpecies
USE moduleList
use moduleRefParam
USE moduleRandom
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
INTEGER:: nPart !Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision
INTEGER:: rnd !random index
TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision
INTEGER:: ij, k
REAL(8):: sigmaVrelMaxNew
TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
nPart = self%listPart_in%amount
!Computes iterations if there is more than one particle in the cell
IF (nPart > 1) THEN
!Probability of collision
pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
!Increases the collisional fraction of the cell
self%collFrac = self%collFrac + REAL(nPart)*pMax*0.5D0
!Number of collisions in the cell
self%nColl = FLOOR(self%collFrac)
IF (self%nColl > 0) THEN
!Converts the list of particles to an array for easy access
partTemp = self%listPart_in%convert2Array()
END IF
DO n = 1, self%nColl
!Select random numbers
rnd = random(1, nPart)
part_i => partTemp(rnd)%part
rnd = random(1, nPart)
part_j => partTemp(rnd)%part
ij = interactionIndex(part_i%sp, part_j%sp)
sigmaVrelMaxNew = 0.D0
DO k = 1, interactionMatrix(ij)%amount
CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
END DO
!Update maximum cross section*v_rel per each collision
IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
self%sigmaVrelMax = sigmaVrelMaxNew
END IF
!Removes one collision from the collisional fraction
self%collFrac = self%collFrac - 1.D0
END DO
END IF
END SUBROUTINE collision
SUBROUTINE printOutputGmsh(self, t)
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring
time = DBLE(t)*tauMin*ti_ref
DO i = 1, nSpecies
WRITE(tstring, iterationFormat) t
fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Density ' // species(i)%obj%name // ' (m^-3)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Velocity ' // species(i)%obj%name // ' (m/s)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 3
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Pressure ' // species(i)%obj%name // ' (Pa)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
END DO
WRITE(60, "(A)") '$EndNodeData'
WRITE(60, "(A)") '$NodeData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Temperature ' // species(i)%obj%name // ' (K)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
WRITE(60, "(A)") '$EndNodeData'
CLOSE (60)
END DO
END SUBROUTINE printOutputGmsh
SUBROUTINE printCollGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleCollisions
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring
IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, iterationFormat) t
fileName='OUTPUT_' // tstring// '_Collisions.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (60, file = path // folder // '/' // fileName)
WRITE(60, "(A)") '$MeshFormat'
WRITE(60, "(A)") '2.2 0 8'
WRITE(60, "(A)") '$EndMeshFormat'
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
END DO
WRITE(60, "(A)") '$EndElementData'
CLOSE(60)
END IF
END SUBROUTINE printCollGmsh
SUBROUTINE printEMGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=iterationDigits):: tstring
REAL(8):: xi(1:3)
xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, iterationFormat) t
fileName='OUTPUT_' // tstring// '_EMField.msh'
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (20, file = path // folder // '/' // fileName)
WRITE(20, "(A)") '$MeshFormat'
WRITE(20, "(A)") '2.2 0 8'
WRITE(20, "(A)") '$EndMeshFormat'
WRITE(20, "(A)") '$NodeData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Potential (V)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 1
WRITE(20, *) self%numNodes
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO
WRITE(20, "(A)") '$EndNodeData'
WRITE(20, "(A)") '$ElementData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Electric Field (V/m)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 3
WRITE(20, *) self%numVols
DO e=1, self%numVols
WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
END DO
WRITE(20, "(A)") '$EndElementData'
CLOSE(20)
END IF
END SUBROUTINE printEMGmsh
END MODULE moduleMesh